X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=db69823050fc191d5ebddb0cbd3baeba0a04e841;hb=9fe87bd08ab29b6fe53364f387d7a2e3a8d39994;hp=6f0125d97eb5e287aae5f90fc8d9a359375bd1d1;hpb=be32c14cd8e48fe0a207cd7030cb9cd46f894678;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 6f0125d..db69823 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,9 +20,35 @@ */ package jalview.analysis; +import jalview.datamodel.AlignedCodon; +import jalview.datamodel.AlignedCodonFrame; +import jalview.datamodel.Alignment; +import jalview.datamodel.AlignmentAnnotation; +import jalview.datamodel.AlignmentI; +import jalview.datamodel.DBRefEntry; +import jalview.datamodel.DBRefSource; +import jalview.datamodel.FeatureProperties; +import jalview.datamodel.IncompleteCodonException; +import jalview.datamodel.Mapping; +import jalview.datamodel.SearchResults; +import jalview.datamodel.Sequence; +import jalview.datamodel.SequenceFeature; +import jalview.datamodel.SequenceGroup; +import jalview.datamodel.SequenceI; +import jalview.io.gff.SequenceOntologyFactory; +import jalview.io.gff.SequenceOntologyI; +import jalview.schemes.ResidueProperties; +import jalview.util.Comparison; +import jalview.util.DBRefUtils; +import jalview.util.MapList; +import jalview.util.MappingUtils; +import jalview.util.StringUtils; + import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; +import java.util.Collections; +import java.util.Comparator; import java.util.HashMap; import java.util.HashSet; import java.util.Iterator; @@ -30,21 +56,10 @@ import java.util.LinkedHashMap; import java.util.List; import java.util.Map; import java.util.Map.Entry; +import java.util.NoSuchElementException; import java.util.Set; import java.util.TreeMap; -import jalview.datamodel.AlignedCodon; -import jalview.datamodel.AlignedCodonFrame; -import jalview.datamodel.AlignmentAnnotation; -import jalview.datamodel.AlignmentI; -import jalview.datamodel.Mapping; -import jalview.datamodel.SearchResults; -import jalview.datamodel.Sequence; -import jalview.datamodel.SequenceGroup; -import jalview.datamodel.SequenceI; -import jalview.schemes.ResidueProperties; -import jalview.util.MapList; - /** * grab bag of useful alignment manipulation operations Expect these to be * refactored elsewhere at some point. @@ -70,18 +85,22 @@ public class AlignmentUtils for (SequenceI s : core.getSequences()) { SequenceI newSeq = s.deriveSequence(); - if (newSeq.getStart() > maxoffset + final int newSeqStart = newSeq.getStart() - 1; + if (newSeqStart > maxoffset && newSeq.getDatasetSequence().getStart() < s.getStart()) { - maxoffset = newSeq.getStart(); + maxoffset = newSeqStart; } sq.add(newSeq); } if (flankSize > -1) { - maxoffset = flankSize; + maxoffset = Math.min(maxoffset, flankSize); } - // now add offset to create a new expanded alignment + + /* + * now add offset left and right to create an expanded alignment + */ for (SequenceI s : sq) { SequenceI ds = s; @@ -91,8 +110,8 @@ public class AlignmentUtils } int s_end = s.findPosition(s.getStart() + s.getLength()); // find available flanking residues for sequence - int ustream_ds = s.getStart() - ds.getStart(), dstream_ds = ds - .getEnd() - s_end; + int ustream_ds = s.getStart() - ds.getStart(); + int dstream_ds = ds.getEnd() - s_end; // build new flanked sequence @@ -108,27 +127,27 @@ public class AlignmentUtils offset = maxoffset - flankSize; ustream_ds = flankSize; } - if (flankSize < dstream_ds) + if (flankSize <= dstream_ds) { - dstream_ds = flankSize; + dstream_ds = flankSize - 1; } } + // TODO use Character.toLowerCase to avoid creating String objects? char[] upstream = new String(ds.getSequence(s.getStart() - 1 - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray(); - char[] downstream = new String(ds.getSequence(s_end - 1, s_end + 1 + char[] downstream = new String(ds.getSequence(s_end - 1, s_end + dstream_ds)).toLowerCase().toCharArray(); char[] coreseq = s.getSequence(); char[] nseq = new char[offset + upstream.length + downstream.length + coreseq.length]; char c = core.getGapCharacter(); - // TODO could lowercase the flanking regions + int p = 0; for (; p < offset; p++) { nseq[p] = c; } - // s.setSequence(new String(upstream).toLowerCase()+new String(coreseq) + - // new String(downstream).toLowerCase()); + System.arraycopy(upstream, 0, nseq, p, upstream.length); System.arraycopy(coreseq, 0, nseq, p + upstream.length, coreseq.length); @@ -146,6 +165,7 @@ public class AlignmentUtils { for (AlignmentAnnotation aa : s.getAnnotation()) { + aa.adjustForAlignment(); // JAL-1712 fix newAl.addAnnotation(aa); } } @@ -216,9 +236,8 @@ public class AlignmentUtils * @param cdnaAlignment * @return */ - public static boolean mapProteinToCdna( - final AlignmentI proteinAlignment, - final AlignmentI cdnaAlignment) + public static boolean mapProteinAlignmentToCdna( + final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment) { if (proteinAlignment == null || cdnaAlignment == null) { @@ -264,7 +283,7 @@ public class AlignmentUtils final AlignmentI cdnaAlignment, Set mappedDna, Set mappedProtein, boolean xrefsOnly) { - boolean mappingPerformed = false; + boolean mappingExistsOrAdded = false; List thisSeqs = proteinAlignment.getSequences(); for (SequenceI aaSeq : thisSeqs) { @@ -282,7 +301,7 @@ public class AlignmentUtils * many-to-many, as that would risk mixing species with similar cDNA * sequences. */ - if (xrefsOnly && !CrossRef.haveCrossRef(aaSeq, cdnaSeq)) + if (xrefsOnly && !AlignmentUtils.haveCrossRef(aaSeq, cdnaSeq)) { continue; } @@ -297,14 +316,18 @@ public class AlignmentUtils { continue; } - if (!mappingExists(proteinAlignment.getCodonFrames(), + if (mappingExists(proteinAlignment.getCodonFrames(), aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence())) { - MapList map = mapProteinToCdna(aaSeq, cdnaSeq); + mappingExistsOrAdded = true; + } + else + { + MapList map = mapCdnaToProtein(aaSeq, cdnaSeq); if (map != null) { acf.addMap(cdnaSeq, aaSeq, map); - mappingPerformed = true; + mappingExistsOrAdded = true; proteinMapped = true; mappedDna.add(cdnaSeq); mappedProtein.add(aaSeq); @@ -316,19 +339,19 @@ public class AlignmentUtils proteinAlignment.addCodonFrame(acf); } } - return mappingPerformed; + return mappingExistsOrAdded; } /** * Answers true if the mappings include one between the given (dataset) * sequences. */ - public static boolean mappingExists(Set set, + public static boolean mappingExists(List mappings, SequenceI aaSeq, SequenceI cdnaSeq) { - if (set != null) + if (mappings != null) { - for (AlignedCodonFrame acf : set) + for (AlignedCodonFrame acf : mappings) { if (cdnaSeq == acf.getDnaForAaSeq(aaSeq)) { @@ -340,16 +363,22 @@ public class AlignmentUtils } /** - * Build a mapping (if possible) of a protein to a cDNA sequence. The cDNA - * must be three times the length of the protein, possibly after ignoring - * start and/or stop codons, and must translate to the protein. Returns null - * if no mapping is determined. + * Builds a mapping (if possible) of a cDNA to a protein sequence. + *
    + *
  • first checks if the cdna translates exactly to the protein sequence
    • + *
  • else checks for translation after removing a STOP codon
    • + *
  • else checks for translation after removing a START codon
    • + *
  • if that fails, inspect CDS features on the cDNA sequence
    • + *
+ * Returns null if no mapping is determined. * - * @param proteinSeqs + * @param proteinSeq + * the aligned protein sequence * @param cdnaSeq + * the aligned cdna sequence * @return */ - public static MapList mapProteinToCdna(SequenceI proteinSeq, + public static MapList mapCdnaToProtein(SequenceI proteinSeq, SequenceI cdnaSeq) { /* @@ -373,13 +402,13 @@ public class AlignmentUtils */ final int mappedLength = 3 * aaSeqChars.length; int cdnaLength = cdnaSeqChars.length; - int cdnaStart = 1; - int cdnaEnd = cdnaLength; - final int proteinStart = 1; - final int proteinEnd = aaSeqChars.length; + int cdnaStart = cdnaSeq.getStart(); + int cdnaEnd = cdnaSeq.getEnd(); + final int proteinStart = proteinSeq.getStart(); + final int proteinEnd = proteinSeq.getEnd(); /* - * If lengths don't match, try ignoring stop codon. + * If lengths don't match, try ignoring stop codon (if present) */ if (cdnaLength != mappedLength && cdnaLength > 2) { @@ -399,28 +428,31 @@ public class AlignmentUtils /* * If lengths still don't match, try ignoring start codon. */ + int startOffset = 0; if (cdnaLength != mappedLength && cdnaLength > 2 && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase() - .equals( - ResidueProperties.START)) + .equals(ResidueProperties.START)) { + startOffset += 3; cdnaStart += 3; cdnaLength -= 3; } - if (cdnaLength != mappedLength) - { - return null; - } - if (!translatesAs(cdnaSeqChars, cdnaStart - 1, aaSeqChars)) + if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars)) { - return null; + /* + * protein is translation of dna (+/- start/stop codons) + */ + MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] + { proteinStart, proteinEnd }, 3, 1); + return map; } - MapList map = new MapList(new int[] - { cdnaStart, cdnaEnd }, new int[] - { proteinStart, proteinEnd }, 3, 1); - return map; + + /* + * translation failed - try mapping CDS annotated regions of dna + */ + return mapCdsToProtein(cdnaSeq, proteinSeq); } /** @@ -436,23 +468,28 @@ public class AlignmentUtils protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart, char[] aaSeqChars) { - int aaResidue = 0; - for (int i = cdnaStart; i < cdnaSeqChars.length - 2 - && aaResidue < aaSeqChars.length; i += 3, aaResidue++) + if (cdnaSeqChars == null || aaSeqChars == null) + { + return false; + } + + int aaPos = 0; + int dnaPos = cdnaStart; + for (; dnaPos < cdnaSeqChars.length - 2 + && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++) { - String codon = String.valueOf(cdnaSeqChars, i, 3); - final String translated = ResidueProperties.codonTranslate( - codon); + String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + final String translated = ResidueProperties.codonTranslate(codon); + /* - * ? allow X in protein to match untranslatable in dna ? + * allow * in protein to match untranslatable in dna */ - final char aaRes = aaSeqChars[aaResidue]; - if ((translated == null || "STOP".equals(translated)) && aaRes == 'X') + final char aaRes = aaSeqChars[aaPos]; + if ((translated == null || "STOP".equals(translated)) && aaRes == '*') { continue; } - if (translated == null - || !(aaRes == translated.charAt(0))) + if (translated == null || !(aaRes == translated.charAt(0))) { // debug // System.out.println(("Mismatch at " + i + "/" + aaResidue + ": " @@ -460,8 +497,32 @@ public class AlignmentUtils return false; } } - // fail if we didn't match all of the aa sequence - return (aaResidue == aaSeqChars.length); + + /* + * check we matched all of the protein sequence + */ + if (aaPos != aaSeqChars.length) + { + return false; + } + + /* + * check we matched all of the dna except + * for optional trailing STOP codon + */ + if (dnaPos == cdnaSeqChars.length) + { + return true; + } + if (dnaPos == cdnaSeqChars.length - 3) + { + String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + if ("STOP".equals(ResidueProperties.codonTranslate(codon))) + { + return true; + } + } + return false; } /** @@ -483,7 +544,8 @@ public class AlignmentUtils boolean preserveUnmappedGaps) { /* - * Get any mappings from the source alignment to the target (dataset) sequence. + * Get any mappings from the source alignment to the target (dataset) + * sequence. */ // TODO there may be one AlignedCodonFrame per dataset sequence, or one with // all mappings. Would it help to constrain this? @@ -492,11 +554,11 @@ public class AlignmentUtils { return false; } - + /* * Locate the aligned source sequence whose dataset sequence is mapped. We - * just take the first match here (as we can't align cDNA like more than one - * protein sequence). + * just take the first match here (as we can't align like more than one + * sequence). */ SequenceI alignFrom = null; AlignedCodonFrame mapping = null; @@ -509,7 +571,7 @@ public class AlignmentUtils break; } } - + if (alignFrom == null) { return false; @@ -522,8 +584,8 @@ public class AlignmentUtils /** * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to * match residues and codons. Flags control whether existing gaps in unmapped - * (intron) and mapped (exon) regions are preserved or not. Gaps linking intro - * and exon are only retained if both flags are set. + * (intron) and mapped (exon) regions are preserved or not. Gaps between + * intron and exon are only retained if both flags are set. * * @param alignTo * @param alignFrom @@ -534,15 +596,12 @@ public class AlignmentUtils * @param preserveMappedGaps */ public static void alignSequenceAs(SequenceI alignTo, - SequenceI alignFrom, - AlignedCodonFrame mapping, String myGap, char sourceGap, - boolean preserveMappedGaps, boolean preserveUnmappedGaps) + SequenceI alignFrom, AlignedCodonFrame mapping, String myGap, + char sourceGap, boolean preserveMappedGaps, + boolean preserveUnmappedGaps) { // TODO generalise to work for Protein-Protein, dna-dna, dna-protein - final char[] thisSeq = alignTo.getSequence(); - final char[] thatAligned = alignFrom.getSequence(); - StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length); - + // aligned and dataset sequence positions, all base zero int thisSeqPos = 0; int sourceDsPos = 0; @@ -551,11 +610,17 @@ public class AlignmentUtils char myGapChar = myGap.charAt(0); int ratio = myGap.length(); - /* - * Traverse the aligned protein sequence. - */ + int fromOffset = alignFrom.getStart() - 1; + int toOffset = alignTo.getStart() - 1; int sourceGapMappedLength = 0; boolean inExon = false; + final char[] thisSeq = alignTo.getSequence(); + final char[] thatAligned = alignFrom.getSequence(); + StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length); + + /* + * Traverse the 'model' aligned sequence + */ for (char sourceChar : thatAligned) { if (sourceChar == sourceGap) @@ -565,20 +630,22 @@ public class AlignmentUtils } /* - * Found a residue. Locate its mapped codon (start) position. + * Found a non-gap character. Locate its mapped region if any. */ sourceDsPos++; // Note mapping positions are base 1, our sequence positions base 0 int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom, - sourceDsPos); + sourceDsPos + fromOffset); if (mappedPos == null) { /* - * Abort realignment if unmapped protein. Or could ignore it?? + * unmapped position; treat like a gap */ - System.err.println("Can't align: no codon mapping to residue " - + sourceDsPos + "(" + sourceChar + ")"); - return; + sourceGapMappedLength += ratio; + // System.err.println("Can't align: no codon mapping to residue " + // + sourceDsPos + "(" + sourceChar + ")"); + // return; + continue; } int mappedCodonStart = mappedPos[0]; // position (1...) of codon start @@ -594,14 +661,15 @@ public class AlignmentUtils * But then 'align dna as protein' doesn't make much sense otherwise. */ int intronLength = 0; - while (basesWritten < mappedCodonEnd && thisSeqPos < thisSeq.length) + while (basesWritten + toOffset < mappedCodonEnd + && thisSeqPos < thisSeq.length) { final char c = thisSeq[thisSeqPos++]; if (c != myGapChar) { basesWritten++; - - if (basesWritten < mappedCodonStart) + int sourcePosition = basesWritten + toOffset; + if (sourcePosition < mappedCodonStart) { /* * Found an unmapped (intron) base. First add in any preceding gaps @@ -618,7 +686,7 @@ public class AlignmentUtils } else { - final boolean startOfCodon = basesWritten == mappedCodonStart; + final boolean startOfCodon = sourcePosition == mappedCodonStart; int gapsToAdd = calculateGapsToInsert(preserveMappedGaps, preserveUnmappedGaps, sourceGapMappedLength, inExon, trailingCopiedGap.length(), intronLength, startOfCodon); @@ -647,8 +715,8 @@ public class AlignmentUtils } /* - * At end of protein sequence. Copy any remaining dna sequence, optionally - * including (intron) gaps. We do not copy trailing gaps in protein. + * At end of model aligned sequence. Copy any remaining target sequence, optionally + * including (intron) gaps. */ while (thisSeqPos < thisSeq.length) { @@ -657,6 +725,20 @@ public class AlignmentUtils { thisAligned.append(c); } + sourceGapMappedLength--; + } + + /* + * finally add gaps to pad for any trailing source gaps or + * unmapped characters + */ + if (preserveUnmappedGaps) + { + while (sourceGapMappedLength > 0) + { + thisAligned.append(myGapChar); + sourceGapMappedLength--; + } } /* @@ -679,8 +761,8 @@ public class AlignmentUtils */ protected static int calculateGapsToInsert(boolean preserveMappedGaps, boolean preserveUnmappedGaps, int sourceGapMappedLength, - boolean inExon, int trailingGapLength, - int intronLength, final boolean startOfCodon) + boolean inExon, int trailingGapLength, int intronLength, + final boolean startOfCodon) { int gapsToAdd = 0; if (startOfCodon) @@ -804,8 +886,8 @@ public class AlignmentUtils // mapping is from protein to nucleotide toDna = true; // should ideally get gap count ratio from mapping - gap = String.valueOf(new char[] - { gapCharacter, gapCharacter, gapCharacter }); + gap = String.valueOf(new char[] { gapCharacter, gapCharacter, + gapCharacter }); } else { @@ -850,7 +932,7 @@ public class AlignmentUtils { mapping.markMappedRegion(seq, pos, sr); } - newseq.append(sr.toString()); + newseq.append(sr.getCharacters()); if (first) { first = false; @@ -884,30 +966,153 @@ public class AlignmentUtils */ public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) { - Set mappings = protein.getCodonFrames(); + List unmappedProtein = new ArrayList(); + Map> alignedCodons = buildCodonColumnsMap( + protein, dna, unmappedProtein); + return alignProteinAs(protein, alignedCodons, unmappedProtein); + } + + /** + * Builds a map whose key is an aligned codon position (3 alignment column + * numbers base 0), and whose value is a map from protein sequence to each + * protein's peptide residue for that codon. The map generates an ordering of + * the codons, and allows us to read off the peptides at each position in + * order to assemble 'aligned' protein sequences. + * + * @param protein + * the protein alignment + * @param dna + * the coding dna alignment + * @param unmappedProtein + * any unmapped proteins are added to this list + * @return + */ + protected static Map> buildCodonColumnsMap( + AlignmentI protein, AlignmentI dna, + List unmappedProtein) + { + /* + * maintain a list of any proteins with no mappings - these will be + * rendered 'as is' in the protein alignment as we can't align them + */ + unmappedProtein.addAll(protein.getSequences()); + + List mappings = protein.getCodonFrames(); /* * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of * {dnaSequence, {proteinSequence, codonProduct}} at that position. The * comparator keeps the codon positions ordered. */ - Map> alignedCodons = new TreeMap>( + Map> alignedCodons = new TreeMap>( new CodonComparator()); + for (SequenceI dnaSeq : dna.getSequences()) { for (AlignedCodonFrame mapping : mappings) { - Mapping seqMap = mapping.getMappingForSequence(dnaSeq); SequenceI prot = mapping.findAlignedSequence( dnaSeq.getDatasetSequence(), protein); if (prot != null) { + Mapping seqMap = mapping.getMappingForSequence(dnaSeq); addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), seqMap, alignedCodons); + unmappedProtein.remove(prot); + } + } + } + + /* + * Finally add any unmapped peptide start residues (e.g. for incomplete + * codons) as if at the codon position before the second residue + */ + int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); + addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); + + return alignedCodons; + } + + /** + * Scans for any protein mapped from position 2 (meaning unmapped start + * position e.g. an incomplete codon), and synthesizes a 'codon' for it at the + * preceding position in the alignment + * + * @param alignedCodons + * the codon-to-peptide map + * @param mappedSequenceCount + * the number of distinct sequences in the map + */ + protected static void addUnmappedPeptideStarts( + Map> alignedCodons, + int mappedSequenceCount) + { + // TODO delete this ugly hack once JAL-2022 is resolved + // i.e. we can model startPhase > 0 (incomplete start codon) + + List sequencesChecked = new ArrayList(); + AlignedCodon lastCodon = null; + Map toAdd = new HashMap(); + + for (Entry> entry : alignedCodons + .entrySet()) + { + for (Entry sequenceCodon : entry.getValue() + .entrySet()) + { + SequenceI seq = sequenceCodon.getKey(); + if (sequencesChecked.contains(seq)) + { + continue; + } + sequencesChecked.add(seq); + AlignedCodon codon = sequenceCodon.getValue(); + if (codon.peptideCol > 1) + { + System.err + .println("Problem mapping protein with >1 unmapped start positions: " + + seq.getName()); + } + else if (codon.peptideCol == 1) + { + /* + * first position (peptideCol == 0) was unmapped - add it + */ + if (lastCodon != null) + { + AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1, + lastCodon.pos2, lastCodon.pos3, String.valueOf(seq + .getCharAt(0)), 0); + toAdd.put(seq, firstPeptide); + } + else + { + /* + * unmapped residue at start of alignment (no prior column) - + * 'insert' at nominal codon [0, 0, 0] + */ + AlignedCodon firstPeptide = new AlignedCodon(0, 0, 0, + String.valueOf(seq.getCharAt(0)), 0); + toAdd.put(seq, firstPeptide); + } + } + if (sequencesChecked.size() == mappedSequenceCount) + { + // no need to check past first mapped position in all sequences + break; } } + lastCodon = entry.getKey(); + } + + /* + * add any new codons safely after iterating over the map + */ + for (Entry startCodon : toAdd.entrySet()) + { + addCodonToMap(alignedCodons, startCodon.getValue(), + startCodon.getKey()); } - return alignProteinAs(protein, alignedCodons); } /** @@ -918,10 +1123,12 @@ public class AlignmentUtils * @param alignedCodons * an ordered map of codon positions (columns), with sequence/peptide * values present in each column + * @param unmappedProtein * @return */ protected static int alignProteinAs(AlignmentI protein, - Map> alignedCodons) + Map> alignedCodons, + List unmappedProtein) { /* * Prefill aligned sequences with gaps before inserting aligned protein @@ -933,18 +1140,22 @@ public class AlignmentUtils String allGaps = String.valueOf(gaps); for (SequenceI seq : protein.getSequences()) { - seq.setSequence(allGaps); + if (!unmappedProtein.contains(seq)) + { + seq.setSequence(allGaps); + } } int column = 0; for (AlignedCodon codon : alignedCodons.keySet()) { - final Map columnResidues = alignedCodons.get(codon); - for (Entry entry : columnResidues - .entrySet()) + final Map columnResidues = alignedCodons + .get(codon); + for (Entry entry : columnResidues.entrySet()) { // place translated codon at its column position in sequence - entry.getKey().getSequence()[column] = entry.getValue().charAt(0); + entry.getKey().getSequence()[column] = entry.getValue().product + .charAt(0); } column++; } @@ -968,25 +1179,52 @@ public class AlignmentUtils * the map we are building up */ static void addCodonPositions(SequenceI dna, SequenceI protein, - char gapChar, - Mapping seqMap, - Map> alignedCodons) + char gapChar, Mapping seqMap, + Map> alignedCodons) { Iterator codons = seqMap.getCodonIterator(dna, gapChar); + + /* + * add codon positions, and their peptide translations, to the alignment + * map, while remembering the first codon mapped + */ while (codons.hasNext()) { - AlignedCodon codon = codons.next(); - Map seqProduct = alignedCodons.get(codon); - if (seqProduct == null) + try + { + AlignedCodon codon = codons.next(); + addCodonToMap(alignedCodons, codon, protein); + } catch (IncompleteCodonException e) { - seqProduct = new HashMap(); - alignedCodons.put(codon, seqProduct); + // possible incomplete trailing codon - ignore + } catch (NoSuchElementException e) + { + // possibly peptide lacking STOP } - seqProduct.put(protein, codon.product); } } /** + * Helper method to add a codon-to-peptide entry to the aligned codons map + * + * @param alignedCodons + * @param codon + * @param protein + */ + protected static void addCodonToMap( + Map> alignedCodons, + AlignedCodon codon, SequenceI protein) + { + Map seqProduct = alignedCodons.get(codon); + if (seqProduct == null) + { + seqProduct = new HashMap(); + alignedCodons.put(codon, seqProduct); + } + seqProduct.put(protein, codon); + } + + /** * Returns true if a cDNA/Protein mapping either exists, or could be made, * between at least one pair of sequences in the two alignments. Currently, * the logic is: @@ -1004,6 +1242,11 @@ public class AlignmentUtils */ public static boolean isMappable(AlignmentI al1, AlignmentI al2) { + if (al1 == null || al2 == null) + { + return false; + } + /* * Require one nucleotide and one protein */ @@ -1013,7 +1256,7 @@ public class AlignmentUtils } AlignmentI dna = al1.isNucleotide() ? al1 : al2; AlignmentI protein = dna == al1 ? al2 : al1; - Set mappings = protein.getCodonFrames(); + List mappings = protein.getCodonFrames(); for (SequenceI dnaSeq : dna.getSequences()) { for (SequenceI proteinSeq : protein.getSequences()) @@ -1036,18 +1279,26 @@ public class AlignmentUtils * @param mappings * @return */ - public static boolean isMappable(SequenceI dnaSeq, SequenceI proteinSeq, - Set mappings) + protected static boolean isMappable(SequenceI dnaSeq, + SequenceI proteinSeq, List mappings) { - SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq.getDatasetSequence(); + if (dnaSeq == null || proteinSeq == null) + { + return false; + } + + SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq + .getDatasetSequence(); SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq : proteinSeq.getDatasetSequence(); - - /* - * Already mapped? - */ - for (AlignedCodonFrame mapping : mappings) { - if ( proteinDs == mapping.getAaForDnaSeq(dnaDs)) { + + for (AlignedCodonFrame mapping : mappings) + { + if (proteinDs == mapping.getAaForDnaSeq(dnaDs)) + { + /* + * already mapped + */ return true; } } @@ -1056,7 +1307,7 @@ public class AlignmentUtils * Just try to make a mapping (it is not yet stored), test whether * successful. */ - return mapProteinToCdna(proteinDs, dnaDs) != null; + return mapCdnaToProtein(proteinDs, dnaDs) != null; } /** @@ -1075,7 +1326,8 @@ public class AlignmentUtils * the alignment to check for presence of annotations */ public static void findAddableReferenceAnnotations( - List sequenceScope, Map labelForCalcId, + List sequenceScope, + Map labelForCalcId, final Map> candidates, AlignmentI al) { @@ -1083,7 +1335,7 @@ public class AlignmentUtils { return; } - + /* * For each sequence in scope, make a list of any annotations on the * underlying dataset sequence which are not already on the alignment. @@ -1111,8 +1363,7 @@ public class AlignmentUtils * sequence. */ final Iterable matchedAlignmentAnnotations = al - .findAnnotations(seq, dsann.getCalcId(), - dsann.label); + .findAnnotations(seq, dsann.getCalcId(), dsann.label); if (!matchedAlignmentAnnotations.iterator().hasNext()) { result.add(dsann); @@ -1160,7 +1411,7 @@ public class AlignmentUtils endRes = selectionGroup.getEndRes(); } copyAnn.restrict(startRes, endRes); - + /* * Add to the sequence (sets copyAnn.datasetSequence), unless the * original annotation is already on the sequence. @@ -1198,8 +1449,7 @@ public class AlignmentUtils Collection types, List forSequences, boolean anyType, boolean doShow) { - for (AlignmentAnnotation aa : al - .getAlignmentAnnotation()) + for (AlignmentAnnotation aa : al.getAlignmentAnnotation()) { if (anyType || types.contains(aa.label)) { @@ -1212,4 +1462,853 @@ public class AlignmentUtils } } } + + /** + * Returns true if either sequence has a cross-reference to the other + * + * @param seq1 + * @param seq2 + * @return + */ + public static boolean haveCrossRef(SequenceI seq1, SequenceI seq2) + { + // Note: moved here from class CrossRef as the latter class has dependencies + // not availability to the applet's classpath + return hasCrossRef(seq1, seq2) || hasCrossRef(seq2, seq1); + } + + /** + * Returns true if seq1 has a cross-reference to seq2. Currently this assumes + * that sequence name is structured as Source|AccessionId. + * + * @param seq1 + * @param seq2 + * @return + */ + public static boolean hasCrossRef(SequenceI seq1, SequenceI seq2) + { + if (seq1 == null || seq2 == null) + { + return false; + } + String name = seq2.getName(); + final DBRefEntry[] xrefs = seq1.getDBRefs(); + if (xrefs != null) + { + for (DBRefEntry xref : xrefs) + { + String xrefName = xref.getSource() + "|" + xref.getAccessionId(); + // case-insensitive test, consistent with DBRefEntry.equalRef() + if (xrefName.equalsIgnoreCase(name)) + { + return true; + } + } + } + return false; + } + + /** + * Constructs an alignment consisting of the mapped (CDS) regions in the given + * nucleotide sequences, and updates mappings to match. The new sequences are + * aligned as per the original sequence, with entirely gapped columns (codon + * interrupted by intron) omitted. + * + * @param dna + * aligned dna sequences + * @param mappings + * from dna to protein; these are replaced with new mappings + * @param al + * @return an alignment whose sequences are the cds-only parts of the dna + * sequences (or null if no mappings are found) + */ + public static AlignmentI makeCdsAlignment(SequenceI[] dna, + List mappings, AlignmentI al) + { + List cdsColumns = findCdsColumns(dna); + + /* + * create CDS sequences and new mappings + * (from cdna to cds, and cds to peptide) + */ + List newMappings = new ArrayList(); + List cdsSequences = new ArrayList(); + char gap = al.getGapCharacter(); + + for (SequenceI dnaSeq : dna) + { + final SequenceI ds = dnaSeq.getDatasetSequence(); + List seqMappings = MappingUtils + .findMappingsForSequence(ds, mappings); + for (AlignedCodonFrame acf : seqMappings) + { + AlignedCodonFrame newMapping = new AlignedCodonFrame(); + final List mappedCds = makeCdsSequences(dnaSeq, acf, + cdsColumns, newMapping, gap); + if (!mappedCds.isEmpty()) + { + cdsSequences.addAll(mappedCds); + newMappings.add(newMapping); + } + } + } + AlignmentI newAl = new Alignment( + cdsSequences.toArray(new SequenceI[cdsSequences.size()])); + + /* + * add new sequences to the shared dataset, set it on the new alignment + */ + List dsseqs = al.getDataset().getSequences(); + for (SequenceI seq : newAl.getSequences()) + { + if (!dsseqs.contains(seq.getDatasetSequence())) + { + dsseqs.add(seq.getDatasetSequence()); + } + } + newAl.setDataset(al.getDataset()); + + /* + * Replace the old mappings with the new ones + */ + mappings.clear(); + mappings.addAll(newMappings); + + return newAl; + } + + /** + * Returns a consolidated list of column ranges where at least one sequence + * has a CDS feature. This assumes CDS features are on genomic sequence i.e. + * are for contiguous CDS ranges (no gaps). + * + * @param seqs + * @return + */ + public static List findCdsColumns(SequenceI[] seqs) + { + // TODO use refactored code from AlignViewController + // markColumnsContainingFeatures, not reinvent the wheel! + + List result = new ArrayList(); + for (SequenceI seq : seqs) + { + result.addAll(findCdsColumns(seq)); + } + + /* + * sort and compact the list into ascending, non-overlapping ranges + */ + Collections.sort(result, new Comparator() + { + @Override + public int compare(int[] o1, int[] o2) + { + return Integer.compare(o1[0], o2[0]); + } + }); + result = MapList.coalesceRanges(result); + + return result; + } + + public static List findCdsColumns(SequenceI seq) + { + List result = new ArrayList(); + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + SequenceFeature[] sfs = seq.getSequenceFeatures(); + if (sfs != null) + { + for (SequenceFeature sf : sfs) + { + if (so.isA(sf.getType(), SequenceOntologyI.CDS)) + { + int colStart = seq.findIndex(sf.getBegin()); + int colEnd = seq.findIndex(sf.getEnd()); + result.add(new int[] { colStart, colEnd }); + } + } + } + return result; + } + + /** + * Answers true if all sequences have a gap at (or do not extend to) the + * specified column position (base 1) + * + * @param seqs + * @param col + * @return + */ + public static boolean isGappedColumn(List seqs, int col) + { + if (seqs != null) + { + for (SequenceI seq : seqs) + { + if (!Comparison.isGap(seq.getCharAt(col - 1))) + { + return false; + } + } + } + return true; + } + + /** + * Returns the column ranges (base 1) of each aligned sequence that are + * involved in any mapping. This is a helper method for aligning protein + * products of aligned transcripts. + * + * @param mappedSequences + * (possibly gapped) dna sequences + * @param mappings + * @return + */ + protected static List> getMappedColumns( + List mappedSequences, List mappings) + { + List> result = new ArrayList>(); + for (SequenceI seq : mappedSequences) + { + List columns = new ArrayList(); + List seqMappings = MappingUtils + .findMappingsForSequence(seq, mappings); + for (AlignedCodonFrame mapping : seqMappings) + { + List maps = mapping.getMappingsForSequence(seq); + for (Mapping map : maps) + { + /* + * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc } + * Find and add the overall aligned column range for each + */ + for (int[] cdsRange : map.getMap().getFromRanges()) + { + int startPos = cdsRange[0]; + int endPos = cdsRange[1]; + int startCol = seq.findIndex(startPos); + int endCol = seq.findIndex(endPos); + columns.add(new int[] { startCol, endCol }); + } + } + } + result.add(columns); + } + return result; + } + + /** + * Helper method to make cds-only sequences and populate their mappings to + * protein products + *

+ * For example, if ggCCaTTcGAg has mappings [3, 4, 6, 7, 9, 10] to protein + * then generate a sequence CCTTGA with mapping [1, 6] to the same protein + * residues + *

+ * Typically eukaryotic dna will include cds encoding for a single peptide + * sequence i.e. return a single result. Bacterial dna may have overlapping + * cds mappings coding for multiple peptides so return multiple results + * (example EMBL KF591215). + * + * @param dnaSeq + * a dna aligned sequence + * @param mapping + * containing one or more mappings of the sequence to protein + * @param ungappedCdsColumns + * @param newMappings + * the new mapping to populate, from the cds-only sequences to their + * mapped protein sequences + * @return + */ + protected static List makeCdsSequences(SequenceI dnaSeq, + AlignedCodonFrame mapping, List ungappedCdsColumns, + AlignedCodonFrame newMappings, char gapChar) + { + List cdsSequences = new ArrayList(); + List seqMappings = mapping.getMappingsForSequence(dnaSeq); + + for (Mapping seqMapping : seqMappings) + { + SequenceI cds = makeCdsSequence(dnaSeq, seqMapping, + ungappedCdsColumns, gapChar); + cds.createDatasetSequence(); + cdsSequences.add(cds); + + /* + * add new mappings, from dna to cds, and from cds to peptide + */ + MapList dnaToCds = addCdsMappings(dnaSeq.getDatasetSequence(), cds, + seqMapping, newMappings); + + /* + * transfer any features on dna that overlap the CDS + */ + transferFeatures(dnaSeq, cds, dnaToCds, null, SequenceOntologyI.CDS); + } + return cdsSequences; + } + + /** + * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the + * feature start/end ranges, optionally omitting specified feature types. + * Returns the number of features copied. + * + * @param fromSeq + * @param toSeq + * @param select + * if not null, only features of this type are copied (including + * subtypes in the Sequence Ontology) + * @param mapping + * the mapping from 'fromSeq' to 'toSeq' + * @param omitting + */ + public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, + MapList mapping, String select, String... omitting) + { + SequenceI copyTo = toSeq; + while (copyTo.getDatasetSequence() != null) + { + copyTo = copyTo.getDatasetSequence(); + } + + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + int count = 0; + SequenceFeature[] sfs = fromSeq.getSequenceFeatures(); + if (sfs != null) + { + for (SequenceFeature sf : sfs) + { + String type = sf.getType(); + if (select != null && !so.isA(type, select)) + { + continue; + } + boolean omit = false; + for (String toOmit : omitting) + { + if (type.equals(toOmit)) + { + omit = true; + } + } + if (omit) + { + continue; + } + + /* + * locate the mapped range - null if either start or end is + * not mapped (no partial overlaps are calculated) + */ + int start = sf.getBegin(); + int end = sf.getEnd(); + int[] mappedTo = mapping.locateInTo(start, end); + /* + * if whole exon range doesn't map, try interpreting it + * as 5' or 3' exon overlapping the CDS range + */ + if (mappedTo == null) + { + mappedTo = mapping.locateInTo(end, end); + if (mappedTo != null) + { + /* + * end of exon is in CDS range - 5' overlap + * to a range from the start of the peptide + */ + mappedTo[0] = 1; + } + } + if (mappedTo == null) + { + mappedTo = mapping.locateInTo(start, start); + if (mappedTo != null) + { + /* + * start of exon is in CDS range - 3' overlap + * to a range up to the end of the peptide + */ + mappedTo[1] = toSeq.getLength(); + } + } + if (mappedTo != null) + { + SequenceFeature copy = new SequenceFeature(sf); + copy.setBegin(Math.min(mappedTo[0], mappedTo[1])); + copy.setEnd(Math.max(mappedTo[0], mappedTo[1])); + copyTo.addSequenceFeature(copy); + count++; + } + } + } + return count; + } + + /** + * Creates and adds mappings + *

    + *
  • from cds to peptide
    • + *
  • from dna to cds
    • + *
+ * and returns the dna-to-cds mapping + * + * @param dnaSeq + * @param cdsSeq + * @param dnaMapping + * @param newMappings + * @return + */ + protected static MapList addCdsMappings(SequenceI dnaSeq, + SequenceI cdsSeq, Mapping dnaMapping, + AlignedCodonFrame newMappings) + { + cdsSeq.createDatasetSequence(); + + /* + * CDS to peptide is just a contiguous 3:1 mapping, with + * the peptide ranges taken unchanged from the dna mapping + */ + List cdsRanges = new ArrayList(); + SequenceI cdsDataset = cdsSeq.getDatasetSequence(); + cdsRanges.add(new int[] { 1, cdsDataset.getLength() }); + MapList cdsToPeptide = new MapList(cdsRanges, dnaMapping.getMap() + .getToRanges(), 3, 1); + newMappings.addMap(cdsDataset, dnaMapping.getTo(), cdsToPeptide); + + /* + * dna 'from' ranges map 1:1 to the contiguous extracted CDS + */ + MapList dnaToCds = new MapList(dnaMapping.getMap().getFromRanges(), + cdsRanges, 1, 1); + newMappings.addMap(dnaSeq, cdsDataset, dnaToCds); + return dnaToCds; + } + + /** + * Makes and returns a CDS-only sequence, where the CDS regions are identified + * as the 'from' ranges of the mapping on the dna. + * + * @param dnaSeq + * nucleotide sequence + * @param seqMapping + * mappings from CDS regions of nucleotide + * @param ungappedCdsColumns + * @return + */ + protected static SequenceI makeCdsSequence(SequenceI dnaSeq, + Mapping seqMapping, List ungappedCdsColumns, char gapChar) + { + int cdsWidth = MappingUtils.getLength(ungappedCdsColumns); + + /* + * populate CDS columns with the aligned + * column character if that column is mapped (which may be a gap + * if an intron interrupts a codon), else with a gap + */ + List fromRanges = seqMapping.getMap().getFromRanges(); + char[] cdsChars = new char[cdsWidth]; + int pos = 0; + for (int[] columns : ungappedCdsColumns) + { + for (int i = columns[0]; i <= columns[1]; i++) + { + char dnaChar = dnaSeq.getCharAt(i - 1); + if (Comparison.isGap(dnaChar)) + { + cdsChars[pos] = gapChar; + } + else + { + int seqPos = dnaSeq.findPosition(i - 1); + if (MappingUtils.contains(fromRanges, seqPos)) + { + cdsChars[pos] = dnaChar; + } + else + { + cdsChars[pos] = gapChar; + } + } + pos++; + } + } + SequenceI cdsSequence = new Sequence(dnaSeq.getName(), + String.valueOf(cdsChars)); + + transferDbRefs(seqMapping.getTo(), cdsSequence); + + return cdsSequence; + } + + /** + * Locate any xrefs to CDS databases on the protein product and attach to the + * CDS sequence. Also add as a sub-token of the sequence name. + * + * @param from + * @param to + */ + protected static void transferDbRefs(SequenceI from, SequenceI to) + { + String cdsAccId = FeatureProperties.getCodingFeature(DBRefSource.EMBL); + DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(from.getDBRefs(), + DBRefSource.CODINGDBS); + if (cdsRefs != null) + { + for (DBRefEntry cdsRef : cdsRefs) + { + to.addDBRef(new DBRefEntry(cdsRef)); + cdsAccId = cdsRef.getAccessionId(); + } + } + if (!to.getName().contains(cdsAccId)) + { + to.setName(to.getName() + "|" + cdsAccId); + } + } + + /** + * Returns a mapping from dna to protein by inspecting sequence features of + * type "CDS" on the dna. + * + * @param dnaSeq + * @param proteinSeq + * @return + */ + public static MapList mapCdsToProtein(SequenceI dnaSeq, + SequenceI proteinSeq) + { + List ranges = findCdsPositions(dnaSeq); + int mappedDnaLength = MappingUtils.getLength(ranges); + + int proteinLength = proteinSeq.getLength(); + int proteinStart = proteinSeq.getStart(); + int proteinEnd = proteinSeq.getEnd(); + + /* + * incomplete start codon may mean X at start of peptide + * we ignore both for mapping purposes + */ + if (proteinSeq.getCharAt(0) == 'X') + { + // todo JAL-2022 support startPhase > 0 + proteinStart++; + proteinLength--; + } + List proteinRange = new ArrayList(); + + /* + * dna length should map to protein (or protein plus stop codon) + */ + int codesForResidues = mappedDnaLength / 3; + if (codesForResidues == (proteinLength + 1)) + { + // assuming extra codon is for STOP and not in peptide + codesForResidues--; + } + if (codesForResidues == proteinLength) + { + proteinRange.add(new int[] { proteinStart, proteinEnd }); + return new MapList(ranges, proteinRange, 3, 1); + } + return null; + } + + /** + * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of + * start/end positions of sequence features of type "CDS" (or a sub-type of + * CDS in the Sequence Ontology) + * + * @param dnaSeq + * @return + */ + public static List findCdsPositions(SequenceI dnaSeq) + { + List result = new ArrayList(); + SequenceFeature[] sfs = dnaSeq.getSequenceFeatures(); + if (sfs == null) + { + return result; + } + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + for (SequenceFeature sf : sfs) + { + /* + * process a CDS feature (or a sub-type of CDS) + */ + if (so.isA(sf.getType(), SequenceOntologyI.CDS)) + { + int phase = 0; + try { + phase = Integer.parseInt(sf.getPhase()); + } catch (NumberFormatException e) + { + // ignore + } + /* + * phase > 0 on first codon means 5' incomplete - skip to the start + * of the next codon; example ENST00000496384 + */ + int begin = sf.getBegin(); + int end = sf.getEnd(); + if (result.isEmpty()) + { + // TODO JAL-2022 support start phase > 0 + begin += phase; + if (begin > end) + { + continue; // shouldn't happen? + } + } + result.add(new int[] { begin, end }); + } + } + return result; + } + + /** + * Maps exon features from dna to protein, and computes variants in peptide + * product generated by variants in dna, and adds them as sequence_variant + * features on the protein sequence. Returns the number of variant features + * added. + * + * @param dnaSeq + * @param peptide + * @param dnaToProtein + */ + public static int computeProteinFeatures(SequenceI dnaSeq, + SequenceI peptide, MapList dnaToProtein) + { + while (dnaSeq.getDatasetSequence() != null) + { + dnaSeq = dnaSeq.getDatasetSequence(); + } + while (peptide.getDatasetSequence() != null) + { + peptide = peptide.getDatasetSequence(); + } + + transferFeatures(dnaSeq, peptide, dnaToProtein, + SequenceOntologyI.EXON); + + LinkedHashMap variants = buildDnaVariantsMap( + dnaSeq, dnaToProtein); + + /* + * scan codon variations, compute peptide variants and add to peptide sequence + */ + int count = 0; + for (Entry variant : variants.entrySet()) + { + int peptidePos = variant.getKey(); + String[][] codonVariants = variant.getValue(); + String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); // 0-based + List peptideVariants = computePeptideVariants(codonVariants, + residue); + if (!peptideVariants.isEmpty()) + { + String desc = StringUtils.listToDelimitedString(peptideVariants, + ", "); + SequenceFeature sf = new SequenceFeature( + SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, + peptidePos, 0f, null); + peptide.addSequenceFeature(sf); + count++; + } + } + + /* + * ugly sort to get sequence features in start position order + * - would be better to store in Sequence as a TreeSet instead? + */ + Arrays.sort(peptide.getSequenceFeatures(), + new Comparator() + { + @Override + public int compare(SequenceFeature o1, SequenceFeature o2) + { + int c = Integer.compare(o1.getBegin(), o2.getBegin()); + return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd()) + : c; + } + }); + return count; + } + + /** + * Builds a map whose key is position in the protein sequence, and value is an + * array of all variants for the coding codon positions + * + * @param dnaSeq + * @param dnaToProtein + * @return + */ + static LinkedHashMap buildDnaVariantsMap( + SequenceI dnaSeq, MapList dnaToProtein) + { + /* + * map from peptide position to all variant features of the codon for it + * LinkedHashMap ensures we add the peptide features in sequence order + */ + LinkedHashMap variants = new LinkedHashMap(); + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + + SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures(); + if (dnaFeatures == null) + { + return variants; + } + + int dnaStart = dnaSeq.getStart(); + int[] lastCodon = null; + int lastPeptidePostion = 0; + + /* + * build a map of codon variations for peptides + */ + for (SequenceFeature sf : dnaFeatures) + { + int dnaCol = sf.getBegin(); + if (dnaCol != sf.getEnd()) + { + // not handling multi-locus variant features + continue; + } + if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT)) + { + int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol); + if (mapsTo == null) + { + // feature doesn't lie within coding region + continue; + } + int peptidePosition = mapsTo[0]; + String[][] codonVariants = variants.get(peptidePosition); + if (codonVariants == null) + { + codonVariants = new String[3][]; + variants.put(peptidePosition, codonVariants); + } + + /* + * extract dna variants to a string array + */ + String alls = (String) sf.getValue("alleles"); + if (alls == null) + { + continue; + } + String[] alleles = alls.toUpperCase().split(","); + int i = 0; + for (String allele : alleles) + { + alleles[i++] = allele.trim(); // lose any space characters "A, G" + } + + /* + * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10] + */ + int[] codon = peptidePosition == lastPeptidePostion ? lastCodon + : MappingUtils.flattenRanges(dnaToProtein.locateInFrom( + peptidePosition, peptidePosition)); + lastPeptidePostion = peptidePosition; + lastCodon = codon; + + /* + * save nucleotide (and this variant) for each codon position + */ + for (int codonPos = 0; codonPos < 3; codonPos++) + { + String nucleotide = String.valueOf( + dnaSeq.getCharAt(codon[codonPos] - dnaStart)) + .toUpperCase(); + if (codonVariants[codonPos] == null) + { + /* + * record current dna base + */ + codonVariants[codonPos] = new String[] { nucleotide }; + } + if (codon[codonPos] == dnaCol) + { + /* + * add alleles to dna base (and any previously found alleles) + */ + String[] known = codonVariants[codonPos]; + String[] dnaVariants = new String[alleles.length + known.length]; + System.arraycopy(known, 0, dnaVariants, 0, known.length); + System.arraycopy(alleles, 0, dnaVariants, known.length, + alleles.length); + codonVariants[codonPos] = dnaVariants; + } + } + } + } + return variants; + } + + /** + * Returns a sorted, non-redundant list of all peptide translations generated + * by the given dna variants, excluding the current residue value + * + * @param codonVariants + * an array of base values (acgtACGT) for codon positions 1, 2, 3 + * @param residue + * the current residue translation + * @return + */ + static List computePeptideVariants( + String[][] codonVariants, String residue) + { + List result = new ArrayList(); + for (String base1 : codonVariants[0]) + { + for (String base2 : codonVariants[1]) + { + for (String base3 : codonVariants[2]) + { + String codon = base1 + base2 + base3; + /* + * get peptide translation of codon e.g. GAT -> D + * note that variants which are not single alleles, + * e.g. multibase variants or HGMD_MUTATION etc + * are ignored here + */ + String peptide = codon.contains("-") ? "-" + : (codon.length() > 3 ? null : ResidueProperties + .codonTranslate(codon)); + if (peptide != null && !result.contains(peptide) + && !peptide.equalsIgnoreCase(residue)) + { + result.add(peptide); + } + } + } + } + + /* + * sort alphabetically with STOP at the end + */ + Collections.sort(result, new Comparator() + { + + @Override + public int compare(String o1, String o2) + { + if ("STOP".equals(o1)) + { + return 1; + } + else if ("STOP".equals(o2)) + { + return -1; + } + else + { + return o1.compareTo(o2); + } + } + }); + return result; + } }