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;
}