use of com.milaboratory.core.Target in project mixcr by milaboratory.
the class ActionAlign method go.
@Override
@SuppressWarnings("unchecked")
public void go(ActionHelper helper) throws Exception {
// FIXME remove in 2.2
if (actionParameters.printNonFunctionalWarnings())
System.out.println("WARNING: -wf / --non-functional-warnings option is deprecated, will be removed in 2.2 " + "release. Use --verbose instead.");
// Saving initial timestamp
long beginTimestamp = System.currentTimeMillis();
// Getting aligner parameters
VDJCAlignerParameters alignerParameters = actionParameters.getAlignerParameters();
// FIXME remove in 2.3
if (actionParameters.getSaveOriginalReads()) {
System.out.println("WARNING: -g / --save-reads option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
alignerParameters.setSaveOriginalReads(true);
}
// FIXME remove in 2.3
if (actionParameters.getSaveReadDescription()) {
System.out.println("WARNING: -a / --save-description option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
alignerParameters.setSaveOriginalReads(true);
}
if (!actionParameters.overrides.isEmpty()) {
// Perform parameters overriding
alignerParameters = JsonOverrider.override(alignerParameters, VDJCAlignerParameters.class, actionParameters.overrides);
if (alignerParameters == null)
throw new ProcessException("Failed to override some parameter.");
}
// FIXME remove in 2.2
if (actionParameters.allowDifferentVJLoci != null && actionParameters.allowDifferentVJLoci) {
System.out.println("Warning: usage of --diff-loci is deprecated. Use -OallowChimeras=true instead.");
alignerParameters.setAllowChimeras(true);
}
// Creating aligner
VDJCAligner aligner = VDJCAligner.createAligner(alignerParameters, actionParameters.isInputPaired(), !actionParameters.getNoMerge());
// Detect if automatic featureToAlign correction is required
int totalV = 0, totalVErrors = 0, hasVRegion = 0;
GeneFeature correctingFeature = alignerParameters.getVAlignerParameters().getGeneFeatureToAlign().hasReversedRegions() ? GeneFeature.VRegionWithP : GeneFeature.VRegion;
VDJCLibrary library = VDJCLibraryRegistry.getDefault().getLibrary(actionParameters.library, actionParameters.species);
System.out.println("Reference library: " + library.getLibraryId());
// Printing library level warnings, if specified for the library
if (!library.getWarnings().isEmpty()) {
System.out.println("Library warnings:");
for (String l : library.getWarnings()) System.out.println(l);
}
// Printing citation notice, if specified for the library
if (!library.getCitations().isEmpty()) {
System.out.println("Please cite:");
for (String l : library.getCitations()) System.out.println(l);
}
for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
if (gene.getGeneType() == GeneType.Variable) {
totalV++;
if (!alignerParameters.containsRequiredFeature(gene)) {
totalVErrors++;
if (gene.getPartitioning().isAvailable(correctingFeature))
hasVRegion++;
}
}
}
// Performing V featureToAlign correction if needed
if (totalVErrors > totalV * 0.9 && hasVRegion > totalVErrors * 0.8) {
System.out.println("WARNING: forcing -OvParameters.geneFeatureToAlign=" + GeneFeature.encode(correctingFeature) + " since current gene feature (" + GeneFeature.encode(alignerParameters.getVAlignerParameters().getGeneFeatureToAlign()) + ") is absent in " + Util.PERCENT_FORMAT.format(100.0 * totalVErrors / totalV) + "% of V genes.");
alignerParameters.getVAlignerParameters().setGeneFeatureToAlign(correctingFeature);
}
int numberOfExcludedNFGenes = 0;
int numberOfExcludedFGenes = 0;
for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
NucleotideSequence featureSequence = alignerParameters.extractFeatureToAlign(gene);
// exclusionReason is null ==> gene is not excluded
String exclusionReason = null;
if (featureSequence == null)
exclusionReason = "absent " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
else if (featureSequence.containsWildcards())
exclusionReason = "wildcard symbols in " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
if (exclusionReason == null)
// If there are no reasons to exclude the gene, adding it to aligner
aligner.addGene(gene);
else {
if (gene.isFunctional()) {
++numberOfExcludedFGenes;
if (actionParameters.verbose())
System.out.println("WARNING: Functional gene " + gene.getName() + " excluded due to " + exclusionReason);
} else
++numberOfExcludedNFGenes;
}
}
if (actionParameters.printWarnings() && numberOfExcludedFGenes > 0)
System.out.println("WARNING: " + numberOfExcludedFGenes + " functional genes were excluded, re-run " + "with --verbose option to see the list of excluded genes and exclusion reason.");
if (actionParameters.verbose() && numberOfExcludedNFGenes > 0)
System.out.println("WARNING: " + numberOfExcludedNFGenes + " non-functional genes excluded.");
if (aligner.getVGenesToAlign().isEmpty())
throw new ProcessException("No V genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
if (aligner.getJGenesToAlign().isEmpty())
throw new ProcessException("No J genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
AlignerReport report = new AlignerReport();
report.setStartMillis(beginTimestamp);
report.setInputFiles(actionParameters.getInputsForReport());
report.setOutputFiles(actionParameters.getOutputsForReport());
report.setCommandLine(helper.getCommandLineArguments());
// Attaching report to aligner
aligner.setEventsListener(report);
try (SequenceReaderCloseable<? extends SequenceRead> reader = actionParameters.createReader();
VDJCAlignmentsWriter writer = actionParameters.getOutputName().equals(".") ? null : new VDJCAlignmentsWriter(actionParameters.getOutputName());
SequenceWriter notAlignedWriter = actionParameters.failedReadsR1 == null ? null : (actionParameters.isInputPaired() ? new PairedFastqWriter(actionParameters.failedReadsR1, actionParameters.failedReadsR2) : new SingleFastqWriter(actionParameters.failedReadsR1))) {
if (writer != null)
writer.header(aligner);
OutputPort<? extends SequenceRead> sReads = reader;
CanReportProgress progress = (CanReportProgress) reader;
if (actionParameters.limit != 0) {
sReads = new CountLimitingOutputPort<>(sReads, actionParameters.limit);
progress = SmartProgressReporter.extractProgress((CountLimitingOutputPort<?>) sReads);
}
final boolean writeAllResults = actionParameters.getWriteAllResults();
EnumMap<GeneType, VDJCHit[]> emptyHits = new EnumMap<>(GeneType.class);
for (GeneType gt : GeneType.values()) if (alignerParameters.getGeneAlignerParameters(gt) != null)
emptyHits.put(gt, new VDJCHit[0]);
final PairedEndReadsLayout readsLayout = alignerParameters.getReadsLayout();
SmartProgressReporter.startProgressReport("Alignment", progress);
OutputPort<Chunk<? extends SequenceRead>> mainInputReads = CUtils.buffered((OutputPort) chunked(sReads, 64), 16);
OutputPort<VDJCAlignmentResult> alignments = unchunked(new ParallelProcessor(mainInputReads, chunked(aligner), actionParameters.threads));
for (VDJCAlignmentResult result : CUtils.it(new OrderedOutputPort<>(alignments, new Indexer<VDJCAlignmentResult>() {
@Override
public long getIndex(VDJCAlignmentResult o) {
return o.read.getId();
}
}))) {
VDJCAlignments alignment = result.alignment;
SequenceRead read = result.read;
if (alignment == null) {
if (writeAllResults) // Creating empty alignment object if alignment for current read failed
{
Target target = readsLayout.createTargets(read)[0];
alignment = new VDJCAlignments(emptyHits, target.targets, SequenceHistory.RawSequence.of(read.getId(), target), alignerParameters.isSaveOriginalReads() ? new SequenceRead[] { read } : null);
} else {
if (notAlignedWriter != null)
notAlignedWriter.write(result.read);
continue;
}
}
if (alignment.isChimera())
report.onChimera();
if (writer != null)
writer.write(alignment);
}
if (writer != null)
writer.setNumberOfProcessedReads(reader.getNumberOfReads());
}
report.setFinishMillis(System.currentTimeMillis());
// Writing report to stout
System.out.println("============= Report ==============");
Util.writeReportToStdout(report);
if (actionParameters.report != null)
Util.writeReport(actionParameters.report, report);
if (actionParameters.jsonReport != null)
Util.writeJsonReport(actionParameters.jsonReport, report);
}
use of com.milaboratory.core.Target in project mixcr by milaboratory.
the class VDJCAlignerPVFirst method process0.
@Override
protected VDJCAlignmentResult<PairedRead> process0(final PairedRead input) {
Target[] targets = getTargets(input);
// Creates helper classes for each PTarget
PAlignmentHelper[] helpers = createInitialHelpers(targets);
VDJCAlignmentResult<PairedRead> result = parameters.getAllowPartialAlignments() ? processPartial(input, helpers) : processStrict(input, helpers);
// if sAligner == null (which means --no-merge option), no merge will be performed
if (result.alignment != null && sAligner != null) {
final VDJCAlignments alignment = result.alignment;
final TargetMerger.TargetMergingResult mergeResult = alignmentsMerger.merge(new AlignedTarget(alignment, 0), new AlignedTarget(alignment, 1), false);
if (mergeResult.failedDueInconsistentAlignments()) {
GeneType geneType = mergeResult.getFailedMergedGeneType();
int removeId = alignment.getBestHit(geneType).getAlignment(0).getScore() > alignment.getBestHit(geneType).getAlignment(1).getScore() ? 1 : 0;
if (listener != null)
listener.onTopHitSequenceConflict(input, alignment, geneType);
return new VDJCAlignmentResult<>(input, alignment.removeBestHitAlignment(geneType, removeId));
} else if (mergeResult.isSuccessful()) {
assert mergeResult.isUsingAlignments();
NSequenceWithQuality alignedTarget = mergeResult.getResult().getTarget();
SingleRead sRead = new SingleReadImpl(input.getId(), alignedTarget, "");
VDJCAlignmentResult<SingleRead> sResult = sAligner.process0(sRead);
if (sResult.alignment == null)
return result;
VDJCAlignments sAlignment = sResult.alignment.setHistory(new SequenceHistory[] { new SequenceHistory.Merge(SequenceHistory.OverlapType.AlignmentOverlap, result.alignment.getHistory(0), result.alignment.getHistory(1), mergeResult.getOffset(), mergeResult.getMismatched()) }, new SequenceRead[] { input });
if (listener != null)
listener.onSuccessfulAlignmentOverlap(input, sAlignment);
return new VDJCAlignmentResult<>(input, sAlignment);
}
}
return result;
}
use of com.milaboratory.core.Target in project mixcr by milaboratory.
the class VDJCAlignerS method processStrict.
private VDJCAlignmentResult<SingleRead> processStrict(SingleRead input) {
Target[] targets = parameters.getReadsLayout().createTargets(input);
// Algorithm below relies on this fact
assert targets.length <= 2;
boolean anyIsFull = false, anyHasJ = false, anyHasV = false;
KVJResultsForSingle[] results = new KVJResultsForSingle[targets.length];
for (int i = 0; i < results.length; i++) {
results[i] = align(targets[i]);
results[i].performChainFilteringIfNeeded();
anyIsFull |= results[i].hasKVAndJHits();
anyHasJ |= results[i].hasKJHits();
anyHasV |= results[i].hasKVHits();
}
if (!anyIsFull) {
if (!anyHasJ && !anyHasV)
onFailedAlignment(input, VDJCAlignmentFailCause.NoHits);
else if (!anyHasJ)
onFailedAlignment(input, VDJCAlignmentFailCause.NoJHits);
else
onFailedAlignment(input, VDJCAlignmentFailCause.NoVHits);
return new VDJCAlignmentResult<>(input);
}
KVJResultsForSingle topResult = null;
// Calculating best result
for (KVJResultsForSingle result : results) if (result.hasKVAndJHits())
topResult = topResult != null ? null : result;
if (topResult == null) {
// Finalizing alignment for both results to determine who is the best
for (KVJResultsForSingle result : results) {
result.alignDC();
if (topResult == null || topResult.sumScore() < result.sumScore())
topResult = result;
}
} else
// Align C and D genes only for best result
topResult.alignDC();
// Checking minimal sum score
if (topResult.sumScore() < parameters.getMinSumScore()) {
onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
return new VDJCAlignmentResult<>(input);
}
// Filtering hits basing on minSumScore
topResult.calculateHits(parameters.getMinSumScore(), parameters.getMaxHits());
if (topResult.hasVAndJHits()) {
VDJCAlignments alignment = topResult.toVDJCAlignments(input.getId(), input);
onSuccessfulAlignment(input, alignment);
return new VDJCAlignmentResult<>(input, alignment);
} else {
onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
return new VDJCAlignmentResult<>(input);
}
}
use of com.milaboratory.core.Target in project mixcr by milaboratory.
the class VDJCAlignerS method processPartial.
private VDJCAlignmentResult<SingleRead> processPartial(SingleRead input) {
Target[] targets = parameters.getReadsLayout().createTargets(input);
KVJResultsForSingle[] results = new KVJResultsForSingle[targets.length];
for (int i = 0; i < results.length; i++) {
results[i] = align(targets[i]);
results[i].performChainFilteringIfNeeded();
}
KVJResultsForSingle topResult = null;
for (KVJResultsForSingle result : results) {
if (!result.hasNoKVNorKJHits())
result.alignDC();
if (topResult == null || topResult.sumScore() < result.sumScore())
topResult = result;
}
if (topResult.hasNoKVNorKJHits()) {
onFailedAlignment(input, VDJCAlignmentFailCause.NoHits);
return new VDJCAlignmentResult<>(input);
}
// Checking minimal sum score
if (topResult.sumScore() < parameters.getMinSumScore()) {
onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
return new VDJCAlignmentResult<>(input);
}
// Filtering hits basing on minSumScore
topResult.calculateHits(parameters.getMinSumScore(), parameters.getMaxHits());
VDJCAlignments alignment = topResult.toVDJCAlignments(input.getId(), input);
// Final check
if (!parameters.getAllowNoCDR3PartAlignments()) {
// CDR3 Begin / End
if (!alignment.getPartitionedTarget(0).getPartitioning().isAvailable(reqPointL) && !alignment.getPartitionedTarget(0).getPartitioning().isAvailable(reqPointR)) {
onFailedAlignment(input, VDJCAlignmentFailCause.NoCDR3Parts);
return new VDJCAlignmentResult<>(input);
}
}
// Read successfully aligned
onSuccessfulAlignment(input, alignment);
return new VDJCAlignmentResult<>(input, alignment);
}
Aggregations