import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
+import jalview.api.DBRefEntryI;
import jalview.datamodel.AlignedCodon;
import jalview.datamodel.AlignedCodonFrame;
import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping;
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;
*/
public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
{
+ if (protein.isNucleotide() || !dna.isNucleotide())
+ {
+ System.err.println("Wrong alignment type in alignProteinAsDna");
+ return 0;
+ }
List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
protein, dna, unmappedProtein);
}
/**
+ * Realigns the given dna to match the alignment of the protein, using codon
+ * mappings to translate aligned peptide positions to codons.
+ *
+ * @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<AlignedCodonFrame> mappings = protein.getCodonFrames();
+ int alignedCount = 0;
+ for (SequenceI dnaSeq : dna.getSequences())
+ {
+ if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings))
+ {
+ alignedCount++;
+ }
+ }
+ return alignedCount;
+ }
+
+ /**
+ * Helper method to align (if possible) the dna sequence (cds only) to match
+ * the alignment of a mapped protein sequence
+ *
+ * @param dnaSeq
+ * @param protein
+ * @param mappings
+ * @return
+ */
+ static boolean alignCdsSequenceAsProtein(SequenceI dnaSeq,
+ AlignmentI protein, List<AlignedCodonFrame> mappings)
+ {
+ List<AlignedCodonFrame> dnaMappings = MappingUtils
+ .findMappingsForSequence(dnaSeq, mappings);
+ for (AlignedCodonFrame mapping : dnaMappings)
+ {
+ SequenceI peptide = mapping.findAlignedSequence(dnaSeq, protein);
+ if (peptide != null)
+ {
+ Mapping map = mapping.getMappingBetween(dnaSeq, peptide);
+ if (map != null)
+ {
+ /*
+ * traverse peptide adding gaps or codons to new cds sequence
+ */
+ MapList mapList = map.getMap();
+ }
+ }
+ }
+ 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
* 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
* sequences (or null if no mappings are found)
*/
public static AlignmentI makeCdsAlignment(SequenceI[] dna,
- AlignmentI dataset, AlignmentI products)
+ 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<SequenceI> foundSeqs = new ArrayList<SequenceI>();
List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
HashSet<SequenceI> productSeqs = null;
if (products != null)
{
productSeqs = new HashSet<SequenceI>();
- for (SequenceI seq : products.getSequences())
+ 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<AlignedCodonFrame> seqMappings = MappingUtils
- .findMappingsForSequence(seq, mappings);
+ .findMappingsForSequence(dnaSeq, mappings);
for (AlignedCodonFrame mapping : seqMappings)
{
- List<Mapping> mappingsFromSequence = mapping.getMappingsFromSequence(seq);
+ List<Mapping> 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)
}
/*
- * 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 direct lookup instead of double loops to find mappings
-
SequenceI proteinProduct = aMapping.getTo();
-
- /*
- * skip if not mapped to one of a specified set of proteins
- */
if (productSeqs != null && !productSeqs.contains(proteinProduct))
{
continue;
}
- for (AlignedCodonFrame acf : MappingUtils
- .findMappingsForSequence(proteinProduct, mappings))
+ /*
+ * 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)
{
- for (SequenceToSequenceMapping map : acf.getMappings())
+ if (!foundSeqs.contains(cdsSeq))
{
- if (map.getMapping().getMap().getFromRatio() == 3
- && proteinProduct == map.getMapping().getTo()
- && seqDss != map.getFromSeq())
+ foundSeqs.add(cdsSeq);
+ SequenceI derivedSequence = cdsSeq.deriveSequence();
+ cdsSeqs.add(derivedSequence);
+ if (!dataset.getSequences().contains(cdsSeq))
{
- /*
- * 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();
+ dataset.addSequence(cdsSeq);
}
}
- }
- 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;
}
* didn't find mapped CDS sequence - construct it and add
* its dataset sequence to the dataset
*/
- cdsSeq = makeCdsSequence(seq.getDatasetSequence(), aMapping);
+ cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping);
SequenceI cdsSeqDss = cdsSeq.createDatasetSequence();
cdsSeqs.add(cdsSeq);
if (!dataset.getSequences().contains(cdsSeqDss))
*/
List<int[]> 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(cdsSeq, proteinProduct, cdsToProteinMap);
/*
* guard against duplicating the mapping if repeating this action
}
/*
+ * copy protein's dbrefs to CDS sequence
+ * this enables Get Cross-References from CDS alignment
+ */
+ DBRefEntry[] proteinRefs = DBRefUtils.selectDbRefs(false,
+ proteinProduct.getDBRefs());
+ if (proteinRefs != null)
+ {
+ for (DBRefEntry ref : proteinRefs)
+ {
+ DBRefEntry cdsToProteinRef = new DBRefEntry(ref);
+ cdsToProteinRef.setMap(new Mapping(proteinProduct,
+ cdsToProteinMap));
+ cdsSeqDss.addDBRef(cdsToProteinRef);
+ }
+ }
+
+ /*
* add another mapping from original 'from' range to CDS
*/
- AlignedCodonFrame dnaToProteinMapping = new AlignedCodonFrame();
- map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
+ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
+ MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+ cdsRange, 1,
1);
- dnaToProteinMapping.addMap(seq.getDatasetSequence(), cdsSeq, map);
- if (!mappings.contains(dnaToProteinMapping))
+ dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq,
+ 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
+ */
+ DBRefEntryI dnaRef = dnaDss.getSourceDBRef();
+ if (dnaRef != null)
+ {
+ // assuming cds version same as dna ?!?
+ DBRefEntry proteinToCdsRef = new DBRefEntry(dnaRef.getSource(),
+ dnaRef.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, cdsToProteinMap, null,
+ SequenceOntologyI.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<AlignedCodonFrame> mappings,
+ SequenceI dnaSeq, List<AlignedCodonFrame> 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 - 3)
+ {
+ return seqDss;
+ }
+
+ /*
+ * looks like we found the dna-to-protein mapping; search for the
+ * corresponding cds-to-protein mapping
+ */
+ List<AlignedCodonFrame> mappingsToPeptide = MappingUtils
+ .findMappingsForSequence(proteinProduct, mappings);
+ for (AlignedCodonFrame acf : mappingsToPeptide)
+ {
+ for (SequenceToSequenceMapping map : acf.getMappings())
+ {
+ Mapping mapping = map.getMapping();
+ if (mapping != aMapping && mapping.getMap().getFromRatio() == 3
+ && 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<AlignedCodonFrame> 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).
}
}
- 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);
+ // newSeq.setDescription(mapFromId);
+
return newSeq;
}
import static org.testng.AssertJUnit.assertEquals;
import static org.testng.AssertJUnit.assertFalse;
+import static org.testng.AssertJUnit.assertNotNull;
import static org.testng.AssertJUnit.assertNull;
import static org.testng.AssertJUnit.assertSame;
import static org.testng.AssertJUnit.assertTrue;
@Test(groups = { "Functional" })
public void testMakeCdsAlignment()
{
+ /*
+ * scenario:
+ * dna1 --> [4, 6] [10,12] --> pep1
+ * dna2 --> [1, 3] [7, 9] [13,15] --> pep1
+ */
SequenceI dna1 = new Sequence("dna1", "aaaGGGcccTTTaaa");
SequenceI dna2 = new Sequence("dna2", "GGGcccTTTaaaCCC");
SequenceI pep1 = new Sequence("pep1", "GF");
SequenceI pep2 = new Sequence("pep2", "GFP");
+ pep1.addDBRef(new DBRefEntry("UNIPROT", "0", "pep1"));
+ pep2.addDBRef(new DBRefEntry("UNIPROT", "0", "pep2"));
dna1.createDatasetSequence();
dna2.createDatasetSequence();
pep1.createDatasetSequence();
AlignmentI dna = new Alignment(new SequenceI[] { dna1, dna2 });
dna.setDataset(null);
+ /*
+ * need a sourceDbRef if we are to construct dbrefs to the CDS
+ * sequence
+ */
+ DBRefEntry dbref = new DBRefEntry("ENSEMBL", "0", "dna1");
+ dna1.getDatasetSequence().setSourceDBRef(dbref);
+ dbref = new DBRefEntry("ENSEMBL", "0", "dna2");
+ dna2.getDatasetSequence().setSourceDBRef(dbref);
+
+ /*
+ * CDS sequences are 'discovered' from dna-to-protein mappings on the alignment
+ * dataset (e.g. added from dbrefs by CrossRef.findXrefSequences)
+ */
MapList map = new MapList(new int[] { 4, 6, 10, 12 },
new int[] { 1, 2 }, 3, 1);
AlignedCodonFrame acf = new AlignedCodonFrame();
AlignmentI cds = AlignmentUtils.makeCdsAlignment(new SequenceI[] {
dna1, dna2 }, dna.getDataset(), null);
+ /*
+ * verify cds sequences
+ */
assertEquals(2, cds.getSequences().size());
assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());
assertEquals("GGGTTTCCC", cds.getSequenceAt(1).getSequenceAsString());
* verify shared, extended alignment dataset
*/
assertSame(dna.getDataset(), cds.getDataset());
- assertTrue(dna.getDataset().getSequences()
- .contains(cds.getSequenceAt(0).getDatasetSequence()));
- assertTrue(dna.getDataset().getSequences()
- .contains(cds.getSequenceAt(1).getDatasetSequence()));
+ SequenceI cds1Dss = cds.getSequenceAt(0).getDatasetSequence();
+ SequenceI cds2Dss = cds.getSequenceAt(1).getDatasetSequence();
+ assertTrue(dna.getDataset().getSequences().contains(cds1Dss));
+ assertTrue(dna.getDataset().getSequences().contains(cds2Dss));
+
+ /*
+ * verify CDS has a dbref with mapping to peptide
+ */
+ assertNotNull(cds1Dss.getDBRefs());
+ assertEquals(1, cds1Dss.getDBRefs().length);
+ dbref = cds1Dss.getDBRefs()[0];
+ assertEquals("UNIPROT", dbref.getSource());
+ assertEquals("0", dbref.getVersion());
+ assertEquals("pep1", dbref.getAccessionId());
+ assertNotNull(dbref.getMap());
+ assertSame(pep1.getDatasetSequence(), dbref.getMap().getTo());
+ MapList cdsMapping = new MapList(new int[] { 1, 6 },
+ new int[] { 1, 2 }, 3, 1);
+ assertEquals(cdsMapping, dbref.getMap().getMap());
/*
- * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
- * the mappings are on the shared alignment dataset
+ * verify peptide has added a dbref with reverse mapping to CDS
*/
- List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
+ assertNotNull(pep1.getDBRefs());
+ assertEquals(2, pep1.getDBRefs().length);
+ dbref = pep1.getDBRefs()[1];
+ assertEquals("ENSEMBL", dbref.getSource());
+ assertEquals("0", dbref.getVersion());
+ assertEquals("CDS|dna1", dbref.getAccessionId());
+ assertNotNull(dbref.getMap());
+ assertSame(cds1Dss, dbref.getMap().getTo());
+ assertEquals(cdsMapping.getInverse(), dbref.getMap().getMap());
+
/*
+ * Verify mappings from CDS to peptide, cDNA to CDS, and cDNA to peptide
+ * the mappings are on the shared alignment dataset
* 6 mappings, 2*(DNA->CDS), 2*(DNA->Pep), 2*(CDS->Pep)
*/
+ List<AlignedCodonFrame> cdsMappings = cds.getDataset().getCodonFrames();
assertEquals(6, cdsMappings.size());
/*
SearchResults sr = MappingUtils.buildSearchResults(pep1, 1, mappings);
assertEquals(1, sr.getResults().size());
Match m = sr.getResults().get(0);
- assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
+ assertSame(cds1Dss, m.getSequence());
assertEquals(1, m.getStart());
assertEquals(3, m.getEnd());
// map F to TTT
sr = MappingUtils.buildSearchResults(pep1, 2, mappings);
m = sr.getResults().get(0);
- assertSame(cds.getSequenceAt(0).getDatasetSequence(), m.getSequence());
+ assertSame(cds1Dss, m.getSequence());
assertEquals(4, m.getStart());
assertEquals(6, m.getEnd());
sr = MappingUtils.buildSearchResults(pep2, 1, mappings);
assertEquals(1, sr.getResults().size());
m = sr.getResults().get(0);
- assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+ assertSame(cds2Dss, m.getSequence());
assertEquals(1, m.getStart());
assertEquals(3, m.getEnd());
// map F to TTT
sr = MappingUtils.buildSearchResults(pep2, 2, mappings);
m = sr.getResults().get(0);
- assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+ assertSame(cds2Dss, m.getSequence());
assertEquals(4, m.getStart());
assertEquals(6, m.getEnd());
// map P to CCC
sr = MappingUtils.buildSearchResults(pep2, 3, mappings);
m = sr.getResults().get(0);
- assertSame(cds.getSequenceAt(1).getDatasetSequence(), m.getSequence());
+ assertSame(cds2Dss, m.getSequence());
assertEquals(7, m.getStart());
assertEquals(9, m.getEnd());
}
SequenceI cdsSeq = cds.get(0);
assertEquals("GGGTTT", cdsSeq.getSequenceAsString());
// assertEquals("dna1|A12345", cdsSeq.getName());
- assertEquals("dna1|pep1", cdsSeq.getName());
+ assertEquals("CDS|dna1", cdsSeq.getName());
// assertEquals(1, cdsSeq.getDBRefs().length);
// DBRefEntry cdsRef = cdsSeq.getDBRefs()[0];
// assertEquals("EMBLCDS", cdsRef.getSource());
cdsSeq = cds.get(1);
assertEquals("aaaccc", cdsSeq.getSequenceAsString());
// assertEquals("dna1|A12346", cdsSeq.getName());
- assertEquals("dna1|pep2", cdsSeq.getName());
+ assertEquals("CDS|dna1", cdsSeq.getName());
// assertEquals(1, cdsSeq.getDBRefs().length);
// cdsRef = cdsSeq.getDBRefs()[0];
// assertEquals("EMBLCDS", cdsRef.getSource());
cdsSeq = cds.get(2);
assertEquals("aaaTTT", cdsSeq.getSequenceAsString());
// assertEquals("dna1|A12347", cdsSeq.getName());
- assertEquals("dna1|pep3", cdsSeq.getName());
+ assertEquals("CDS|dna1", cdsSeq.getName());
// assertEquals(1, cdsSeq.getDBRefs().length);
// cdsRef = cdsSeq.getDBRefs()[0];
// assertEquals("EMBLCDS", cdsRef.getSource());
* execute method under test to find CDS for EMBL peptides only
*/
AlignmentI cds = AlignmentUtils.makeCdsAlignment(new SequenceI[] {
- dna1, dna2 }, dna.getDataset(), emblPeptides);
+ dna1, dna2 }, dna.getDataset(), emblPeptides.getSequencesArray());
assertEquals(2, cds.getSequences().size());
assertEquals("GGGTTT", cds.getSequenceAt(0).getSequenceAsString());