cdsRange = Collections.singletonList(new int[] { 1,
+ cdsSeq.getLength() });
+ 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
+ */
+ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
+ MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+ cdsRange, 1,
+ 1);
+ dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq,
+ dnaToCdsMap);
+ if (!mappings.contains(dnaToCdsMapping))
+ {
+ 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(dnaSeq, cdsSeq, cdsToProteinMap, null,
+ SequenceOntologyI.CDS);
}
}
}
- AlignmentI al = new Alignment(
- cdsSequences.toArray(new SequenceI[cdsSequences.size()]));
- al.setDataset(null);
- /*
- * Replace the old mappings with the new ones
- */
- mappings.clear();
- mappings.addAll(newMappings);
+ AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
+ .size()]));
+ cds.setDataset(dataset);
- return al;
+ return cds;
}
/**
- * Helper method to make cds-only sequences and populate their mappings to
- * protein products
- *
- * For example, if ggCCaTTcGAg has mappings [3, 4, 6, 7, 9, 10] to protein
- * then generate a sequence CCTTGA with mapping [1, 6] to the same protein
- * residues
- *
- * Typically eukaryotic dna will include cds encoding for a single peptide
- * sequence i.e. return a single result. Bacterial dna may have overlapping
- * cds mappings coding for multiple peptides so return multiple results
- * (example EMBL KF591215).
+ * 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 dataset sequence
- * @param mapping
- * containing one or more mappings of the sequence to protein
- * @param newMappings
- * the new mapping to populate, from the cds-only sequences to their
- * mapped protein sequences
+ * 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
*/
- protected static List makeCdsSequences(SequenceI dnaSeq,
- AlignedCodonFrame mapping, AlignedCodonFrame newMappings)
+ static SequenceI findCdsForProtein(List mappings,
+ SequenceI dnaSeq, List seqMappings,
+ Mapping aMapping)
{
- List cdsSequences = new ArrayList();
- List seqMappings = mapping.getMappingsForSequence(dnaSeq);
+ /*
+ * 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();
- for (Mapping seqMapping : seqMappings)
+ /*
+ * 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)
{
- SequenceI cds = makeCdsSequence(dnaSeq, seqMapping);
- cdsSequences.add(cds);
-
- /*
- * add new mappings, from dna to cds, and from cds to peptide
- */
- MapList dnaToCds = addCdsMappings(dnaSeq, cds, seqMapping,
- newMappings);
-
- /*
- * transfer any features on dna that overlap the CDS
- */
- transferFeatures(dnaSeq, cds, dnaToCds, "CDS" /* SequenceOntology.CDS */);
+ return seqDss;
}
- return cdsSequences;
+
+ /*
+ * 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() == 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 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
+ * @return CDS sequence (as a dataset sequence)
+ */
+ static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping)
+ {
+ char[] seqChars = seq.getSequence();
+ List fromRanges = mapping.getMap().getFromRanges();
+ int cdsWidth = MappingUtils.getLength(fromRanges);
+ char[] newSeqChars = new char[cdsWidth];
+
+ int newPos = 0;
+ for (int[] range : fromRanges)
+ {
+ if (range[0] <= range[1])
+ {
+ // forward strand mapping - just copy the range
+ int length = range[1] - range[0] + 1;
+ System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
+ length);
+ newPos += length;
+ }
+ else
+ {
+ // reverse strand mapping - copy and complement one by one
+ for (int i = range[0]; i >= range[1]; i--)
+ {
+ newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
+ }
+ }
+ }
+
+ /*
+ * 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;
}
/**
- * Transfers any co-located features on 'fromSeq' to 'toSeq', adjusting the
+ * 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.
*
* @param fromSeq
* @param toSeq
+ * @param select
+ * if not null, only features of this type are copied (including
+ * subtypes in the Sequence Ontology)
* @param mapping
* the mapping from 'fromSeq' to 'toSeq'
* @param omitting
*/
- protected static void transferFeatures(SequenceI fromSeq,
- SequenceI toSeq, MapList mapping, String... omitting)
+ public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
+ MapList mapping, String select, String... omitting)
{
SequenceI copyTo = toSeq;
while (copyTo.getDatasetSequence() != null)
@@ -1430,12 +1907,18 @@ public class AlignmentUtils
copyTo = copyTo.getDatasetSequence();
}
+ SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+ int count = 0;
SequenceFeature[] sfs = fromSeq.getSequenceFeatures();
if (sfs != null)
{
for (SequenceFeature sf : sfs)
{
String type = sf.getType();
+ if (select != null && !so.isA(type, select))
+ {
+ continue;
+ }
boolean omit = false;
for (String toOmit : omitting)
{
@@ -1453,119 +1936,863 @@ public class AlignmentUtils
* locate the mapped range - null if either start or end is
* not mapped (no partial overlaps are calculated)
*/
- int[] mappedTo = mapping.locateInTo(sf.getBegin(), sf.getEnd());
+ int start = sf.getBegin();
+ int end = sf.getEnd();
+ int[] mappedTo = mapping.locateInTo(start, end);
+ /*
+ * if whole exon range doesn't map, try interpreting it
+ * as 5' or 3' exon overlapping the CDS range
+ */
+ if (mappedTo == null)
+ {
+ mappedTo = mapping.locateInTo(end, end);
+ if (mappedTo != null)
+ {
+ /*
+ * end of exon is in CDS range - 5' overlap
+ * to a range from the start of the peptide
+ */
+ mappedTo[0] = 1;
+ }
+ }
+ if (mappedTo == null)
+ {
+ mappedTo = mapping.locateInTo(start, start);
+ if (mappedTo != null)
+ {
+ /*
+ * start of exon is in CDS range - 3' overlap
+ * to a range up to the end of the peptide
+ */
+ mappedTo[1] = toSeq.getLength();
+ }
+ }
if (mappedTo != null)
{
SequenceFeature copy = new SequenceFeature(sf);
copy.setBegin(Math.min(mappedTo[0], mappedTo[1]));
copy.setEnd(Math.max(mappedTo[0], mappedTo[1]));
copyTo.addSequenceFeature(copy);
+ count++;
}
}
}
+ return count;
}
/**
- * Creates and adds mappings
- *
- * - from cds to peptide
- * - from dna to cds
- *
- * and returns the dna-to-cds mapping
+ * Returns a mapping from dna to protein by inspecting sequence features of
+ * type "CDS" on the dna.
*
* @param dnaSeq
- * @param cdsSeq
- * @param dnaMapping
- * @param newMappings
+ * @param proteinSeq
* @return
*/
- protected static MapList addCdsMappings(SequenceI dnaSeq,
- SequenceI cdsSeq,
- Mapping dnaMapping, AlignedCodonFrame newMappings)
+ public static MapList mapCdsToProtein(SequenceI dnaSeq,
+ SequenceI proteinSeq)
{
- cdsSeq.createDatasetSequence();
+ List ranges = findCdsPositions(dnaSeq);
+ int mappedDnaLength = MappingUtils.getLength(ranges);
+
+ int proteinLength = proteinSeq.getLength();
+ int proteinStart = proteinSeq.getStart();
+ int proteinEnd = proteinSeq.getEnd();
/*
- * CDS to peptide is just a contiguous 3:1 mapping, with
- * the peptide ranges taken unchanged from the dna mapping
+ * incomplete start codon may mean X at start of peptide
+ * we ignore both for mapping purposes
*/
- List cdsRanges = new ArrayList();
- cdsRanges.add(new int[] { 1, cdsSeq.getLength() });
- MapList cdsToPeptide = new MapList(cdsRanges, dnaMapping.getMap()
- .getToRanges(), 3, 1);
- newMappings.addMap(cdsSeq.getDatasetSequence(), dnaMapping.getTo(),
- cdsToPeptide);
+ if (proteinSeq.getCharAt(0) == 'X')
+ {
+ // todo JAL-2022 support startPhase > 0
+ proteinStart++;
+ proteinLength--;
+ }
+ List proteinRange = new ArrayList();
/*
- * dna 'from' ranges map 1:1 to the contiguous extracted CDS
+ * dna length should map to protein (or protein plus stop codon)
*/
- MapList dnaToCds = new MapList(
- dnaMapping.getMap().getFromRanges(), cdsRanges, 1, 1);
- newMappings.addMap(dnaSeq, cdsSeq.getDatasetSequence(), dnaToCds);
- return dnaToCds;
+ int codesForResidues = mappedDnaLength / 3;
+ if (codesForResidues == (proteinLength + 1))
+ {
+ // assuming extra codon is for STOP and not in peptide
+ codesForResidues--;
+ }
+ if (codesForResidues == proteinLength)
+ {
+ proteinRange.add(new int[] { proteinStart, proteinEnd });
+ return new MapList(ranges, proteinRange, 3, 1);
+ }
+ return null;
}
/**
- * Makes and returns a CDS-only sequence, where the CDS regions are identified
- * as the 'from' ranges of the mapping on the dna.
+ * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of
+ * start/end positions of sequence features of type "CDS" (or a sub-type of
+ * CDS in the Sequence Ontology). The ranges are sorted into ascending start
+ * position order, so this method is only valid for linear CDS in the same
+ * sense as the protein product.
*
* @param dnaSeq
- * nucleotide sequence
- * @param seqMapping
- * mappings from CDS regions of nucleotide
* @return
*/
- protected static SequenceI makeCdsSequence(SequenceI dnaSeq,
- Mapping seqMapping)
+ public static List findCdsPositions(SequenceI dnaSeq)
+ {
+ List result = new ArrayList();
+ SequenceFeature[] sfs = dnaSeq.getSequenceFeatures();
+ if (sfs == null)
+ {
+ return result;
+ }
+
+ SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+ int startPhase = 0;
+
+ for (SequenceFeature sf : sfs)
+ {
+ /*
+ * process a CDS feature (or a sub-type of CDS)
+ */
+ if (so.isA(sf.getType(), SequenceOntologyI.CDS))
+ {
+ int phase = 0;
+ try
+ {
+ phase = Integer.parseInt(sf.getPhase());
+ } catch (NumberFormatException e)
+ {
+ // ignore
+ }
+ /*
+ * phase > 0 on first codon means 5' incomplete - skip to the start
+ * of the next codon; example ENST00000496384
+ */
+ int begin = sf.getBegin();
+ int end = sf.getEnd();
+ if (result.isEmpty())
+ {
+ begin += phase;
+ if (begin > end)
+ {
+ // shouldn't happen!
+ System.err
+ .println("Error: start phase extends beyond start CDS in "
+ + dnaSeq.getName());
+ }
+ }
+ result.add(new int[] { begin, end });
+ }
+ }
+
+ /*
+ * remove 'startPhase' positions (usually 0) from the first range
+ * so we begin at the start of a complete codon
+ */
+ if (!result.isEmpty())
+ {
+ // TODO JAL-2022 correctly model start phase > 0
+ result.get(0)[0] += startPhase;
+ }
+
+ /*
+ * Finally sort ranges by start position. This avoids a dependency on
+ * keeping features in order on the sequence (if they are in order anyway,
+ * the sort will have almost no work to do). The implicit assumption is CDS
+ * ranges are assembled in order. Other cases should not use this method,
+ * but instead construct an explicit mapping for CDS (e.g. EMBL parsing).
+ */
+ Collections.sort(result, new Comparator()
+ {
+ @Override
+ public int compare(int[] o1, int[] o2)
+ {
+ return Integer.compare(o1[0], o2[0]);
+ }
+ });
+ return result;
+ }
+
+ /**
+ * Maps exon features from dna to protein, and computes variants in peptide
+ * product generated by variants in dna, and adds them as sequence_variant
+ * features on the protein sequence. Returns the number of variant features
+ * added.
+ *
+ * @param dnaSeq
+ * @param peptide
+ * @param dnaToProtein
+ */
+ public static int computeProteinFeatures(SequenceI dnaSeq,
+ SequenceI peptide, MapList dnaToProtein)
+ {
+ while (dnaSeq.getDatasetSequence() != null)
+ {
+ dnaSeq = dnaSeq.getDatasetSequence();
+ }
+ while (peptide.getDatasetSequence() != null)
+ {
+ peptide = peptide.getDatasetSequence();
+ }
+
+ transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON);
+
+ /*
+ * compute protein variants from dna variants and codon mappings;
+ * NB - alternatively we could retrieve this using the REST service e.g.
+ * http://rest.ensembl.org/overlap/translation
+ * /ENSP00000288602?feature=transcript_variation;content-type=text/xml
+ * which would be a bit slower but possibly more reliable
+ */
+
+ /*
+ * build a map with codon variations for each potentially varying peptide
+ */
+ LinkedHashMap[]> variants = buildDnaVariantsMap(
+ dnaSeq, dnaToProtein);
+
+ /*
+ * scan codon variations, compute peptide variants and add to peptide sequence
+ */
+ int count = 0;
+ for (Entry[]> variant : variants.entrySet())
+ {
+ int peptidePos = variant.getKey();
+ List[] codonVariants = variant.getValue();
+ count += computePeptideVariants(peptide, peptidePos, codonVariants);
+ }
+
+ /*
+ * sort to get sequence features in start position order
+ * - would be better to store in Sequence as a TreeSet or NCList?
+ */
+ if (peptide.getSequenceFeatures() != null)
+ {
+ Arrays.sort(peptide.getSequenceFeatures(),
+ new Comparator()
+ {
+ @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;
+ }
+
+ /**
+ * Computes non-synonymous peptide variants from codon variants and adds them
+ * as sequence_variant features on the protein sequence (one feature per
+ * allele variant). Selected attributes (variant id, clinical significance)
+ * are copied over to the new features.
+ *
+ * @param peptide
+ * the protein sequence
+ * @param peptidePos
+ * the position to compute peptide variants for
+ * @param codonVariants
+ * a list of dna variants per codon position
+ * @return the number of features added
+ */
+ static int computePeptideVariants(SequenceI peptide, int peptidePos,
+ List[] codonVariants)
{
- StringBuilder newSequence = new StringBuilder(dnaSeq.getLength());
- final char[] dna = dnaSeq.getSequence();
- int offset = dnaSeq.getStart() - 1;
+ String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
+ int count = 0;
+ String base1 = codonVariants[0].get(0).base;
+ String base2 = codonVariants[1].get(0).base;
+ String base3 = codonVariants[2].get(0).base;
/*
- * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc }
+ * variants in first codon base
*/
- final List dnaCdsRanges = seqMapping.getMap().getFromRanges();
- for (int[] range : dnaCdsRanges)
+ for (DnaVariant var : codonVariants[0])
{
- // TODO handle reverse mapping as well (range[1] < range[0])
- for (int pos = range[0]; pos <= range[1]; pos++)
+ if (var.variant != null)
{
- newSequence.append(dna[pos - offset - 1]);
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base + base2 + base3;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
}
}
- SequenceI cds = new Sequence(dnaSeq.getName(),
- newSequence.toString());
+ /*
+ * variants in second codon base
+ */
+ for (DnaVariant var : codonVariants[1])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base1 + base + base3;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
- transferDbRefs(seqMapping.getTo(), cds);
+ /*
+ * variants in third codon base
+ */
+ for (DnaVariant var : codonVariants[2])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base1 + base2 + base;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
- return cds;
+ return count;
}
/**
- * Locate any xrefs to CDS databases on the protein product and attach to the
- * CDS sequence. Also add as a sub-token of the sequence name.
+ * Helper method that adds a peptide variant feature, provided the given codon
+ * translates to a value different to the current residue (is a non-synonymous
+ * variant). ID and clinical_significance attributes of the dna variant (if
+ * present) are copied to the new feature.
*
- * @param from
- * @param to
+ * @param peptide
+ * @param peptidePos
+ * @param residue
+ * @param var
+ * @param codon
+ * @return true if a feature was added, else false
*/
- protected static void transferDbRefs(SequenceI from, SequenceI to)
+ static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
+ String residue, DnaVariant var, String codon)
{
- String cdsAccId = FeatureProperties.getCodingFeature(DBRefSource.EMBL);
- DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(from.getDBRef(),
- DBRefSource.CODINGDBS);
- if (cdsRefs != null)
+ /*
+ * get peptide translation of codon e.g. GAT -> D
+ * note that variants which are not single alleles,
+ * e.g. multibase variants or HGMD_MUTATION etc
+ * are currently ignored here
+ */
+ String trans = codon.contains("-") ? "-"
+ : (codon.length() > 3 ? null : ResidueProperties
+ .codonTranslate(codon));
+ if (trans != null && !trans.equals(residue))
{
- for (DBRefEntry cdsRef : cdsRefs)
+ String residue3Char = StringUtils
+ .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
+ String trans3Char = StringUtils
+ .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
+ String desc = "p." + residue3Char + peptidePos + trans3Char;
+ // set score to 0f so 'graduated colour' option is offered! JAL-2060
+ SequenceFeature sf = new SequenceFeature(
+ SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
+ peptidePos, 0f, "Jalview");
+ StringBuilder attributes = new StringBuilder(32);
+ String id = (String) var.variant.getValue(ID);
+ if (id != null)
+ {
+ if (id.startsWith(SEQUENCE_VARIANT))
+ {
+ id = id.substring(SEQUENCE_VARIANT.length());
+ }
+ sf.setValue(ID, id);
+ attributes.append(ID).append("=").append(id);
+ // TODO handle other species variants
+ StringBuilder link = new StringBuilder(32);
+ try
+ {
+ 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());
+ } catch (UnsupportedEncodingException e)
+ {
+ // as if
+ }
+ }
+ String clinSig = (String) var.variant
+ .getValue(CLINICAL_SIGNIFICANCE);
+ if (clinSig != null)
+ {
+ sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
+ attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
+ .append(clinSig);
+ }
+ peptide.addSequenceFeature(sf);
+ if (attributes.length() > 0)
{
- to.addDBRef(new DBRefEntry(cdsRef));
- cdsAccId = cdsRef.getAccessionId();
+ sf.setAttributes(attributes.toString());
}
+ return true;
}
- if (!to.getName().contains(cdsAccId))
+ return false;
+ }
+
+ /**
+ * Builds a map whose key is position in the protein sequence, and value is a
+ * list of the base and all variants for each corresponding codon position
+ *
+ * @param dnaSeq
+ * @param dnaToProtein
+ * @return
+ */
+ static LinkedHashMap[]> buildDnaVariantsMap(
+ SequenceI dnaSeq, MapList dnaToProtein)
+ {
+ /*
+ * map from peptide position to all variants of the codon which codes for it
+ * LinkedHashMap ensures we keep the peptide features in sequence order
+ */
+ LinkedHashMap[]> variants = new LinkedHashMap[]>();
+ SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+
+ SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
+ if (dnaFeatures == null)
{
- to.setName(to.getName() + "|" + cdsAccId);
+ return variants;
}
+
+ int dnaStart = dnaSeq.getStart();
+ int[] lastCodon = null;
+ int lastPeptidePostion = 0;
+
+ /*
+ * build a map of codon variations for peptides
+ */
+ for (SequenceFeature sf : dnaFeatures)
+ {
+ int dnaCol = sf.getBegin();
+ if (dnaCol != sf.getEnd())
+ {
+ // not handling multi-locus variant features
+ continue;
+ }
+ if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT))
+ {
+ int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
+ if (mapsTo == null)
+ {
+ // feature doesn't lie within coding region
+ continue;
+ }
+ int peptidePosition = mapsTo[0];
+ List[] codonVariants = variants.get(peptidePosition);
+ if (codonVariants == null)
+ {
+ codonVariants = new ArrayList[3];
+ codonVariants[0] = new ArrayList();
+ codonVariants[1] = new ArrayList();
+ codonVariants[2] = new ArrayList();
+ variants.put(peptidePosition, codonVariants);
+ }
+
+ /*
+ * extract dna variants to a string array
+ */
+ String alls = (String) sf.getValue("alleles");
+ if (alls == null)
+ {
+ continue;
+ }
+ String[] alleles = alls.toUpperCase().split(",");
+ int i = 0;
+ for (String allele : alleles)
+ {
+ alleles[i++] = allele.trim(); // lose any space characters "A, G"
+ }
+
+ /*
+ * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
+ */
+ int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
+ : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
+ peptidePosition, peptidePosition));
+ lastPeptidePostion = peptidePosition;
+ lastCodon = codon;
+
+ /*
+ * save nucleotide (and any variant) for each codon position
+ */
+ for (int codonPos = 0; codonPos < 3; codonPos++)
+ {
+ String nucleotide = String.valueOf(
+ dnaSeq.getCharAt(codon[codonPos] - dnaStart))
+ .toUpperCase();
+ List codonVariant = codonVariants[codonPos];
+ if (codon[codonPos] == dnaCol)
+ {
+ if (!codonVariant.isEmpty()
+ && codonVariant.get(0).variant == null)
+ {
+ /*
+ * already recorded base value, add this variant
+ */
+ codonVariant.get(0).variant = sf;
+ }
+ else
+ {
+ /*
+ * add variant with base value
+ */
+ codonVariant.add(new DnaVariant(nucleotide, sf));
+ }
+ }
+ else if (codonVariant.isEmpty())
+ {
+ /*
+ * record (possibly non-varying) base value
+ */
+ codonVariant.add(new DnaVariant(nucleotide));
+ }
+ }
+ }
+ }
+ return variants;
+ }
+
+ /**
+ * Makes an alignment with a copy of the given sequences, adding in any
+ * non-redundant sequences which are mapped to by the cross-referenced
+ * sequences.
+ *
+ * @param seqs
+ * @param xrefs
+ * @param dataset
+ * the alignment dataset shared by the new copy
+ * @return
+ */
+ public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
+ SequenceI[] xrefs, AlignmentI dataset)
+ {
+ AlignmentI copy = new Alignment(new Alignment(seqs));
+ copy.setDataset(dataset);
+
+ SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
+ if (xrefs != null)
+ {
+ for (SequenceI xref : xrefs)
+ {
+ DBRefEntry[] dbrefs = xref.getDBRefs();
+ if (dbrefs != null)
+ {
+ for (DBRefEntry dbref : dbrefs)
+ {
+ if (dbref.getMap() == null || dbref.getMap().getTo() == null)
+ {
+ continue;
+ }
+ SequenceI mappedTo = dbref.getMap().getTo();
+ SequenceI match = matcher.findIdMatch(mappedTo);
+ if (match == null)
+ {
+ matcher.add(mappedTo);
+ copy.addSequence(mappedTo);
+ }
+ }
+ }
+ }
+ }
+ return copy;
+ }
+
+ /**
+ * Try to align sequences in 'unaligned' to match the alignment of their
+ * mapped regions in 'aligned'. For example, could use this to align CDS
+ * sequences which are mapped to their parent cDNA sequences.
+ *
+ * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For
+ * dna-to-protein or protein-to-dna use alternative methods.
+ *
+ * @param unaligned
+ * sequences to be aligned
+ * @param aligned
+ * holds aligned sequences and their mappings
+ * @return
+ */
+ 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;
+
+ for (SequenceI seq : unaligned.getSequences())
+ {
+ if (!unmapped.contains(seq))
+ {
+ char[] newSeq = new char[width];
+ Arrays.fill(newSeq, gap);
+ int newCol = 0;
+ int lastCol = 0;
+
+ /*
+ * traverse the map to find columns populated
+ * by our sequence
+ */
+ for (Integer column : columnMap.keySet())
+ {
+ Character c = columnMap.get(column).get(seq);
+ if (c != null)
+ {
+ /*
+ * sequence has a character at this position
+ *
+ */
+ newSeq[newCol] = c;
+ lastCol = newCol;
+ }
+ newCol++;
+ }
+
+ /*
+ * trim trailing gaps
+ */
+ if (lastCol < width)
+ {
+ char[] tmp = new char[lastCol + 1];
+ System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
+ newSeq = tmp;
+ }
+ seq.setSequence(String.valueOf(newSeq));
+ realignedCount++;
+ }
+ }
+ return realignedCount;
+ }
+
+ /**
+ * 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 alignedDatasets = new HashMap();
+ 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.
+ *
+ * @param unaligned
+ * @param aligned
+ * @param unmapped
+ * @return
+ */
+ static Map> buildMappedColumnsMap(
+ AlignmentI unaligned, AlignmentI aligned, List unmapped)
+ {
+ /*
+ * Map will hold, for each aligned column position, a map of
+ * {unalignedSequence, characterPerSequence} at that position.
+ * TreeMap keeps the entries in ascending column order.
+ */
+ Map> map = new TreeMap>();
+
+ /*
+ * record any sequences that have no mapping so can't be realigned
+ */
+ unmapped.addAll(unaligned.getSequences());
+
+ List mappings = aligned.getCodonFrames();
+
+ for (SequenceI seq : unaligned.getSequences())
+ {
+ for (AlignedCodonFrame mapping : mappings)
+ {
+ SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned);
+ if (fromSeq != null)
+ {
+ Mapping seqMap = mapping.getMappingBetween(fromSeq, seq);
+ if (addMappedPositions(seq, fromSeq, seqMap, map))
+ {
+ unmapped.remove(seq);
+ }
+ }
+ }
+ }
+ return map;
+ }
+
+ /**
+ * Helper method that adds to a map the mapped column positions of a sequence.
+ * For example if aaTT-Tg-gAAA is mapped to TTTAAA then the map should record
+ * that columns 3,4,6,10,11,12 map to characters T,T,T,A,A,A of the mapped to
+ * sequence.
+ *
+ * @param seq
+ * the sequence whose column positions we are recording
+ * @param fromSeq
+ * a sequence that is mapped to the first sequence
+ * @param seqMap
+ * the mapping from 'fromSeq' to 'seq'
+ * @param map
+ * a map to add the column positions (in fromSeq) of the mapped
+ * positions of seq
+ * @return
+ */
+ static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq,
+ Mapping seqMap, Map> map)
+ {
+ if (seqMap == null)
+ {
+ 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();
+
+ /*
+ * traverse [start, end, start, end...] ranges in fromSeq
+ */
+ for (int[] fromRange : seqMap.getMap().getFromRanges())
+ {
+ for (int i = 0; i < fromRange.length - 1; i += 2)
+ {
+ boolean forward = fromRange[i + 1] >= fromRange[i];
+
+ /*
+ * find the range mapped to (sequence positions base 1)
+ */
+ int[] range = seqMap.locateMappedRange(fromRange[i],
+ fromRange[i + 1]);
+ if (range == null)
+ {
+ System.err.println("Error in mapping " + seqMap + " from "
+ + fromSeq.getName());
+ return false;
+ }
+ int fromCol = fromSeq.findIndex(fromRange[i]);
+ int mappedCharPos = range[0];
+
+ /*
+ * walk over the 'from' aligned sequence in forward or reverse
+ * direction; when a non-gap is found, record the column position
+ * of the next character of the mapped-to sequence; stop when all
+ * the characters of the range have been counted
+ */
+ while (mappedCharPos <= range[1] && fromCol <= fromChars.length
+ && fromCol >= 0)
+ {
+ if (!Comparison.isGap(fromChars[fromCol - 1]))
+ {
+ /*
+ * mapped from sequence has a character in this column
+ * record the column position for the mapped to character
+ */
+ Map seqsMap = map.get(fromCol);
+ if (seqsMap == null)
+ {
+ seqsMap = new HashMap();
+ map.put(fromCol, seqsMap);
+ }
+ seqsMap.put(seq, toChars[mappedCharPos - toStart]);
+ mappedCharPos++;
+ }
+ fromCol += (forward ? 1 : -1);
+ }
+ }
+ }
+ return true;
+ }
+
+ // strictly temporary hack until proper criteria for aligning protein to cds
+ // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot
+ public static boolean looksLikeEnsembl(AlignmentI alignment)
+ {
+ for (SequenceI seq : alignment.getSequences())
+ {
+ String name = seq.getName();
+ if (!name.startsWith("ENSG") && !name.startsWith("ENST"))
+ {
+ return false;
+ }
+ }
+ return true;
}
}