use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters 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.VDJCAlignerParameters in project mixcr by milaboratory.
the class ActionAssembleContigs method go.
@Override
public void go(ActionHelper helper) throws Exception {
// TODO FIX!!!!!!!!!!!!!
if (parameters.threads > 1)
throw new ParameterException("Multithreaded processing is not supported yet.");
long beginTimestamp = System.currentTimeMillis();
FullSeqAssemblerParameters p = FullSeqAssemblerParameters.getByName("default");
if (!parameters.overrides.isEmpty()) {
// Perform parameters overriding
p = JsonOverrider.override(p, FullSeqAssemblerParameters.class, parameters.overrides);
if (p == null)
throw new ProcessException("Failed to override some parameter.");
}
final FullSeqAssemblerReport report = new FullSeqAssemblerReport();
FullSeqAssemblerParameters assemblerParameters = p;
int totalClonesCount = 0;
List<VDJCGene> genes;
VDJCAlignerParameters alignerParameters;
CloneAssemblerParameters cloneAssemblerParameters;
try (ClnAReader reader = new ClnAReader(parameters.getInputFileName(), VDJCLibraryRegistry.getDefault());
PrimitivO tmpOut = new PrimitivO(new BufferedOutputStream(new FileOutputStream(parameters.getOutputFileName())))) {
final CloneFactory cloneFactory = new CloneFactory(reader.getAssemblerParameters().getCloneFactoryParameters(), reader.getAssemblingFeatures(), reader.getGenes(), reader.getAlignerParameters().getFeaturesToAlignMap());
alignerParameters = reader.getAlignerParameters();
cloneAssemblerParameters = reader.getAssemblerParameters();
genes = reader.getGenes();
IOUtil.registerGeneReferences(tmpOut, genes, alignerParameters);
ClnAReader.CloneAlignmentsPort cloneAlignmentsPort = reader.clonesAndAlignments();
SmartProgressReporter.startProgressReport("Assembling", cloneAlignmentsPort);
OutputPort<Clone[]> parallelProcessor = new ParallelProcessor<>(cloneAlignmentsPort, cloneAlignments -> {
FullSeqAssembler fullSeqAssembler = new FullSeqAssembler(cloneFactory, assemblerParameters, cloneAlignments.clone, alignerParameters);
fullSeqAssembler.setReport(report);
FullSeqAssembler.RawVariantsData rawVariantsData = fullSeqAssembler.calculateRawData(() -> {
try {
return cloneAlignments.alignments();
} catch (IOException e) {
throw new RuntimeException(e);
}
});
return fullSeqAssembler.callVariants(rawVariantsData);
}, parameters.threads);
for (Clone[] clones : CUtils.it(parallelProcessor)) {
totalClonesCount += clones.length;
for (Clone cl : clones) tmpOut.writeObject(cl);
}
assert report.getInitialCloneCount() == reader.numberOfClones();
}
assert report.getFinalCloneCount() == totalClonesCount;
assert report.getFinalCloneCount() >= report.getInitialCloneCount();
Clone[] clones = new Clone[totalClonesCount];
try (PrimitivI tmpIn = new PrimitivI(new BufferedInputStream(new FileInputStream(parameters.getOutputFileName())))) {
IOUtil.registerGeneReferences(tmpIn, genes, alignerParameters);
int i = 0;
for (Clone clone : CUtils.it(new PipeDataInputReader<>(Clone.class, tmpIn, totalClonesCount))) clones[i++] = clone;
}
Arrays.sort(clones, Comparator.comparingDouble(c -> -c.getCount()));
for (int i = 0; i < clones.length; i++) clones[i] = clones[i].setId(i);
CloneSet cloneSet = new CloneSet(Arrays.asList(clones), genes, alignerParameters.getFeaturesToAlignMap(), alignerParameters, cloneAssemblerParameters);
try (CloneSetIO.CloneSetWriter writer = new CloneSetIO.CloneSetWriter(cloneSet, parameters.getOutputFileName())) {
SmartProgressReporter.startProgressReport(writer);
writer.write();
}
ReportWrapper reportWrapper = new ReportWrapper(command(), report);
reportWrapper.setStartMillis(beginTimestamp);
reportWrapper.setInputFiles(parameters.getInputFileName());
reportWrapper.setOutputFiles(parameters.getOutputFileName());
reportWrapper.setCommandLine(helper.getCommandLineArguments());
reportWrapper.setFinishMillis(System.currentTimeMillis());
// Writing report to stout
System.out.println("============= Report ==============");
Util.writeReportToStdout(report);
if (parameters.report != null)
Util.writeReport(parameters.report, reportWrapper);
if (parameters.jsonReport != null)
Util.writeJsonReport(parameters.jsonReport, reportWrapper);
}
use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.
the class CloneSetIO method readClns.
public static CloneSet readClns(InputStream inputStream, VDJCLibraryRegistry libraryRegistry) {
PrimitivI input = new PrimitivI(inputStream);
// Registering custom serializer
input.getSerializersManager().registerCustomSerializer(GeneFeature.class, new GeneFeatureSerializer(true));
byte[] magicBytes = new byte[MAGIC_LENGTH];
input.readFully(magicBytes);
String magicString = new String(magicBytes);
SerializersManager serializersManager = input.getSerializersManager();
switch(magicString) {
case MAGIC_V5:
serializersManager.registerCustomSerializer(Clone.class, new CompatibilityIO.CloneSerializerV5());
break;
case MAGIC:
break;
default:
throw new RuntimeException("Unsupported file format; .clns file of version " + magicString + " while you are running MiXCR " + MAGIC);
}
String versionInfo = input.readUTF();
if (!magicString.equals(MAGIC))
// Dropping this field for v5 files
input.readObject(GeneFeature[].class);
VDJCAlignerParameters alignerParameters;
CloneAssemblerParameters assemblerParameters;
if (magicString.equals(MAGIC)) {
alignerParameters = input.readObject(VDJCAlignerParameters.class);
assemblerParameters = input.readObject(CloneAssemblerParameters.class);
} else {
alignerParameters = null;
assemblerParameters = null;
}
EnumMap<GeneType, GeneFeature> alignedFeatures = IO.readGF2GTMap(input);
List<VDJCGene> genes = IOUtil.readAndRegisterGeneReferences(input, libraryRegistry, new GT2GFAdapter(alignedFeatures));
int count = input.readInt();
List<Clone> clones = new ArrayList<>(count);
for (int i = 0; i < count; i++) clones.add(input.readObject(Clone.class));
CloneSet cloneSet = new CloneSet(clones, genes, alignedFeatures, alignerParameters, assemblerParameters);
cloneSet.versionInfo = versionInfo;
return cloneSet;
}
use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.
the class VDJCObjectExtenderTest method testTripleRead.
@Test
public void testTripleRead() throws Exception {
final boolean print = true;
Well44497b rg = new Well44497b(12312);
final VDJCAlignerParameters rnaSeqParams = VDJCParametersPresets.getByName("rna-seq");
final PartialAlignmentsAssemblerAligner aligner = new PartialAlignmentsAssemblerAligner(rnaSeqParams);
final VDJCLibrary lib = VDJCLibraryRegistry.getDefault().getLibrary("default", "hs");
for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").getGenes()) if (gene.isFunctional())
aligner.addGene(gene);
TargetBuilder.VDJCGenes genes = new TargetBuilder.VDJCGenes(lib, "TRBV12-3*00", "TRBD1*00", "TRBJ1-3*00", "TRBC2*00");
// | 310 | 338 | 438
// 250V + 60CDR3 (20V 7N 10D 3N 20J) + 28J + 100C + 100N
final NucleotideSequence baseSeq = TargetBuilder.generateSequence(genes, "{CDR3Begin(-250)}V*270 NNNNNNN {DBegin(0)}D*10 NNN {CDR3End(-20):FR4End} {CBegin}C*100 N*100", rg);
F4 goAssert = new F4() {
@Override
public VDJCAlignments go(int len, int offset1, int offset2, int offset3) {
NucleotideSequence seq1 = baseSeq.getRange(offset1, offset1 + len);
NucleotideSequence seq2 = baseSeq.getRange(offset2, offset2 + len);
NucleotideSequence seq3 = offset3 == -1 ? null : baseSeq.getRange(offset3, offset3 + len);
VDJCAlignmentResult<VDJCMultiRead> alignment = offset3 == -1 ? aligner.process(MiXCRTestUtils.createMultiRead(seq1, seq2)) : aligner.process(MiXCRTestUtils.createMultiRead(seq1, seq2, seq3));
VDJCAlignments al = alignment.alignment;
Assert.assertNotNull(al);
VDJCObjectExtender<VDJCAlignments> extender = new VDJCObjectExtender<>(Chains.TCR, (byte) 35, rnaSeqParams.getVAlignerParameters().getParameters().getScoring(), rnaSeqParams.getJAlignerParameters().getParameters().getScoring(), 100, 70, ReferencePoint.CDR3Begin, ReferencePoint.CDR3End);
MiXCRTestUtils.assertAlignments(al);
if (print) {
MiXCRTestUtils.printAlignment(al);
System.out.println();
System.out.println("-------------------------------------------");
System.out.println();
}
VDJCAlignments processed = extender.process(al);
if (print) {
MiXCRTestUtils.printAlignment(processed);
System.out.println();
System.out.println("===========================================");
System.out.println();
System.out.println();
}
MiXCRTestUtils.assertAlignments(processed);
Assert.assertEquals(al.getFeature(GeneFeature.VDJunction), processed.getFeature(GeneFeature.VDJunction));
Assert.assertEquals(al.getFeature(GeneFeature.DJJunction), processed.getFeature(GeneFeature.DJJunction));
Assert.assertEquals(al.getFeature(GeneFeature.VJJunction), processed.getFeature(GeneFeature.VJJunction));
return processed;
}
};
VDJCAlignments a1 = goAssert.go(60, 245, 307, 450);
Assert.assertEquals(2, a1.numberOfTargets());
VDJCAlignments a2 = goAssert.go(60, 245, 315, 450);
Assert.assertEquals(3, a2.numberOfTargets());
VDJCAlignments a3 = goAssert.go(60, 245, 315, -1);
Assert.assertEquals(2, a3.numberOfTargets());
VDJCAlignments a4 = goAssert.go(60, 245, 307, -1);
Assert.assertEquals(1, a4.numberOfTargets());
VDJCAlignments a5 = goAssert.go(53, 252, 307, -1);
Assert.assertEquals(1, a5.numberOfTargets());
VDJCAlignments a6 = goAssert.go(53, 252, 307, 450);
Assert.assertEquals(2, a6.numberOfTargets());
}
use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.
the class PartialAlignmentsAssembler method searchOverlaps.
public void searchOverlaps(VDJCAlignmentsReader reader) {
final VDJCAlignerParameters alignerParameters = reader.getParameters();
PartialAlignmentsAssemblerAligner aligner = new PartialAlignmentsAssemblerAligner(alignerParameters);
targetMerger.setAlignerParameters(alignerParameters);
for (VDJCGene gene : reader.getUsedGenes()) aligner.addGene(gene);
for (final VDJCAlignments alignment : CUtils.it(reader)) {
total.incrementAndGet();
if (alignment.getFeature(GeneFeature.CDR3) != null) {
containsCDR3.incrementAndGet();
if (!overlappedOnly) {
totalWritten.incrementAndGet();
writer.write(alignment);
}
continue;
}
if (alreadyMergedIds.contains(alignment.getAlignmentsIndex()))
continue;
final OverlapSearchResult searchResult = searchOverlaps(alignment, alignerParameters.isAllowChimeras());
// Common procedure to cancel processing of current input alignment if it fails to pass some good
// overlap filtering criterion
final Runnable cancelCurrentResult = new Runnable() {
@Override
public void run() {
if (searchResult != null)
searchResult.cancel();
if (writePartial && !overlappedOnly && notInLeftIndexIds.contains(alignment.getAlignmentsIndex())) {
totalWritten.incrementAndGet();
partialAsIs.incrementAndGet();
writer.write(alignment);
}
}
};
if (searchResult == null) {
cancelCurrentResult.run();
continue;
}
List<AlignedTarget> mergedTargets = searchResult.result;
VDJCMultiRead mRead = new VDJCMultiRead(mergedTargets);
final VDJCAlignments mAlignment = aligner.process(mRead).alignment;
// Checking number of overlapped non-template (NRegion) letters
int overlapTargetId = -1;
Range overlapRange = null;
for (int i = 0; i < mergedTargets.size(); i++) {
overlapRange = AlignedTarget.getOverlapRange(mergedTargets.get(i));
if (overlapRange != null) {
overlapTargetId = i;
break;
}
}
if (overlapTargetId == -1) {
// No alignments for Best V Hit and Best J Hit in central (overlapped) target
cancelCurrentResult.run();
continue;
}
int targetLength = mergedTargets.get(overlapTargetId).getTarget().size();
VDJCHit bestVHit = mAlignment.getBestHit(GeneType.Variable), bestJHit = mAlignment.getBestHit(GeneType.Joining);
if (bestVHit == null || bestJHit == null || bestVHit.getAlignment(overlapTargetId) == null || bestJHit.getAlignment(overlapTargetId) == null) {
cancelCurrentResult.run();
continue;
}
int ndnRegionBegin = 0;
int ndnRegionEnd = targetLength;
Alignment<NucleotideSequence> vAlignment = bestVHit.getAlignment(overlapTargetId);
if (vAlignment != null)
ndnRegionBegin = vAlignment.getSequence2Range().getTo();
Alignment<NucleotideSequence> jAlignment = bestJHit.getAlignment(overlapTargetId);
if (jAlignment != null)
ndnRegionEnd = jAlignment.getSequence2Range().getFrom();
RangeSet nRegion = ndnRegionBegin >= ndnRegionEnd ? RangeSet.EMPTY : RangeSet.create(ndnRegionBegin, ndnRegionEnd);
Range dRange = mAlignment.getPartitionedTarget(overlapTargetId).getPartitioning().getRange(GeneFeature.DRegionTrimmed);
if (dRange != null)
nRegion = nRegion.subtract(dRange);
RangeSet nRegionInOverlap = nRegion.intersection(overlapRange);
int actualNRegionLength = nRegion.totalLength();
int minimalN = Math.min(minimalNOverlap, actualNRegionLength);
if (nRegionInOverlap.totalLength() < minimalN) {
droppedSmallOverlapNRegion.incrementAndGet();
cancelCurrentResult.run();
continue;
}
// DDDDDDDDDJJJJJJJJJJJJJJJ
if (minimalN == 0 && (!overlapRange.contains(ndnRegionBegin - 1) || !overlapRange.contains(ndnRegionEnd))) {
droppedNoNRegion.incrementAndGet();
cancelCurrentResult.run();
continue;
}
overlapped.incrementAndGet();
totalWritten.incrementAndGet();
writer.write(mAlignment);
// Saving alignment that where merge to prevent it's use as left part
alreadyMergedIds.add(alignment.getAlignmentsIndex());
alreadyMergedIds.add(searchResult.KMerInfo.alignments.getAlignmentsIndex());
}
if (writePartial && !overlappedOnly)
for (List<KMerInfo> kMerInfos : kToIndexLeft.valueCollection()) for (KMerInfo kMerInfo : kMerInfos) if (alreadyMergedIds.add(kMerInfo.alignments.getAlignmentsIndex())) {
totalWritten.incrementAndGet();
partialAsIs.incrementAndGet();
writer.write(kMerInfo.getAlignments());
}
writer.setNumberOfProcessedReads(reader.getNumberOfReads() - overlapped.get());
}
Aggregations