X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=34fe221a98d9099fe4e072f403a7bf12b61d580e;hb=37de9310bec3501cbc6381e0c3dcb282fcaad812;hp=6a9dc7bdd1e06ba500d842436d05acabd4873220;hpb=dafd78a3307dd856d328139887036c156b44a203;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 6a9dc7b..34fe221 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -39,6 +39,7 @@ 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; @@ -71,22 +72,26 @@ import java.util.TreeMap; public class AlignmentUtils { + private static final int CODON_LENGTH = 3; + private static final String SEQUENCE_VARIANT = "sequence_variant:"; + private static final String ID = "ID"; /** * A data model to hold the 'normal' base value at a position, and an optional * sequence variant feature */ - static class DnaVariant + static final class DnaVariant { - String base; + final String base; SequenceFeature variant; DnaVariant(String nuc) { base = nuc; + variant = null; } DnaVariant(String nuc, SequenceFeature var) @@ -94,6 +99,11 @@ public class AlignmentUtils base = nuc; variant = var; } + + public String getSource() + { + return variant == null ? null : variant.getFeatureGroup(); + } } /** @@ -426,7 +436,7 @@ public class AlignmentUtils /* * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping) */ - final int mappedLength = 3 * aaSeqChars.length; + final int mappedLength = CODON_LENGTH * aaSeqChars.length; int cdnaLength = cdnaSeqChars.length; int cdnaStart = cdnaSeq.getStart(); int cdnaEnd = cdnaSeq.getEnd(); @@ -438,14 +448,14 @@ public class AlignmentUtils */ if (cdnaLength != mappedLength && cdnaLength > 2) { - String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3) - .toUpperCase(); + String lastCodon = String.valueOf(cdnaSeqChars, + cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(); for (String stop : ResidueProperties.STOP) { if (lastCodon.equals(stop)) { - cdnaEnd -= 3; - cdnaLength -= 3; + cdnaEnd -= CODON_LENGTH; + cdnaLength -= CODON_LENGTH; break; } } @@ -457,12 +467,12 @@ public class AlignmentUtils int startOffset = 0; if (cdnaLength != mappedLength && cdnaLength > 2 - && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase() + && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase() .equals(ResidueProperties.START)) { - startOffset += 3; - cdnaStart += 3; - cdnaLength -= 3; + startOffset += CODON_LENGTH; + cdnaStart += CODON_LENGTH; + cdnaLength -= CODON_LENGTH; } if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars)) @@ -471,7 +481,7 @@ public class AlignmentUtils * protein is translation of dna (+/- start/stop codons) */ MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] - { proteinStart, proteinEnd }, 3, 1); + { proteinStart, proteinEnd }, CODON_LENGTH, 1); return map; } @@ -501,10 +511,9 @@ public class AlignmentUtils int aaPos = 0; int dnaPos = cdnaStart; - for (; dnaPos < cdnaSeqChars.length - 2 - && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++) + for (; dnaPos < cdnaSeqChars.length - 2 && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++) { - String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); final String translated = ResidueProperties.codonTranslate(codon); /* @@ -540,9 +549,9 @@ public class AlignmentUtils { return true; } - if (dnaPos == cdnaSeqChars.length - 3) + if (dnaPos == cdnaSeqChars.length - CODON_LENGTH) { - String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); + String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); if ("STOP".equals(ResidueProperties.codonTranslate(codon))) { return true; @@ -850,6 +859,11 @@ public class AlignmentUtils */ public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) { + if (protein.isNucleotide() || !dna.isNucleotide()) + { + System.err.println("Wrong alignment type in alignProteinAsDna"); + return 0; + } List unmappedProtein = new ArrayList(); Map> alignedCodons = buildCodonColumnsMap( protein, dna, unmappedProtein); @@ -857,6 +871,178 @@ public class AlignmentUtils } /** + * Realigns the given dna to match the alignment of the protein, using codon + * mappings to translate aligned peptide positions to codons. + * + * Always produces a padded CDS alignment. + * + * @param dna + * the alignment whose sequences are realigned by this method + * @param protein + * the protein alignment whose alignment we are 'copying' + * @return the number of sequences that were realigned + */ + public static int alignCdsAsProtein(AlignmentI dna, AlignmentI protein) + { + if (protein.isNucleotide() || !dna.isNucleotide()) + { + System.err.println("Wrong alignment type in alignProteinAsDna"); + return 0; + } + // todo: implement this + List mappings = protein.getCodonFrames(); + int alignedCount = 0; + int width = 0; // alignment width for padding CDS + for (SequenceI dnaSeq : dna.getSequences()) + { + if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings, + dna.getGapCharacter())) + { + alignedCount++; + } + width = Math.max(dnaSeq.getLength(), width); + } + int oldwidth; + int diff; + for (SequenceI dnaSeq : dna.getSequences()) + { + oldwidth = dnaSeq.getLength(); + diff = width - oldwidth; + if (diff > 0) + { + dnaSeq.insertCharAt(oldwidth, diff, dna.getGapCharacter()); + } + } + return alignedCount; + } + + /** + * Helper method to align (if possible) the dna sequence to match the + * alignment of a mapped protein sequence. This is currently limited to + * handling coding sequence only. + * + * @param cdsSeq + * @param protein + * @param mappings + * @param gapChar + * @return + */ + static boolean alignCdsSequenceAsProtein(SequenceI cdsSeq, + AlignmentI protein, List mappings, char gapChar) + { + SequenceI cdsDss = cdsSeq.getDatasetSequence(); + if (cdsDss == null) + { + System.err + .println("alignCdsSequenceAsProtein needs aligned sequence!"); + return false; + } + + List dnaMappings = MappingUtils + .findMappingsForSequence(cdsSeq, mappings); + for (AlignedCodonFrame mapping : dnaMappings) + { + SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); + if (peptide != null) + { + int peptideLength = peptide.getLength(); + Mapping map = mapping.getMappingBetween(cdsSeq, peptide); + if (map != null) + { + MapList mapList = map.getMap(); + if (map.getTo() == peptide.getDatasetSequence()) + { + mapList = mapList.getInverse(); + } + int cdsLength = cdsDss.getLength(); + int mappedFromLength = MappingUtils.getLength(mapList + .getFromRanges()); + int mappedToLength = MappingUtils + .getLength(mapList.getToRanges()); + boolean addStopCodon = (cdsLength == mappedFromLength + * CODON_LENGTH + CODON_LENGTH) + || (peptide.getDatasetSequence().getLength() == mappedFromLength - 1); + if (cdsLength != mappedToLength && !addStopCodon) + { + System.err + .println(String + .format("Can't align cds as protein (length mismatch %d/%d): %s", + cdsLength, mappedToLength, + cdsSeq.getName())); + } + + /* + * pre-fill the aligned cds sequence with gaps + */ + char[] alignedCds = new char[peptideLength * CODON_LENGTH + + (addStopCodon ? CODON_LENGTH : 0)]; + Arrays.fill(alignedCds, gapChar); + + /* + * walk over the aligned peptide sequence and insert mapped + * codons for residues in the aligned cds sequence + */ + char[] alignedPeptide = peptide.getSequence(); + char[] nucleotides = cdsDss.getSequence(); + int copiedBases = 0; + int cdsStart = cdsDss.getStart(); + int proteinPos = peptide.getStart() - 1; + int cdsCol = 0; + for (char residue : alignedPeptide) + { + if (Comparison.isGap(residue)) + { + cdsCol += CODON_LENGTH; + } + else + { + proteinPos++; + int[] codon = mapList.locateInTo(proteinPos, proteinPos); + if (codon == null) + { + // e.g. incomplete start codon, X in peptide + cdsCol += CODON_LENGTH; + } + else + { + for (int j = codon[0]; j <= codon[1]; j++) + { + char mappedBase = nucleotides[j - cdsStart]; + alignedCds[cdsCol++] = mappedBase; + copiedBases++; + } + } + } + } + + /* + * append stop codon if not mapped from protein, + * closing it up to the end of the mapped sequence + */ + if (copiedBases == nucleotides.length - CODON_LENGTH) + { + for (int i = alignedCds.length - 1; i >= 0; i--) + { + if (!Comparison.isGap(alignedCds[i])) + { + cdsCol = i + 1; // gap just after end of sequence + break; + } + } + for (int i = nucleotides.length - CODON_LENGTH; i < nucleotides.length; i++) + { + alignedCds[cdsCol++] = nucleotides[i]; + } + } + cdsSeq.setSequence(new String(alignedCds)); + return true; + } + } + } + return false; + } + + /** * 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 @@ -913,7 +1099,7 @@ public class AlignmentUtils // TODO resolve JAL-2022 so this fudge can be removed int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); - + return alignedCodons; } @@ -1333,15 +1519,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; + } } } } @@ -1400,43 +1590,64 @@ public class AlignmentUtils * added to the alignment dataset. * * @param dna - * aligned dna sequences + * aligned nucleotide (dna or cds) sequences * @param dataset - * - throws error if not given a dataset + * the alignment dataset the sequences belong to + * @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, - AlignmentI dataset) + AlignmentI dataset, SequenceI[] products) { - if (dataset.getDataset() != null) + if (dataset == null || dataset.getDataset() != null) { - throw new Error( + throw new IllegalArgumentException( "IMPLEMENTATION ERROR: dataset.getDataset() must be null!"); } + List foundSeqs = new ArrayList(); List cdsSeqs = new ArrayList(); List mappings = dataset.getCodonFrames(); - + HashSet productSeqs = null; + if (products != null) + { + productSeqs = new HashSet(); + for (SequenceI seq : products) + { + productSeqs.add(seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence()); + } + } /* - * 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 :-( + * Construct CDS sequences from mappings on the alignment dataset. + * The logic is: + * - find the protein product(s) mapped to from each dna sequence + * - if the mapping covers the whole dna sequence (give or take start/stop + * codon), take the dna as the CDS sequence + * - else search dataset mappings for a suitable dna sequence, i.e. one + * whose whole sequence is mapped to the protein + * - if no sequence found, construct one from the dna sequence and mapping + * (and add it to dataset so it is found if this is repeated) */ - for (SequenceI seq : dna) + for (SequenceI dnaSeq : dna) { - SequenceI seqDss = seq.getDatasetSequence() == null ? seq : seq - .getDatasetSequence(); + SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + List seqMappings = MappingUtils - .findMappingsForSequence(seq, mappings); + .findMappingsForSequence(dnaSeq, mappings); for (AlignedCodonFrame mapping : seqMappings) { - List mappingsFromSequence = mapping.getMappingsFromSequence(seq); + List mappingsFromSequence = mapping + .getMappingsFromSequence(dnaSeq); for (Mapping aMapping : mappingsFromSequence) { - if (aMapping.getMap().getFromRatio() == 1) + MapList mapList = aMapping.getMap(); + if (mapList.getFromRatio() == 1) { /* * not a dna-to-protein mapping (likely dna-to-cds) @@ -1445,43 +1656,32 @@ public class AlignmentUtils } /* - * check for an existing CDS sequence i.e. a 3:1 mapping to - * the dna mapping's product + * skip if mapping is not to one of the target set of proteins */ - SequenceI cdsSeq = null; - // TODO better mappings collection data model so we can do - // a table lookup instead of double loops to find mappings SequenceI proteinProduct = aMapping.getTo(); - for (AlignedCodonFrame acf : MappingUtils - .findMappingsForSequence(proteinProduct, mappings)) + if (productSeqs != null && !productSeqs.contains(proteinProduct)) { - 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(); - } - } + continue; } + + /* + * try to locate the CDS from the dataset mappings; + * guard against duplicate results (for the case that protein has + * dbrefs to both dna and cds sequences) + */ + SequenceI cdsSeq = findCdsForProtein(mappings, dnaSeq, + seqMappings, aMapping); 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)) + if (!foundSeqs.contains(cdsSeq)) { - dataset.addSequence(cdsSeq); + foundSeqs.add(cdsSeq); + SequenceI derivedSequence = cdsSeq.deriveSequence(); + cdsSeqs.add(derivedSequence); + if (!dataset.getSequences().contains(cdsSeq)) + { + dataset.addSequence(cdsSeq); + } } continue; } @@ -1490,11 +1690,20 @@ public class AlignmentUtils * 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(); + cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping, + dataset).deriveSequence(); + // cdsSeq has a name constructed as CDS| + // will be either the accession for the coding sequence, + // marked in the /via/ dbref to the protein product accession + // or it will be the original nucleotide accession. + SequenceI cdsSeqDss = cdsSeq.getDatasetSequence(); + cdsSeqs.add(cdsSeq); + if (!dataset.getSequences().contains(cdsSeqDss)) { + // check if this sequence is a newly created one + // so needs adding to the dataset dataset.addSequence(cdsSeqDss); } @@ -1503,11 +1712,12 @@ public class AlignmentUtils */ List cdsRange = Collections.singletonList(new int[] { 1, cdsSeq.getLength() }); - MapList map = new MapList(cdsRange, aMapping.getMap() - .getToRanges(), aMapping.getMap().getFromRatio(), - aMapping.getMap().getToRatio()); + MapList cdsToProteinMap = new MapList(cdsRange, + mapList.getToRanges(), mapList.getFromRatio(), + mapList.getToRatio()); AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); - cdsToProteinMapping.addMap(cdsSeq, proteinProduct, map); + cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct, + cdsToProteinMap); /* * guard against duplicating the mapping if repeating this action @@ -1517,44 +1727,174 @@ public class AlignmentUtils mappings.add(cdsToProteinMapping); } + propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(), + proteinProduct, aMapping); /* * 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)) + AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); + MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), + cdsRange, 1, 1); + dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, + dnaToCdsMap); + if (!mappings.contains(dnaToCdsMapping)) { - mappings.add(dnaToProteinMapping); + mappings.add(dnaToCdsMapping); } + /* + * add DBRef with mapping from protein to CDS + * (this enables Get Cross-References from protein alignment) + * This is tricky because we can't have two DBRefs with the + * same source and accession, so need a different accession for + * the CDS from the dna sequence + */ + + // specific use case: + // Genomic contig ENSCHR:1, contains coding regions for ENSG01, + // ENSG02, ENSG03, with transcripts and products similarly named. + // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01 + + // JBPNote: ?? can't actually create an example that demonstrates we + // need to + // synthesize an xref. + + for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs()) + { + // creates a complementary cross-reference to the source sequence's + // primary reference. + + DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(), + primRef.getSource() + ":" + primRef.getVersion(), + primRef.getAccessionId()); + cdsCrossRef + .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap))); + cdsSeqDss.addDBRef(cdsCrossRef); + + // problem here is that the cross-reference is synthesized - + // cdsSeq.getName() may be like 'CDS|dnaaccession' or + // 'CDS|emblcdsacc' + // assuming cds version same as dna ?!? + + DBRefEntry proteinToCdsRef = new DBRefEntry( + primRef.getSource(), primRef.getVersion(), + cdsSeq.getName()); + // + proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap + .getInverse())); + proteinProduct.addDBRef(proteinToCdsRef); + } /* * transfer any features on dna that overlap the CDS */ - transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS); + transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null, + SequenceOntologyI.CDS); } } } AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs .size()])); - cds.setDataset((Alignment) dataset); + cds.setDataset(dataset); return cds; } /** + * A helper method that finds a CDS sequence in the alignment dataset that is + * mapped to the given protein sequence, and either is, or has a mapping from, + * the given dna sequence. + * + * @param mappings + * set of all mappings on the dataset + * @param dnaSeq + * a dna (or cds) sequence we are searching from + * @param seqMappings + * the set of mappings involving dnaSeq + * @param aMapping + * an initial candidate from seqMappings + * @return + */ + static SequenceI findCdsForProtein(List mappings, + SequenceI dnaSeq, List seqMappings, + Mapping aMapping) + { + /* + * TODO a better dna-cds-protein mapping data representation to allow easy + * navigation; until then this clunky looping around lists of mappings + */ + SequenceI seqDss = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + SequenceI proteinProduct = aMapping.getTo(); + + /* + * is this mapping from the whole dna sequence (i.e. CDS)? + * allowing for possible stop codon on dna but not peptide + */ + int mappedFromLength = MappingUtils.getLength(aMapping.getMap() + .getFromRanges()); + int dnaLength = seqDss.getLength(); + if (mappedFromLength == dnaLength + || mappedFromLength == dnaLength - CODON_LENGTH) + { + return seqDss; + } + + /* + * looks like we found the dna-to-protein mapping; search for the + * corresponding cds-to-protein mapping + */ + List mappingsToPeptide = MappingUtils + .findMappingsForSequence(proteinProduct, mappings); + for (AlignedCodonFrame acf : mappingsToPeptide) + { + for (SequenceToSequenceMapping map : acf.getMappings()) + { + Mapping mapping = map.getMapping(); + if (mapping != aMapping + && mapping.getMap().getFromRatio() == CODON_LENGTH + && proteinProduct == mapping.getTo() + && seqDss != map.getFromSeq()) + { + mappedFromLength = MappingUtils.getLength(mapping.getMap() + .getFromRanges()); + if (mappedFromLength == map.getFromSeq().getLength()) + { + /* + * found a 3:1 mapping to the protein product which covers + * the whole dna sequence i.e. is from CDS; finally check it + * is from the dna start sequence + */ + SequenceI cdsSeq = map.getFromSeq(); + List dnaToCdsMaps = MappingUtils + .findMappingsForSequence(cdsSeq, seqMappings); + if (!dnaToCdsMaps.isEmpty()) + { + return cdsSeq; + } + } + } + } + } + return null; + } + + /** * 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 seq * @param mapping + * @param dataset + * - existing dataset. We check for sequences that look like the CDS + * we are about to construct, if one exists already, then we will + * just return that one. * @return CDS sequence (as a dataset sequence) */ - static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) + static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, + AlignmentI dataset) { char[] seqChars = seq.getSequence(); List fromRanges = mapping.getMap().getFromRanges(); @@ -1582,12 +1922,131 @@ public class AlignmentUtils } } - SequenceI newSeq = new Sequence(seq.getName() + "|" - + mapping.getTo().getName(), newSeqChars, 1, newPos); + /* + * assign 'from id' held in the mapping if set (e.g. EMBL protein_id), + * else generate a sequence name + */ + String mapFromId = mapping.getMappedFromId(); + String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName()); + SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos); + if (dataset != null) + { + SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName()); + if (matches != null) + { + boolean matched = false; + for (SequenceI mtch : matches) + { + if (mtch.getStart() != newSeq.getStart()) + { + continue; + } + if (mtch.getEnd() != newSeq.getEnd()) + { + continue; + } + if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence())) + { + continue; + } + if (!matched) + { + matched = true; + newSeq = mtch; + } + else + { + System.err + .println("JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + + mtch.toString()); + } + } + } + } + // newSeq.setDescription(mapFromId); + return newSeq; } /** + * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to + * the given mapping. + * + * @param cdsSeq + * @param contig + * @param mapping + * @return list of DBRefEntrys added. + */ + public static List propagateDBRefsToCDS(SequenceI cdsSeq, + SequenceI contig, SequenceI proteinProduct, Mapping mapping) + { + + // gather direct refs from contig congrent with mapping + List direct = new ArrayList(); + HashSet directSources = new HashSet(); + if (contig.getDBRefs() != null) + { + for (DBRefEntry dbr : contig.getDBRefs()) + { + if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap()) + { + MapList map = dbr.getMap().getMap(); + // check if map is the CDS mapping + if (mapping.getMap().equals(map)) + { + direct.add(dbr); + directSources.add(dbr.getSource()); + } + } + } + } + DBRefEntry[] onSource = DBRefUtils.selectRefs( + proteinProduct.getDBRefs(), + directSources.toArray(new String[0])); + List propagated = new ArrayList(); + + // and generate appropriate mappings + for (DBRefEntry cdsref : direct) + { + // clone maplist and mapping + MapList cdsposmap = new MapList(Arrays.asList(new int[][] { new int[] + { cdsSeq.getStart(), cdsSeq.getEnd() } }), cdsref.getMap().getMap() + .getToRanges(), 3, 1); + Mapping cdsmap = new Mapping(cdsref.getMap().getTo(), cdsref.getMap() + .getMap()); + + // create dbref + DBRefEntry newref = new DBRefEntry(cdsref.getSource(), + cdsref.getVersion(), cdsref.getAccessionId(), new Mapping( + cdsmap.getTo(), cdsposmap)); + + // and see if we can map to the protein product for this mapping. + // onSource is the filtered set of accessions on protein that we are + // tranferring, so we assume accession is the same. + if (cdsmap.getTo() == null && onSource != null) + { + List sourceRefs = DBRefUtils.searchRefs(onSource, + cdsref.getAccessionId()); + if (sourceRefs != null) + { + for (DBRefEntry srcref : sourceRefs) + { + if (srcref.getSource().equalsIgnoreCase(cdsref.getSource())) + { + // we have found a complementary dbref on the protein product, so + // update mapping's getTo + newref.getMap().setTo(proteinProduct); + } + } + } + } + cdsSeq.addDBRef(newref); + propagated.add(newref); + } + return propagated; + } + + /** * 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. @@ -1716,7 +2175,7 @@ public class AlignmentUtils /* * dna length should map to protein (or protein plus stop codon) */ - int codesForResidues = mappedDnaLength / 3; + int codesForResidues = mappedDnaLength / CODON_LENGTH; if (codesForResidues == (proteinLength + 1)) { // assuming extra codon is for STOP and not in peptide @@ -1725,7 +2184,7 @@ public class AlignmentUtils if (codesForResidues == proteinLength) { proteinRange.add(new int[] { proteinStart, proteinEnd }); - return new MapList(ranges, proteinRange, 3, 1); + return new MapList(ranges, proteinRange, CODON_LENGTH, 1); } return null; } @@ -2001,7 +2460,7 @@ public class AlignmentUtils * are currently ignored here */ String trans = codon.contains("-") ? "-" - : (codon.length() > 3 ? null : ResidueProperties + : (codon.length() > CODON_LENGTH ? null : ResidueProperties .codonTranslate(codon)); if (trans != null && !trans.equals(residue)) { @@ -2013,7 +2472,7 @@ public class AlignmentUtils // set score to 0f so 'graduated colour' option is offered! JAL-2060 SequenceFeature sf = new SequenceFeature( SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, - peptidePos, 0f, "Jalview"); + peptidePos, 0f, var.getSource()); StringBuilder attributes = new StringBuilder(32); String id = (String) var.variant.getValue(ID); if (id != null) @@ -2024,11 +2483,13 @@ public class AlignmentUtils } sf.setValue(ID, id); attributes.append(ID).append("=").append(id); - // TODO handle other species variants + // TODO handle other species variants JAL-2064 StringBuilder link = new StringBuilder(32); try { - link.append(desc).append(" ").append(id) + link.append(desc) + .append(" ") + .append(id) .append("|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=") .append(URLEncoder.encode(id, "UTF-8")); sf.addLink(link.toString()); @@ -2037,8 +2498,7 @@ public class AlignmentUtils // as if } } - String clinSig = (String) var.variant - .getValue(CLINICAL_SIGNIFICANCE); + String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE); if (clinSig != null) { sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); @@ -2063,6 +2523,7 @@ public class AlignmentUtils * @param dnaToProtein * @return */ + @SuppressWarnings("unchecked") static LinkedHashMap[]> buildDnaVariantsMap( SequenceI dnaSeq, MapList dnaToProtein) { @@ -2106,7 +2567,7 @@ public class AlignmentUtils List[] codonVariants = variants.get(peptidePosition); if (codonVariants == null) { - codonVariants = new ArrayList[3]; + codonVariants = new ArrayList[CODON_LENGTH]; codonVariants[0] = new ArrayList(); codonVariants[1] = new ArrayList(); codonVariants[2] = new ArrayList(); @@ -2140,7 +2601,7 @@ public class AlignmentUtils /* * save nucleotide (and any variant) for each codon position */ - for (int codonPos = 0; codonPos < 3; codonPos++) + for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++) { String nucleotide = String.valueOf( dnaSeq.getCharAt(codon[codonPos] - dnaStart)) @@ -2184,33 +2645,16 @@ public class AlignmentUtils * * @param seqs * @param xrefs + * @param dataset + * the alignment dataset shared by the new copy * @return */ public static AlignmentI makeCopyAlignment(SequenceI[] seqs, - SequenceI[] xrefs) + SequenceI[] xrefs, AlignmentI dataset) { 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); - } - + copy.setDataset(dataset); + boolean isProtein = !copy.isNucleotide(); SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) { @@ -2221,7 +2665,8 @@ public class AlignmentUtils { for (DBRefEntry dbref : dbrefs) { - if (dbref.getMap() == null || dbref.getMap().getTo() == null) + if (dbref.getMap() == null || dbref.getMap().getTo() == null + || dbref.getMap().getTo().isProtein() != isProtein) { continue; } @@ -2255,19 +2700,32 @@ public class AlignmentUtils */ public static int alignAs(AlignmentI unaligned, AlignmentI aligned) { + /* + * easy case - aligning a copy of aligned sequences + */ + if (alignAsSameSequences(unaligned, aligned)) + { + return unaligned.getHeight(); + } + + /* + * fancy case - aligning via mappings between sequences + */ List unmapped = new ArrayList(); Map> columnMap = buildMappedColumnsMap( unaligned, aligned, unmapped); int width = columnMap.size(); char gap = unaligned.getGapCharacter(); int realignedCount = 0; + // TODO: verify this loop scales sensibly for very wide/high alignments for (SequenceI seq : unaligned.getSequences()) { if (!unmapped.contains(seq)) { char[] newSeq = new char[width]; - Arrays.fill(newSeq, gap); + Arrays.fill(newSeq, gap); // JBPComment - doubt this is faster than the + // Integer iteration below int newCol = 0; int lastCol = 0; @@ -2289,7 +2747,7 @@ public class AlignmentUtils } newCol++; } - + /* * trim trailing gaps */ @@ -2299,6 +2757,7 @@ public class AlignmentUtils System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1); newSeq = tmp; } + // TODO: optimise SequenceI to avoid char[]->String->char[] seq.setSequence(String.valueOf(newSeq)); realignedCount++; } @@ -2307,6 +2766,72 @@ public class AlignmentUtils } /** + * If unaligned and aligned sequences share the same dataset sequences, then + * simply copies the aligned sequences to the unaligned sequences and returns + * true; else returns false + * + * @param unaligned + * - sequences to be aligned based on aligned + * @param aligned + * - 'guide' alignment containing sequences derived from same dataset + * as unaligned + * @return + */ + static boolean alignAsSameSequences(AlignmentI unaligned, + AlignmentI aligned) + { + if (aligned.getDataset() == null || unaligned.getDataset() == null) + { + return false; // should only pass alignments with datasets here + } + + // map from dataset sequence to alignment sequence(s) + Map> alignedDatasets = new HashMap>(); + for (SequenceI seq : aligned.getSequences()) + { + SequenceI ds = seq.getDatasetSequence(); + if (alignedDatasets.get(ds) == null) + { + alignedDatasets.put(ds, new ArrayList()); + } + alignedDatasets.get(ds).add(seq); + } + + /* + * first pass - check whether all sequences to be aligned share a dataset + * sequence with an aligned sequence + */ + for (SequenceI seq : unaligned.getSequences()) + { + if (!alignedDatasets.containsKey(seq.getDatasetSequence())) + { + return false; + } + } + + /* + * second pass - copy aligned sequences; + * heuristic rule: pair off sequences in order for the case where + * more than one shares the same dataset sequence + */ + for (SequenceI seq : unaligned.getSequences()) + { + List alignedSequences = alignedDatasets.get(seq + .getDatasetSequence()); + // TODO: getSequenceAsString() will be deprecated in the future + // TODO: need to leave to SequenceI implementor to update gaps + seq.setSequence(alignedSequences.get(0).getSequenceAsString()); + if (alignedSequences.size() > 0) + { + // pop off aligned sequences (except the last one) + alignedSequences.remove(0); + } + } + + return true; + } + + /** * Returns a map whose key is alignment column number (base 1), and whose * values are a map of sequence characters in that column. * @@ -2320,13 +2845,13 @@ public class AlignmentUtils { /* * Map will hold, for each aligned column position, a map of - * {unalignedSequence, sequenceCharacter} at that position. + * {unalignedSequence, characterPerSequence} 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 + * record any sequences that have no mapping so can't be realigned */ unmapped.addAll(unaligned.getSequences()); @@ -2375,6 +2900,15 @@ public class AlignmentUtils 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(); @@ -2408,7 +2942,8 @@ public class AlignmentUtils * of the next character of the mapped-to sequence; stop when all * the characters of the range have been counted */ - while (mappedCharPos <= range[1]) + while (mappedCharPos <= range[1] && fromCol <= fromChars.length + && fromCol >= 0) { if (!Comparison.isGap(fromChars[fromCol - 1])) {