More reporting around viral consensus generation

Author: bbimber

SequenceAnalysis/resources/external/viral_consensus.sh

@@ -23,39 +23,43 @@ OUT_IUPAC=`basename $BAM .bam`".consensus.iupac.fasta"
 # call variants
 echo 'Calling variants'
 $BCFTOOLS mpileup -Ou -d 20000 -f $FASTA $BAM | $BCFTOOLS call --ploidy 1 --threads $THREADS -mv -Oz -o $VCF_CALLS
-$BCFTOOLS index $VCF_CALLS
+$BCFTOOLS index -t $VCF_CALLS
 COUNT1=`$BCFTOOLS view -H $VCF_CALLS | wc -l`
 echo 'Variants called: '$COUNT1
 echo -e 'VariantsCalled\t'$COUNT1 >> $REPORT
 
 # normalize indels
 echo 'Normalize indels'
-VCF_NORM=calls.norm.bcf
-$BCFTOOLS norm -f $FASTA --threads $THREADS -Ob -o $VCF_NORM $VCF_CALLS
-$BCFTOOLS index -f $VCF_NORM
+VCF_NORM=calls.norm.vcf.gz
+$BCFTOOLS norm -f $FASTA --threads $THREADS -Oz -o $VCF_NORM $VCF_CALLS
+$BCFTOOLS index -t -f $VCF_NORM
 COUNT2=`$BCFTOOLS view -H $VCF_NORM | wc -l`
 echo 'Variants remaining: '$COUNT2
 echo -e 'VariantsAfterNorm\t'$COUNT2 >> $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.bcf
-$BCFTOOLS filter --IndelGap 5 -Ob -o $VCF_INDEL_FILTER $VCF_NORM
-$BCFTOOLS index -f $VCF_INDEL_FILTER
+VCF_INDEL_FILTER=calls.norm.flt-indels.vcf.gz
+$BCFTOOLS filter --IndelGap 5 -Oz -o $VCF_INDEL_FILTER $VCF_NORM
+$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 user the filtered version:
-VCF_FOR_CONSENSUS=$VCF_NORM
+#At the moment, do not use the 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
+
 $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}.csi
+rm ${VCF_CALLS}.tbi
 
-rm $VCF_NORM
-rm ${VCF_NORM}.csi
+#Previously moved:
+#rm $VCF_NORM
+#rm ${VCF_NORM}.tbi
 
 rm $VCF_INDEL_FILTER
-rm ${VCF_INDEL_FILTER}.csi
+rm ${VCF_INDEL_FILTER}.tbi

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

@@ -5,6 +5,7 @@
 import htsjdk.samtools.util.CloseableIterator;
 import htsjdk.samtools.util.IOUtil;
 import htsjdk.samtools.util.Interval;
+import htsjdk.variant.variantcontext.Allele;
 import htsjdk.variant.variantcontext.VariantContext;
 import htsjdk.variant.vcf.VCFFileReader;
 import org.apache.log4j.Logger;
@@ -37,7 +38,12 @@
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Collections;
+import java.util.HashMap;
+import java.util.HashSet;
 import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import java.util.stream.Collectors;
 
 public class LofreqAnalysis extends AbstractCommandPipelineStep<LofreqAnalysis.LofreqWrapper> implements AnalysisStep
 {
@@ -118,6 +124,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         int totalVariants = 0;
         int totalGT1 = 0;
         int totalGT50 = 0;
+        Map<String, Double> alleleToAF = new HashMap<>();
         int totalIndelGT1 = 0;
         try (VCFFileReader reader = new VCFFileReader(outputVcfSnpEff);CloseableIterator<VariantContext> it = reader.iterator())
         {
@@ -137,6 +144,10 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
                 if (vc.hasAttribute("AF") && vc.getAttributeAsDouble("AF", 0.0) > 0.5)
                 {
                     totalGT50++;
+                    String key = getHashKey(vc);
+                    Double af = alleleToAF.getOrDefault(key, 0.0);
+                    af = Math.max(af, vc.getAttributeAsDouble("AF", 0.0));
+                    alleleToAF.put(key, af);
                 }
             }
         }
@@ -228,10 +239,51 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
         getPipelineCtx().getLogger().info("Total intervals of these gaps: " + gapIntervals);
 
         consensusWrapper.execute(Arrays.asList("/bin/bash", script.getPath(), inputBam.getPath(), referenceGenome.getWorkingFastaFile().getPath(), mask.getPath()));
+        File calls = new File(inputBam.getParentFile(), FileUtil.getBaseName(inputBam) + ".calls.vcf.gz");
+
+        Set<VariantContext> variantsBcftoolsOnly = new HashSet<>();
+        try (VCFFileReader reader = new VCFFileReader(calls);CloseableIterator<VariantContext> it = reader.iterator())
+        {
+            while (it.hasNext())
+            {
+                VariantContext vc = it.next();
+                String key = getHashKey(vc);
+                if (alleleToAF.containsKey(key))
+                {
+                    //Variant shared
+                    alleleToAF.remove(key);
+                }
+                else
+                {
+                    variantsBcftoolsOnly.add(vc);
+                }
+            }
+        }
 
         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);
+
+        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)));
+
+            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));
+
+            description += "\n" + "WARNING: " + alleleToAF.size() + " variants detected in lofreq and not bcftools";
+        }
+
+
         output.addSequenceOutput(outputVcfSnpEff, "LoFreq: " + rs.getName(), CATEGORY, rs.getReadsetId(), null, referenceGenome.getGenomeId(), description);
-        output.addSequenceOutput(coverageOut, "Depth of Coverage: " + rs.getName(), "Depth of Coverage", rs.getReadsetId(), null, referenceGenome.getGenomeId(), description);
+        output.addSequenceOutput(coverageOut, "Depth of Coverage: " + rs.getName(), "Depth of Coverage", rs.getReadsetId(), null, referenceGenome.getGenomeId(), null);
+
+        File consensusFasta = new File(inputBam.getParentFile(), FileUtil.getBaseName(inputBam.getName()) + ".consensus.fasta");
+        output.addSequenceOutput(consensusFasta, "Consensus: " + rs.getName(), "Viral Consensus Sequence", rs.getReadsetId(), null, referenceGenome.getGenomeId(), description);
 
         return output;
     }
@@ -242,6 +294,11 @@ public Output performAnalysisPerSampleLocal(AnalysisModel model, File inputBam,
         return null;
     }
 
+    private String getHashKey(VariantContext vc)
+    {
+        return vc.getContig() + "<>" + vc.getStart() + vc.getAlternateAlleles().stream().map(Allele::getBaseString).collect(Collectors.joining(";"));
+    }
+
     public static class LofreqWrapper extends AbstractCommandWrapper
     {
         public LofreqWrapper(Logger log)