Skip to content
Snippets Groups Projects
Commit 488f3a25 authored by Brandi Cantarel's avatar Brandi Cantarel
Browse files

sv filtering

parent 8d9c724f
Branches
Tags
No related merge requests found
...@@ -15,7 +15,7 @@ W1:while (my $line = <IN>) { ...@@ -15,7 +15,7 @@ W1:while (my $line = <IN>) {
chomp($line); chomp($line);
if ($line =~ m/^#/) { if ($line =~ m/^#/) {
if ($line =~ m/^#CHROM/) { if ($line =~ m/^#CHROM/) {
print OUT "##INFO=<ID=AF,Number=A,Type=Integer,Description=\"Alternate allele observation frequency\">\n"; print VCFOUT "##INFO=<ID=AF,Number=A,Type=Integer,Description=\"Alternate allele observation frequency\">\n";
my @header = split(/\t/,$line); my @header = split(/\t/,$line);
($chrom, $pos,$id,$ref,$alt,$score, ($chrom, $pos,$id,$ref,$alt,$score,
$filter,$info,$format,@gtheader) = split(/\t/, $line); $filter,$info,$format,@gtheader) = split(/\t/, $line);
...@@ -29,6 +29,7 @@ W1:while (my $line = <IN>) { ...@@ -29,6 +29,7 @@ W1:while (my $line = <IN>) {
} }
} }
print VCFOUT $line,"\n"; print VCFOUT $line,"\n";
next;
} }
my ($chrom, $pos,$id,$ref,$alt,$score, my ($chrom, $pos,$id,$ref,$alt,$score,
$filter,$annot,$format,@gts) = split(/\t/, $line); $filter,$annot,$format,@gts) = split(/\t/, $line);
......
...@@ -15,7 +15,7 @@ W1:while (my $line = <IN>) { ...@@ -15,7 +15,7 @@ W1:while (my $line = <IN>) {
chomp($line); chomp($line);
if ($line =~ m/^#/) { if ($line =~ m/^#/) {
if ($line =~ m/^#CHROM/) { if ($line =~ m/^#CHROM/) {
print OUT "##INFO=<ID=AF,Number=A,Type=Integer,Description=\"Alternate allele observation frequency\">\n"; print VCFOUT "##INFO=<ID=AF,Number=A,Type=Integer,Description=\"Alternate allele observation frequency\">\n";
my @header = split(/\t/,$line); my @header = split(/\t/,$line);
($chrom, $pos,$id,$ref,$alt,$score, ($chrom, $pos,$id,$ref,$alt,$score,
$filter,$info,$format,@gtheader) = split(/\t/, $line); $filter,$info,$format,@gtheader) = split(/\t/, $line);
...@@ -29,6 +29,7 @@ W1:while (my $line = <IN>) { ...@@ -29,6 +29,7 @@ W1:while (my $line = <IN>) {
} }
} }
print VCFOUT $line,"\n"; print VCFOUT $line,"\n";
next;
} }
my ($chrom, $pos,$id,$ref,$alt,$score, my ($chrom, $pos,$id,$ref,$alt,$score,
$filter,$annot,$format,@gts) = split(/\t/, $line); $filter,$annot,$format,@gts) = split(/\t/, $line);
......
...@@ -77,7 +77,6 @@ then ...@@ -77,7 +77,6 @@ then
delly2 call -t INS -o ${pair_id}.delly_insertion.bcf -q 30 -g ${reffa} ${sbam} ${normal} delly2 call -t INS -o ${pair_id}.delly_insertion.bcf -q 30 -g ${reffa} ${sbam} ${normal}
#delly2 filter -o ${pair_id}.delly_tra.bcf -f somatic -s samples.tsv ${pair_id}.delly_translocations.bcf #delly2 filter -o ${pair_id}.delly_tra.bcf -f somatic -s samples.tsv ${pair_id}.delly_translocations.bcf
else else
$tid=
#RUN DELLY #RUN DELLY
delly2 call -t BND -o ${pair_id}.delly_translocations.bcf -q 30 -g ${reffa} ${sbam} delly2 call -t BND -o ${pair_id}.delly_translocations.bcf -q 30 -g ${reffa} ${sbam}
delly2 call -t DUP -o ${pair_id}.delly_duplications.bcf -q 30 -g ${reffa} ${sbam} delly2 call -t DUP -o ${pair_id}.delly_duplications.bcf -q 30 -g ${reffa} ${sbam}
...@@ -87,19 +86,23 @@ then ...@@ -87,19 +86,23 @@ then
#delly2 filter -o ${pair_id}.delly_tra.bcf -f germline ${pair_id}.delly_translocations.bcf #delly2 filter -o ${pair_id}.delly_tra.bcf -f germline ${pair_id}.delly_translocations.bcf
fi fi
#MERGE DELLY AND MAKE BED #MERGE DELLY AND MAKE BED
bcftools concat -a -O v ${pair_id}.delly_duplications.bcf ${pair_id}.delly_inversions.bcf ${pair_id}.delly_translocations.bcf ${pair_id}.delly_deletion.bcf ${pair_id}.delly_insertion.bcf | vcf-sort -t temp | bgzip > ${pair_id}.delly.svar.vcf.gz bcftools concat -a -O v ${pair_id}.delly_duplications.bcf ${pair_id}.delly_inversions.bcf ${pair_id}.delly_translocations.bcf ${pair_id}.delly_deletion.bcf ${pair_id}.delly_insertion.bcf | vcf-sort -t temp | bgzip > ${pair_id}.delly.svar.vcf.gz
bash $baseDir/norm_annot.sh -r ${index_path} -p ${pair_id}.delly.sv -v ${pair_id}.delly.svar.vcf.gz -s bash $baseDir/norm_annot.sh -r ${index_path} -p ${pair_id}.delly.sv -v ${pair_id}.delly.svar.vcf.gz -s
java -jar $SNPEFF_HOME/SnpSift.jar filter "( GEN[*].DP >= 20 )" ${pair_id}.delly.sv.norm.vcf.gz | java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config ${snpeffgeno} - | bgzip > ${pair_id}.delly.vcf.gz java -jar $SNPEFF_HOME/SnpSift.jar filter "( GEN[*].DP >= 20 )" ${pair_id}.delly.sv.norm.vcf.gz | java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config ${snpeffgeno} - | bgzip > ${pair_id}.delly.vcf.gz
if [[ $filter == 1 ]] if [[ $filter == 1 ]]
then then
if [[ -z $tid ]]
then
tid=`samtools view -H ${sbam} |grep '^@RG' |perl -pi -e 's/\t/\n/g' |grep ID |cut -f 2 -d ':'`
fi
zgrep '#CHROM' ${pair_id}.delly.vcf.gz > ${pair_id}.delly.genefusion.txt zgrep '#CHROM' ${pair_id}.delly.vcf.gz > ${pair_id}.delly.genefusion.txt
zcat ${pair_id}.delly.vcf.gz | $SNPEFF_HOME/scripts/vcfEffOnePerLine.pl |java -jar $SNPEFF_HOME/SnpSift.jar extractFields - CHROM POS CHR2 END ANN[*].EFFECT ANN[*].GENE ANN[*].BIOTYPE FILTER FORMAT GEN[*] |grep -E 'gene_fusion|feature_fusion' | sort -u >> ${pair_id}.delly.genefusion.txt zcat ${pair_id}.delly.vcf.gz | $SNPEFF_HOME/scripts/vcfEffOnePerLine.pl |java -jar $SNPEFF_HOME/SnpSift.jar extractFields - CHROM POS CHR2 END ANN[*].EFFECT ANN[*].GENE ANN[*].BIOTYPE FILTER FORMAT GEN[*] |grep -E 'gene_fusion|feature_fusion' | sort -u >> ${pair_id}.dgf.txt
mv ${pair_id}.delly.vcf.gz ${pair_id}.delly.ori.vcf.gz mv ${pair_id}.delly.vcf.gz ${pair_id}.delly.ori.vcf.gz
bash $baseDir/filter_delly.pl -t $tid -p $pair_id -i ${pair_id}.delly.ori.vcf.gz perl $baseDir/filter_delly.pl -t $tid -p $pair_id -i ${pair_id}.delly.ori.vcf.gz
bgzip ${pair_id}.delly.vcf bgzip ${pair_id}.delly.vcf
cat ${pair_id}.potentialfusion.txt >> ${pair_id}.delly.genefusion.txt cat ${pair_id}.delly.potentialfusion.txt ${pair_id}.dgf.txt |sort -u > ${pair_id}.delly.genefusion.txt
fi fi
elif [[ $method == 'svaba' ]] elif [[ $method == 'svaba' ]]
then then
if [[ -n ${normal} ]] if [[ -n ${normal} ]]
...@@ -120,11 +123,11 @@ then ...@@ -120,11 +123,11 @@ then
if [[ $filter == 1 ]] if [[ $filter == 1 ]]
then then
zgrep '#CHROM' ${pair_id}.svaba.sv.vcf.gz > ${pair_id}.svaba.genefusion.txt zgrep '#CHROM' ${pair_id}.svaba.sv.vcf.gz > ${pair_id}.svaba.genefusion.txt
zcat ${pair_id}.svaba.sv.vcf.gz | $SNPEFF_HOME/scripts/vcfEffOnePerLine.pl |java -jar $SNPEFF_HOME/SnpSift.jar extractFields - CHROM POS ALT ID ANN[*].EFFECT ANN[*].GENE ANN[*].BIOTYPE FILTER FORMAT GEN[*] |grep -E 'gene_fusion|feature_fusion' | sort -u >> ${pair_id}.svaba.genefusion.txt zcat ${pair_id}.svaba.sv.vcf.gz | $SNPEFF_HOME/scripts/vcfEffOnePerLine.pl |java -jar $SNPEFF_HOME/SnpSift.jar extractFields - CHROM POS ALT ID ANN[*].EFFECT ANN[*].GENE ANN[*].BIOTYPE FILTER FORMAT GEN[*] |grep -E 'gene_fusion|feature_fusion' | sort -u >> ${pair_id}.sgf.txt
mv ${pair_id}.svaba.vcf.gz ${pair_id}.svaba.ori.vcf.gz mv ${pair_id}.svaba.vcf.gz ${pair_id}.svaba.ori.vcf.gz
bash $baseDir/filter_svaba.pl -t $tid -p $pair_id -i ${pair_id}.svaba.ori.vcf.gz perl $baseDir/filter_svaba.pl -t $tid -p $pair_id -i ${pair_id}.svaba.ori.vcf.gz -v
bgzip ${pair_id}.svaba.vcf bgzip ${pair_id}.svaba.vcf
cat ${pair_id}.potentialfusion.txt >> ${pair_id}.svaba.genefusion.txt cat ${pair_id}.svaba.potentialfusion.txt ${pair_id}.sgf.txt | sort -u >> ${pair_id}.svaba.genefusion.txt
fi fi
elif [[ $method == 'lumpy' ]] elif [[ $method == 'lumpy' ]]
then then
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment