use of com.milaboratory.core.sequence.NucleotideSequence in project repseqio by repseqio.
the class GenerateClonesAction method go.
@Override
public void go(ActionHelper helper) throws Exception {
GCloneModel model = GModels.getGCloneModelByName(params.getModelName());
GCloneGenerator generator = model.create(new Well19937c(params.getSeed()), VDJCLibraryRegistry.getDefault());
VDJCLibrary library = VDJCLibraryRegistry.getDefault().getLibrary(model.libraryId());
try (BufferedOutputStream s = new BufferedOutputStream(params.getOutput().equals(".") ? System.out : new FileOutputStream(params.getOutput()), 128 * 1024)) {
s.write(GlobalObjectMappers.toOneLine(model.libraryId()).getBytes());
s.write('\n');
ObjectWriter writer = GlobalObjectMappers.ONE_LINE.writerFor(new TypeReference<GClone>() {
}).withAttribute(VDJCGene.JSON_CURRENT_LIBRARY_ATTRIBUTE_KEY, library);
OUTER: for (int i = 0; i < params.numberOfClones; i++) {
GClone clone = generator.sample();
for (GGene g : clone.genes.values()) {
NucleotideSequence cdr3 = g.getFeature(GeneFeature.CDR3);
if (params.isInFrame())
if (cdr3.size() % 3 != 0) {
--i;
continue OUTER;
}
if (params.isNoStops())
if (AminoAcidSequence.translateFromCenter(cdr3).containStops()) {
--i;
continue OUTER;
}
}
writer.writeValue(new CloseShieldOutputStream(s), clone);
s.write('\n');
}
}
}
use of com.milaboratory.core.sequence.NucleotideSequence in project repseqio by repseqio.
the class InferAnchorPointsAction method go.
public void go(GeneFeature geneFeature, String inputFile, String outputFile) throws Exception {
VDJCLibraryRegistry.resetDefaultRegistry();
VDJCLibraryRegistry reg = VDJCLibraryRegistry.getDefault();
// Map of library names
Map<String, String> libraryNameToAddress = new HashMap<>();
// Registering reference library
int i = 0;
for (String refAddress : params.getReference()) {
String name = REFERENCE_LIBRARY_PREFIX + (i++);
reg.registerLibraries(refAddress, name);
libraryNameToAddress.put(name, refAddress);
}
if (params.getReference().isEmpty()) {
reg.loadAllLibraries("default");
for (VDJCLibrary library : reg.getLoadedLibraries()) libraryNameToAddress.put(library.getName(), "built-in");
}
// Registering target library
reg.registerLibraries(inputFile, TARGET_LIBRARY_NAME);
// Compile gene filter
Pattern namePattern = params.name == null ? null : Pattern.compile(params.name);
List<VDJCLibrary> refLibraries = reg.getLoadedLibrariesByNamePattern(REFERENCE_LIBRARY_PATTERN);
SimpleBatchAlignerParameters<AminoAcidSequence> aParams = new SimpleBatchAlignerParameters<>(5, 0.4f, params.getAbsoluteMinScore(geneFeature), true, AffineGapAlignmentScoring.getAminoAcidBLASTScoring(BLASTMatrix.BLOSUM62, -10, -1));
SimpleBatchAligner<AminoAcidSequence, Ref> aligner = new SimpleBatchAligner<>(aParams);
int dbSize = 0;
for (VDJCLibrary lib : refLibraries) {
for (VDJCGene gene : lib.getGenes()) {
NucleotideSequence nSeq = gene.getFeature(geneFeature);
if (nSeq == null)
continue;
ReferencePoint frameReference = GeneFeature.getFrameReference(geneFeature);
ReferencePoints partitioning = gene.getPartitioning();
if (frameReference == null)
continue;
int relativePosition = partitioning.getRelativePosition(geneFeature, frameReference);
if (relativePosition < 0)
continue;
TranslationParameters frame = withIncompleteCodon(relativePosition);
AminoAcidSequence aaSequence = AminoAcidSequence.translate(nSeq, frame);
aligner.addReference(aaSequence, new Ref(gene, frame, nSeq.size()));
++dbSize;
}
}
System.out.println("DB size: " + dbSize);
System.out.println();
// Checking that db is not empty
if (dbSize == 0)
throw new RuntimeException("No reference genes.");
ArrayList<VDJCLibraryData> result = new ArrayList<>();
ByteArrayOutputStream bos = new ByteArrayOutputStream();
PrintStream bufferPS = new PrintStream(bos);
// Iteration over target genes
for (VDJCLibrary lib : reg.getLoadedLibrariesByName(TARGET_LIBRARY_NAME)) {
ArrayList<VDJCGeneData> genes = new ArrayList<>();
for (VDJCGene targetGene : lib.getGenes()) {
bos.reset();
PrintStream ps = params.outputOnlyModified() ? bufferPS : System.out;
if (namePattern != null && !namePattern.matcher(targetGene.getName()).matches()) {
if (!params.outputOnlyModified())
genes.add(targetGene.getData());
continue;
}
ps.println("Processing: " + targetGene.getName() + " (" + (targetGene.isFunctional() ? "F" : "P") + ") " + targetGene.getChains());
// Getting gene feature sequence from target gene
NucleotideSequence nSeq = targetGene.getFeature(geneFeature);
if (nSeq == null) {
ps.println("Failed to extract " + GeneFeature.encode(geneFeature));
ps.println("================");
ps.println();
if (!params.outputOnlyModified())
genes.add(targetGene.getData());
continue;
}
// Alignment result
AlignmentResult<AlignmentHit<AminoAcidSequence, Ref>> bestAResult = null;
TranslationParameters bestFrame = null;
// Searching for best alignment
for (TranslationParameters frame : TRANSLATION_PARAMETERS) {
AminoAcidSequence aaSeq = AminoAcidSequence.translate(nSeq, frame);
AlignmentResult<AlignmentHit<AminoAcidSequence, Ref>> r = aligner.align(aaSeq);
if (r != null && r.hasHits() && (bestAResult == null || bestAResult.getBestHit().getAlignment().getScore() < r.getBestHit().getAlignment().getScore())) {
bestAResult = r;
bestFrame = frame;
}
}
if (bestFrame == null) {
ps.println("No alignments found.");
if (!params.outputOnlyModified())
genes.add(targetGene.getData());
continue;
}
List<AlignmentHit<AminoAcidSequence, Ref>> hits = bestAResult.getHits();
VDJCGeneData targetGeneData = targetGene.getData().clone();
boolean anyPointChanged = false;
for (int ai = 0; ai < hits.size(); ai++) {
// Accumulate output
ByteArrayOutputStream localBos = new ByteArrayOutputStream();
PrintStream localPS = new PrintStream(localBos);
Alignment<AminoAcidSequence> bestAlignment = hits.get(ai).getAlignment();
Ref bestRef = hits.get(ai).getRecordPayload();
VDJCGene bestReferenceGene = bestRef.gene;
localPS.println("Aligned with " + bestReferenceGene.getName() + " from " + libraryNameToAddress.get(bestReferenceGene.getParentLibrary().getName()) + " ; Score = " + bestAlignment.getScore());
AlignmentHelper alignmentHelper = bestAlignment.getAlignmentHelper();
for (AlignmentHelper h : alignmentHelper.split(150)) localPS.println(h + "\n");
ReferencePoints targetPartitioning = targetGene.getPartitioning();
ReferencePoints referencePartitioning = bestReferenceGene.getPartitioning();
for (GeneFeature.ReferenceRange range : geneFeature) for (ReferencePoint point : range.getIntermediatePoints()) {
localPS.print(point + ": ");
boolean isAvailable = targetPartitioning.isAvailable(point);
if (isAvailable) {
localPS.println("already set");
}
if (!referencePartitioning.isAvailable(point)) {
if (!isAvailable)
localPS.println("not set in reference gene");
continue;
}
int ntPositionInReference = referencePartitioning.getRelativePosition(geneFeature, point);
// Projecting position
AminoAcidSequence.AminoAcidSequencePosition aaPositionInReferece = AminoAcidSequence.convertNtPositionToAA(ntPositionInReference, bestRef.ntSeqLength, bestRef.frame);
if (aaPositionInReferece == null) {
if (!isAvailable)
localPS.println("failed to convert to aa position in ref");
continue;
}
int aaPositionInTarget = Alignment.aabs(bestAlignment.convertToSeq2Position(aaPositionInReferece.aminoAcidPosition));
if (aaPositionInTarget == -1) {
if (!isAvailable)
localPS.println("failed to project using alignment");
continue;
}
int ntPositionInTarget = AminoAcidSequence.convertAAPositionToNt(aaPositionInTarget, nSeq.size(), bestFrame);
if (ntPositionInTarget == -1) {
if (!isAvailable)
localPS.println("failed");
continue;
}
ntPositionInTarget += aaPositionInReferece.positionInTriplet;
ntPositionInTarget = targetPartitioning.getAbsolutePosition(geneFeature, ntPositionInTarget);
if (ntPositionInTarget == -1) {
if (!isAvailable)
localPS.println("failed");
continue;
}
if (isAvailable) {
int existingPosition = targetPartitioning.getPosition(point);
if (existingPosition != ntPositionInTarget) {
localPS.println("inferred position differs from existing value. existing: " + existingPosition + ", inferred: " + ntPositionInTarget);
}
continue;
}
localPS.println(ntPositionInTarget);
targetGeneData.getAnchorPoints().put(point, (long) ntPositionInTarget);
anyPointChanged = true;
}
if (!anyPointChanged) {
if (!params.outputOnlyModified() && ai == 0)
ps.write(localBos.toByteArray());
} else {
ps.write(localBos.toByteArray());
break;
}
}
ps.println("================");
ps.println();
if (anyPointChanged && params.outputOnlyModified())
System.out.write(bos.toByteArray());
genes.add(targetGeneData);
}
result.add(new VDJCLibraryData(lib.getData(), genes));
}
VDJCDataUtils.writeToFile(result, outputFile, false);
}
use of com.milaboratory.core.sequence.NucleotideSequence in project repseqio by repseqio.
the class TsvAction method go.
@Override
public void go(ActionHelper helper) throws Exception {
VDJCLibraryRegistry reg = VDJCLibraryRegistry.getDefault();
if (!"default".equals(params.getInput()))
reg.registerLibraries(params.getInput());
else
reg.loadAllLibraries("default");
Pattern chainPattern = params.chain == null ? null : Pattern.compile(params.chain);
Pattern namePattern = params.name == null ? null : Pattern.compile(params.name);
Long taxonFilter = params.taxonId;
if (taxonFilter == null && params.species != null)
taxonFilter = reg.resolveSpecies(params.species);
try (BufferedWriter writer = new BufferedWriter(new OutputStreamWriter(params.getOutputStream(), StandardCharsets.UTF_8))) {
writer.write("Gene\tChains\tFeature\tStart\tStop\tSource\tSequence\n");
for (VDJCLibrary lib : reg.getLoadedLibraries()) {
if (taxonFilter != null && taxonFilter != lib.getTaxonId())
continue;
for (VDJCGene gene : lib.getGenes()) {
if (chainPattern != null) {
boolean y = false;
for (String s : gene.getChains()) if (y |= chainPattern.matcher(s).matches())
break;
if (!y)
continue;
}
if (namePattern != null && !namePattern.matcher(gene.getName()).matches())
continue;
for (GeneFeatureWithOriginalName feature : params.features) {
GeneFeature geneFeature = feature.feature;
NucleotideSequence featureSequence = gene.getFeature(geneFeature);
if (featureSequence == null)
continue;
// Don't output start and end positions for composite gene features
Long start = geneFeature.isComposite() ? null : gene.getData().getAnchorPoints().get(geneFeature.getFirstPoint());
Long end = geneFeature.isComposite() ? null : gene.getData().getAnchorPoints().get(geneFeature.getLastPoint());
NucleotideSequence nSequence = gene.getFeature(geneFeature);
List<String> tokens = Arrays.asList(gene.getGeneName(), gene.getChains().toString(), feature.originalName, // (so essentially 1-based inclusive). Report both as 1-based.
(start == null ? "" : params.isOneBased() ? String.valueOf(start + 1) : String.valueOf(start)), (end == null ? "" : String.valueOf(end)), gene.getData().getBaseSequence().getOrigin().toString(), nSequence.toString());
String delim = "";
for (String t : tokens) {
writer.write(delim);
writer.write(t);
delim = "\t";
}
writer.write('\n');
}
}
}
}
}
use of com.milaboratory.core.sequence.NucleotideSequence in project repseqio by repseqio.
the class DebugAction method go.
@Override
public void go(ActionHelper helper) throws Exception {
VDJCLibraryRegistry reg = VDJCLibraryRegistry.getDefault();
reg.registerLibraries(params.getInput());
Pattern namePattern = params.name == null ? null : Pattern.compile(params.name);
GeneFeature cdr3FirstTriplet = new GeneFeature(ReferencePoint.CDR3Begin, 0, 3);
GeneFeature cdr3LastTriplet = new GeneFeature(ReferencePoint.CDR3End, -3, 0);
GeneFeature vIntronDonor = new GeneFeature(ReferencePoint.VIntronBegin, 0, 2);
GeneFeature vIntronAcceptor = new GeneFeature(ReferencePoint.VIntronEnd, -2, 0);
for (VDJCLibrary lib : reg.getLoadedLibraries()) {
for (VDJCGene gene : lib.getGenes()) {
if (namePattern != null && !namePattern.matcher(gene.getName()).matches())
continue;
// First generate list of warning messages
List<String> warnings = new ArrayList<>();
if (gene.isFunctional() || params.getCheckAll()) {
NucleotideSequence l3;
switch(gene.getGeneType()) {
case Variable:
// Flag AA residues flanking CDR3
l3 = gene.getFeature(cdr3FirstTriplet);
if (l3 == null)
warnings.add("Unable to find CDR3 start");
else if (l3.size() != 3)
warnings.add("Unable to translate sequence: " + gene.getName() + " / " + l3);
else if (AminoAcidSequence.translate(l3).codeAt(0) != AminoAcidAlphabet.C)
warnings.add("CDR3 does not start with C, was: " + l3.toString() + " / " + AminoAcidSequence.translate(l3).toString() + " / CDR3Begin: " + gene.getData().getAnchorPoints().get(ReferencePoint.CDR3Begin));
// Flag suspicious exon borders
// https://schneider.ncifcrf.gov/gallery/SequenceLogoSculpture.gif
NucleotideSequence vIntronDonorSeq = gene.getFeature(vIntronDonor);
if (vIntronDonorSeq != null && !vIntronDonorSeq.toString().equals("GT") && !vIntronDonorSeq.toString().equals("GC"))
warnings.add("Expected VIntron sequence to start with GT, was: " + vIntronDonorSeq.toString());
NucleotideSequence vIntronAcceptorSeq = gene.getFeature(vIntronAcceptor);
if (vIntronAcceptorSeq != null && !vIntronAcceptorSeq.toString().equals("AG"))
warnings.add("Expected VIntron sequence to end with AG, was: " + vIntronAcceptorSeq.toString());
ReferencePoints partitioning = gene.getPartitioning();
if (partitioning.isAvailable(GeneFeature.VTranscriptWithout5UTR)) {
// Iterating over all reading-frame bound anchor points of V gene
for (ReferencePoint anchorPoint : ReferencePoint.DefaultReferencePoints) {
if (anchorPoint.getGeneType() != GeneType.Variable || !anchorPoint.isTripletBoundary())
continue;
// And checking that they are in the same frame as Start (L1Begin)
int relativePosition = partitioning.getRelativePosition(GeneFeature.VTranscriptWithout5UTR, anchorPoint);
if (// Point is defined
relativePosition >= 0 && relativePosition % 3 != 0)
warnings.add("Expected " + anchorPoint + " to have position dividable by three inside VTranscriptWithout5UTR. " + "This may indicate an error in the L2 boundaries.");
}
}
break;
case Joining:
// Flag AA residues flanking CDR3
l3 = gene.getFeature(cdr3LastTriplet);
if (l3 == null)
warnings.add("Unable to find CDR3 end");
else if (l3.size() != 3)
warnings.add("Unable to translate sequence: " + gene.getName() + " / " + l3);
else if (AminoAcidSequence.translate(l3).codeAt(0) != AminoAcidAlphabet.W && AminoAcidSequence.translate(l3).codeAt(0) != AminoAcidAlphabet.F)
warnings.add("CDR3 does not end with W or F, was: " + l3.toString() + " / " + AminoAcidSequence.translate(l3).toString() + " / CDR3End: " + gene.getData().getAnchorPoints().get(ReferencePoint.CDR3End));
break;
}
for (GeneFeature geneFeature : aaGeneFeatures.get(gene.getGeneType())) {
AminoAcidSequence aaSequence = getAminoAcidSequence(gene, geneFeature, gene.getFeature(geneFeature));
if (aaSequence != null) {
// Flag if contains stop codon
if (aaSequence.numberOfStops() > 0)
warnings.add(GeneFeature.encode(geneFeature) + " contains a stop codon");
}
}
}
if (params.getProblemOnly() && warnings.isEmpty())
continue;
System.out.println(gene.getName() + " (" + (gene.isFunctional() ? "F" : "P") + ") " + gene.getChains() + " : " + lib.getTaxonId());
if (!warnings.isEmpty()) {
System.out.println();
System.out.println("WARNINGS: ");
for (String warning : warnings) {
System.out.println(warning);
}
System.out.println();
}
for (GeneFeature geneFeature : geneFeatures.get(gene.getGeneType())) {
System.out.println();
System.out.println(GeneFeature.encode(geneFeature));
NucleotideSequence nSequence = gene.getFeature(geneFeature);
AminoAcidSequence aaSequence = getAminoAcidSequence(gene, geneFeature, nSequence);
System.out.print("N: ");
if (nSequence == null)
System.out.println("Not Available");
else
System.out.println(nSequence);
if (GeneFeature.getFrameReference(geneFeature) != null) {
System.out.print("AA: ");
if (aaSequence == null)
System.out.println("Not Available");
else
System.out.println(aaSequence);
}
}
System.out.println("=========");
System.out.println();
}
}
}
use of com.milaboratory.core.sequence.NucleotideSequence in project repseqio by repseqio.
the class FastaAction method go.
@Override
public void go(ActionHelper helper) throws Exception {
GeneFeature geneFeature = params.getGeneFeature();
VDJCLibraryRegistry reg = VDJCLibraryRegistry.getDefault();
if (!"default".equals(params.getInput()))
reg.registerLibraries(params.getInput());
else
reg.loadAllLibraries("default");
Pattern chainPattern = params.chain == null ? null : Pattern.compile(params.chain);
Pattern namePattern = params.name == null ? null : Pattern.compile(params.name);
Long taxonFilter = params.taxonId;
if (taxonFilter == null && params.species != null)
taxonFilter = reg.resolveSpecies(params.species);
try (FastaWriter<NucleotideSequence> writer = CLIUtils.createSingleFastaWriter(params.getOutput())) {
for (VDJCLibrary lib : reg.getLoadedLibraries()) {
if (taxonFilter != null && taxonFilter != lib.getTaxonId())
continue;
for (VDJCGene gene : lib.getGenes()) {
if (chainPattern != null) {
boolean y = false;
for (String s : gene.getChains()) if (y |= chainPattern.matcher(s).matches())
break;
if (!y)
continue;
}
if (namePattern != null && !namePattern.matcher(gene.getName()).matches())
continue;
NucleotideSequence featureSequence = gene.getFeature(geneFeature);
if (featureSequence == null)
continue;
writer.write(gene.getName() + "|" + (gene.isFunctional() ? "F" : "P") + "|taxonId=" + gene.getParentLibrary().getTaxonId(), featureSequence);
}
}
}
}
Aggregations