Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SW implementation is command line argument in FilterAlignmentArtifacts #7105

Merged
merged 4 commits into from
Oct 26, 2021
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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.*;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we want to use explicit imports

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

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