Skip to content

Commit

Permalink
SW implementation is command line argument in FilterAlignmentArtifacts (
Browse files Browse the repository at this point in the history
  • Loading branch information
davidbenjamin committed Oct 26, 2021
1 parent 64a47cd commit f74cf29
Showing 1 changed file with 9 additions and 7 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,7 @@
import htsjdk.variant.variantcontext.writer.VariantContextWriter;
import htsjdk.variant.vcf.VCFHeader;
import htsjdk.variant.vcf.VCFHeaderLine;
import org.broadinstitute.barclay.argparser.Argument;
import org.broadinstitute.barclay.argparser.ArgumentCollection;
import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
import org.broadinstitute.barclay.argparser.ExperimentalFeature;
import org.broadinstitute.barclay.argparser.*;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.gatk.nativebindings.smithwaterman.SWParameters;
import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
Expand Down Expand Up @@ -112,7 +109,6 @@ public class FilterAlignmentArtifacts extends MultiVariantWalkerGroupedOnStart {
public static final int DEFAULT_MAX_GROUPED_SPAN = 10_000;
private static final int MIN_UNITIG_LENGTH = 30;
private static final int ASSEMBLY_PADDING = 50;
private static final SmithWatermanAligner ALIGNER = SmithWatermanAligner.getAligner(SmithWatermanAligner.Implementation.FASTEST_AVAILABLE);

@Argument(fullName = StandardArgumentDefinitions.OUTPUT_LONG_NAME, shortName = StandardArgumentDefinitions.OUTPUT_SHORT_NAME,
doc="The output filtered VCF file", optional=false)
Expand All @@ -136,10 +132,15 @@ public class FilterAlignmentArtifacts extends MultiVariantWalkerGroupedOnStart {
@Argument(fullName= AssemblyBasedCallerArgumentCollection.BAM_OUTPUT_LONG_NAME, shortName= AssemblyBasedCallerArgumentCollection.BAM_OUTPUT_SHORT_NAME, doc="File to which assembled haplotypes should be written", optional = true)
public String bamOutputPath = null;

@Advanced
@Argument(fullName = AssemblyBasedCallerArgumentCollection.SMITH_WATERMAN_LONG_NAME, doc = "Smith-Waterman implementation to use, generally FASTEST_AVAILABLE is the right choice", optional = true)
public SmithWatermanAligner.Implementation smithWatermanImplementation = SmithWatermanAligner.Implementation.FASTEST_AVAILABLE;


@ArgumentCollection
protected RealignmentArgumentCollection realignmentArgumentCollection = new RealignmentArgumentCollection();

private SmithWatermanAligner smithWatermanAligner;
private VariantContextWriter vcfWriter;
private RealignmentEngine realignmentEngine;
private SAMFileHeader bamHeader;
Expand Down Expand Up @@ -174,6 +175,7 @@ protected int defaultMaxGroupedSpan() {

@Override
public void onTraversalStart() {
smithWatermanAligner = SmithWatermanAligner.getAligner(smithWatermanImplementation);
realignmentEngine = new RealignmentEngine(realignmentArgumentCollection);
vcfWriter = createVCFWriter(outputVcf);

Expand Down Expand Up @@ -211,15 +213,15 @@ public void apply(List<VariantContext> variantContexts, ReferenceContext referen


// TODO: give this tool M2 Assembler args to allow override default M2ArgumentCollection?
final AssemblyResultSet assemblyResult = AssemblyBasedCallerUtils.assembleReads(assemblyRegion, Collections.emptyList(), MTAC, bamHeader, samplesList, logger, referenceReader, assemblyEngine, ALIGNER, false);
final AssemblyResultSet assemblyResult = AssemblyBasedCallerUtils.assembleReads(assemblyRegion, Collections.emptyList(), MTAC, bamHeader, samplesList, logger, referenceReader, assemblyEngine, smithWatermanAligner, false);
final AssemblyRegion regionForGenotyping = assemblyResult.getRegionForGenotyping();

final Map<String,List<GATKRead>> reads = AssemblyBasedCallerUtils.splitReadsBySample(samplesList, bamHeader, regionForGenotyping.getReads());

final AlleleLikelihoods<GATKRead, Haplotype> readLikelihoods = likelihoodCalculationEngine.computeReadLikelihoods(assemblyResult,samplesList,reads);
readLikelihoods.switchToNaturalLog();
final SWParameters readToHaplotypeSWParameters = MTAC.getReadToHaplotypeSWParameters();
final Map<GATKRead,GATKRead> readRealignments = AssemblyBasedCallerUtils.realignReadsToTheirBestHaplotype(readLikelihoods, assemblyResult.getReferenceHaplotype(), assemblyResult.getPaddedReferenceLoc(), ALIGNER, readToHaplotypeSWParameters);
final Map<GATKRead,GATKRead> readRealignments = AssemblyBasedCallerUtils.realignReadsToTheirBestHaplotype(readLikelihoods, assemblyResult.getReferenceHaplotype(), assemblyResult.getPaddedReferenceLoc(), smithWatermanAligner, readToHaplotypeSWParameters);
readLikelihoods.changeEvidence(readRealignments);
writeBamOutput(assemblyResult, readLikelihoods, new HashSet<>(readLikelihoods.alleles()), regionForGenotyping.getSpan());

Expand Down

0 comments on commit f74cf29

Please sign in to comment.