*/
package jalview.analysis;
+import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
+
import jalview.datamodel.AlignedCodon;
import jalview.datamodel.AlignedCodonFrame;
import jalview.datamodel.Alignment;
import jalview.datamodel.AlignmentAnnotation;
import jalview.datamodel.AlignmentI;
import jalview.datamodel.DBRefEntry;
-import jalview.datamodel.DBRefSource;
-import jalview.datamodel.FeatureProperties;
import jalview.datamodel.IncompleteCodonException;
import jalview.datamodel.Mapping;
-import jalview.datamodel.SearchResults;
import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceFeature;
import jalview.datamodel.SequenceGroup;
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;
+import java.io.UnsupportedEncodingException;
+import java.net.URLEncoder;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
public class AlignmentUtils
{
+ private static final String SEQUENCE_VARIANT = "sequence_variant:";
+ private static final String ID = "ID";
+
+ /**
+ * A data model to hold the 'normal' base value at a position, and an optional
+ * sequence variant feature
+ */
+ static class DnaVariant
+ {
+ String base;
+
+ SequenceFeature variant;
+
+ DnaVariant(String nuc)
+ {
+ base = nuc;
+ }
+
+ DnaVariant(String nuc, SequenceFeature var)
+ {
+ base = nuc;
+ variant = var;
+ }
+ }
+
/**
* given an existing alignment, create a new alignment including all, or up to
* flankSize additional symbols from each sequence's dataset sequence
AlignedCodonFrame mapping = null;
for (AlignedCodonFrame mp : mappings)
{
- alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al);
+ alignFrom = mp.findAlignedSequence(seq, al);
if (alignFrom != null)
{
mapping = mp;
}
/**
- * Returns a list of sequences mapped from the given sequences and aligned
- * (gapped) in the same way. For example, the cDNA for aligned protein, where
- * a single gap in protein generates three gaps in cDNA.
- *
- * @param sequences
- * @param gapCharacter
- * @param mappings
- * @return
- */
- public static List<SequenceI> getAlignedTranslation(
- List<SequenceI> sequences, char gapCharacter,
- Set<AlignedCodonFrame> mappings)
- {
- List<SequenceI> alignedSeqs = new ArrayList<SequenceI>();
-
- for (SequenceI seq : sequences)
- {
- List<SequenceI> mapped = getAlignedTranslation(seq, gapCharacter,
- mappings);
- alignedSeqs.addAll(mapped);
- }
- return alignedSeqs;
- }
-
- /**
- * Returns sequences aligned 'like' the source sequence, as mapped by the
- * given mappings. Normally we expect zero or one 'mapped' sequences, but this
- * will support 1-to-many as well.
- *
- * @param seq
- * @param gapCharacter
- * @param mappings
- * @return
- */
- protected static List<SequenceI> getAlignedTranslation(SequenceI seq,
- char gapCharacter, Set<AlignedCodonFrame> mappings)
- {
- List<SequenceI> result = new ArrayList<SequenceI>();
- for (AlignedCodonFrame mapping : mappings)
- {
- if (mapping.involvesSequence(seq))
- {
- SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping);
- if (mapped != null)
- {
- result.add(mapped);
- }
- }
- }
- return result;
- }
-
- /**
- * Returns the translation of 'seq' (as held in the mapping) with
- * corresponding alignment (gaps).
- *
- * @param seq
- * @param gapCharacter
- * @param mapping
- * @return
- */
- protected static SequenceI getAlignedTranslation(SequenceI seq,
- char gapCharacter, AlignedCodonFrame mapping)
- {
- String gap = String.valueOf(gapCharacter);
- boolean toDna = false;
- int fromRatio = 1;
- SequenceI mapTo = mapping.getDnaForAaSeq(seq);
- if (mapTo != null)
- {
- // mapping is from protein to nucleotide
- toDna = true;
- // should ideally get gap count ratio from mapping
- gap = String.valueOf(new char[] { gapCharacter, gapCharacter,
- gapCharacter });
- }
- else
- {
- // mapping is from nucleotide to protein
- mapTo = mapping.getAaForDnaSeq(seq);
- fromRatio = 3;
- }
- StringBuilder newseq = new StringBuilder(seq.getLength()
- * (toDna ? 3 : 1));
-
- int residueNo = 0; // in seq, base 1
- int[] phrase = new int[fromRatio];
- int phraseOffset = 0;
- int gapWidth = 0;
- boolean first = true;
- final Sequence alignedSeq = new Sequence("", "");
-
- for (char c : seq.getSequence())
- {
- if (c == gapCharacter)
- {
- gapWidth++;
- if (gapWidth >= fromRatio)
- {
- newseq.append(gap);
- gapWidth = 0;
- }
- }
- else
- {
- phrase[phraseOffset++] = residueNo + 1;
- if (phraseOffset == fromRatio)
- {
- /*
- * Have read a whole codon (or protein residue), now translate: map
- * source phrase to positions in target sequence add characters at
- * these positions to newseq Note mapping positions are base 1, our
- * sequence positions base 0.
- */
- SearchResults sr = new SearchResults();
- for (int pos : phrase)
- {
- mapping.markMappedRegion(seq, pos, sr);
- }
- newseq.append(sr.getCharacters());
- if (first)
- {
- first = false;
- // Hack: Copy sequence dataset, name and description from
- // SearchResults.match[0].sequence
- // TODO? carry over sequence names from original 'complement'
- // alignment
- SequenceI mappedTo = sr.getResultSequence(0);
- alignedSeq.setName(mappedTo.getName());
- alignedSeq.setDescription(mappedTo.getDescription());
- alignedSeq.setDatasetSequence(mappedTo);
- }
- phraseOffset = 0;
- }
- residueNo++;
- }
- }
- alignedSeq.setSequence(newseq.toString());
- return alignedSeq;
- }
-
- /**
* Realigns the given protein to match the alignment of the dna, using codon
* mappings to translate aligned codon positions to protein residues.
*
{
for (AlignedCodonFrame mapping : mappings)
{
- SequenceI prot = mapping.findAlignedSequence(
- dnaSeq.getDatasetSequence(), protein);
+ SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein);
if (prot != null)
{
Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
* Finally add any unmapped peptide start residues (e.g. for incomplete
* codons) as if at the codon position before the second residue
*/
+ // TODO resolve JAL-2022 so this fudge can be removed
int mappedSequenceCount = protein.getHeight() - unmappedProtein.size();
addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount);
/**
* Constructs an alignment consisting of the mapped (CDS) regions in the given
- * nucleotide sequences, and updates mappings to match. The new sequences are
- * aligned as per the original sequence, with entirely gapped columns (codon
- * interrupted by intron) omitted.
+ * nucleotide sequences, and updates mappings to match. The CDS sequences are
+ * added to the original alignment's dataset, which is shared by the new
+ * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are
+ * added to the alignment dataset.
*
* @param dna
* aligned dna sequences
* @param mappings
- * from dna to protein; these are replaced with new mappings
+ * from dna to protein
* @param al
* @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)
{
- List<int[]> cdsColumns = findCdsColumns(dna);
-
- /*
- * create CDS sequences and new mappings
- * (from cdna to cds, and cds to peptide)
- */
- List<AlignedCodonFrame> newMappings = new ArrayList<AlignedCodonFrame>();
- List<SequenceI> cdsSequences = new ArrayList<SequenceI>();
- char gap = al.getGapCharacter();
-
- for (SequenceI dnaSeq : dna)
+ List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
+
+ for (SequenceI seq : dna)
{
- final SequenceI ds = dnaSeq.getDatasetSequence();
+ AlignedCodonFrame cdsMappings = new AlignedCodonFrame();
List<AlignedCodonFrame> seqMappings = MappingUtils
- .findMappingsForSequence(ds, mappings);
- for (AlignedCodonFrame acf : seqMappings)
+ .findMappingsForSequence(seq, mappings);
+ List<AlignedCodonFrame> alignmentMappings = al.getCodonFrames();
+ for (AlignedCodonFrame mapping : seqMappings)
{
- AlignedCodonFrame newMapping = new AlignedCodonFrame();
- final List<SequenceI> mappedCds = makeCdsSequences(dnaSeq, acf,
- cdsColumns, newMapping, gap);
- if (!mappedCds.isEmpty())
+ for (Mapping aMapping : mapping.getMappingsFromSequence(seq))
{
- cdsSequences.addAll(mappedCds);
- newMappings.add(newMapping);
+ SequenceI cdsSeq = makeCdsSequence(seq.getDatasetSequence(),
+ aMapping);
+ cdsSeqs.add(cdsSeq);
+
+ /*
+ * 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);
+
+ /*
+ * add another mapping from original 'from' range to CDS
+ */
+ map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
+ 1);
+ cdsMappings.addMap(seq.getDatasetSequence(), cdsSeq, map);
+
+ alignmentMappings.add(cdsMappings);
+
+ /*
+ * transfer any features on dna that overlap the CDS
+ */
+ transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS);
}
}
}
- AlignmentI newAl = new Alignment(
- cdsSequences.toArray(new SequenceI[cdsSequences.size()]));
/*
- * add new sequences to the shared dataset, set it on the new alignment
+ * add CDS seqs to shared dataset
*/
- List<SequenceI> dsseqs = al.getDataset().getSequences();
- for (SequenceI seq : newAl.getSequences())
+ Alignment dataset = al.getDataset();
+ for (SequenceI seq : cdsSeqs)
{
- if (!dsseqs.contains(seq.getDatasetSequence()))
+ if (!dataset.getSequences().contains(seq.getDatasetSequence()))
{
- dsseqs.add(seq.getDatasetSequence());
+ dataset.addSequence(seq.getDatasetSequence());
}
}
- newAl.setDataset(al.getDataset());
-
- /*
- * 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 newAl;
+ return cds;
}
/**
- * Returns a consolidated list of column ranges where at least one sequence
- * has a CDS feature. This assumes CDS features are on genomic sequence i.e.
- * are for contiguous CDS ranges (no gaps).
+ * 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 seqs
+ * @param seq
+ * @param mapping
* @return
*/
- public static List<int[]> findCdsColumns(SequenceI[] seqs)
+ static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping)
{
- // TODO use refactored code from AlignViewController
- // markColumnsContainingFeatures, not reinvent the wheel!
-
- List<int[]> result = new ArrayList<int[]>();
- for (SequenceI seq : seqs)
- {
- result.addAll(findCdsColumns(seq));
- }
+ char[] seqChars = seq.getSequence();
+ List<int[]> fromRanges = mapping.getMap().getFromRanges();
+ int cdsWidth = MappingUtils.getLength(fromRanges);
+ char[] newSeqChars = new char[cdsWidth];
- /*
- * sort and compact the list into ascending, non-overlapping ranges
- */
- Collections.sort(result, new Comparator<int[]>()
+ int newPos = 0;
+ for (int[] range : fromRanges)
{
- @Override
- public int compare(int[] o1, int[] o2)
+ if (range[0] <= range[1])
{
- return Integer.compare(o1[0], o2[0]);
+ // 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;
}
- });
- result = MapList.coalesceRanges(result);
-
- return result;
- }
-
- public static List<int[]> findCdsColumns(SequenceI seq)
- {
- List<int[]> result = new ArrayList<int[]>();
- SequenceOntologyI so = SequenceOntologyFactory.getInstance();
- SequenceFeature[] sfs = seq.getSequenceFeatures();
- if (sfs != null)
- {
- for (SequenceFeature sf : sfs)
- {
- if (so.isA(sf.getType(), SequenceOntologyI.CDS))
- {
- int colStart = seq.findIndex(sf.getBegin());
- int colEnd = seq.findIndex(sf.getEnd());
- result.add(new int[] { colStart, colEnd });
- }
- }
- }
- return result;
- }
-
- /**
- * Answers true if all sequences have a gap at (or do not extend to) the
- * specified column position (base 1)
- *
- * @param seqs
- * @param col
- * @return
- */
- public static boolean isGappedColumn(List<SequenceI> seqs, int col)
- {
- if (seqs != null)
- {
- for (SequenceI seq : seqs)
- {
- if (!Comparison.isGap(seq.getCharAt(col - 1)))
- {
- return false;
- }
- }
- }
- return true;
- }
-
- /**
- * Returns the column ranges (base 1) of each aligned sequence that are
- * involved in any mapping. This is a helper method for aligning protein
- * products of aligned transcripts.
- *
- * @param mappedSequences
- * (possibly gapped) dna sequences
- * @param mappings
- * @return
- */
- protected static List<List<int[]>> getMappedColumns(
- List<SequenceI> mappedSequences, List<AlignedCodonFrame> mappings)
- {
- List<List<int[]>> result = new ArrayList<List<int[]>>();
- for (SequenceI seq : mappedSequences)
- {
- List<int[]> columns = new ArrayList<int[]>();
- List<AlignedCodonFrame> seqMappings = MappingUtils
- .findMappingsForSequence(seq, mappings);
- for (AlignedCodonFrame mapping : seqMappings)
+ else
{
- List<Mapping> maps = mapping.getMappingsForSequence(seq);
- for (Mapping map : maps)
+ // reverse strand mapping - copy and complement one by one
+ for (int i = range[0]; i >= range[1]; i--)
{
- /*
- * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc }
- * Find and add the overall aligned column range for each
- */
- for (int[] cdsRange : map.getMap().getFromRanges())
- {
- int startPos = cdsRange[0];
- int endPos = cdsRange[1];
- int startCol = seq.findIndex(startPos);
- int endCol = seq.findIndex(endPos);
- columns.add(new int[] { startCol, endCol });
- }
+ newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
}
}
- result.add(columns);
}
- return result;
- }
-
- /**
- * Helper method to make cds-only sequences and populate their mappings to
- * protein products
- * <p>
- * 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
- * <p>
- * 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).
- *
- * @param dnaSeq
- * a dna aligned sequence
- * @param mapping
- * containing one or more mappings of the sequence to protein
- * @param ungappedCdsColumns
- * @param newMappings
- * the new mapping to populate, from the cds-only sequences to their
- * mapped protein sequences
- * @return
- */
- protected static List<SequenceI> makeCdsSequences(SequenceI dnaSeq,
- AlignedCodonFrame mapping, List<int[]> ungappedCdsColumns,
- AlignedCodonFrame newMappings, char gapChar)
- {
- List<SequenceI> cdsSequences = new ArrayList<SequenceI>();
- List<Mapping> seqMappings = mapping.getMappingsForSequence(dnaSeq);
-
- for (Mapping seqMapping : seqMappings)
- {
- SequenceI cds = makeCdsSequence(dnaSeq, seqMapping,
- ungappedCdsColumns, gapChar);
- cds.createDatasetSequence();
- cdsSequences.add(cds);
-
- /*
- * add new mappings, from dna to cds, and from cds to peptide
- */
- MapList dnaToCds = addCdsMappings(dnaSeq.getDatasetSequence(), cds,
- seqMapping, newMappings);
- /*
- * transfer any features on dna that overlap the CDS
- */
- transferFeatures(dnaSeq, cds, dnaToCds, null, SequenceOntologyI.CDS);
- }
- return cdsSequences;
+ SequenceI newSeq = new Sequence(seq.getName() + "|"
+ + mapping.getTo().getName(), newSeqChars, 1, newPos);
+ newSeq.createDatasetSequence();
+ return newSeq;
}
/**
public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
MapList mapping, String select, String... omitting)
{
+ Map<String, Integer> featureCounts = new HashMap<String, Integer>();
+
SequenceI copyTo = toSeq;
while (copyTo.getDatasetSequence() != null)
{
copy.setBegin(Math.min(mappedTo[0], mappedTo[1]));
copy.setEnd(Math.max(mappedTo[0], mappedTo[1]));
copyTo.addSequenceFeature(copy);
+ int featureNumber = featureCounts.containsKey(type) ? featureCounts
+ .get(type) : -1;
+ featureNumber++;
+ featureCounts.put(type, featureNumber);
+ copy.setFeatureNumber(featureNumber);
count++;
}
}
}
/**
- * Creates and adds mappings
- * <ul>
- * <li>from cds to peptide</li>
- * <li>from dna to cds</li>
- * </ul>
- * and returns the dna-to-cds mapping
- *
- * @param dnaSeq
- * @param cdsSeq
- * @param dnaMapping
- * @param newMappings
- * @return
- */
- protected static MapList addCdsMappings(SequenceI dnaSeq,
- SequenceI cdsSeq, Mapping dnaMapping,
- AlignedCodonFrame newMappings)
- {
- cdsSeq.createDatasetSequence();
-
- /*
- * CDS to peptide is just a contiguous 3:1 mapping, with
- * the peptide ranges taken unchanged from the dna mapping
- */
- List<int[]> cdsRanges = new ArrayList<int[]>();
- SequenceI cdsDataset = cdsSeq.getDatasetSequence();
- cdsRanges.add(new int[] { 1, cdsDataset.getLength() });
- MapList cdsToPeptide = new MapList(cdsRanges, dnaMapping.getMap()
- .getToRanges(), 3, 1);
- newMappings.addMap(cdsDataset, dnaMapping.getTo(), cdsToPeptide);
-
- /*
- * dna 'from' ranges map 1:1 to the contiguous extracted CDS
- */
- MapList dnaToCds = new MapList(dnaMapping.getMap().getFromRanges(),
- cdsRanges, 1, 1);
- newMappings.addMap(dnaSeq, cdsDataset, dnaToCds);
- return dnaToCds;
- }
-
- /**
- * Makes and returns a CDS-only sequence, where the CDS regions are identified
- * as the 'from' ranges of the mapping on the dna.
- *
- * @param dnaSeq
- * nucleotide sequence
- * @param seqMapping
- * mappings from CDS regions of nucleotide
- * @param ungappedCdsColumns
- * @return
- */
- protected static SequenceI makeCdsSequence(SequenceI dnaSeq,
- Mapping seqMapping, List<int[]> ungappedCdsColumns, char gapChar)
- {
- int cdsWidth = MappingUtils.getLength(ungappedCdsColumns);
-
- /*
- * populate CDS columns with the aligned
- * column character if that column is mapped (which may be a gap
- * if an intron interrupts a codon), else with a gap
- */
- List<int[]> fromRanges = seqMapping.getMap().getFromRanges();
- char[] cdsChars = new char[cdsWidth];
- int pos = 0;
- for (int[] columns : ungappedCdsColumns)
- {
- for (int i = columns[0]; i <= columns[1]; i++)
- {
- char dnaChar = dnaSeq.getCharAt(i - 1);
- if (Comparison.isGap(dnaChar))
- {
- cdsChars[pos] = gapChar;
- }
- else
- {
- int seqPos = dnaSeq.findPosition(i - 1);
- if (MappingUtils.contains(fromRanges, seqPos))
- {
- cdsChars[pos] = dnaChar;
- }
- else
- {
- cdsChars[pos] = gapChar;
- }
- }
- pos++;
- }
- }
- SequenceI cdsSequence = new Sequence(dnaSeq.getName(),
- String.valueOf(cdsChars));
-
- transferDbRefs(seqMapping.getTo(), cdsSequence);
-
- return cdsSequence;
- }
-
- /**
- * 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.
- *
- * @param from
- * @param to
- */
- protected static void transferDbRefs(SequenceI from, SequenceI to)
- {
- String cdsAccId = FeatureProperties.getCodingFeature(DBRefSource.EMBL);
- DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(from.getDBRefs(),
- DBRefSource.CODINGDBS);
- if (cdsRefs != null)
- {
- for (DBRefEntry cdsRef : cdsRefs)
- {
- to.addDBRef(new DBRefEntry(cdsRef));
- cdsAccId = cdsRef.getAccessionId();
- }
- }
- if (!to.getName().contains(cdsAccId))
- {
- to.setName(to.getName() + "|" + cdsAccId);
- }
- }
-
- /**
* Returns a mapping from dna to protein by inspecting sequence features of
* type "CDS" on the dna.
*
{
List<int[]> ranges = findCdsPositions(dnaSeq);
int mappedDnaLength = MappingUtils.getLength(ranges);
-
+
int proteinLength = proteinSeq.getLength();
int proteinStart = proteinSeq.getStart();
int proteinEnd = proteinSeq.getEnd();
-
+
/*
* incomplete start codon may mean X at start of peptide
* we ignore both for mapping purposes
proteinLength--;
}
List<int[]> proteinRange = new ArrayList<int[]>();
-
+
/*
* dna length should map to protein (or protein plus stop codon)
*/
/**
* 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)
+ * 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
* @return
{
return result;
}
+
SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+ int startPhase = 0;
+
for (SequenceFeature sf : sfs)
{
/*
if (so.isA(sf.getType(), SequenceOntologyI.CDS))
{
int phase = 0;
- try {
+ try
+ {
phase = Integer.parseInt(sf.getPhase());
} catch (NumberFormatException e)
{
int end = sf.getEnd();
if (result.isEmpty())
{
- // TODO JAL-2022 support start phase > 0
begin += phase;
if (begin > end)
{
- continue; // shouldn't happen?
+ // 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<int[]>()
+ {
+ @Override
+ public int compare(int[] o1, int[] o2)
+ {
+ return Integer.compare(o1[0], o2[0]);
+ }
+ });
return result;
}
{
peptide = peptide.getDatasetSequence();
}
-
- transferFeatures(dnaSeq, peptide, dnaToProtein,
- SequenceOntologyI.EXON);
-
- LinkedHashMap<Integer, String[][]> variants = buildDnaVariantsMap(
+
+ 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<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
dnaSeq, dnaToProtein);
-
+
/*
* scan codon variations, compute peptide variants and add to peptide sequence
*/
int count = 0;
- for (Entry<Integer, String[][]> variant : variants.entrySet())
+ for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
{
int peptidePos = variant.getKey();
- String[][] codonVariants = variant.getValue();
- String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); // 0-based
- List<String> peptideVariants = computePeptideVariants(codonVariants,
- residue);
- if (!peptideVariants.isEmpty())
- {
- String desc = StringUtils.listToDelimitedString(peptideVariants,
- ", ");
- SequenceFeature sf = new SequenceFeature(
- SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
- peptidePos, 0f, null);
- peptide.addSequenceFeature(sf);
- count++;
- }
+ List<DnaVariant>[] codonVariants = variant.getValue();
+ count += computePeptideVariants(peptide, peptidePos, codonVariants);
}
-
+
/*
- * ugly sort to get sequence features in start position order
- * - would be better to store in Sequence as a TreeSet instead?
+ * 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>()
}
/**
- * Builds a map whose key is position in the protein sequence, and value is an
- * array of all variants for the coding codon positions
+ * 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<DnaVariant>[] codonVariants)
+ {
+ 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;
+
+ /*
+ * variants in first codon base
+ */
+ for (DnaVariant var : codonVariants[0])
+ {
+ if (var.variant != null)
+ {
+ 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++;
+ }
+ }
+ }
+ }
+ }
+
+ /*
+ * 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++;
+ }
+ }
+ }
+ }
+ }
+
+ /*
+ * 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 count;
+ }
+
+ /**
+ * 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 peptide
+ * @param peptidePos
+ * @param residue
+ * @param var
+ * @param codon
+ * @return true if a feature was added, else false
+ */
+ static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
+ String residue, DnaVariant var, String codon)
+ {
+ /*
+ * 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))
+ {
+ String desc = residue + "->" + trans;
+ // set score to 0f so 'graduated colour' option is offered! JAL-2060
+ SequenceFeature sf = new SequenceFeature(
+ SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
+ peptidePos, 0f, null);
+ 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)
+ {
+ sf.setAttributes(attributes.toString());
+ }
+ return true;
+ }
+ 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<Integer, String[][]> buildDnaVariantsMap(
+ static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
SequenceI dnaSeq, MapList dnaToProtein)
{
/*
- * map from peptide position to all variant features of the codon for it
- * LinkedHashMap ensures we add the peptide features in sequence order
+ * 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<Integer, String[][]> variants = new LinkedHashMap<Integer, String[][]>();
+ LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<Integer, List<DnaVariant>[]>();
SequenceOntologyI so = SequenceOntologyFactory.getInstance();
-
+
SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
if (dnaFeatures == null)
{
return variants;
}
-
+
int dnaStart = dnaSeq.getStart();
int[] lastCodon = null;
int lastPeptidePostion = 0;
-
+
/*
* build a map of codon variations for peptides
*/
continue;
}
int peptidePosition = mapsTo[0];
- String[][] codonVariants = variants.get(peptidePosition);
+ List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
if (codonVariants == null)
{
- codonVariants = new String[3][];
+ codonVariants = new ArrayList[3];
+ codonVariants[0] = new ArrayList<DnaVariant>();
+ codonVariants[1] = new ArrayList<DnaVariant>();
+ codonVariants[2] = new ArrayList<DnaVariant>();
variants.put(peptidePosition, codonVariants);
}
-
+
/*
* extract dna variants to a string array
*/
{
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]
*/
peptidePosition, peptidePosition));
lastPeptidePostion = peptidePosition;
lastCodon = codon;
-
+
/*
- * save nucleotide (and this variant) for each codon position
+ * 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();
- if (codonVariants[codonPos] == null)
+ List<DnaVariant> codonVariant = codonVariants[codonPos];
+ if (codon[codonPos] == dnaCol)
{
- /*
- * record current dna base
- */
- codonVariants[codonPos] = new String[] { nucleotide };
+ 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));
+ }
}
- if (codon[codonPos] == dnaCol)
+ else if (codonVariant.isEmpty())
{
/*
- * add alleles to dna base (and any previously found alleles)
+ * record (possibly non-varying) base value
*/
- String[] known = codonVariants[codonPos];
- String[] dnaVariants = new String[alleles.length + known.length];
- System.arraycopy(known, 0, dnaVariants, 0, known.length);
- System.arraycopy(alleles, 0, dnaVariants, known.length,
- alleles.length);
- codonVariants[codonPos] = dnaVariants;
+ codonVariant.add(new DnaVariant(nucleotide));
}
}
}
}
/**
- * Returns a sorted, non-redundant list of all peptide translations generated
- * by the given dna variants, excluding the current residue value
+ * 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 codonVariants
- * an array of base values (acgtACGT) for codon positions 1, 2, 3
- * @param residue
- * the current residue translation
+ * @param seqs
+ * @param xrefs
* @return
*/
- static List<String> computePeptideVariants(
- String[][] codonVariants, String residue)
+ public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
+ SequenceI[] xrefs)
{
- List<String> result = new ArrayList<String>();
- for (String base1 : codonVariants[0])
+ AlignmentI copy = new Alignment(new Alignment(seqs));
+
+ SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
+ if (xrefs != null)
{
- for (String base2 : codonVariants[1])
+ for (SequenceI xref : xrefs)
{
- for (String base3 : codonVariants[2])
+ DBRefEntry[] dbrefs = xref.getDBRefs();
+ if (dbrefs != null)
{
- String codon = base1 + base2 + base3;
- /*
- * 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 ignored here
- */
- String peptide = codon.contains("-") ? "-"
- : (codon.length() > 3 ? null : ResidueProperties
- .codonTranslate(codon));
- if (peptide != null && !result.contains(peptide)
- && !peptide.equalsIgnoreCase(residue))
+ for (DBRefEntry dbref : dbrefs)
{
- result.add(peptide);
+ 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)
+ {
+ List<SequenceI> unmapped = new ArrayList<SequenceI>();
+ Map<Integer, Map<SequenceI, Character>> 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;
+ }
+
+ /**
+ * 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<Integer, Map<SequenceI, Character>> buildMappedColumnsMap(
+ AlignmentI unaligned, AlignmentI aligned, List<SequenceI> unmapped)
+ {
+ /*
+ * Map will hold, for each aligned column position, a map of
+ * {unalignedSequence, sequenceCharacter} at that position.
+ * TreeMap keeps the entries in ascending column order.
+ */
+ Map<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
+
/*
- * sort alphabetically with STOP at the end
+ * r any sequences that have no mapping so can't be realigned
*/
- Collections.sort(result, new Comparator<String>()
+ unmapped.addAll(unaligned.getSequences());
+
+ List<AlignedCodonFrame> mappings = aligned.getCodonFrames();
+
+ for (SequenceI seq : unaligned.getSequences())
{
-
- @Override
- public int compare(String o1, String o2)
+ for (AlignedCodonFrame mapping : mappings)
{
- if ("STOP".equals(o1))
+ SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned);
+ if (fromSeq != null)
{
- return 1;
+ Mapping seqMap = mapping.getMappingBetween(fromSeq, seq);
+ if (addMappedPositions(seq, fromSeq, seqMap, map))
+ {
+ unmapped.remove(seq);
+ }
}
- else if ("STOP".equals(o2))
+ }
+ }
+ return map;
+ }
+
+ /**
+ * Helper method that adds to a map the mapped column positions of a sequence. <br>
+ * 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<Integer, Map<SequenceI, Character>> map)
+ {
+ if (seqMap == null)
+ {
+ return false;
+ }
+
+ 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)
{
- return -1;
+ System.err.println("Error in mapping " + seqMap + " from "
+ + fromSeq.getName());
+ return false;
}
- else
+ 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])
{
- return o1.compareTo(o2);
+ 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<SequenceI, Character> seqsMap = map.get(fromCol);
+ if (seqsMap == null)
+ {
+ seqsMap = new HashMap<SequenceI, Character>();
+ map.put(fromCol, seqsMap);
+ }
+ seqsMap.put(seq, toChars[mappedCharPos - toStart]);
+ mappedCharPos++;
+ }
+ fromCol += (forward ? 1 : -1);
}
}
- });
- return result;
+ }
+ 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;
}
}
git://source.jalview.org / jalview.git / blobdiff