From e6be1c74e88acd3605a09f96630ac0093016b9b1 Mon Sep 17 00:00:00 2001
From: Brandi Cantarel <brandi.cantarel@utsouthwestern.edu>
Date: Tue, 18 Feb 2020 11:55:57 -0600
Subject: [PATCH] updates for tmp directory, indels, svs and remove hard links
to files
---
alignment/bamqc.sh | 6 +++---
alignment/markdups.sh | 8 +++++---
variants/gatkrunner.sh | 6 +++---
variants/parse_pindel.pl | 5 ++++-
variants/svcalling.sh | 17 ++++++++++++-----
variants/unionvcf.pl | 2 +-
6 files changed, 28 insertions(+), 16 deletions(-)
diff --git a/alignment/bamqc.sh b/alignment/bamqc.sh
index a464836..cee7b80 100644
--- a/alignment/bamqc.sh
+++ b/alignment/bamqc.sh
@@ -59,12 +59,12 @@ if [[ $nuctype == 'dna' ]]; then
samtools view -@ $NPROC -b -L ${bed} -o ${pair_id}.ontarget.bam ${sbam}
samtools index -@ $NPROC ${pair_id}.ontarget.bam
samtools flagstat ${pair_id}.ontarget.bam > ${pair_id}.ontarget.flagstat.txt
- java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectAlignmentSummaryMetrics R=${index_path}/genome.fa I=${pair_id}.ontarget.bam OUTPUT=${pair_id}.alignmentsummarymetrics.txt --TMP_DIR=${tmpdir}
- java -Xmx64g -Djava.io.tmpdir=${tmpdir} -XX:ParallelGCThreads=$NPROC -jar $PICARD/picard.jar EstimateLibraryComplexity I=${sbam} OUTPUT=${pair_id}.libcomplex.txt --TMP_DIR=${tmpdir}
+ java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectAlignmentSummaryMetrics R=${index_path}/genome.fa I=${pair_id}.ontarget.bam OUTPUT=${pair_id}.alignmentsummarymetrics.txt TMP_DIR=${tmpdir}
+ java -Xmx64g -Djava.io.tmpdir=${tmpdir} -XX:ParallelGCThreads=$NPROC -jar $PICARD/picard.jar EstimateLibraryComplexity I=${sbam} OUTPUT=${pair_id}.libcomplex.txt TMP_DIR=${tmpdir}
samtools view -@ $NPROC -b -q 1 ${sbam} | bedtools coverage -sorted -hist -g ${index_path}/genomefile.txt -b stdin -a ${bed} > ${pair_id}.mapqualcov.txt
samtools view -@ $NPROC ${sbam} | awk '{sum+=$5} END { print "Mean MAPQ =",sum/NR}' > ${pair_id}.meanmap.txt
fi
- java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectInsertSizeMetrics INPUT=${sbam} HISTOGRAM_FILE=${pair_id}.hist.ps REFERENCE_SEQUENCE=${index_path}/genome.fa OUTPUT=${pair_id}.hist.txt --TMP_DIR=${tmpdir}
+ java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectInsertSizeMetrics INPUT=${sbam} HISTOGRAM_FILE=${pair_id}.hist.ps REFERENCE_SEQUENCE=${index_path}/genome.fa OUTPUT=${pair_id}.hist.txt TMP_DIR=${tmpdir}
bedtools coverage -sorted -g ${index_path}/genomefile.txt -a ${bed} -b ${sbam} -hist > ${pair_id}.covhist.txt
grep ^all ${pair_id}.covhist.txt > ${pair_id}.genomecov.txt
perl $baseDir/calculate_depthcov.pl ${pair_id}.covhist.txt
diff --git a/alignment/markdups.sh b/alignment/markdups.sh
index 306343b..dbe1b8e 100644
--- a/alignment/markdups.sh
+++ b/alignment/markdups.sh
@@ -71,9 +71,9 @@ elif [ $algo == 'fgbio_umi' ]
then
module load fgbio bwakit/0.7.15 bwa/intel/0.7.17 samtools/gcc/1.8
samtools index -@ $NPROC ${sbam}
- fgbio --tmp-dir ./ GroupReadsByUmi -s identity -i ${sbam} -o ${pair_id}.group.bam --family-size-histogram ${pair_id}.umihist.txt -e 0 -m 0
- fgbio --tmp-dir ./ CallMolecularConsensusReads -i ${pair_id}.group.bam -p consensus -M 1 -o ${pair_id}.consensus.bam -S ':none:'
- samtools index ${pair_id}.consensus.bam
+ fgbio --tmp-dir ./ GroupReadsByUmi -s identity -i ${sbam} -o group.bam --family-size-histogram ${pair_id}.umihist.txt -e 0 -m 0
+ fgbio --tmp-dir ./ CallMolecularConsensusReads -i group.bam -p consensus -M 1 -o ${pair_id}.consensus.bam -S ':none:'
+ samtools index -@ $NPROC ${pair_id}.consensus.bam
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
@@ -85,6 +85,8 @@ then
samtools view -1 out.sam > ${pair_id}.consensus.bam
fi
samtools sort --threads $NPROC -o ${pair_id}.dedup.bam ${pair_id}.consensus.bam
+ samtools sort --threads $NPROC -o ${pair_id}.group.bam group.bam
+ samtools index -@ $NPROC ${pair_id}.group.bam
else
cp ${sbam} ${pair_id}.dedup.bam
fi
diff --git a/variants/gatkrunner.sh b/variants/gatkrunner.sh
index 2bb2273..df24cc7 100755
--- a/variants/gatkrunner.sh
+++ b/variants/gatkrunner.sh
@@ -60,7 +60,7 @@ fi
source /etc/profile.d/modules.sh
module load gatk/4.1.4.0 samtools/gcc/1.8
which samtools
-/cm/shared/apps/samtools/gcc/1.8/bin/samtools index -@ $NPROC ${sbam}
+samtools index -@ $NPROC ${sbam}
if [[ $algo == 'gatkbam_rna' ]]
then
@@ -72,12 +72,12 @@ then
gatk SplitNCigarReads -R ${reffa} -I ${pair_id}.rg_added_sorted.bam -O ${pair_id}.split.bam
gatk --java-options "-Xmx32g" BaseRecalibrator -I ${pair_id}.split.bam --known-sites ${index_path}/dbSnp.gatk4.vcf.gz -R ${reffa} -O ${pair_id}.recal_data.table --use-original-qualities
gatk --java-options "-Xmx32g" ApplyBQSR -I ${pair_id}.split.bam -R ${reffa} -O ${pair_id}.final.bam --use-original-qualities -bqsr ${pair_id}.recal_data.table
- /cm/shared/apps/samtools/gcc/1.8/bin/samtools index -@ $NPROC ${pair_id}.final.bam
+ samtools index -@ $NPROC ${pair_id}.final.bam
elif [[ $algo == 'gatkbam' ]]
then
gatk --java-options "-Xmx32g" BaseRecalibrator -I ${sbam} --known-sites ${index_path}/dbSnp.gatk4.vcf.gz -R ${reffa} -O ${pair_id}.recal_data.table --use-original-qualities
gatk --java-options "-Xmx32g" ApplyBQSR -I ${sbam} -R ${reffa} -O ${pair_id}.final.bam --use-original-qualities -bqsr ${pair_id}.recal_data.table
- /cm/shared/apps/samtools/gcc/1.8/bin/samtools index -@ $NPROC ${pair_id}.final.bam
+ samtools index -@ $NPROC ${pair_id}.final.bam
elif [[ $algo == 'abra2' ]]
then
diff --git a/variants/parse_pindel.pl b/variants/parse_pindel.pl
index 3f00316..ecc47d2 100755
--- a/variants/parse_pindel.pl
+++ b/variants/parse_pindel.pl
@@ -79,12 +79,15 @@ while (my $line = <VCF>) {
if ($hash{SVTYPE} eq 'DUP:TANDEM') {
print DUP join("\t",$chrom,$pos,$id,$ref,$alt,$score,$filter,$annot,$newformat,@newgts),"\n";
}elsif ($hash{SVTYPE} eq 'DEL' || $hash{SVTYPE} eq 'INS') {
- if (abs($hash{SVLEN}) < 20) {
+ if (abs($hash{SVLEN}) < 100) {
print SI join("\t",$chrom,$pos,$id,$ref,$alt,$score,$filter,$annot,$newformat,@newgts),"\n";
}else {
$newalt = "<".$hash{SVTYPE}.">";
print SV join("\t",$chrom,$pos,$id,'N',$newalt,$score,$filter,$annot,$newformat,@newgts),"\n";
}
+ }else {
+ $newalt = "<".$hash{SVTYPE}.">";
+ print SV join("\t",$chrom,$pos,$id,'N',$newalt,$score,$filter,$annot,$newformat,@newgts),"\n";
}
}
close VCF;
diff --git a/variants/svcalling.sh b/variants/svcalling.sh
index 75a062d..17830e7 100755
--- a/variants/svcalling.sh
+++ b/variants/svcalling.sh
@@ -59,10 +59,10 @@ fi
source /etc/profile.d/modules.sh
module load htslib/gcc/1.8 samtools/gcc/1.8 bcftools/gcc/1.8 bedtools/2.26.0 snpeff/4.3q vcftools/0.1.14
export PATH=/project/shared/bicf_workflow_ref/seqprg/bin:$PATH
+mkdir temp
if [[ $method == 'delly' ]]
then
- mkdir temp
module load delly2/v0.7.7-multi
if [[ -n ${normal} ]]
then
@@ -105,10 +105,17 @@ then
else
svaba run -p $NPROC -G ${reffa} -t ${sbam} -a ${pair_id}
fi
- vcf-concat ${pair_id}.svaba.unfiltered*sv.vcf | vcf-sort -t temp > svaba.unfiltered.vcf
- bash $baseDir/norm_annot.sh -r ${index_path} -p svaba -v svaba.unfiltered.vcf -s
- java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config ${snpeffgeno} svaba.norm.vcf | bgzip > svaba.vcf.gz
-
+ vcf-concat ${pair_id}.svaba.unfiltered*sv.vcf | vcf-sort -t temp > svaba.sv.vcf
+ cat svaba.sv.vcf | java -jar $SNPEFF_HOME/SnpSift.jar filter "( GEN[*].AD[0] >= 20 )" | bgzip > svaba.vcf.gz
+ tabix svaba.sv.vcf.gz
+ bash $baseDir/norm_annot.sh -r ${index_path} -p svaba.sv -v svaba.sv.vcf.gz -s
+ java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config ${snpeffgeno} svaba.sv.norm.vcf | bgzip > ${pair_id}.svaba.sv.vcf.gz
+ vcf-concat ${pair_id}.svaba.unfiltered*indel.vcf | vcf-sort -t temp > svaba.indel.vcf
+ cat svaba.indel.vcf | java -jar $SNPEFF_HOME/SnpSift.jar filter "( GEN[*].AD[0] >= 20 )" | bgzip > svaba.indel.vcf.gz
+ tabix svaba.indel.vcf.gz
+ bash $baseDir/norm_annot.sh -r ${index_path} -p svaba.indel -v svaba.indel.vcf.gz -s
+ java -Xmx10g -jar $SNPEFF_HOME/snpEff.jar -no-intergenic -lof -c $SNPEFF_HOME/snpEff.config ${snpeffgeno} svaba.indel.norm.vcf | bgzip > ${pair_id}.svaba.vcf.gz
+
fi
if [[ $method == 'lumpy' ]]
diff --git a/variants/unionvcf.pl b/variants/unionvcf.pl
index c63fbae..84eedd5 100755
--- a/variants/unionvcf.pl
+++ b/variants/unionvcf.pl
@@ -3,7 +3,7 @@
my $headerfile = shift @ARGV;
my @vcffiles = @ARGV;
-my @callers = ('fb','mutect','gatk','platypus','strelka2','shimmer','virmid');
+my @callers = ('fb','mutect','gatk','platypus','strelka2','shimmer','pindel','svaba');
%algos = map {$_=>1} @callers;
open HEADER, "<$headerfile" or die $!;
--
GitLab