fully working exome pipeline

Author: Dmitry Pushkarev

apply_bed_filter/apply_bed_filter.pl

@@ -0,0 +1,52 @@
+#!/usr/bin/perl
+#
+# Applies BED filter to a VCF file
+$filter=$ARGV[0];
+die "Filter file is not specified" unless -e $filter;
+die "Incorrect filter format" unless $filter =~/\.(bed|BED)/;
+
+
+#hash based coarse filter
+open(FL,"$filter")||die;
+while (<FL>)
+{
+    @a=split("\t");
+    $a[1]=int($a[1]/10000);
+    $a[2]=int($a[2]/10000);
+
+    for ($i=$a[1]-1;$i<=$a[2]+1;$i++)
+    {
+	$OK{"$a[0]-$i"}=1;
+    }
+}
+$tmp=int(rand()*1000000)."bed";
+$tmp="out.bed";
+
+open(OUT,">$tmp");
+
+while (<STDIN>)
+{
+    if (/^\#/)
+    {$header.=$_;next;}
+    @a=split("\t");
+
+    $red=int($a[1]/10000);
+
+    next unless $OK{"$a[0]-$red"};
+    print OUT "$a[0]\t$a[1]\t$a[1]\t$_";
+}
+
+close OUT;
+
+print $header;
+open(FL,"bedtools intersect -u -a $tmp -b $filter|") ||die;
+while (<FL>)
+{
+    @a=split("\t");
+    for ($i=3;$i<scalar(@a);$i++)
+    {
+	print "\t" unless $i==3;
+	print $a[$i];
+    }
+
+}

apply_bed_filter/apply_bed_filter.sh

@@ -0,0 +1,71 @@
+#!/bin/bash
+
+set -e -x
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+exit
+#!/usr/bin/perl
+#
+# Applies BED filter to a VCF file
+$filter=$ARGV[0];
+die "Filter file is not specified" unless -e $filter;
+die "Incorrect filter format" unless $filter =~/\.(bed|BED)/;
+
+
+#hash based coarse filter
+open(FL,"$filter")||die;
+while (<FL>)
+{
+    @a=split("\t");
+    $a[1]=int($a[1]/10000);
+    $a[2]=int($a[2]/10000);
+
+    for ($i=$a[1]-1;$i<=$a[2]+1;$i++)
+    {
+	$OK{"$a[0]-$i"}=1;
+    }
+}
+$tmp=int(rand()*1000000)."bed";
+$tmp="out.bed";
+
+open(OUT,">$tmp");
+
+while (<STDIN>)
+{
+    if (/^\#/)
+    {$header.=$_;next;}
+    @a=split("\t");
+
+    $red=int($a[1]/10000);
+
+    next unless $OK{"$a[0]-$red"};
+    print OUT "$a[0]\t$a[1]\t$a[1]\t$_";
+}
+
+close OUT;
+
+print $header;
+open(FL,"bedtools intersect -u -a $tmp -b $filter|") ||die;
+while (<FL>)
+{
+    @a=split("\t");
+    for ($i=3;$i<scalar(@a);$i++)
+    {
+	print "\t" unless $i==3;
+	print $a[$i];
+    }
+
+}

download_reference.sh

@@ -21,5 +21,6 @@ if [ ! -e /tmp/gcat_set_025_2.fastq.gz ]
 then
 wget	-O /tmp/gcat_set_025_1.fastq.gz http://gapp-west.s3.amazonaws.com/sample_exome/025_Bioplanet_GCAT_30x/gcat_set_025_1.fastq.gz
 wget	-O /tmp/gcat_set_025_2.fastq.gz http://gapp-west.s3.amazonaws.com/sample_exome/025_Bioplanet_GCAT_30x/gcat_set_025_2.fastq.gz
+wget	-O /tmp/gcat_set_025.bed.gz     http://gapp-west.s3.amazonaws.com/sample_exome/025_Bioplanet_GCAT_30x/gcat_set_025.bed.gz
 
 fi

gcat_set_025.bed.gz

@@ -127,9 +127,7 @@ do
 	-recalFile indels.recal.tmp \
 	-tranchesFile indels.tranches \
 	-rscriptFile indels.plots.R \
-	-nt $cpu_cores 
-
-	mv indels.recal.tmp indels.recal
+	-nt $cpu_cores && mv indels.recal.tmp indels.recal
     fi
 
 done
@@ -154,7 +152,7 @@ java -Xmx7g -jar $gatk_jar \
     -T CombineVariants  \
     -R $gatk_data/hg19/ucsc.hg19.fasta \
     --variant snps.recalibrated.filtered.vcf \
-    --variant indels.recalibrated.filtered.vcf \ 
+    --variant indels.recalibrated.filtered.vcf \
     -o recalibrated.filtered.vcf 
 
 bgzip recalibrated.filtered.vcf