X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=34fe221a98d9099fe4e072f403a7bf12b61d580e;hb=37de9310bec3501cbc6381e0c3dcb282fcaad812;hp=cc803844163915dd2414396ef3c6c95d24606dc1;hpb=3410919969a9a70e54738d65489ab732abb9010c;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index cc80384..34fe221 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -72,22 +72,26 @@ import java.util.TreeMap; public class AlignmentUtils { + private static final int CODON_LENGTH = 3; + private static final String SEQUENCE_VARIANT = "sequence_variant:"; + private static final String ID = "ID"; /** * A data model to hold the 'normal' base value at a position, and an optional * sequence variant feature */ - static class DnaVariant + static final class DnaVariant { - String base; + final String base; SequenceFeature variant; DnaVariant(String nuc) { base = nuc; + variant = null; } DnaVariant(String nuc, SequenceFeature var) @@ -95,6 +99,11 @@ public class AlignmentUtils base = nuc; variant = var; } + + public String getSource() + { + return variant == null ? null : variant.getFeatureGroup(); + } } /** @@ -427,7 +436,7 @@ public class AlignmentUtils /* * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping) */ - final int mappedLength = 3 * aaSeqChars.length; + final int mappedLength = CODON_LENGTH * aaSeqChars.length; int cdnaLength = cdnaSeqChars.length; int cdnaStart = cdnaSeq.getStart(); int cdnaEnd = cdnaSeq.getEnd(); @@ -439,14 +448,14 @@ public class AlignmentUtils */ if (cdnaLength != mappedLength && cdnaLength > 2) { - String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3) - .toUpperCase(); + String lastCodon = String.valueOf(cdnaSeqChars, + cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(); for (String stop : ResidueProperties.STOP) { if (lastCodon.equals(stop)) { - cdnaEnd -= 3; - cdnaLength -= 3; + cdnaEnd -= CODON_LENGTH; + cdnaLength -= CODON_LENGTH; break; } } @@ -458,12 +467,12 @@ public class AlignmentUtils int startOffset = 0; if (cdnaLength != mappedLength && cdnaLength > 2 - && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase() + && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase() .equals(ResidueProperties.START)) { - startOffset += 3; - cdnaStart += 3; - cdnaLength -= 3; + startOffset += CODON_LENGTH; + cdnaStart += CODON_LENGTH; + cdnaLength -= CODON_LENGTH; } if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars)) @@ -472,7 +481,7 @@ public class AlignmentUtils * protein is translation of dna (+/- start/stop codons) */ MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] - { proteinStart, proteinEnd }, 3, 1); + { proteinStart, proteinEnd }, CODON_LENGTH, 1); return map; } @@ -502,10 +511,9 @@ public class AlignmentUtils int aaPos = 0; int dnaPos = cdnaStart; - for (; dnaPos < cdnaSeqChars.length - 2 - && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++) + for (; dnaPos < cdnaSeqChars.length - 2 && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++) { - String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); final String translated = ResidueProperties.codonTranslate(codon); /* @@ -541,9 +549,9 @@ public class AlignmentUtils { return true; } - if (dnaPos == cdnaSeqChars.length - 3) + if (dnaPos == cdnaSeqChars.length - CODON_LENGTH) { - String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); if ("STOP".equals(ResidueProperties.codonTranslate(codon))) { return true; @@ -894,7 +902,8 @@ public class AlignmentUtils } width = Math.max(dnaSeq.getLength(), width); } - int oldwidth, diff; + int oldwidth; + int diff; for (SequenceI dnaSeq : dna.getSequences()) { oldwidth = dnaSeq.getLength(); @@ -928,15 +937,15 @@ public class AlignmentUtils .println("alignCdsSequenceAsProtein needs aligned sequence!"); return false; } - + List dnaMappings = MappingUtils .findMappingsForSequence(cdsSeq, mappings); for (AlignedCodonFrame mapping : dnaMappings) { SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); - int peptideLength = peptide.getLength(); if (peptide != null) { + int peptideLength = peptide.getLength(); Mapping map = mapping.getMappingBetween(cdsSeq, peptide); if (map != null) { @@ -950,7 +959,8 @@ public class AlignmentUtils .getFromRanges()); int mappedToLength = MappingUtils .getLength(mapList.getToRanges()); - boolean addStopCodon = (cdsLength == mappedFromLength * 3 + 3) + boolean addStopCodon = (cdsLength == mappedFromLength + * CODON_LENGTH + CODON_LENGTH) || (peptide.getDatasetSequence().getLength() == mappedFromLength - 1); if (cdsLength != mappedToLength && !addStopCodon) { @@ -964,8 +974,8 @@ public class AlignmentUtils /* * pre-fill the aligned cds sequence with gaps */ - char[] alignedCds = new char[peptideLength * 3 - + (addStopCodon ? 3 : 0)]; + char[] alignedCds = new char[peptideLength * CODON_LENGTH + + (addStopCodon ? CODON_LENGTH : 0)]; Arrays.fill(alignedCds, gapChar); /* @@ -982,7 +992,7 @@ public class AlignmentUtils { if (Comparison.isGap(residue)) { - cdsCol += 3; + cdsCol += CODON_LENGTH; } else { @@ -991,7 +1001,7 @@ public class AlignmentUtils if (codon == null) { // e.g. incomplete start codon, X in peptide - cdsCol += 3; + cdsCol += CODON_LENGTH; } else { @@ -1009,7 +1019,7 @@ public class AlignmentUtils * append stop codon if not mapped from protein, * closing it up to the end of the mapped sequence */ - if (copiedBases == nucleotides.length - 3) + if (copiedBases == nucleotides.length - CODON_LENGTH) { for (int i = alignedCds.length - 1; i >= 0; i--) { @@ -1019,7 +1029,7 @@ public class AlignmentUtils break; } } - for (int i = nucleotides.length - 3; i < nucleotides.length; i++) + for (int i = nucleotides.length - CODON_LENGTH; i < nucleotides.length; i++) { alignedCds[cdsCol++] = nucleotides[i]; } @@ -1089,7 +1099,7 @@ public class AlignmentUtils // TODO resolve JAL-2022 so this fudge can be removed int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); - + return alignedCodons; } @@ -1680,15 +1690,20 @@ public class AlignmentUtils * didn't find mapped CDS sequence - construct it and add * its dataset sequence to the dataset */ - cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping); + cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping, + dataset).deriveSequence(); // cdsSeq has a name constructed as CDS| // will be either the accession for the coding sequence, // marked in the /via/ dbref to the protein product accession // or it will be the original nucleotide accession. - SequenceI cdsSeqDss = cdsSeq.createDatasetSequence(); + SequenceI cdsSeqDss = cdsSeq.getDatasetSequence(); + cdsSeqs.add(cdsSeq); + if (!dataset.getSequences().contains(cdsSeqDss)) { + // check if this sequence is a newly created one + // so needs adding to the dataset dataset.addSequence(cdsSeqDss); } @@ -1697,8 +1712,9 @@ public class AlignmentUtils */ List cdsRange = Collections.singletonList(new int[] { 1, cdsSeq.getLength() }); - MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(), - mapList.getFromRatio(), mapList.getToRatio()); + MapList cdsToProteinMap = new MapList(cdsRange, + mapList.getToRanges(), mapList.getFromRatio(), + mapList.getToRatio()); AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct, cdsToProteinMap); @@ -1718,8 +1734,7 @@ public class AlignmentUtils */ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), - cdsRange, 1, - 1); + cdsRange, 1, 1); dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, dnaToCdsMap); if (!mappings.contains(dnaToCdsMapping)) @@ -1734,16 +1749,16 @@ public class AlignmentUtils * same source and accession, so need a different accession for * the CDS from the dna sequence */ - + // specific use case: // Genomic contig ENSCHR:1, contains coding regions for ENSG01, // ENSG02, ENSG03, with transcripts and products similarly named. // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01 - + // JBPNote: ?? can't actually create an example that demonstrates we // need to // synthesize an xref. - + for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs()) { // creates a complementary cross-reference to the source sequence's @@ -1773,7 +1788,7 @@ public class AlignmentUtils /* * transfer any features on dna that overlap the CDS */ - transferFeatures(dnaSeq, cdsSeq, cdsToProteinMap, null, + transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null, SequenceOntologyI.CDS); } } @@ -1820,7 +1835,8 @@ public class AlignmentUtils int mappedFromLength = MappingUtils.getLength(aMapping.getMap() .getFromRanges()); int dnaLength = seqDss.getLength(); - if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - 3) + if (mappedFromLength == dnaLength + || mappedFromLength == dnaLength - CODON_LENGTH) { return seqDss; } @@ -1836,7 +1852,8 @@ public class AlignmentUtils for (SequenceToSequenceMapping map : acf.getMappings()) { Mapping mapping = map.getMapping(); - if (mapping != aMapping && mapping.getMap().getFromRatio() == 3 + if (mapping != aMapping + && mapping.getMap().getFromRatio() == CODON_LENGTH && proteinProduct == mapping.getTo() && seqDss != map.getFromSeq()) { @@ -1870,9 +1887,14 @@ public class AlignmentUtils * * @param seq * @param mapping + * @param dataset + * - existing dataset. We check for sequences that look like the CDS + * we are about to construct, if one exists already, then we will + * just return that one. * @return CDS sequence (as a dataset sequence) */ - static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) + static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, + AlignmentI dataset) { char[] seqChars = seq.getSequence(); List fromRanges = mapping.getMap().getFromRanges(); @@ -1899,7 +1921,7 @@ public class AlignmentUtils } } } - + /* * assign 'from id' held in the mapping if set (e.g. EMBL protein_id), * else generate a sequence name @@ -1907,6 +1929,40 @@ public class AlignmentUtils String mapFromId = mapping.getMappedFromId(); String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName()); SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos); + if (dataset != null) + { + SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName()); + if (matches != null) + { + boolean matched = false; + for (SequenceI mtch : matches) + { + if (mtch.getStart() != newSeq.getStart()) + { + continue; + } + if (mtch.getEnd() != newSeq.getEnd()) + { + continue; + } + if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence())) + { + continue; + } + if (!matched) + { + matched = true; + newSeq = mtch; + } + else + { + System.err + .println("JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + + mtch.toString()); + } + } + } + } // newSeq.setDescription(mapFromId); return newSeq; @@ -2119,7 +2175,7 @@ public class AlignmentUtils /* * dna length should map to protein (or protein plus stop codon) */ - int codesForResidues = mappedDnaLength / 3; + int codesForResidues = mappedDnaLength / CODON_LENGTH; if (codesForResidues == (proteinLength + 1)) { // assuming extra codon is for STOP and not in peptide @@ -2128,7 +2184,7 @@ public class AlignmentUtils if (codesForResidues == proteinLength) { proteinRange.add(new int[] { proteinStart, proteinEnd }); - return new MapList(ranges, proteinRange, 3, 1); + return new MapList(ranges, proteinRange, CODON_LENGTH, 1); } return null; } @@ -2404,7 +2460,7 @@ public class AlignmentUtils * are currently ignored here */ String trans = codon.contains("-") ? "-" - : (codon.length() > 3 ? null : ResidueProperties + : (codon.length() > CODON_LENGTH ? null : ResidueProperties .codonTranslate(codon)); if (trans != null && !trans.equals(residue)) { @@ -2416,7 +2472,7 @@ public class AlignmentUtils // set score to 0f so 'graduated colour' option is offered! JAL-2060 SequenceFeature sf = new SequenceFeature( SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, - peptidePos, 0f, "Jalview"); + peptidePos, 0f, var.getSource()); StringBuilder attributes = new StringBuilder(32); String id = (String) var.variant.getValue(ID); if (id != null) @@ -2427,11 +2483,13 @@ public class AlignmentUtils } sf.setValue(ID, id); attributes.append(ID).append("=").append(id); - // TODO handle other species variants + // TODO handle other species variants JAL-2064 StringBuilder link = new StringBuilder(32); try { - link.append(desc).append(" ").append(id) + link.append(desc) + .append(" ") + .append(id) .append("|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=") .append(URLEncoder.encode(id, "UTF-8")); sf.addLink(link.toString()); @@ -2440,8 +2498,7 @@ public class AlignmentUtils // as if } } - String clinSig = (String) var.variant - .getValue(CLINICAL_SIGNIFICANCE); + String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE); if (clinSig != null) { sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); @@ -2466,6 +2523,7 @@ public class AlignmentUtils * @param dnaToProtein * @return */ + @SuppressWarnings("unchecked") static LinkedHashMap[]> buildDnaVariantsMap( SequenceI dnaSeq, MapList dnaToProtein) { @@ -2509,7 +2567,7 @@ public class AlignmentUtils List[] codonVariants = variants.get(peptidePosition); if (codonVariants == null) { - codonVariants = new ArrayList[3]; + codonVariants = new ArrayList[CODON_LENGTH]; codonVariants[0] = new ArrayList(); codonVariants[1] = new ArrayList(); codonVariants[2] = new ArrayList(); @@ -2543,7 +2601,7 @@ public class AlignmentUtils /* * save nucleotide (and any variant) for each codon position */ - for (int codonPos = 0; codonPos < 3; codonPos++) + for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++) { String nucleotide = String.valueOf( dnaSeq.getCharAt(codon[codonPos] - dnaStart)) @@ -2596,7 +2654,7 @@ public class AlignmentUtils { AlignmentI copy = new Alignment(new Alignment(seqs)); copy.setDataset(dataset); - + boolean isProtein = !copy.isNucleotide(); SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) { @@ -2607,7 +2665,8 @@ public class AlignmentUtils { for (DBRefEntry dbref : dbrefs) { - if (dbref.getMap() == null || dbref.getMap().getTo() == null) + if (dbref.getMap() == null || dbref.getMap().getTo() == null + || dbref.getMap().getTo().isProtein() != isProtein) { continue; } @@ -2688,7 +2747,7 @@ public class AlignmentUtils } newCol++; } - + /* * trim trailing gaps */