X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=aa7cb181ee45b97937bda7c5fee15f2503783381;hb=768a4920722ec46f7c16e35780391e39a76d02ff;hp=db69823050fc191d5ebddb0cbd3baeba0a04e841;hpb=6c52cc0b81ae3abdc3c5f6f88a23364a0246351a;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index db69823..aa7cb18 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,17 +20,17 @@ */ package jalview.analysis; +import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE; + import jalview.datamodel.AlignedCodon; import jalview.datamodel.AlignedCodonFrame; +import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping; 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; @@ -39,11 +39,12 @@ 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.io.UnsupportedEncodingException; +import java.net.URLEncoder; import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; @@ -70,6 +71,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 @@ -564,7 +590,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; @@ -813,148 +839,6 @@ 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. - * - * @param sequences - * @param gapCharacter - * @param mappings - * @return - */ - public static List getAlignedTranslation( - List sequences, char gapCharacter, - Set mappings) - { - List alignedSeqs = new ArrayList(); - - for (SequenceI seq : sequences) - { - List mapped = getAlignedTranslation(seq, gapCharacter, - mappings); - alignedSeqs.addAll(mapped); - } - return alignedSeqs; - } - - /** - * 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. - * - * @param seq - * @param gapCharacter - * @param mappings - * @return - */ - protected static List getAlignedTranslation(SequenceI seq, - char gapCharacter, Set mappings) - { - List result = new ArrayList(); - for (AlignedCodonFrame mapping : mappings) - { - if (mapping.involvesSequence(seq)) - { - SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping); - if (mapped != null) - { - result.add(mapped); - } - } - } - return result; - } - - /** - * Returns the translation of 'seq' (as held in the mapping) with - * corresponding alignment (gaps). - * - * @param seq - * @param gapCharacter - * @param mapping - * @return - */ - protected static SequenceI getAlignedTranslation(SequenceI seq, - char gapCharacter, AlignedCodonFrame mapping) - { - 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)); - - 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("", ""); - - for (char c : seq.getSequence()) - { - if (c == gapCharacter) - { - gapWidth++; - if (gapWidth >= fromRatio) - { - newseq.append(gap); - gapWidth = 0; - } - } - else - { - phrase[phraseOffset++] = residueNo + 1; - if (phraseOffset == fromRatio) - { - /* - * 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. - */ - SearchResults sr = new SearchResults(); - for (int pos : phrase) - { - mapping.markMappedRegion(seq, pos, sr); - } - newseq.append(sr.getCharacters()); - if (first) - { - 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); - } - phraseOffset = 0; - } - residueNo++; - } - } - 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. * @@ -1011,8 +895,7 @@ public class AlignmentUtils { for (AlignedCodonFrame mapping : mappings) { - SequenceI prot = mapping.findAlignedSequence( - dnaSeq.getDatasetSequence(), protein); + SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein); if (prot != null) { Mapping seqMap = mapping.getMappingForSequence(dnaSeq); @@ -1027,6 +910,7 @@ public class AlignmentUtils * 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); @@ -1449,15 +1333,19 @@ public class AlignmentUtils Collection types, List forSequences, boolean anyType, boolean doShow) { - for (AlignmentAnnotation aa : al.getAlignmentAnnotation()) + AlignmentAnnotation[] anns = al.getAlignmentAnnotation(); + if (anns != null) { - if (anyType || types.contains(aa.label)) + for (AlignmentAnnotation aa : anns) { - if ((aa.sequenceRef != null) - && (forSequences == null || forSequences - .contains(aa.sequenceRef))) + if (anyType || types.contains(aa.label)) { - aa.visible = doShow; + if ((aa.sequenceRef != null) + && (forSequences == null || forSequences + .contains(aa.sequenceRef))) + { + aa.visible = doShow; + } } } } @@ -1510,243 +1398,220 @@ public class AlignmentUtils /** * 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. + * 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 - * @param al + * @param dataset + * - throws error if not given a dataset + * @param products + * (optional) to restrict results to CDS that map to specified + * protein products * @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) + AlignmentI dataset, AlignmentI products) { - 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) + if (dataset.getDataset() != null) { - final SequenceI ds = dnaSeq.getDatasetSequence(); - List seqMappings = MappingUtils - .findMappingsForSequence(ds, mappings); - for (AlignedCodonFrame acf : seqMappings) + throw new Error( + "IMPLEMENTATION ERROR: dataset.getDataset() must be null!"); + } + List cdsSeqs = new ArrayList(); + List mappings = dataset.getCodonFrames(); + HashSet productSeqs = null; + if (products != null) + { + productSeqs = new HashSet(); + for (SequenceI seq : products.getSequences()) { - AlignedCodonFrame newMapping = new AlignedCodonFrame(); - final List mappedCds = makeCdsSequences(dnaSeq, acf, - cdsColumns, newMapping, gap); - if (!mappedCds.isEmpty()) - { - cdsSequences.addAll(mappedCds); - newMappings.add(newMapping); - } + productSeqs.add(seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence()); } } - AlignmentI newAl = new Alignment( - cdsSequences.toArray(new SequenceI[cdsSequences.size()])); /* - * add new sequences to the shared dataset, set it on the new alignment + * 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 :-( */ - List dsseqs = al.getDataset().getSequences(); - for (SequenceI seq : newAl.getSequences()) + for (SequenceI seq : dna) { - if (!dsseqs.contains(seq.getDatasetSequence())) + SequenceI seqDss = seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence(); + List seqMappings = MappingUtils + .findMappingsForSequence(seq, mappings); + for (AlignedCodonFrame mapping : seqMappings) { - dsseqs.add(seq.getDatasetSequence()); - } - } - newAl.setDataset(al.getDataset()); + List mappingsFromSequence = mapping.getMappingsFromSequence(seq); - /* - * Replace the old mappings with the new ones - */ - mappings.clear(); - mappings.addAll(newMappings); + for (Mapping aMapping : mappingsFromSequence) + { + if (aMapping.getMap().getFromRatio() == 1) + { + /* + * not a dna-to-protein mapping (likely dna-to-cds) + */ + continue; + } - return newAl; - } + /* + * check for an existing CDS sequence i.e. a 3:1 mapping to + * the dna mapping's product + */ + SequenceI cdsSeq = null; - /** - * 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! + // TODO better mappings collection data model so we can do + // a direct lookup instead of double loops to find mappings - List result = new ArrayList(); - for (SequenceI seq : seqs) - { - result.addAll(findCdsColumns(seq)); - } + SequenceI proteinProduct = aMapping.getTo(); - /* - * 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); + /* + * skip if not mapped to one of a specified set of proteins + */ + if (productSeqs != null && !productSeqs.contains(proteinProduct)) + { + continue; + } - return result; - } + for (AlignedCodonFrame acf : MappingUtils + .findMappingsForSequence(proteinProduct, mappings)) + { + for (SequenceToSequenceMapping map : acf.getMappings()) + { + if (map.getMapping().getMap().getFromRatio() == 3 + && proteinProduct == map.getMapping().getTo() + && seqDss != map.getFromSeq()) + { + /* + * found a 3:1 mapping to the protein product which is not + * from the dna sequence...assume it is from the CDS sequence + * TODO mappings data model that brings together related + * dna-cds-protein mappings in one object + */ + cdsSeq = map.getFromSeq(); + } + } + } + if (cdsSeq != null) + { + /* + * mappings are always to dataset sequences so create an aligned + * sequence to own it; add the dataset sequence to the dataset + */ + SequenceI derivedSequence = cdsSeq.deriveSequence(); + cdsSeqs.add(derivedSequence); + if (!dataset.getSequences().contains(cdsSeq)) + { + dataset.addSequence(cdsSeq); + } + continue; + } - 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; - } + /* + * didn't find mapped CDS sequence - construct it and add + * its dataset sequence to the dataset + */ + cdsSeq = makeCdsSequence(seq.getDatasetSequence(), aMapping); + SequenceI cdsSeqDss = cdsSeq.createDatasetSequence(); + cdsSeqs.add(cdsSeq); + if (!dataset.getSequences().contains(cdsSeqDss)) + { + dataset.addSequence(cdsSeqDss); + } - /** - * 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; - } + /* + * 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()); + AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); + cdsToProteinMapping.addMap(cdsSeq, proteinProduct, map); - /** - * 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 + * guard against duplicating the mapping if repeating this action */ - for (int[] cdsRange : map.getMap().getFromRanges()) + if (!mappings.contains(cdsToProteinMapping)) { - int startPos = cdsRange[0]; - int endPos = cdsRange[1]; - int startCol = seq.findIndex(startPos); - int endCol = seq.findIndex(endPos); - columns.add(new int[] { startCol, endCol }); + mappings.add(cdsToProteinMapping); } + + /* + * add another mapping from original 'from' range to CDS + */ + AlignedCodonFrame dnaToProteinMapping = new AlignedCodonFrame(); + map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1, + 1); + dnaToProteinMapping.addMap(seq.getDatasetSequence(), cdsSeq, map); + if (!mappings.contains(dnaToProteinMapping)) + { + mappings.add(dnaToProteinMapping); + } + + + /* + * transfer any features on dna that overlap the CDS + */ + transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS); } } - result.add(columns); } - return result; + + AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs + .size()])); + cds.setDataset(dataset); + + 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 aligned sequence + * @param seq * @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 + * @return CDS sequence (as a dataset sequence) */ - protected static List makeCdsSequences(SequenceI dnaSeq, - AlignedCodonFrame mapping, List ungappedCdsColumns, - AlignedCodonFrame newMappings, char gapChar) + 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, - 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); + 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); + return newSeq; } /** @@ -1846,128 +1711,6 @@ public class AlignmentUtils } /** - * 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. * @@ -1980,11 +1723,11 @@ public class AlignmentUtils { 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 @@ -1996,7 +1739,7 @@ public class AlignmentUtils proteinLength--; } List proteinRange = new ArrayList(); - + /* * dna length should map to protein (or protein plus stop codon) */ @@ -2017,7 +1760,9 @@ public class AlignmentUtils /** * 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) + * 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 * @return @@ -2030,7 +1775,10 @@ public class AlignmentUtils { return result; } + SequenceOntologyI so = SequenceOntologyFactory.getInstance(); + int startPhase = 0; + for (SequenceFeature sf : sfs) { /* @@ -2039,7 +1787,8 @@ public class AlignmentUtils if (so.isA(sf.getType(), SequenceOntologyI.CDS)) { int phase = 0; - try { + try + { phase = Integer.parseInt(sf.getPhase()); } catch (NumberFormatException e) { @@ -2053,16 +1802,44 @@ public class AlignmentUtils int end = sf.getEnd(); if (result.isEmpty()) { - // TODO JAL-2022 support start phase > 0 begin += phase; if (begin > end) { - continue; // shouldn't happen? + // shouldn't happen! + System.err + .println("Error: start phase extends beyond start CDS in " + + dnaSeq.getName()); } } result.add(new int[] { begin, end }); } } + + /* + * remove 'startPhase' positions (usually 0) from the first range + * so we begin at the start of a complete codon + */ + if (!result.isEmpty()) + { + // 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) + { + return Integer.compare(o1[0], o2[0]); + } + }); return result; } @@ -2087,82 +1864,252 @@ public class AlignmentUtils { peptide = peptide.getDatasetSequence(); } - - transferFeatures(dnaSeq, peptide, dnaToProtein, - SequenceOntologyI.EXON); - - LinkedHashMap variants = buildDnaVariantsMap( + + transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON); + + /* + * 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 + */ + + /* + * 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()) + 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()) + 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 desc = StringUtils.listToDelimitedString(peptideVariants, - ", "); - SequenceFeature sf = new SequenceFeature( - SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, - peptidePos, 0f, null); - peptide.addSequenceFeature(sf); - count++; + 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++; + } + } + } } } - + /* - * ugly sort to get sequence features in start position order - * - would be better to store in Sequence as a TreeSet instead? + * variants in second codon base */ - Arrays.sort(peptide.getSequenceFeatures(), - new Comparator() + 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)) { - @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; - } - }); + 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; } /** - * Builds a map whose key is position in the protein sequence, and value is an - * array of all variants for the coding codon positions + * 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 peptide + * @param peptidePos + * @param residue + * @param var + * @param codon + * @return true if a feature was added, else false + */ + static boolean addPeptideVariant(SequenceI peptide, int peptidePos, + String residue, DnaVariant var, String codon) + { + /* + * 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)) + { + 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) + { + sf.setAttributes(attributes.toString()); + } + return true; + } + 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( + 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 + * 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(); + 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 */ @@ -2183,13 +2130,16 @@ public class AlignmentUtils continue; } int peptidePosition = mapsTo[0]; - String[][] codonVariants = variants.get(peptidePosition); + List[] codonVariants = variants.get(peptidePosition); if (codonVariants == null) { - codonVariants = new String[3][]; + 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 */ @@ -2204,7 +2154,7 @@ public class AlignmentUtils { 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] */ @@ -2213,33 +2163,40 @@ public class AlignmentUtils peptidePosition, peptidePosition)); lastPeptidePostion = peptidePosition; lastCodon = codon; - + /* - * save nucleotide (and this variant) for each codon position + * 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(); - if (codonVariants[codonPos] == null) + List codonVariant = codonVariants[codonPos]; + if (codon[codonPos] == dnaCol) { - /* - * record current dna base - */ - codonVariants[codonPos] = new String[] { nucleotide }; + 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)); + } } - if (codon[codonPos] == dnaCol) + else if (codonVariant.isEmpty()) { /* - * add alleles to dna base (and any previously found alleles) + * record (possibly non-varying) base value */ - 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; + codonVariant.add(new DnaVariant(nucleotide)); } } } @@ -2248,67 +2205,265 @@ public class AlignmentUtils } /** - * Returns a sorted, non-redundant list of all peptide translations generated - * by the given dna variants, excluding the current residue value + * 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 codonVariants - * an array of base values (acgtACGT) for codon positions 1, 2, 3 - * @param residue - * the current residue translation + * @param seqs + * @param xrefs + * @param dataset + * the alignment dataset shared by the new copy * @return */ - static List computePeptideVariants( - String[][] codonVariants, String residue) + public static AlignmentI makeCopyAlignment(SequenceI[] seqs, + SequenceI[] xrefs, AlignmentI dataset) { - List result = new ArrayList(); - for (String base1 : codonVariants[0]) + AlignmentI copy = new Alignment(new Alignment(seqs)); + copy.setDataset(dataset); + + SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); + if (xrefs != null) { - for (String base2 : codonVariants[1]) + for (SequenceI xref : xrefs) { - for (String base3 : codonVariants[2]) + DBRefEntry[] dbrefs = xref.getDBRefs(); + if (dbrefs != null) { - 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)) + for (DBRefEntry dbref : dbrefs) { - result.add(peptide); + 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) + { /* - * sort alphabetically with STOP at the end + * Map will hold, for each aligned column position, a map of + * {unalignedSequence, characterPerSequence} at that position. + * TreeMap keeps the entries in ascending column order. */ - Collections.sort(result, new Comparator() + Map> map = new TreeMap>(); + + /* + * record any sequences that have no mapping so can't be realigned + */ + unmapped.addAll(unaligned.getSequences()); + + List mappings = aligned.getCodonFrames(); + + for (SequenceI seq : unaligned.getSequences()) { - - @Override - public int compare(String o1, String o2) + for (AlignedCodonFrame mapping : mappings) { - if ("STOP".equals(o1)) + SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned); + if (fromSeq != null) { - return 1; + Mapping seqMap = mapping.getMappingBetween(fromSeq, seq); + if (addMappedPositions(seq, fromSeq, seqMap, map)) + { + unmapped.remove(seq); + } } - else if ("STOP".equals(o2)) + } + } + 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; + } + + /* + * invert mapping if it is from unaligned to aligned sequence + */ + if (seqMap.getTo() == fromSeq.getDatasetSequence()) + { + seqMap = new Mapping(seq.getDatasetSequence(), seqMap.getMap() + .getInverse()); + } + + 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) { - return -1; + System.err.println("Error in mapping " + seqMap + " from " + + fromSeq.getName()); + return false; } - else + 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] && fromCol <= fromChars.length + && fromCol >= 0) { - return o1.compareTo(o2); + 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 result; + } + 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; } }