X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=949c47a69d9f8c559feb9f344df8a9033f36f654;hb=9c39e96af6b84257604da448101505361dced686;hp=41538eb431646924fe2c2c0f645ac2d0e3dc1af5;hpb=f6d5cf647066635879824367816e4717e07bbb9d;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 41538eb..949c47a 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,29 +20,35 @@ */ package jalview.analysis; +import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE; + 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.SequenceOntology; +import jalview.io.gff.SequenceOntologyFactory; +import jalview.io.gff.SequenceOntologyI; import jalview.schemes.ResidueProperties; -import jalview.util.DBRefUtils; +import jalview.util.Comparison; import jalview.util.MapList; import jalview.util.MappingUtils; +import jalview.util.StringUtils; +import java.io.UnsupportedEncodingException; +import java.net.URLEncoder; 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; @@ -50,6 +56,7 @@ 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; @@ -63,6 +70,31 @@ import java.util.TreeMap; public class AlignmentUtils { + 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 + { + String base; + + SequenceFeature variant; + + DnaVariant(String nuc) + { + base = nuc; + } + + DnaVariant(String nuc, SequenceFeature var) + { + base = nuc; + variant = var; + } + } + /** * given an existing alignment, create a new alignment including all, or up to * flankSize additional symbols from each sequence's dataset sequence @@ -316,7 +348,7 @@ public class AlignmentUtils } else { - MapList map = mapProteinSequenceToCdna(aaSeq, cdnaSeq); + MapList map = mapCdnaToProtein(aaSeq, cdnaSeq); if (map != null) { acf.addMap(cdnaSeq, aaSeq, map); @@ -356,16 +388,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. + * + * 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 mapProteinSequenceToCdna(SequenceI proteinSeq, + public static MapList mapCdnaToProtein(SequenceI proteinSeq, SequenceI cdnaSeq) { /* @@ -395,7 +433,7 @@ public class AlignmentUtils 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) { @@ -426,17 +464,20 @@ public class AlignmentUtils cdnaLength -= 3; } - if (cdnaLength != mappedLength) - { - return null; - } - if (!translatesAs(cdnaSeqChars, startOffset, 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); } /** @@ -457,16 +498,18 @@ public class AlignmentUtils return false; } - int aaResidue = 0; - for (int i = cdnaStart; i < cdnaSeqChars.length - 2 - && aaResidue < aaSeqChars.length; i += 3, aaResidue++) + int aaPos = 0; + int dnaPos = cdnaStart; + for (; dnaPos < cdnaSeqChars.length - 2 + && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++) { - String codon = String.valueOf(cdnaSeqChars, i, 3); + String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); final String translated = ResidueProperties.codonTranslate(codon); + /* * allow * in protein to match untranslatable in dna */ - final char aaRes = aaSeqChars[aaResidue]; + final char aaRes = aaSeqChars[aaPos]; if ((translated == null || "STOP".equals(translated)) && aaRes == '*') { continue; @@ -479,8 +522,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; } /** @@ -522,7 +589,7 @@ public class AlignmentUtils AlignedCodonFrame mapping = null; for (AlignedCodonFrame mp : mappings) { - alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al); + alignFrom = mp.findAlignedSequence(seq, al); if (alignFrom != null) { mapping = mp; @@ -771,187 +838,164 @@ public class AlignmentUtils } /** - * Returns a list of sequences mapped from the given sequences and aligned - * (gapped) in the same way. For example, the cDNA for aligned protein, where - * a single gap in protein generates three gaps in cDNA. + * Realigns the given protein to match the alignment of the dna, using codon + * mappings to translate aligned codon positions to protein residues. * - * @param sequences - * @param gapCharacter - * @param mappings - * @return + * @param protein + * the alignment whose sequences are realigned by this method + * @param dna + * the dna alignment whose alignment we are 'copying' + * @return the number of sequences that were realigned */ - public static List getAlignedTranslation( - List sequences, char gapCharacter, - Set mappings) + public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) { - List alignedSeqs = new ArrayList(); - - for (SequenceI seq : sequences) - { - List mapped = getAlignedTranslation(seq, gapCharacter, - mappings); - alignedSeqs.addAll(mapped); - } - return alignedSeqs; + List unmappedProtein = new ArrayList(); + Map> alignedCodons = buildCodonColumnsMap( + protein, dna, unmappedProtein); + return alignProteinAs(protein, alignedCodons, unmappedProtein); } /** - * Returns sequences aligned 'like' the source sequence, as mapped by the - * given mappings. Normally we expect zero or one 'mapped' sequences, but this - * will support 1-to-many as well. + * 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 seq - * @param gapCharacter - * @param mappings + * @param protein + * the protein alignment + * @param dna + * the coding dna alignment + * @param unmappedProtein + * any unmapped proteins are added to this list * @return */ - protected static List getAlignedTranslation(SequenceI seq, - char gapCharacter, Set mappings) + protected static Map> buildCodonColumnsMap( + AlignmentI protein, AlignmentI dna, + List unmappedProtein) { - List result = new ArrayList(); - for (AlignedCodonFrame mapping : mappings) + /* + * 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>( + new CodonComparator()); + + for (SequenceI dnaSeq : dna.getSequences()) { - if (mapping.involvesSequence(seq)) + for (AlignedCodonFrame mapping : mappings) { - SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping); - if (mapped != null) + SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein); + if (prot != null) { - result.add(mapped); + Mapping seqMap = mapping.getMappingForSequence(dnaSeq); + addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), + seqMap, alignedCodons); + unmappedProtein.remove(prot); } } } - return result; + + /* + * Finally add any unmapped peptide start residues (e.g. for incomplete + * codons) as if at the codon position before the second residue + */ + // TODO resolve JAL-2022 so this fudge can be removed + int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); + addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); + + return alignedCodons; } /** - * Returns the translation of 'seq' (as held in the mapping) with - * corresponding alignment (gaps). + * 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 seq - * @param gapCharacter - * @param mapping - * @return + * @param alignedCodons + * the codon-to-peptide map + * @param mappedSequenceCount + * the number of distinct sequences in the map */ - protected static SequenceI getAlignedTranslation(SequenceI seq, - char gapCharacter, AlignedCodonFrame mapping) + protected static void addUnmappedPeptideStarts( + Map> alignedCodons, + int mappedSequenceCount) { - String gap = String.valueOf(gapCharacter); - boolean toDna = false; - int fromRatio = 1; - SequenceI mapTo = mapping.getDnaForAaSeq(seq); - if (mapTo != null) - { - // mapping is from protein to nucleotide - toDna = true; - // should ideally get gap count ratio from mapping - gap = String.valueOf(new char[] { gapCharacter, gapCharacter, - gapCharacter }); - } - else - { - // mapping is from nucleotide to protein - mapTo = mapping.getAaForDnaSeq(seq); - fromRatio = 3; - } - StringBuilder newseq = new StringBuilder(seq.getLength() - * (toDna ? 3 : 1)); + // TODO delete this ugly hack once JAL-2022 is resolved + // i.e. we can model startPhase > 0 (incomplete start codon) - int residueNo = 0; // in seq, base 1 - int[] phrase = new int[fromRatio]; - int phraseOffset = 0; - int gapWidth = 0; - boolean first = true; - final Sequence alignedSeq = new Sequence("", ""); + List sequencesChecked = new ArrayList(); + AlignedCodon lastCodon = null; + Map toAdd = new HashMap(); - for (char c : seq.getSequence()) + for (Entry> entry : alignedCodons + .entrySet()) { - if (c == gapCharacter) + for (Entry sequenceCodon : entry.getValue() + .entrySet()) { - gapWidth++; - if (gapWidth >= fromRatio) + SequenceI seq = sequenceCodon.getKey(); + if (sequencesChecked.contains(seq)) { - newseq.append(gap); - gapWidth = 0; + continue; } - } - else - { - phrase[phraseOffset++] = residueNo + 1; - if (phraseOffset == fromRatio) + 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) { /* - * Have read a whole codon (or protein residue), now translate: map - * source phrase to positions in target sequence add characters at - * these positions to newseq Note mapping positions are base 1, our - * sequence positions base 0. + * first position (peptideCol == 0) was unmapped - add it */ - SearchResults sr = new SearchResults(); - for (int pos : phrase) + if (lastCodon != null) { - mapping.markMappedRegion(seq, pos, sr); + AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1, + lastCodon.pos2, lastCodon.pos3, String.valueOf(seq + .getCharAt(0)), 0); + toAdd.put(seq, firstPeptide); } - newseq.append(sr.getCharacters()); - if (first) + else { - first = false; - // Hack: Copy sequence dataset, name and description from - // SearchResults.match[0].sequence - // TODO? carry over sequence names from original 'complement' - // alignment - SequenceI mappedTo = sr.getResultSequence(0); - alignedSeq.setName(mappedTo.getName()); - alignedSeq.setDescription(mappedTo.getDescription()); - alignedSeq.setDatasetSequence(mappedTo); + /* + * 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); } - phraseOffset = 0; } - residueNo++; + if (sequencesChecked.size() == mappedSequenceCount) + { + // no need to check past first mapped position in all sequences + break; + } } + lastCodon = entry.getKey(); } - alignedSeq.setSequence(newseq.toString()); - return alignedSeq; - } - - /** - * Realigns the given protein to match the alignment of the dna, using codon - * mappings to translate aligned codon positions to protein residues. - * - * @param protein - * the alignment whose sequences are realigned by this method - * @param dna - * the dna alignment whose alignment we are 'copying' - * @return the number of sequences that were realigned - */ - public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) - { - List unmappedProtein = new ArrayList(); - 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. + * add any new codons safely after iterating over the map */ - Map> alignedCodons = new TreeMap>( - new CodonComparator()); - for (SequenceI dnaSeq : dna.getSequences()) + for (Entry startCodon : toAdd.entrySet()) { - for (AlignedCodonFrame mapping : mappings) - { - Mapping seqMap = mapping.getMappingForSequence(dnaSeq); - SequenceI prot = mapping.findAlignedSequence( - dnaSeq.getDatasetSequence(), protein); - if (prot != null) - { - addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), - seqMap, alignedCodons); - unmappedProtein.remove(prot); - } - } + addCodonToMap(alignedCodons, startCodon.getValue(), + startCodon.getKey()); } - return alignProteinAs(protein, alignedCodons, unmappedProtein); } /** @@ -966,7 +1010,7 @@ public class AlignmentUtils * @return */ protected static int alignProteinAs(AlignmentI protein, - Map> alignedCodons, + Map> alignedCodons, List unmappedProtein) { /* @@ -988,12 +1032,13 @@ public class AlignmentUtils int column = 0; for (AlignedCodon codon : alignedCodons.keySet()) { - final Map columnResidues = alignedCodons + final Map columnResidues = alignedCodons .get(codon); - for (Entry entry : columnResidues.entrySet()) + 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++; } @@ -1018,23 +1063,51 @@ public class AlignmentUtils */ static void addCodonPositions(SequenceI dna, SequenceI protein, char gapChar, Mapping seqMap, - Map> alignedCodons) + 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) + { + // possible incomplete trailing codon - ignore + } catch (NoSuchElementException e) { - seqProduct = new HashMap(); - alignedCodons.put(codon, seqProduct); + // 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: @@ -1117,7 +1190,7 @@ public class AlignmentUtils * Just try to make a mapping (it is not yet stored), test whether * successful. */ - return mapProteinSequenceToCdna(proteinDs, dnaDs) != null; + return mapCdnaToProtein(proteinDs, dnaDs) != null; } /** @@ -1319,97 +1392,131 @@ public class AlignmentUtils } /** - * Constructs an alignment consisting of the mapped cds regions in the given - * nucleotide sequences, and updates mappings to match. + * Constructs an alignment consisting of the mapped (CDS) regions in the given + * nucleotide sequences, and updates mappings to match. The CDS sequences are + * added to the original alignment's dataset, which is shared by the new + * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are + * added to the alignment dataset. * * @param dna * aligned dna sequences * @param mappings - * from dna to protein; these are replaced with new mappings + * from dna to protein + * @param al * @return an alignment whose sequences are the cds-only parts of the dna - * sequences (or null if no cds are found) + * sequences (or null if no mappings are found) */ public static AlignmentI makeCdsAlignment(SequenceI[] dna, - List mappings) + List mappings, AlignmentI al) { - List newMappings = new ArrayList(); - List cdsSequences = new ArrayList(); - - for (SequenceI dnaSeq : dna) + List cdsSeqs = new ArrayList(); + + /* + * construct CDS sequences from the (cds-to-protein) mappings made earlier; + * this makes it possible to model multiple products from dna (e.g. EMBL); + * however it does mean we don't have the EMBL protein_id (a property on + * the CDS features) in order to make the CDS sequence name :-( + */ + for (SequenceI seq : dna) { - final SequenceI ds = dnaSeq.getDatasetSequence(); + AlignedCodonFrame cdsMappings = new AlignedCodonFrame(); List seqMappings = MappingUtils - .findMappingsForSequence(ds, mappings); - for (AlignedCodonFrame acf : seqMappings) + .findMappingsForSequence(seq, mappings); + List alignmentMappings = al.getCodonFrames(); + for (AlignedCodonFrame mapping : seqMappings) { - AlignedCodonFrame newMapping = new AlignedCodonFrame(); - final List mappedCds = makeCdsSequences(ds, acf, - newMapping); - if (!mappedCds.isEmpty()) + for (Mapping aMapping : mapping.getMappingsFromSequence(seq)) { - cdsSequences.addAll(mappedCds); - newMappings.add(newMapping); + SequenceI cdsSeq = makeCdsSequence(seq.getDatasetSequence(), + aMapping); + cdsSeqs.add(cdsSeq); + + /* + * add a mapping from CDS to the (unchanged) mapped to range + */ + List cdsRange = Collections.singletonList(new int[] { 1, + cdsSeq.getLength() }); + MapList map = new MapList(cdsRange, aMapping.getMap() + .getToRanges(), aMapping.getMap().getFromRatio(), + aMapping.getMap().getToRatio()); + cdsMappings.addMap(cdsSeq, aMapping.getTo(), map); + + /* + * add another mapping from original 'from' range to CDS + */ + map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1, + 1); + cdsMappings.addMap(seq.getDatasetSequence(), cdsSeq, map); + + alignmentMappings.add(cdsMappings); + + /* + * transfer any features on dna that overlap the CDS + */ + transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS); } } } - AlignmentI al = new Alignment( - cdsSequences.toArray(new SequenceI[cdsSequences.size()])); - al.setDataset(null); /* - * Replace the old mappings with the new ones + * add CDS seqs to shared dataset */ - mappings.clear(); - mappings.addAll(newMappings); + Alignment dataset = al.getDataset(); + for (SequenceI seq : cdsSeqs) + { + if (!dataset.getSequences().contains(seq.getDatasetSequence())) + { + dataset.addSequence(seq.getDatasetSequence()); + } + } + AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs + .size()])); + cds.setDataset(dataset); - return al; + return cds; } /** - * 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). + * Helper method that makes a CDS sequence as defined by the mappings from the + * given sequence i.e. extracts the 'mapped from' ranges (which may be on + * forward or reverse strand). * - * @param dnaSeq - * a dna dataset sequence + * @param seq * @param mapping - * containing one or more mappings of the sequence to protein - * @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, AlignedCodonFrame newMappings) + static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) { - List cdsSequences = new ArrayList(); - List seqMappings = mapping.getMappingsForSequence(dnaSeq); + char[] seqChars = seq.getSequence(); + List fromRanges = mapping.getMap().getFromRanges(); + int cdsWidth = MappingUtils.getLength(fromRanges); + char[] newSeqChars = new char[cdsWidth]; - for (Mapping seqMapping : seqMappings) + int newPos = 0; + for (int[] range : fromRanges) { - SequenceI cds = makeCdsSequence(dnaSeq, seqMapping); - cdsSequences.add(cds); - - /* - * add new mappings, from dna to cds, and from cds to peptide - */ - MapList dnaToCds = addCdsMappings(dnaSeq, cds, seqMapping, - newMappings); - - /* - * transfer any features on dna that overlap the CDS - */ - transferFeatures(dnaSeq, cds, dnaToCds, null, "CDS" /* SequenceOntology.CDS */); + if (range[0] <= range[1]) + { + // forward strand mapping - just copy the range + int length = range[1] - range[0] + 1; + System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos, + length); + newPos += length; + } + else + { + // reverse strand mapping - copy and complement one by one + for (int i = range[0]; i >= range[1]; i--) + { + newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); + } + } } - return cdsSequences; + + SequenceI newSeq = new Sequence(seq.getName() + "|" + + mapping.getTo().getName(), newSeqChars, 1, newPos); + newSeq.createDatasetSequence(); + return newSeq; } /** @@ -1435,7 +1542,7 @@ public class AlignmentUtils copyTo = copyTo.getDatasetSequence(); } - SequenceOntology so = SequenceOntology.getInstance(); + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); int count = 0; SequenceFeature[] sfs = fromSeq.getSequenceFeatures(); if (sfs != null) @@ -1509,106 +1616,766 @@ public class AlignmentUtils } /** - * Creates and adds mappings - *

    - *
  • from cds to peptide
  • - *
  • from dna to cds
  • - *
- * and returns the dna-to-cds mapping + * Returns a mapping from dna to protein by inspecting sequence features of + * type "CDS" on the dna. * * @param dnaSeq - * @param cdsSeq - * @param dnaMapping - * @param newMappings + * @param proteinSeq * @return */ - protected static MapList addCdsMappings(SequenceI dnaSeq, - SequenceI cdsSeq, - Mapping dnaMapping, AlignedCodonFrame newMappings) + public static MapList mapCdsToProtein(SequenceI dnaSeq, + SequenceI proteinSeq) { - cdsSeq.createDatasetSequence(); + List ranges = findCdsPositions(dnaSeq); + int mappedDnaLength = MappingUtils.getLength(ranges); + + int proteinLength = proteinSeq.getLength(); + int proteinStart = proteinSeq.getStart(); + int proteinEnd = proteinSeq.getEnd(); /* - * CDS to peptide is just a contiguous 3:1 mapping, with - * the peptide ranges taken unchanged from the dna mapping + * incomplete start codon may mean X at start of peptide + * we ignore both for mapping purposes */ - List cdsRanges = new ArrayList(); - cdsRanges.add(new int[] { 1, cdsSeq.getLength() }); - MapList cdsToPeptide = new MapList(cdsRanges, dnaMapping.getMap() - .getToRanges(), 3, 1); - newMappings.addMap(cdsSeq.getDatasetSequence(), dnaMapping.getTo(), - cdsToPeptide); + if (proteinSeq.getCharAt(0) == 'X') + { + // todo JAL-2022 support startPhase > 0 + proteinStart++; + proteinLength--; + } + List proteinRange = new ArrayList(); /* - * dna 'from' ranges map 1:1 to the contiguous extracted CDS + * dna length should map to protein (or protein plus stop codon) */ - MapList dnaToCds = new MapList( - dnaMapping.getMap().getFromRanges(), cdsRanges, 1, 1); - newMappings.addMap(dnaSeq, cdsSeq.getDatasetSequence(), dnaToCds); - return dnaToCds; + 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; } /** - * Makes and returns a CDS-only sequence, where the CDS regions are identified - * as the 'from' ranges of the mapping on the dna. + * 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). The ranges are sorted into ascending start + * position order, so this method is only valid for linear CDS in the same + * sense as the protein product. * * @param dnaSeq - * nucleotide sequence - * @param seqMapping - * mappings from CDS regions of nucleotide * @return */ - protected static SequenceI makeCdsSequence(SequenceI dnaSeq, - Mapping seqMapping) + public static List findCdsPositions(SequenceI dnaSeq) { - StringBuilder newSequence = new StringBuilder(dnaSeq.getLength()); - final char[] dna = dnaSeq.getSequence(); - int offset = dnaSeq.getStart() - 1; + List result = new ArrayList(); + SequenceFeature[] sfs = dnaSeq.getSequenceFeatures(); + if (sfs == null) + { + return result; + } + + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + int startPhase = 0; + + 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()) + { + begin += phase; + if (begin > end) + { + // shouldn't happen! + System.err + .println("Error: start phase extends beyond start CDS in " + + dnaSeq.getName()); + } + } + result.add(new int[] { begin, end }); + } + } /* - * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc } + * remove 'startPhase' positions (usually 0) from the first range + * so we begin at the start of a complete codon */ - final List dnaCdsRanges = seqMapping.getMap().getFromRanges(); - for (int[] range : dnaCdsRanges) + if (!result.isEmpty()) { - // TODO handle reverse mapping as well (range[1] < range[0]) - for (int pos = range[0]; pos <= range[1]; pos++) + // TODO JAL-2022 correctly model start phase > 0 + result.get(0)[0] += startPhase; + } + + /* + * Finally sort ranges by start position. This avoids a dependency on + * keeping features in order on the sequence (if they are in order anyway, + * the sort will have almost no work to do). The implicit assumption is CDS + * ranges are assembled in order. Other cases should not use this method, + * but instead construct an explicit mapping for CDS (e.g. EMBL parsing). + */ + Collections.sort(result, new Comparator() + { + @Override + public int compare(int[] o1, int[] o2) { - newSequence.append(dna[pos - offset - 1]); + return Integer.compare(o1[0], o2[0]); } + }); + 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(); } - SequenceI cds = new Sequence(dnaSeq.getName(), - newSequence.toString()); + transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON); - transferDbRefs(seqMapping.getTo(), cds); + /* + * compute protein variants from dna variants and codon mappings; + * NB - alternatively we could retrieve this using the REST service e.g. + * http://rest.ensembl.org/overlap/translation + * /ENSP00000288602?feature=transcript_variation;content-type=text/xml + * which would be a bit slower but possibly more reliable + */ - return cds; + /* + * build a map with codon variations for each potentially varying peptide + */ + 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(); + List[] codonVariants = variant.getValue(); + count += computePeptideVariants(peptide, peptidePos, codonVariants); + } + + /* + * sort to get sequence features in start position order + * - would be better to store in Sequence as a TreeSet or NCList? + */ + if (peptide.getSequenceFeatures() != null) + { + 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; + } + + /** + * Computes non-synonymous peptide variants from codon variants and adds them + * as sequence_variant features on the protein sequence (one feature per + * allele variant). Selected attributes (variant id, clinical significance) + * are copied over to the new features. + * + * @param peptide + * the protein sequence + * @param peptidePos + * the position to compute peptide variants for + * @param codonVariants + * a list of dna variants per codon position + * @return the number of features added + */ + static int computePeptideVariants(SequenceI peptide, int peptidePos, + List[] codonVariants) + { + String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); + int count = 0; + String base1 = codonVariants[0].get(0).base; + String base2 = codonVariants[1].get(0).base; + String base3 = codonVariants[2].get(0).base; + + /* + * variants in first codon base + */ + for (DnaVariant var : codonVariants[0]) + { + if (var.variant != null) + { + String alleles = (String) var.variant.getValue("alleles"); + if (alleles != null) + { + for (String base : alleles.split(",")) + { + String codon = base + base2 + base3; + if (addPeptideVariant(peptide, peptidePos, residue, var, codon)) + { + count++; + } + } + } + } + } + + /* + * variants in second codon base + */ + for (DnaVariant var : codonVariants[1]) + { + if (var.variant != null) + { + String alleles = (String) var.variant.getValue("alleles"); + if (alleles != null) + { + for (String base : alleles.split(",")) + { + String codon = base1 + base + base3; + if (addPeptideVariant(peptide, peptidePos, residue, var, codon)) + { + count++; + } + } + } + } + } + + /* + * variants in third codon base + */ + for (DnaVariant var : codonVariants[2]) + { + if (var.variant != null) + { + String alleles = (String) var.variant.getValue("alleles"); + if (alleles != null) + { + for (String base : alleles.split(",")) + { + String codon = base1 + base2 + base; + if (addPeptideVariant(peptide, peptidePos, residue, var, codon)) + { + count++; + } + } + } + } + } + + return count; } /** - * 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. + * Helper method that adds a peptide variant feature, provided the given codon + * translates to a value different to the current residue (is a non-synonymous + * variant). ID and clinical_significance attributes of the dna variant (if + * present) are copied to the new feature. * - * @param from - * @param to + * @param peptide + * @param peptidePos + * @param residue + * @param var + * @param codon + * @return true if a feature was added, else false */ - protected static void transferDbRefs(SequenceI from, SequenceI to) + static boolean addPeptideVariant(SequenceI peptide, int peptidePos, + String residue, DnaVariant var, String codon) { - String cdsAccId = FeatureProperties.getCodingFeature(DBRefSource.EMBL); - DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(from.getDBRefs(), - DBRefSource.CODINGDBS); - if (cdsRefs != null) + /* + * 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 currently ignored here + */ + String trans = codon.contains("-") ? "-" + : (codon.length() > 3 ? null : ResidueProperties + .codonTranslate(codon)); + if (trans != null && !trans.equals(residue)) { - for (DBRefEntry cdsRef : cdsRefs) + String residue3Char = StringUtils + .toSentenceCase(ResidueProperties.aa2Triplet.get(residue)); + String trans3Char = StringUtils + .toSentenceCase(ResidueProperties.aa2Triplet.get(trans)); + String desc = "p." + residue3Char + peptidePos + trans3Char; + // set score to 0f so 'graduated colour' option is offered! JAL-2060 + SequenceFeature sf = new SequenceFeature( + SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, + peptidePos, 0f, "Jalview"); + StringBuilder attributes = new StringBuilder(32); + String id = (String) var.variant.getValue(ID); + if (id != null) + { + if (id.startsWith(SEQUENCE_VARIANT)) + { + id = id.substring(SEQUENCE_VARIANT.length()); + } + sf.setValue(ID, id); + attributes.append(ID).append("=").append(id); + // TODO handle other species variants + StringBuilder link = new StringBuilder(32); + try + { + 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()); + } catch (UnsupportedEncodingException e) + { + // as if + } + } + String clinSig = (String) var.variant + .getValue(CLINICAL_SIGNIFICANCE); + if (clinSig != null) + { + sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); + attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=") + .append(clinSig); + } + peptide.addSequenceFeature(sf); + if (attributes.length() > 0) { - to.addDBRef(new DBRefEntry(cdsRef)); - cdsAccId = cdsRef.getAccessionId(); + sf.setAttributes(attributes.toString()); } + return true; } - if (!to.getName().contains(cdsAccId)) + return false; + } + + /** + * Builds a map whose key is position in the protein sequence, and value is a + * list of the base and all variants for each corresponding codon position + * + * @param dnaSeq + * @param dnaToProtein + * @return + */ + static LinkedHashMap[]> buildDnaVariantsMap( + SequenceI dnaSeq, MapList dnaToProtein) + { + /* + * map from peptide position to all variants of the codon which codes for it + * LinkedHashMap ensures we keep 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) { - to.setName(to.getName() + "|" + cdsAccId); + 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]; + List[] codonVariants = variants.get(peptidePosition); + if (codonVariants == null) + { + codonVariants = new ArrayList[3]; + codonVariants[0] = new ArrayList(); + codonVariants[1] = new ArrayList(); + codonVariants[2] = new ArrayList(); + 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 any variant) for each codon position + */ + for (int codonPos = 0; codonPos < 3; codonPos++) + { + String nucleotide = String.valueOf( + dnaSeq.getCharAt(codon[codonPos] - dnaStart)) + .toUpperCase(); + List codonVariant = codonVariants[codonPos]; + if (codon[codonPos] == dnaCol) + { + if (!codonVariant.isEmpty() + && codonVariant.get(0).variant == null) + { + /* + * already recorded base value, add this variant + */ + codonVariant.get(0).variant = sf; + } + else + { + /* + * add variant with base value + */ + codonVariant.add(new DnaVariant(nucleotide, sf)); + } + } + else if (codonVariant.isEmpty()) + { + /* + * record (possibly non-varying) base value + */ + codonVariant.add(new DnaVariant(nucleotide)); + } + } + } } + return variants; + } + + /** + * Makes an alignment with a copy of the given sequences, adding in any + * non-redundant sequences which are mapped to by the cross-referenced + * sequences. + * + * @param seqs + * @param xrefs + * @return + */ + public static AlignmentI makeCopyAlignment(SequenceI[] seqs, + SequenceI[] xrefs) + { + AlignmentI copy = new Alignment(new Alignment(seqs)); + + /* + * add mappings between sequences to the new alignment + */ + AlignedCodonFrame mappings = new AlignedCodonFrame(); + copy.addCodonFrame(mappings); + for (int i = 0; i < copy.getHeight(); i++) + { + SequenceI from = seqs[i]; + SequenceI to = copy.getSequenceAt(i); + if (to.getDatasetSequence() != null) + { + to = to.getDatasetSequence(); + } + int start = from.getStart(); + int end = from.getEnd(); + MapList map = new MapList(new int[] { start, end }, new int[] { + start, end }, 1, 1); + mappings.addMap(to, from, map); + } + + SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); + if (xrefs != null) + { + for (SequenceI xref : xrefs) + { + DBRefEntry[] dbrefs = xref.getDBRefs(); + if (dbrefs != null) + { + for (DBRefEntry dbref : dbrefs) + { + if (dbref.getMap() == null || dbref.getMap().getTo() == null) + { + continue; + } + SequenceI mappedTo = dbref.getMap().getTo(); + SequenceI match = matcher.findIdMatch(mappedTo); + if (match == null) + { + matcher.add(mappedTo); + copy.addSequence(mappedTo); + } + } + } + } + } + return copy; + } + + /** + * Try to align sequences in 'unaligned' to match the alignment of their + * mapped regions in 'aligned'. For example, could use this to align CDS + * sequences which are mapped to their parent cDNA sequences. + * + * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For + * dna-to-protein or protein-to-dna use alternative methods. + * + * @param unaligned + * sequences to be aligned + * @param aligned + * holds aligned sequences and their mappings + * @return + */ + public static int alignAs(AlignmentI unaligned, AlignmentI aligned) + { + List unmapped = new ArrayList(); + Map> columnMap = buildMappedColumnsMap( + unaligned, aligned, unmapped); + int width = columnMap.size(); + char gap = unaligned.getGapCharacter(); + int realignedCount = 0; + + for (SequenceI seq : unaligned.getSequences()) + { + if (!unmapped.contains(seq)) + { + char[] newSeq = new char[width]; + Arrays.fill(newSeq, gap); + int newCol = 0; + int lastCol = 0; + + /* + * traverse the map to find columns populated + * by our sequence + */ + for (Integer column : columnMap.keySet()) + { + Character c = columnMap.get(column).get(seq); + if (c != null) + { + /* + * sequence has a character at this position + * + */ + newSeq[newCol] = c; + lastCol = newCol; + } + newCol++; + } + + /* + * trim trailing gaps + */ + if (lastCol < width) + { + char[] tmp = new char[lastCol + 1]; + System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1); + newSeq = tmp; + } + seq.setSequence(String.valueOf(newSeq)); + realignedCount++; + } + } + return realignedCount; + } + + /** + * Returns a map whose key is alignment column number (base 1), and whose + * values are a map of sequence characters in that column. + * + * @param unaligned + * @param aligned + * @param unmapped + * @return + */ + static Map> buildMappedColumnsMap( + AlignmentI unaligned, AlignmentI aligned, List unmapped) + { + /* + * Map will hold, for each aligned column position, a map of + * {unalignedSequence, sequenceCharacter} at that position. + * TreeMap keeps the entries in ascending column order. + */ + Map> map = new TreeMap>(); + + /* + * r any sequences that have no mapping so can't be realigned + */ + unmapped.addAll(unaligned.getSequences()); + + List mappings = aligned.getCodonFrames(); + + for (SequenceI seq : unaligned.getSequences()) + { + for (AlignedCodonFrame mapping : mappings) + { + SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned); + if (fromSeq != null) + { + Mapping seqMap = mapping.getMappingBetween(fromSeq, seq); + if (addMappedPositions(seq, fromSeq, seqMap, map)) + { + unmapped.remove(seq); + } + } + } + } + return map; + } + + /** + * Helper method that adds to a map the mapped column positions of a sequence.
+ * For example if aaTT-Tg-gAAA is mapped to TTTAAA then the map should record + * that columns 3,4,6,10,11,12 map to characters T,T,T,A,A,A of the mapped to + * sequence. + * + * @param seq + * the sequence whose column positions we are recording + * @param fromSeq + * a sequence that is mapped to the first sequence + * @param seqMap + * the mapping from 'fromSeq' to 'seq' + * @param map + * a map to add the column positions (in fromSeq) of the mapped + * positions of seq + * @return + */ + static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq, + Mapping seqMap, Map> map) + { + if (seqMap == null) + { + return false; + } + + char[] fromChars = fromSeq.getSequence(); + int toStart = seq.getStart(); + char[] toChars = seq.getSequence(); + + /* + * traverse [start, end, start, end...] ranges in fromSeq + */ + for (int[] fromRange : seqMap.getMap().getFromRanges()) + { + for (int i = 0; i < fromRange.length - 1; i += 2) + { + boolean forward = fromRange[i + 1] >= fromRange[i]; + + /* + * find the range mapped to (sequence positions base 1) + */ + int[] range = seqMap.locateMappedRange(fromRange[i], + fromRange[i + 1]); + if (range == null) + { + System.err.println("Error in mapping " + seqMap + " from " + + fromSeq.getName()); + return false; + } + int fromCol = fromSeq.findIndex(fromRange[i]); + int mappedCharPos = range[0]; + + /* + * walk over the 'from' aligned sequence in forward or reverse + * direction; when a non-gap is found, record the column position + * of the next character of the mapped-to sequence; stop when all + * the characters of the range have been counted + */ + while (mappedCharPos <= range[1]) + { + if (!Comparison.isGap(fromChars[fromCol - 1])) + { + /* + * mapped from sequence has a character in this column + * record the column position for the mapped to character + */ + Map seqsMap = map.get(fromCol); + if (seqsMap == null) + { + seqsMap = new HashMap(); + map.put(fromCol, seqsMap); + } + seqsMap.put(seq, toChars[mappedCharPos - toStart]); + mappedCharPos++; + } + fromCol += (forward ? 1 : -1); + } + } + } + return true; + } + + // strictly temporary hack until proper criteria for aligning protein to cds + // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot + public static boolean looksLikeEnsembl(AlignmentI alignment) + { + for (SequenceI seq : alignment.getSequences()) + { + String name = seq.getName(); + if (!name.startsWith("ENSG") && !name.startsWith("ENST")) + { + return false; + } + } + return true; } }