update svcalling for somatic

alignment/markdups.sh

@@ -54,11 +54,10 @@ then
     java -Djava.io.tmpdir=./ -Xmx4g  -jar $PICARD/picard.jar MarkDuplicates BARCODE_TAG=RX I=${sbam} O=${pair_id}.dedup.bam M=${pair_id}.dedup.stat.txt
 elif [ $algo == 'fgbio_umi' ]   
 then
+    module load fgbio
     samtools index -@ $SLURM_CPUS_ON_NODE ${sbam}
-    source activate fgbiotools
     fgbio GroupReadsByUmi -s identity -i ${sbam} -o ${pair_id}.group.bam -m 0
     fgbio CallMolecularConsensusReads -i ${pair_id}.group.bam -p consensus -M 1 -o ${pair_id}.consensus.bam -S ':none:'
-    source deactivate
     module load bwa/intel/0.7.15
     samtools index ${pair_id}.consensus.bam
     samtools fastq -1 ${pair_id}.consensus.R1.fastq -2 ${pair_id}.consensus.R2.fastq ${pair_id}.consensus.bam

variants/somatic_vc.sh

@@ -12,12 +12,14 @@ usage(){
   echo "-y --TumorID"
   echo "-i --NormalBAM used for Mantra in the case of UMI consensus"
   echo "-j --TumorBAM used for Mantra in the case of UMI consensus"
+  echo "-b --TargetBed"
+  
   echo "Example: bash somatic_vc.sh -a strelka2 -y ORD1_N_panel1385 -y ORD1_T_panel138 -n ORD1_N_panel1385.final.bam -t ORD1_T_panel1385.final.bam"
   exit 1
 }
 
 OPTIND=1 # Reset OPTIND
-while getopts :n:t:r:x:y:i:j:a:h opt
+while getopts :n:t:r:x:y:i:j:a:b:h opt
 do
     case $opt in 
       r) index_path=$OPTARG;;
@@ -28,6 +30,7 @@ do
       i) mnormal=$OPTARG;;
       j) mtumor=$OPTARG;;
       a) algo=$OPTARG;;
+      b) tbed=$OPTARG;; 
       h) usage;;
     esac
 done
@@ -130,6 +133,6 @@ fi
 if [ $algo == 'lancet' ]
 then
     module load snpeff/4.3q lancet vcftools/0.1.14
-    lancet --tumor ${tumor} --normal ${normal} --ref $reffa --bed $target_panel --num-threads 16 > lancet.vcf
+    lancet --tumor ${tumor} --normal ${normal} --ref $reffa -B $tbed --num-threads 16 > lancet.vcf
     vcf-sort lancet.vcf | vcf-annotate -n --fill-type -n | perl -pe 's/TUMOR/${tid}/' | perl -pe 's/NORMAL/${nid}/g' |bgzip > ${pair_id}.lancet.vcf.gz
 fi

variants/svcalling.sh

@@ -86,12 +86,12 @@ fi
 #MERGE DELLY AND MAKE BED
 
 bcftools concat -a -O v delly_dup.bcf delly_inv.bcf delly_tra.bcf delly_del.bcf delly_ins.bcf| vcf-sort -t temp > delly.vcf
-perl $baseDir/vcf2bed.sv.pl delly.vcf > delly.bed
+perl $baseDir/vcf2bed.sv.pl delly.vcf | sort -V -k 1,1 -k 2,2n > delly.bed
 bgzip delly.vcf
 if [[ -n ${normal} ]]
 then
     tabix delly.vcf.gz
-    bcftools view -O z -o delly.vcf.gz -s ${keepid} ${pair_id}.delly.vcf.gz
+    bcftools view -O z -o ${pair_id}.delly.vcf.gz -s ${keepid} delly.vcf.gz
 else
     mv delly.vcf.gz ${pair_id}.delly.vcf.gz
 fi
@@ -110,19 +110,21 @@ then
     gawk '{ if ($0~"^@") { print; next } else { $10="*"; $11="*"; print } }' OFS="\t" splitters.sam | samtools  view -S -b - | samtools sort -o normal.splitters.bam -
     gawk '{ if ($0~"^@") { print; next } else { $10="*"; $11="*"; print } }' OFS="\t" discordants.sam | samtools  view -S  -b - | samtools sort -o normal.discordants.bam -
     speedseq sv -t $SLURM_CPUS_ON_NODE -o sssv -R ${reffa} -B ${normal},${sbam} -D normal.discordants.bam,discordants.bam -S normal.splitters.bam,splitters.bam -x ${index_path}/exclude_alt.bed
-    bcftools view -O z -o sssv.vcf.gz -s ${keepid} ${pair_id}.sssv.sv.vcf.gz 
+    java -jar $SNPEFF_HOME/SnpSift.jar filter -n "GEN[0].SU > 1" sssv.sv.vcf.gz |bgzip > tumor.sssv.sv.vcf.gz
+    tabix tumor.sssv.sv.vcf.gz
+    bcftools view -O z -o ${pair_id}.sssv.sv.vcf.gz -s ${keepid} tumor.sssv.sv.vcf.gz
 else
     speedseq sv -t $SLURM_CPUS_ON_NODE -o ${pair_id}.sssv -R ${reffa} -B ${sbam} -D discordants.bam -S splitters.bam -x ${index_path}/exclude_alt.bed   
 fi
 
-java -jar $SNPEFF_HOME/SnpSift.jar filter "GEN[0].SU > 2" ${pair_id}.sssv.sv.vcf.gz > lumpy.vcf
+java -jar $SNPEFF_HOME/SnpSift.jar filter "GEN[0].SU > 10" ${pair_id}.sssv.sv.vcf.gz > lumpy.vcf
 perl $baseDir/vcf2bed.sv.pl lumpy.vcf > lumpy.bed
 
 #COMPARE DELLY & LUMPY
 if [[ -n ${normal} ]]
 then
   bedtools multiinter -cluster -header -names novobreak delly lumpy -i novobreak.bed delly.bed lumpy.bed > sv.intersect.bed 
-  grep novobreak sv.intersect.bed |cut -f 1,2,3 |sort -V -k 1,1 -k 2,2n |grep -v start | bedtools intersect -header -b stdin -a ${tid}_${nid}.novobreak.vcf.gz  | perl -p -e 's/SPIKEIN/${tid}/' |bgzip > svt1.vcf.gz
+  grep novobreak sv.intersect.bed |cut -f 1,2,3 |sort -V -k 1,1 -k 2,2n |grep -v start | bedtools intersect -header -b stdin -a ${pair_id}.novobreak.vcf.gz  | perl -p -e 's/SPIKEIN/${tid}/' |bgzip > svt1.vcf.gz
 else
   bedtools multiinter -cluster -header -names delly lumpy -i delly.bed lumpy.bed > sv.intersect.bed 
 fi

variants/unionvcf.pl

@@ -15,7 +15,7 @@ my @sampleorder;
 
 my %headerlines;
 foreach $vcf (@vcffiles) {
-    $caller = (split(/\./,$vcf))[1];
+    $caller = (split(/\./,$vcf))[-3];
     open VCF, "gunzip -c $vcf|" or die $!;
     my @sampleids;
     while (my $line = <VCF>) {