diff --git a/alignment/markdups.sh b/alignment/markdups.sh
index 9042ac32f3a34f6dc8866caba1afa4f090cbe720..55fc2acc6aa03f853a53b16dd8c37fdc16275f5e 100644
--- a/alignment/markdups.sh
+++ b/alignment/markdups.sh
@@ -49,10 +49,10 @@ then
samtools markdup -s --output-fmt BAM -@ $SLURM_CPUS_ON_NODE sort.bam ${pair_id}.dedup.bam
elif [ $algo == 'picard' ]
then
- java -Djava.io.tmpdir=./ -Xmx16g -jar $PICARD/picard.jar MarkDuplicates I=${sbam} O=${pair_id}.dedup.bam M=${pair_id}.dedup.stat.txt
+ java -XX:ParallelGCThreads=$SLURM_CPUS_ON_NODE -Djava.io.tmpdir=./ -Xmx16g -jar $PICARD/picard.jar MarkDuplicates I=${sbam} O=${pair_id}.dedup.bam M=${pair_id}.dedup.stat.txt
elif [ $algo == 'picard_umi' ]
then
- java -Djava.io.tmpdir=./ -Xmx16g -jar $PICARD/picard.jar MarkDuplicates BARCODE_TAG=RX I=${sbam} O=${pair_id}.dedup.bam M=${pair_id}.dedup.stat.txt
+ java -XX:ParallelGCThreads=$SLURM_CPUS_ON_NODE -Djava.io.tmpdir=./ -Xmx16g -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 bwa/intel/0.7.15
@@ -63,8 +63,8 @@ then
samtools fastq -1 ${pair_id}.consensus.R1.fastq -2 ${pair_id}.consensus.R2.fastq ${pair_id}.consensus.bam
gzip ${pair_id}.consensus.R1.fastq
gzip ${pair_id}.consensus.R2.fastq
- bwa mem -M -C -t 2 -R "@RG\tID:${pair_id}\tLB:tx\tPL:illumina\tPU:barcode\tSM:${pair_id}" /project/shared/bicf_workflow_ref/human/GRCh38/genome.fa ${pair_id}.consensus.R1.fastq.gz ${pair_id}.consensus.R2.fastq.gz | samtools view -1 - > ${pair_id}.consensus.bam
- samtools sort --threads 10 -o ${pair_id}.dedup.bam ${pair_id}.consensus.bam
+ bwa mem -M -C -t $SLURM_CPUS_ON_NODE -R "@RG\tID:${pair_id}\tLB:tx\tPL:illumina\tPU:barcode\tSM:${pair_id}" /project/shared/bicf_workflow_ref/human/GRCh38/genome.fa ${pair_id}.consensus.R1.fastq.gz ${pair_id}.consensus.R2.fastq.gz | samtools view -1 - > ${pair_id}.consensus.bam
+ samtools sort --threads $SLURM_CPUS_ON_NODE -o ${pair_id}.dedup.bam ${pair_id}.consensus.bam
else
cp ${sbam} ${pair_id}.dedup.bam
fi
diff --git a/variants/annotvcf.sh b/variants/annotvcf.sh
index 22a44883f21badeb071b4b6afa2e6d3a55a94108..1d841ced55a69dd7922bc16ff9943d0f083250bf 100755
--- a/variants/annotvcf.sh
+++ b/variants/annotvcf.sh
@@ -34,8 +34,10 @@ if [[ $index_path == '/project/shared/bicf_workflow_ref/human/GRCh38' ]]
then
tabix -f ${unionvcf}
bcftools annotate -Oz -a ${index_path}/gnomad.txt.gz -h ${index_path}/gnomad.header -c CHROM,POS,REF,ALT,GNOMAD_HOM,GNOMAD_AF,AF_POPMAX,GNOMAD_HG19_VARIANT -o ${pair_id}.gnomad.vcf.gz ${unionvcf}
- tabix ${pair_id}.gnomad.vcf.gz
- bcftools annotate -Oz -a ${index_path}/repeat_regions.bed.gz -o ${pair_id}.repeat.vcf.gz --columns CHROM,FROM,TO,RepeatType -h /project/shared/bicf_workflow_ref/RepeatType.header ${pair_id}.gnomad.vcf.gz
+ tabix ${pair_id}.gnomad.vcf.gz
+ bcftools annotate -Oz -a ${index_path}/oncokb_hotspot.txt.gz -o ${pair_id}.oncohotspot.vcf.gz -h ${index_path}/oncokb_hotspot.header -c CHROM,FROM,TO,OncoKB_REF,OncoKB_ALT,Gene,OncoKB_ProteinChange,OncoKB_AF,OncoTree_Tissue,OncoTree_MainType,OncoTree_Code,OncoKBHotspot ${pair_id}.gnomad.vcf.gz
+ tabix ${pair_id}.oncohotspot.vcf.gz
+ bcftools annotate -Oz -a ${index_path}/repeat_regions.bed.gz -o ${pair_id}.repeat.vcf.gz --columns CHROM,FROM,TO,RepeatType -h /project/shared/bicf_workflow_ref/RepeatType.header ${pair_id}.oncohotspot.vcf.gz
java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-downstream -no-upstream -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config GRCh38.86 ${pair_id}.repeat.vcf.gz | java -jar $SNPEFF_HOME/SnpSift.jar annotate -id ${index_path}/dbSnp.vcf.gz - | java -jar $SNPEFF_HOME/SnpSift.jar annotate -info CLNSIG,CLNDSDB,CLNDSDBID,CLNDBN,CLNREVSTAT,CLNACC ${index_path}/clinvar.vcf.gz - | java -jar $SNPEFF_HOME/SnpSift.jar annotate -info CNT ${index_path}/cosmic.vcf.gz - | java -Xmx10g -jar $SNPEFF_HOME/SnpSift.jar dbnsfp -v -db ${index_path}/dbNSFP.txt.gz - | bgzip > ${pair_id}.annot.vcf.gz
tabix ${pair_id}.annot.vcf.gz
else
diff --git a/variants/germline_vc.sh b/variants/germline_vc.sh
index 425746bf167c50af151cf47ca95996c479369355..73e0740050cd97bdab0d2be44c77be0dcd604393 100755
--- a/variants/germline_vc.sh
+++ b/variants/germline_vc.sh
@@ -76,7 +76,8 @@ elif [[ $algo == 'hotspot' ]]
then
samtools mpileup -d 99999 -t 'AD,DP,INFO/AD' -uf ${reffa} *.bam > ${pair_id}.mpi
bcftools filter -i "AD[1]/DP > 0.01" ${pair_id}.mpi | bcftools filter -i "DP > 50" | bcftools call -m -A |vcf-annotate -n --fill-type | bcftools norm -c s -f ${reffa} -w 10 -O z -o ${pair_id}.lowfreq.vcf.gz -
- java -jar $SNPEFF_HOME/SnpSift.jar annotate ${index_path}/cosmic.vcf.gz ${pair_id}.lowfreq.vcf.gz | java -jar $SNPEFF_HOME/SnpSift.jar filter "(CNT[*] >0)" - |bgzip > ${pair_id}.hotspot.vcf.gz
+ tabix ${pair_id}.lowfreq.vcf.gz
+ bcftools annotate -Ov -a ${index_path}/oncokb_hotspot.txt.gz -h ${index_path}/oncokb_hotspot.header -c CHROM,FROM,TO,OncoKB_REF,OncoKB_ALT,Gene,OncoKB_ProteinChange,OncoKB_AF,OncoTree_Tissue,OncoTree_MainType,OncoTree_Code,OncoKBHotspot ${pair_id}.lowfreq.vcf.gz | java -jar $SNPEFF_HOME/SnpSift.jar annotate ${index_path}/cosmic.vcf.gz - | grep '#\|CNT\|OncoKBHotspot' | bgzip > ${pair_id}.hotspot.vcf.gz
elif [[ $algo == 'speedseq' ]]
then
module load speedseq/gcc/0.1.2
diff --git a/variants/somatic_vc.sh b/variants/somatic_vc.sh
index f96b9a4f5890da9f9d54642b5f60a93a481ce0cf..4a5ef6e9b1515748381ddceab7c6f54dfcff9441 100755
--- a/variants/somatic_vc.sh
+++ b/variants/somatic_vc.sh
@@ -120,7 +120,7 @@ then
else
awk '{print $1":"$2"-"$3}' ${tbed} | parallel --delay 2 -j 10 "java -Xmx20g -jar \$GATK_JAR -R ${reffa} -D ${dbsnp} -T MuTect2 -stand_call_conf 10 -A FisherStrand -A QualByDepth -A VariantType -A DepthPerAlleleBySample -A HaplotypeScore -A AlleleBalance -I:tumor ${tumor} -I:normal ${normal} --cosmic ${cosmic} -o ${tid}.{}.mutect.vcf -L {}"
fi
- vcf-concat ${tid}*mutect.vcf | vcf-sort | vcf-annotate -n --fill-type | java -jar $SNPEFF_HOME/SnpSift.jar filter -p '((FS <= 60) & GEN[*].DP >= 10)' | perl -pe "s/TUMOR/${tid}/" | perl -pe "s/NORMAL/${nid}/g" |bgzip > ${pair_id}.mutect.vcf.gz
+ vcf-concat ${tid}*mutect.vcf | vcf-sort | vcf-annotate -n --fill-type | java -jar $SNPEFF_HOME/SnpSift.jar filter -p '(GEN[*].DP >= 10)' | perl -pe "s/TUMOR/${tid}/" | perl -pe "s/NORMAL/${nid}/g" |bgzip > ${pair_id}.mutect.vcf.gz
fi
if [ $algo == 'varscan' ]
diff --git a/variants/unionvcf.pl b/variants/unionvcf.pl
index 57c2046753dad942dd4eb6186ee79e5e94b1ac84..1b5ff5d202c6bd49568e6d47aa8e05b4f0d6e6c5 100755
--- a/variants/unionvcf.pl
+++ b/variants/unionvcf.pl
@@ -48,6 +48,7 @@ foreach $vcf (@vcffiles) {
my $newformat = 'GT:DP:AD:AO:RO';
my %newgts;
my %afinfo;
+ my %gtfilt;
my $missingGT = 0;
FG:foreach my $i (0..$#gts) {
my $allele_info = $gts[$i];
@@ -67,6 +68,9 @@ foreach $vcf (@vcffiles) {
$missingGT ++;
next FG;
}
+ if ($gtdata{FT} && $gtdata{FT} =~ m/HighSNVSB/) {
+ $gtfilt{'StrandBias'} = 1;
+ }
if ($gtdata{DP4}) { #varscan uses this
my ($ref_fwd,$ref_rev,$alt_fwd,$alt_rev) = split(',',$gtdata{DP4});
$gtdata{AO} = $alt_fwd+$alt_rev;
@@ -118,6 +122,17 @@ foreach $vcf (@vcffiles) {
push @gtdesc, join(":",$id,$afinfo{$id});
push @newgts, $newgts{$id};
}
+ my @filts = split(";",$filter);
+ my %filterqc = map {$_ => 1} @filts;
+ delete $filterqc{'.'};
+ if ($gtfilt{'StrandBias'} || ($hash{FS} && $hash{FS} > 60) ||
+ ($hash{SAP} && $hash{SAP} > 20)) {
+ $filterqc{strandBias} = 1;
+ }if (scalar(keys %filterqc) > 1) {
+ $filter = join(";",keys %filterqc);
+ }else {
+ $filter = '.';
+ }
$lines{$chrom}{$pos}{$alt}{$caller} = [$chrom,$pos,$id,$ref,$alt,$score,$filter,$annot,$newformat,\@newgts,\@gtdesc];
}
close VCF;