use of com.milaboratory.mixcr.vdjaligners.VDJCAlignmentResult 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.mixcr.vdjaligners.VDJCAlignmentResult in project mixcr by milaboratory.
the class PartialAlignmentsAssemblerAligner method process0.
@Override
@SuppressWarnings("unchecked")
protected VDJCAlignmentResult<VDJCMultiRead> process0(VDJCMultiRead input) {
final int nReads = input.numberOfReads();
EnumMap<GeneType, VDJCHit[]> vdjcHits = new EnumMap<>(GeneType.class);
NSequenceWithQuality[] targets = new NSequenceWithQuality[nReads];
Chains currentChains = Chains.ALL;
// Across all gene types
int lastAlignedTarget = 0;
int firstJTarget = -1;
int lastVTarget = -1;
for (int g = 0; g < GeneType.VJC_REFERENCE.length; g++) {
GeneType gt = GeneType.VJC_REFERENCE[g];
AlignmentHit<NucleotideSequence, VDJCGene>[][] alignmentHits = new AlignmentHit[nReads][];
Arrays.fill(alignmentHits, new AlignmentHit[0]);
for (int targetId = lastAlignedTarget; targetId < nReads; targetId++) {
targets[targetId] = input.getRead(targetId).getData();
final NucleotideSequence sequence = input.getRead(targetId).getData().getSequence();
AlignmentResult<AlignmentHit<NucleotideSequence, VDJCGene>> als;
final BatchAlignerWithBaseWithFilter<NucleotideSequence, VDJCGene, AlignmentHit<NucleotideSequence, VDJCGene>> aligner = getAligner(gt);
if (aligner != null) {
int pointer = 0;
if (g != 0) {
// Not V gene
VDJCHit[] vdjcHits1 = vdjcHits.get(GeneType.VJC_REFERENCE[g - 1]);
Alignment<NucleotideSequence> alignment;
if (vdjcHits1.length != 0 && (alignment = vdjcHits1[0].getAlignment(targetId)) != null)
pointer = alignment.getSequence2Range().getTo();
}
als = aligner.align(sequence, pointer, sequence.size(), getFilter(gt, currentChains));
if (als != null && als.hasHits()) {
lastAlignedTarget = targetId;
if (// V
g == 0)
lastVTarget = targetId;
if (// J
g == 1)
firstJTarget = targetId;
alignmentHits[targetId] = als.getHits().toArray(new AlignmentHit[als.getHits().size()]);
}
}
}
Chains chains = Chains.EMPTY;
for (AlignmentHit<NucleotideSequence, VDJCGene>[] alignmentHit0 : alignmentHits) if (alignmentHit0 != null)
for (AlignmentHit<NucleotideSequence, VDJCGene> hit : alignmentHit0) chains = chains.merge(hit.getRecordPayload().getChains());
currentChains = currentChains.intersection(chains);
vdjcHits.put(gt, combine(parameters.getFeatureToAlign(gt), alignmentHits));
}
boolean fineVAlignmentPerformed = false, fineJAlignmentPerformed = false;
// Additional (fine) alignment step for V gene
VDJCHit[] vHits = vdjcHits.get(GeneType.Variable);
final AlignmentScoring<NucleotideSequence> vScoring = parameters.getVAlignerParameters().getParameters().getScoring();
if (// TODO implement AffineGapAlignmentScoring
vHits != null && vHits.length > 0 && !(vScoring instanceof AffineGapAlignmentScoring) && vdjcHits.get(GeneType.Joining) != null && vdjcHits.get(GeneType.Joining).length > 0) {
int minimalVSpace = getAbsoluteMinScore(parameters.getVAlignerParameters().getParameters()) / vScoring.getMaximalMatchScore();
// Assert
if (firstJTarget == -1)
throw new AssertionError();
for (int targetId = 1; targetId <= firstJTarget; targetId++) {
int vSpace;
final NucleotideSequence sequence2 = targets[targetId].getSequence();
if (vdjcHits.get(GeneType.Joining)[0].getAlignment(targetId) != null && (vSpace = vdjcHits.get(GeneType.Joining)[0].getAlignment(targetId).getSequence2Range().getFrom()) >= minimalVSpace) {
for (int vHitIndex = 0; vHitIndex < vHits.length; vHitIndex++) {
VDJCHit vHit = vHits[vHitIndex];
// Perform fine alignment only if target is not already aligned by fast aligner
if (vHit.getAlignment(targetId) != null)
continue;
Alignment<NucleotideSequence> leftAlignment = vHit.getAlignment(targetId - 1);
if (leftAlignment == null)
continue;
final NucleotideSequence sequence1 = leftAlignment.getSequence1();
final int beginFR3 = vHit.getGene().getPartitioning().getRelativePosition(parameters.getFeatureToAlign(GeneType.Variable), ReferencePoint.FR3Begin);
if (beginFR3 == -1)
continue;
final Alignment alignment = AlignerCustom.alignLinearSemiLocalLeft0((LinearGapAlignmentScoring<NucleotideSequence>) vScoring, sequence1, sequence2, beginFR3, sequence1.size() - beginFR3, 0, vSpace, false, true, NucleotideSequence.ALPHABET, linearMatrixCache.get());
if (alignment.getScore() < getAbsoluteMinScore(parameters.getVAlignerParameters().getParameters()))
continue;
fineVAlignmentPerformed = true;
vHits[vHitIndex] = vHit.setAlignment(targetId, alignment);
}
}
}
}
Arrays.sort(vHits);
vdjcHits.put(GeneType.Variable, cutRelativeScore(vHits, parameters.getVAlignerParameters().getRelativeMinScore(), parameters.getVAlignerParameters().getParameters().getMaxHits()));
// Additional (fine) alignment step for J gene
VDJCHit[] jHits = vdjcHits.get(GeneType.Joining);
final AlignmentScoring<NucleotideSequence> jScoring = parameters.getJAlignerParameters().getParameters().getScoring();
if (// TODO implement AffineGapAlignmentScoring
jHits != null && jHits.length > 0 && !(jScoring instanceof AffineGapAlignmentScoring) && vdjcHits.get(GeneType.Variable) != null && vdjcHits.get(GeneType.Variable).length > 0) {
int minimalJSpace = getAbsoluteMinScore(parameters.getJAlignerParameters().getParameters()) / jScoring.getMaximalMatchScore();
// Assert
if (lastVTarget == -1)
throw new AssertionError();
for (int targetId = lastVTarget; targetId < nReads - 1; targetId++) {
int jSpaceBegin;
final NucleotideSequence sequence2 = targets[targetId].getSequence();
if (vdjcHits.get(GeneType.Variable)[0].getAlignment(targetId) != null && (sequence2.size() - (jSpaceBegin = vdjcHits.get(GeneType.Variable)[0].getAlignment(targetId).getSequence2Range().getTo())) >= minimalJSpace) {
for (int jHitIndex = 0; jHitIndex < jHits.length; jHitIndex++) {
VDJCHit jHit = jHits[jHitIndex];
// Perform fine alignment only if target is not already aligned by fast aligner
if (jHit.getAlignment(targetId) != null)
continue;
Alignment<NucleotideSequence> rightAlignment = jHit.getAlignment(targetId + 1);
if (rightAlignment == null)
continue;
final NucleotideSequence sequence1 = rightAlignment.getSequence1();
final Alignment alignment = AlignerCustom.alignLinearSemiLocalRight0((LinearGapAlignmentScoring) jScoring, sequence1, sequence2, 0, sequence1.size(), jSpaceBegin, sequence2.size() - jSpaceBegin, false, true, NucleotideSequence.ALPHABET, linearMatrixCache.get());
if (alignment.getScore() < getAbsoluteMinScore(parameters.getJAlignerParameters().getParameters()))
continue;
fineJAlignmentPerformed = true;
jHits[jHitIndex] = jHit.setAlignment(targetId, alignment);
}
}
}
}
Arrays.sort(jHits);
vdjcHits.put(GeneType.Joining, cutRelativeScore(jHits, parameters.getJAlignerParameters().getRelativeMinScore(), parameters.getJAlignerParameters().getParameters().getMaxHits()));
int dGeneTarget = -1;
VDJCHit[] vResult = vdjcHits.get(GeneType.Variable);
VDJCHit[] jResult = vdjcHits.get(GeneType.Joining);
if (vResult.length != 0 && jResult.length != 0)
for (int i = 0; i < nReads; i++) if (vResult[0].getAlignment(i) != null && jResult[0].getAlignment(i) != null) {
dGeneTarget = i;
break;
}
// if (fineVAlignmentPerformed && fineJAlignmentPerformed)
// System.out.println("sd");
VDJCHit[] dResult;
if (dGeneTarget == -1)
dResult = new VDJCHit[0];
else {
final Alignment<NucleotideSequence> vAl = vResult[0].getAlignment(dGeneTarget);
final Alignment<NucleotideSequence> jAl = jResult[0].getAlignment(dGeneTarget);
if (vAl == null || jAl == null || singleDAligner == null)
dResult = new VDJCHit[0];
else
dResult = singleDAligner.align(targets[dGeneTarget].getSequence(), getPossibleDLoci(vResult, jResult), vAl.getSequence2Range().getTo(), jAl.getSequence2Range().getFrom(), dGeneTarget, nReads);
}
final VDJCAlignments alignment = new VDJCAlignments(vResult, dResult, jResult, cutRelativeScore(vdjcHits.get(GeneType.Constant), parameters.getCAlignerParameters().getRelativeMinScore(), parameters.getMaxHits()), targets, input.getHistory(), input.getOriginalReads());
return new VDJCAlignmentResult<>(input, alignment);
}
use of com.milaboratory.mixcr.vdjaligners.VDJCAlignmentResult in project mixcr by milaboratory.
the class RunMiXCR method align.
public static AlignResult align(RunMiXCRAnalysis parameters) throws Exception {
VDJCAlignerParameters alignerParameters = parameters.alignerParameters;
VDJCAligner aligner = VDJCAligner.createAligner(alignerParameters, parameters.isInputPaired(), alignerParameters.getMergerParameters() != null);
List<VDJCGene> genes = new ArrayList<>();
for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary(parameters.library, parameters.species).getGenes(parameters.chains)) if (alignerParameters.containsRequiredFeature(gene) && (gene.isFunctional() || !parameters.isFunctionalOnly)) {
genes.add(gene);
aligner.addGene(gene);
}
AlignerReport report = new AlignerReport();
aligner.setEventsListener(report);
try (SequenceReaderCloseable<? extends SequenceRead> reader = parameters.getReader()) {
// start progress reporting
if (reader instanceof CanReportProgress)
SmartProgressReporter.startProgressReport("align", (CanReportProgress) reader);
OutputPort<Chunk<SequenceRead>> mainInputReads = CUtils.buffered((OutputPort) chunked(reader, 64), 16);
OutputPort<VDJCAlignmentResult> alignments = unchunked(new ParallelProcessor(mainInputReads, chunked(aligner), parameters.threads));
List<VDJCAlignments> als = new ArrayList<>();
int ind = 0;
for (VDJCAlignmentResult t : CUtils.it(new OrderedOutputPort<>(alignments, new Indexer<VDJCAlignmentResult>() {
@Override
public long getIndex(VDJCAlignmentResult r) {
return r.read.getId();
}
}))) {
if (t.alignment != null) {
t.alignment.setAlignmentsIndex(ind++);
als.add(t.alignment);
}
}
return new AlignResult(parameters, reader.getNumberOfReads(), report, als, genes, aligner);
}
}
Aggregations