diff --git a/variants/somatic_callers.sh b/variants/somatic_callers.sh new file mode 100644 index 0000000000000000000000000000000000000000..99c6870f833c2efbe0ea46bdbe0ca02bef59decf --- /dev/null +++ b/variants/somatic_callers.sh @@ -0,0 +1,98 @@ +#!/bin/bash +#run_somatic_caller.sh + + +usage(){ + echo "-h --Help documentation for run_somatic_caller.sh" + echo "-a --Somatic Workflow Method: strelka2, virmid, speedseq, mutect, varscan, shimmer" + echo "-n --Normal" + echo "-t --Tumor" + echo "Example: bash somatic_callers.sh -a strelka2 -n ORD1_N_panel1385 -t ORD1_T_panel1385" + exit 1 +} + +OPTIND=1 # Reset OPTIND + +while getopts :n:t:a:h opt +do + case $opt in + n) normal=$OPTARG;; + t) tumor=$OPTARG;; + a) algo=$OPTARG;; + h) usage;; + esac +done + +shift $(($OPTIND -1)) + + +#Check for mandatory options +if [[ -z $normal ]] || [[ -z $tumor ]] || [[ -z $algo ]]; then + usage +fi + +if [[ -z $SLURM_CPUS_ON_NODE ]] + then + SLURM_CPUS_ON_NODE=1 +fi + +index_path=/project/shared/bicf_workflow_ref/GRCh38 + +genome_reference=${index_path}/genome.fa +cosmic_reference=${index_path}/cosmic.vcf.gz +dbSnp_reference=${index_path}/dbSnp.vcf.gz +target_bed=${index_path}/UTSWV2.bed + +if [ $algo == 'strelka2' ] + then + module load strelka/2.8.3 manta/1.2.0 + manta_analysisPath=MantaAnalysisPath + strelka_analysisPath=StrelkaAnalysisPath + mkdir ${manta_analysisPath} + mkdir ${strelka_analysisPath} + /cm/shared/apps/manta/1.2.0/bin/configManta.py --normalBam ${normal}.bam --tumorBam ${tumor}.bam --referenceFasta ${genome_reference} --runDir ${manta_analysisPath} + ${manta_analysisPath}/runWorkflow.py -m local -j 8 + /cm/shared/apps/strelka/2.8.3/bin/configureStrelkaSomaticWorkflow.py --normalBam ${normal}.bam --tumorBam ${tumor}.bam --referenceFasta ${genome_reference} --targeted --indelCandidates ${manta_analysisPath}/results/variants/candidateSmallIndels.vcf.gz --runDir ${strelka_analysisPath} + ${strelka_analysisPath}/runWorkflow.py -m local -j 8 +fi + +if [ $algo == 'virmid' ] + then + module load python/2.7.x-anaconda bedtools/2.25.0 snpeff/4.2 virmid/1.2 vcftools/0.1.14 + virmid -R ${genome_reference} -D ${tumor}.bam -N ${normal}.bam -s ${cosmic_reference} -t $SLURM_CPUS_ON_NODE -M 2000 -c1 10 -c2 10 + perl /project/PHG/PHG_Clinical/clinseq_workflows/scripts/addgt_virmid.pl ${tumor}.final.bam.virmid.som.passed.vcf + perl /project/PHG/PHG_Clinical/clinseq_workflows/scripts/addgt_virmid.pl ${tumor}.final.bam.virmid.loh.passed.vcf + vcf-concat *gt.vcf | vcf-sort | vcf-annotate -n --fill-type -n | java -jar $SNPEFF_HOME/SnpSift.jar filter '((NDP >= 10) & (DDP >= 10))' | perl -pe 's/TUMOR/'${tumor}'/' | perl -pe 's/NORMAL/'${normal}'/g' |bedtools intersect -header -a stdin -b ${target_bed} |bgzip > ${tumor}_${normal}.virmid.vcf.gz +fi + +if [ $algo == 'speedseq' ] + then + module load python/2.7.x-anaconda bedtools/2.25.0 snpeff/4.2 speedseq/20160506 bcftools/intel/1.3 vcftools/0.1.14 + speedseq somatic -q 10 -w ${target_bed} -t $SLURM_CPUS_ON_NODE -o ${tumor}.sssom ${genome_reference} ${normal}.final.bam ${tumor}.final.bam + vcf-annotate -H -n --fill-type ${tumor}.sssom.vcf.gz | java -jar $SNPEFF_HOME/SnpSift.jar filter --pass '((QUAL >= 10) & (GEN[*].DP >= 10))' | perl -pe 's/TUMOR/'${tumor}'/' | perl -pe 's/NORMAL/'${normal}'/g' |bgzip > ${tumor}_${normal}.sspanel.vcf.gz +fi + +if [ $algo == 'mutect' ] + then + module load parallel python/2.7.x-anaconda gatk/3.5 bcftools/intel/1.3 bedtools/2.25.0 snpeff/4.2 vcftools/0.1.14 + cut -f 1 /project/shared/bicf_workflow_ref/GRCh38/genomefile.5M.txt | parallel --delay 2 -j 10 "java -Xmx20g -jar $GATK_JAR -R ${genome_reference} -D ${dbSnp_reference} -T MuTect2 -stand_call_conf 30 -stand_emit_conf 10.0 -A FisherStrand -A QualByDepth -A VariantType -A DepthPerAlleleBySample -A HaplotypeScore -A AlleleBalance -I:tumor ${tumor}.final.bam -I:normal ${normal}.final.bam --cosmic ${cosmic} -o ${tumor}.{}.mutect.vcf -L {}" + vcf-concat ${tumor}*.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/'${tumor}'/' | perl -pe 's/NORMAL/'${normal}'/g' |bgzip > ${tumor}_${normal}.mutect.vcf.gz +fi + +if [ $algo == 'varscan' ] + then + module load python/2.7.x-anaconda bedtools/2.25.0 snpeff/4.2 bcftools/intel/1.3 samtools/intel/1.3 VarScan/2.4.2 speedseq/20160506 vcftools/0.1.14 + sambamba mpileup -L ${target_bed} -t $SLURM_CPUS_ON_NODE ${tumor}.final.bam --samtools "-C 50 -f ${genome_reference}" > t.mpileup + sambamba mpileup -L ${target_bed} -t $SLURM_CPUS_ON_NODE ${normal} --samtools "-C 50 -f ${genome_reference}" > n.mpileup + VarScan somatic n.mpileup t.mpileup ${tumor}.vscan --output-vcf 1 + VarScan copynumber n.mpileup t.mpileup ${tumor}.vscancnv + vcf-concat ${tumor}.vscan*.vcf | vcf-sort | vcf-annotate -n --fill-type -n | java -jar $SNPEFF_HOME/SnpSift.jar filter '((exists SOMATIC) & (GEN[*].DP >= 10))' | perl -pe 's/TUMOR/'${tumor}'/' | perl -pe 's/NORMAL/'${normal}'/g' |bedtools intersect -header -a stdin -b ${target_bed} |bgzip > ${tumor}_${normal}.varscan.vcf.gz +fi + +if [ $algo == 'shimmer' ] + then + module load python/2.7.x-anaconda bedtools/2.25.0 snpeff/4.2 bcftools/intel/1.3 shimmer/0.1.1 vcftools/0.1.14 + shimmer.pl --minqual 25 --ref ${genome_reference} ${normal}.final.bam ${tumor}.final.bam --outdir shimmer 2> shimmer.err + perl /project/PHG/PHG_Clinical/clinseq_workflows/scripts/add_readct_shimmer.pl + vcf-annotate -n --fill-type shimmer/somatic_diffs.readct.vcf | java -jar $SNPEFF_HOME/SnpSift.jar filter '(GEN[*].DP >= 10)' | perl -pe 's/TUMOR/'${tumor}'/' | perl -pe 's/NORMAL/'${normal}'/g' | bedtools intersect -header -a stdin -b ${target_bed} | bgzip > ${tumor}_${normal}.shimmer.vcf.gz +fi