Add depth filtering to bcftools consensus

Author: bbimber

SequenceAnalysis/resources/external/viral_consensus.sh

@@ -6,7 +6,8 @@ set -x
 BAM=$1
 FASTA=$2
 MASK_BED=$3
-BCFTOOLS=${4:-bcftools}
+MIN_DEPTH=$4
+BCFTOOLS=${5:-bcftools}
 
 VCF_CALLS=calls.vcf.gz
 THREADS=1
@@ -37,29 +38,37 @@ COUNT2=`$BCFTOOLS view -H $VCF_NORM | wc -l`
 echo 'Variants remaining: '$COUNT2
 echo -e 'VariantsAfterNorm\t'$COUNT2 >> $REPORT
 
+# filter on depth
+echo 'Filtering indel clusters.  Note: this is not currently used in the consensus.'
+VCF_DEPTH_FILTER=calls.norm.filter.vcf.gz
+$BCFTOOLS filter -e "DP < $MIN_DEPTH" -Oz -o $VCF_DEPTH_FILTER $VCF_NORM
+$BCFTOOLS index -t -f $VCF_DEPTH_FILTER
+COUNT3=`$BCFTOOLS view -H $VCF_DEPTH_FILTER | wc -l`
+echo 'Variants remaining: '$COUNT3
+echo -e 'VariantsAfterDpethFilter\t'$COUNT3 >> $REPORT
+
 # filter adjacent indels within 5bp
 echo 'Filtering indel clusters.  Note: this is not currently used in the consensus.'
 VCF_INDEL_FILTER=calls.norm.flt-indels.vcf.gz
-$BCFTOOLS filter --IndelGap 5 -Oz -o $VCF_INDEL_FILTER $VCF_NORM
+$BCFTOOLS filter --IndelGap 5 -Oz -o $VCF_INDEL_FILTER $VCF_DEPTH_FILTER
 $BCFTOOLS index -t -f $VCF_INDEL_FILTER
 COUNT3=`$BCFTOOLS view -H $VCF_INDEL_FILTER | wc -l`
 echo 'Variants that would remain: '$COUNT3
 echo -e 'VariantsAfterIndelFilter\t'$COUNT3 >> $REPORT
 
-#At the moment, do not use the filtered version:
+#At the moment, do not use the indel-filtered version:
 VCF_FOR_CONSENSUS=`basename $BAM .bam`".calls.vcf.gz"
-mv $VCF_NORM $VCF_FOR_CONSENSUS
-mv ${VCF_NORM}.tbi ${VCF_FOR_CONSENSUS}.tbi
+mv $VCF_DEPTH_FILTER $VCF_FOR_CONSENSUS
+mv ${VCF_DEPTH_FILTER}.tbi ${VCF_FOR_CONSENSUS}.tbi
 
 $BCFTOOLS consensus -f $FASTA -m $MASK_BED -o $OUT $VCF_FOR_CONSENSUS
 $BCFTOOLS consensus -f $FASTA -m $MASK_BED -o $OUT_IUPAC --iupac-codes $VCF_FOR_CONSENSUS
 
 rm $VCF_CALLS
 rm ${VCF_CALLS}.tbi
 
-#Previously moved:
-#rm $VCF_NORM
-#rm ${VCF_NORM}.tbi
+rm $VCF_NORM
+rm ${VCF_NORM}.tbi
 
 rm $VCF_INDEL_FILTER
 rm ${VCF_INDEL_FILTER}.tbi

SequenceAnalysis/src/org/labkey/sequenceanalysis/run/analysis/LofreqAnalysis.java

@@ -34,6 +34,7 @@
 import org.labkey.api.sequenceanalysis.run.AbstractCommandWrapper;
 import org.labkey.api.sequenceanalysis.run.SimpleScriptWrapper;
 import org.labkey.api.util.FileUtil;
+import org.labkey.api.writer.PrintWriters;
 import org.labkey.sequenceanalysis.SequenceAnalysisModule;
 import org.labkey.sequenceanalysis.run.util.DepthOfCoverageWrapper;
 import org.labkey.sequenceanalysis.run.variant.SNPEffStep;
@@ -42,6 +43,7 @@
 import java.io.File;
 import java.io.IOException;
 import java.io.InputStream;
+import java.io.PrintWriter;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Collections;
@@ -247,7 +249,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         getPipelineCtx().getLogger().info("Total positions with coverage below threshold (" + minCoverage + "): " + positionsSkipped);
         getPipelineCtx().getLogger().info("Total intervals of these gaps: " + gapIntervals);
 
-        consensusWrapper.execute(Arrays.asList("/bin/bash", script.getPath(), inputBam.getPath(), referenceGenome.getWorkingFastaFile().getPath(), mask.getPath()));
+        consensusWrapper.execute(Arrays.asList("/bin/bash", script.getPath(), inputBam.getPath(), referenceGenome.getWorkingFastaFile().getPath(), mask.getPath(), String.valueOf(minCoverage)));
         File calls = new File(inputBam.getParentFile(), FileUtil.getBaseName(inputBam) + ".calls.vcf.gz");
 
         Set<VariantContext> variantsBcftoolsOnly = new HashSet<>();
@@ -269,20 +271,20 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
             }
         }
 
-        String description = String.format("Total Variants: %s\nTotal GT 1 PCT: %s\nTotal GT 50 PCT: %s\nTotal Indel GT 1 PCT: %s", totalVariants, totalGT1, totalGT50, totalIndelGT1);
+        String description = String.format("Total Variants: %s\nTotal GT 1 PCT: %s\nTotal GT 50 PCT: %s\nTotal Indel GT 1 PCT: %s\tPositions Below Coverage: %s", totalVariants, totalGT1, totalGT50, totalIndelGT1, positionsSkipped);
 
         if (!variantsBcftoolsOnly.isEmpty())
         {
             getPipelineCtx().getLogger().error("The following variants were in bcftools, but not GT50% in lofreq: ");
-            variantsBcftoolsOnly.forEach(vc -> getPipelineCtx().getLogger().error(getHashKey(vc)));
+            variantsBcftoolsOnly.forEach(vc -> getPipelineCtx().getLogger().error(getHashKey(vc) + ", DP=" + vc.getAttribute("DP")));
 
             description += "\n" + "WARNING: " + variantsBcftoolsOnly.size() + " variants detected in bcftools and not lofreq";
         }
 
         if (!alleleToAF.isEmpty())
         {
             getPipelineCtx().getLogger().error("The following variants were GT50% in lofreq, but not in bcftools: ");
-            alleleToAF.keySet().forEach(vc -> getPipelineCtx().getLogger().error(vc));
+            alleleToAF.keySet().forEach(vc -> getPipelineCtx().getLogger().error(vc + ", AF=" + alleleToAF.get(vc)));
 
             description += "\n" + "WARNING: " + alleleToAF.size() + " variants detected in lofreq and not bcftools";
         }
@@ -293,23 +295,44 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
             throw new PipelineJobException("Expected file not found: " + consensusFasta.getPath());
         }
 
+        DNASequence seq;
         try (InputStream is = IOUtil.openFileForReading(consensusFasta))
         {
             FastaReader<DNASequence, NucleotideCompound> fastaReader = new FastaReader<>(is, new GenericFastaHeaderParser<>(), new DNASequenceCreator(AmbiguityDNACompoundSet.getDNACompoundSet()));
             LinkedHashMap<String, DNASequence> fastaData = fastaReader.process();
 
-            for (String fastaHeader : fastaData.keySet())
-            {
-                AtomicInteger totalN = new AtomicInteger();
-                DNASequence seq = fastaData.get(fastaHeader);
-                seq.forEach(nt -> {
-                    if (nt.getUpperedBase().equals("N")) {
-                        totalN.getAndIncrement();
-                    }
-                });
+            AtomicInteger totalN = new AtomicInteger();
+            seq = fastaData.values().iterator().next();
+            seq.forEach(nt -> {
+                if (nt.getUpperedBase().equals("N")) {
+                    totalN.getAndIncrement();
+                }
+            });
 
-                description += "\nConsensus Ns: " + totalN.get();
+            description += "\nConsensus Ns: " + totalN.get();
+        }
+        catch (IOException e)
+        {
+            throw new PipelineJobException(e);
+        }
+
+        //Replace FASTA header:
+        try (PrintWriter writer = PrintWriters.getPrintWriter(consensusFasta))
+        {
+            StringBuilder header = new StringBuilder();
+            if (rs.getSubjectId() != null)
+            {
+                header.append(rs.getSubjectId()).append("|");
             }
+            else
+            {
+                header.append(rs.getName()).append("|");
+            }
+
+            header.append(rs.getLibraryType() == null ? rs.getApplication() : rs.getLibraryType());
+
+            writer.println(">" + header);
+            writer.println(seq.getSequenceAsString());
         }
         catch (IOException e)
         {