From 88a1f31f0cdddf9dd62496db5183ea3c5facc86a Mon Sep 17 00:00:00 2001
From: naumenko-sa <evolgenomicslab@gmail.com>
Date: Tue, 21 Aug 2018 11:46:28 -0400
Subject: [PATCH] 	modified:   cre.R 	modified:   cre.gemini2txt.sh 
 modified:   cre.gemini_variant_impacts.sh 	modified:   cre.sh 	new
 file:   cre.vcfanno.sh

---
 cre.R                         |  8 ++++----
 cre.gemini2txt.sh             |  2 +-
 cre.gemini_variant_impacts.sh |  2 +-
 cre.sh                        | 13 ++++++++++---
 cre.vcfanno.sh                | 12 ++++++++++++
 5 files changed, 28 insertions(+), 9 deletions(-)
 create mode 100644 cre.vcfanno.sh

diff --git a/cre.R b/cre.R
index 5789d09..284c951 100644
--- a/cre.R
+++ b/cre.R
@@ -166,7 +166,7 @@ create_report = function(family,samples)
     # Column 17 - Ensembl_gene_id
 
     # Column 18 - Gene_description
-    gene_descriptions = read.delim2(paste0(default_tables_path,"/ensembl_w_description.txt"), stringsAsFactors=F)
+    gene_descriptions = read.delim2(paste0(default_tables_path,"/data/ensembl_w_description.txt"), stringsAsFactors=F)
     variants = merge(variants,gene_descriptions,by.x = "Ensembl_gene_id",by.y = "ensembl_gene_id",all.x=T)
     
     # Column 19 - Omim_gene_description - from omim text file
@@ -234,7 +234,7 @@ create_report = function(family,samples)
     # Columns 30-35 - population frequencies
 
     # Columns 36-37, Exac scores
-    exac_scores_file = paste0(default_tables_path,"/exac_scores.txt")
+    exac_scores_file = paste0(default_tables_path,"/data/exac_scores.txt")
     exac_scores = read.delim(exac_scores_file, stringsAsFactors=F)
     variants = merge(variants,exac_scores,all.x=T)
 
@@ -251,12 +251,12 @@ create_report = function(family,samples)
 
     
     # Columns 44,45 - imprinting
-    imprinting_file_name = paste0(default_tables_path,"/imprinting.txt")
+    imprinting_file_name = paste0(default_tables_path,"/data/imprinting.txt")
     imprinting = read.delim(imprinting_file_name, stringsAsFactors=F)
     variants = merge(variants,imprinting,all.x=T)
     
     # Column 46 - pseudoautosomal
-    pseudoautosomal_file_name = paste0(default_tables_path,"/pseudoautosomal.txt")
+    pseudoautosomal_file_name = paste0(default_tables_path,"/data/pseudoautosomal.txt")
     pseudoautosomal = read.delim(pseudoautosomal_file_name, stringsAsFactors=F)
     variants = merge(variants,pseudoautosomal,all.x=T)
     
diff --git a/cre.gemini2txt.sh b/cre.gemini2txt.sh
index 50d194f..a7c6f91 100755
--- a/cre.gemini2txt.sh
+++ b/cre.gemini2txt.sh
@@ -83,7 +83,7 @@ sQuery=$sQuery"v.vep_hgvsc as Nucleotide_change_ensembl,
 	        where 
 	        v.transcript=g.transcript and 
 	        (v.gene=g.gene or g.gene is NULL) "$severity_filter" and 
-	        v.max_aaf_all < 0.01 and 
+	        v.max_aaf_all <= "$max_af" and 
 	        (v.depth >= "$depth_threshold" or v.depth = '' or v.depth is null)"
 
 echo $sQuery
diff --git a/cre.gemini_variant_impacts.sh b/cre.gemini_variant_impacts.sh
index 7264dd0..31e66a4 100755
--- a/cre.gemini_variant_impacts.sh
+++ b/cre.gemini_variant_impacts.sh
@@ -64,7 +64,7 @@ fi
 
 sQuery=$sQuery" from variants v, 
 		variant_impacts i 
-		where "$severity_filter"v.max_aaf_all < 0.01 and 
+		where "$severity_filter"v.max_aaf_all <= "$max_af" and 
 		v.variant_id=i.variant_id and 
 		(v.depth>="$depth_threshold" or v.depth='' or v.depth is null)"
 
diff --git a/cre.sh b/cre.sh
index e62c2ae..05ceb49 100755
--- a/cre.sh
+++ b/cre.sh
@@ -11,6 +11,7 @@
 # 	cleanup= [0|1] default = 0
 # 	make_report=[0|1] default = 1
 # 	type = [ wes.regular (default) | wes.synonymous | wes.fast | rnaseq | wgs | vcf2db (new wes/wgs report with gemini loaded by vcf2db)]
+#	max_af = af filter, default = 0.01
 ####################################################################################################
 
 #PBS -l walltime=20:00:00,nodes=1:ppn=1
@@ -109,8 +110,8 @@ function f_make_report
 	cre.gemini2txt.vcf2db.sh ${family}-ensemble.db $depth_threshold $severity_filter > $family.variants.txt
 	cre.gemini.variant_impacts.vcf2db.sh ${family}-ensemble.db $depth_threshold $severity_filter > $family.variant_impacts.txt
     else
-	cre.gemini2txt.sh ${family}-ensemble.db $depth_threshold $severity_filter
-	cre.gemini_variant_impacts.sh ${family}-ensemble.db $depth_threshold $severity_filter
+	cre.gemini2txt.sh ${family}-ensemble.db $depth_threshold $severity_filter $max_af
+	cre.gemini_variant_impacts.sh ${family}-ensemble.db $depth_threshold $severity_filter $maf_af
     fi
 
     for f in *.vcf.gz;
@@ -125,7 +126,7 @@ function f_make_report
     then
         cat $family.variants.txt | cut -f 26,27 | sed 1d | sed s/chr// | sort -k1,1 -k2,2n > ${family}-ensemble.db.txt.positions
     else
-        cat ${family}-ensemble.db.txt | cut -f 23,24  | sed 1d | sed s/chr// | sort -k1,1 -k2,2n > ${family}-ensemble.db.txt.positions
+        cat ${family}-ensemble.db.txt | cut -f 24,25  | sed 1d | sed s/chr// | sort -k1,1 -k2,2n > ${family}-ensemble.db.txt.positions
     fi
 
     # this may produce duplicate records if two positions from positions file overlap with a variant 
@@ -221,6 +222,12 @@ then
     make_report=1
 fi 
 
+if [ -z $max_af ]
+then
+    max_af=0.01
+fi
+export max_af
+
 if [ $make_report -eq 1 ]
 then
     f_make_report
diff --git a/cre.vcfanno.sh b/cre.vcfanno.sh
new file mode 100644
index 0000000..1e440c5
--- /dev/null
+++ b/cre.vcfanno.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+
+#PBS -l walltime=23:00:00,nodes=1:ppn=1
+#PBS -joe .
+#PBS -d .
+#PBS -l vmem=10g,mem=10g
+
+vcfanno -p 7 -lua /home/naumenko/cre/cre.vcfanno.lua \
+	     -base-path /hpf/largeprojects/ccmbio/naumenko/tools/bcbio/gemini_data \
+	     /home/naumenko/cre/cre.vcfanno.conf \
+	     $vcf | sed -e 's/Number=A/Number=1/g' | bgzip -c > `echo $vcf | sed s/vcf/annotated.vcf/`
+