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.datamodel.Alignment;
import jalview.datamodel.AlignmentAnnotation;
import jalview.datamodel.AlignmentI;
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
Collection<String> types, List<SequenceI> 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;
+ }
}
}
}
* added to the alignment dataset.
*
* @param dna
- * aligned dna sequences
- * @param mappings
- * from dna to protein
- * @param al
+ * aligned nucleotide (dna or cds) sequences
+ * @param 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,
- List<AlignedCodonFrame> mappings, AlignmentI al)
+ AlignmentI dataset, SequenceI[] products)
{
+ if (dataset == null || dataset.getDataset() != null)
+ {
+ throw new IllegalArgumentException(
+ "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
+ }
+ List<SequenceI> foundSeqs = new ArrayList<SequenceI>();
List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
-
- for (SequenceI seq : dna)
+ List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
+ HashSet<SequenceI> productSeqs = null;
+ if (products != null)
{
- AlignedCodonFrame cdsMappings = new AlignedCodonFrame();
+ productSeqs = new HashSet<SequenceI>();
+ for (SequenceI seq : products)
+ {
+ productSeqs.add(seq.getDatasetSequence() == null ? seq : seq
+ .getDatasetSequence());
+ }
+ }
+
+ /*
+ * 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 dnaSeq : dna)
+ {
+ SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
+ : dnaSeq.getDatasetSequence();
+
List<AlignedCodonFrame> seqMappings = MappingUtils
- .findMappingsForSequence(seq, mappings);
- List<AlignedCodonFrame> alignmentMappings = al.getCodonFrames();
+ .findMappingsForSequence(dnaSeq, mappings);
for (AlignedCodonFrame mapping : seqMappings)
{
- for (Mapping aMapping : mapping.getMappingsFromSequence(seq))
+ List<Mapping> mappingsFromSequence = mapping
+ .getMappingsFromSequence(dnaSeq);
+
+ for (Mapping aMapping : mappingsFromSequence)
{
- SequenceI cdsSeq = makeCdsSequence(seq.getDatasetSequence(),
- aMapping);
+ MapList mapList = aMapping.getMap();
+ if (mapList.getFromRatio() == 1)
+ {
+ /*
+ * not a dna-to-protein mapping (likely dna-to-cds)
+ */
+ continue;
+ }
+
+ /*
+ * skip if mapping is not to one of the target set of proteins
+ */
+ SequenceI proteinProduct = aMapping.getTo();
+ if (productSeqs != null && !productSeqs.contains(proteinProduct))
+ {
+ 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)
+ {
+ if (!foundSeqs.contains(cdsSeq))
+ {
+ foundSeqs.add(cdsSeq);
+ 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(dnaSeq.getDatasetSequence(), aMapping);
+ SequenceI cdsSeqDss = cdsSeq.createDatasetSequence();
cdsSeqs.add(cdsSeq);
-
+ if (!dataset.getSequences().contains(cdsSeqDss))
+ {
+ dataset.addSequence(cdsSeqDss);
+ }
+
/*
* add a mapping from CDS to the (unchanged) mapped to range
*/
List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
cdsSeq.getLength() });
- MapList map = new MapList(cdsRange, aMapping.getMap()
- .getToRanges(), aMapping.getMap().getFromRatio(),
- aMapping.getMap().getToRatio());
- cdsMappings.addMap(cdsSeq, aMapping.getTo(), map);
+ MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(),
+ mapList.getFromRatio(), mapList.getToRatio());
+ AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
+ cdsToProteinMapping.addMap(cdsSeq, proteinProduct, cdsToProteinMap);
+
+ /*
+ * guard against duplicating the mapping if repeating this action
+ */
+ if (!mappings.contains(cdsToProteinMapping))
+ {
+ mappings.add(cdsToProteinMapping);
+ }
+
+ /*
+ * 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
*/
- map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
+ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
+ MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+ cdsRange, 1,
1);
- cdsMappings.addMap(seq.getDatasetSequence(), cdsSeq, map);
+ dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq,
+ dnaToCdsMap);
+ if (!mappings.contains(dnaToCdsMapping))
+ {
+ mappings.add(dnaToCdsMapping);
+ }
- alignmentMappings.add(cdsMappings);
+ /*
+ * 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);
}
}
}
+ AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
+ .size()]));
+ 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<AlignedCodonFrame> mappings,
+ SequenceI dnaSeq, List<AlignedCodonFrame> seqMappings,
+ Mapping aMapping)
+ {
/*
- * add CDS seqs to shared dataset
+ * TODO a better dna-cds-protein mapping data representation to allow easy
+ * navigation; until then this clunky looping around lists of mappings
*/
- Alignment dataset = al.getDataset();
- for (SequenceI seq : cdsSeqs)
+ 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)
{
- if (!dataset.getSequences().contains(seq.getDatasetSequence()))
+ 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())
{
- dataset.addSequence(seq.getDatasetSequence());
+ 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;
+ }
+ }
+ }
}
}
- AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
- .size()]));
- cds.setDataset(dataset);
-
- return cds;
+ return null;
}
/**
*
* @param seq
* @param mapping
- * @return
+ * @return CDS sequence (as a dataset sequence)
*/
static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping)
{
}
}
- SequenceI newSeq = new Sequence(seq.getName() + "|"
- + mapping.getTo().getName(), newSeqChars, 1, newPos);
- newSeq.createDatasetSequence();
+ /*
+ * 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;
}
* sort to get sequence features in start position order
* - would be better to store in Sequence as a TreeSet or NCList?
*/
- Arrays.sort(peptide.getSequenceFeatures(),
- new Comparator<SequenceFeature>()
- {
- @Override
- public int compare(SequenceFeature o1, SequenceFeature o2)
+ if (peptide.getSequenceFeatures() != null)
+ {
+ Arrays.sort(peptide.getSequenceFeatures(),
+ new Comparator<SequenceFeature>()
{
- int c = Integer.compare(o1.getBegin(), o2.getBegin());
- return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd())
- : c;
- }
- });
+ @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;
}
*
* @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));
+ 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