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,
+ dna.getGapCharacter()))
+ {
+ alignedCount++;
+ }
+ }
+ 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<AlignedCodonFrame> mappings, char gapChar)
+ {
+ SequenceI cdsDss = cdsSeq.getDatasetSequence();
+ if (cdsDss == null)
+ {
+ System.err
+ .println("alignCdsSequenceAsProtein needs aligned sequence!");
+ return false;
+ }
+
+ List<AlignedCodonFrame> dnaMappings = MappingUtils
+ .findMappingsForSequence(cdsSeq, mappings);
+ for (AlignedCodonFrame mapping : dnaMappings)
+ {
+ SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein);
+ int peptideLength = peptide.getLength();
+ if (peptide != null)
+ {
+ 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 * 3 + 3)
+ || (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 * 3
+ + (addStopCodon ? 3 : 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 += 3;
+ }
+ else
+ {
+ proteinPos++;
+ int[] codon = mapList.locateInTo(proteinPos, proteinPos);
+ if (codon == null)
+ {
+ // e.g. incomplete start codon, X in peptide
+ cdsCol += 3;
+ }
+ 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 - 3)
+ {
+ 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 - 3; 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
* 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;
}
*
* @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);
SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
if (xrefs != null)
*/
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<SequenceI> unmapped = new ArrayList<SequenceI>();
Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
unaligned, aligned, unmapped);
}
/**
+ * 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
+ * @param aligned
+ * @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<SequenceI, SequenceI> alignedDatasets = new HashMap<SequenceI, SequenceI>();
+ for (SequenceI seq : aligned.getSequences())
+ {
+ alignedDatasets.put(seq.getDatasetSequence(), 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
+ */
+ for (SequenceI seq : unaligned.getSequences())
+ {
+ SequenceI alignedSequence = alignedDatasets.get(seq
+ .getDatasetSequence());
+ seq.setSequence(alignedSequence.getSequenceAsString());
+ }
+
+ 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.
*
{
/*
* 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<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
/*
- * 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());
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();
* 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]))
{
git://source.jalview.org / jalview.git / blobdiff