public class AlignmentUtils
{
+ private static final int CODON_LENGTH = 3;
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
+ static final class DnaVariant
{
- String base;
+ final String base;
SequenceFeature variant;
DnaVariant(String nuc)
{
base = nuc;
+ variant = null;
}
DnaVariant(String nuc, SequenceFeature var)
base = nuc;
variant = var;
}
+
+ public String getSource()
+ {
+ return variant == null ? null : variant.getFeatureGroup();
+ }
}
/**
/*
* cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping)
*/
- final int mappedLength = 3 * aaSeqChars.length;
+ final int mappedLength = CODON_LENGTH * aaSeqChars.length;
int cdnaLength = cdnaSeqChars.length;
int cdnaStart = cdnaSeq.getStart();
int cdnaEnd = cdnaSeq.getEnd();
*/
if (cdnaLength != mappedLength && cdnaLength > 2)
{
- String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3)
+ String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - CODON_LENGTH, CODON_LENGTH)
.toUpperCase();
for (String stop : ResidueProperties.STOP)
{
if (lastCodon.equals(stop))
{
- cdnaEnd -= 3;
- cdnaLength -= 3;
+ cdnaEnd -= CODON_LENGTH;
+ cdnaLength -= CODON_LENGTH;
break;
}
}
int startOffset = 0;
if (cdnaLength != mappedLength
&& cdnaLength > 2
- && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase()
+ && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase()
.equals(ResidueProperties.START))
{
- startOffset += 3;
- cdnaStart += 3;
- cdnaLength -= 3;
+ startOffset += CODON_LENGTH;
+ cdnaStart += CODON_LENGTH;
+ cdnaLength -= CODON_LENGTH;
}
if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars))
* protein is translation of dna (+/- start/stop codons)
*/
MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[]
- { proteinStart, proteinEnd }, 3, 1);
+ { proteinStart, proteinEnd }, CODON_LENGTH, 1);
return map;
}
int aaPos = 0;
int dnaPos = cdnaStart;
for (; dnaPos < cdnaSeqChars.length - 2
- && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++)
+ && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++)
{
- String codon = String.valueOf(cdnaSeqChars, dnaPos, 3);
+ String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
final String translated = ResidueProperties.codonTranslate(codon);
/*
{
return true;
}
- if (dnaPos == cdnaSeqChars.length - 3)
+ if (dnaPos == cdnaSeqChars.length - CODON_LENGTH)
{
- String codon = String.valueOf(cdnaSeqChars, dnaPos, 3);
+ String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
if ("STOP".equals(ResidueProperties.codonTranslate(codon)))
{
return true;
* Realigns the given dna to match the alignment of the protein, using codon
* mappings to translate aligned peptide positions to codons.
*
+ * Always produces a padded CDS alignment.
+ *
* @param dna
* the alignment whose sequences are realigned by this method
* @param protein
// todo: implement this
List<AlignedCodonFrame> mappings = protein.getCodonFrames();
int alignedCount = 0;
+ int width = 0; // alignment width for padding CDS
for (SequenceI dnaSeq : dna.getSequences())
{
- if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings))
+ if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings,
+ dna.getGapCharacter()))
{
alignedCount++;
}
+ width = Math.max(dnaSeq.getLength(), width);
+ }
+ int oldwidth;
+ int diff;
+ for (SequenceI dnaSeq : dna.getSequences())
+ {
+ oldwidth = dnaSeq.getLength();
+ diff = width - oldwidth;
+ if (diff > 0)
+ {
+ dnaSeq.insertCharAt(oldwidth, diff, dna.getGapCharacter());
+ }
}
return alignedCount;
}
/**
- * Helper method to align (if possible) the dna sequence (cds only) to match
- * the alignment of a mapped protein sequence
+ * 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 dnaSeq
+ * @param cdsSeq
* @param protein
* @param mappings
+ * @param gapChar
* @return
*/
- static boolean alignCdsSequenceAsProtein(SequenceI dnaSeq,
- AlignmentI protein, List<AlignedCodonFrame> mappings)
+ 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(dnaSeq, mappings);
+ .findMappingsForSequence(cdsSeq, mappings);
for (AlignedCodonFrame mapping : dnaMappings)
{
- SequenceI peptide = mapping.findAlignedSequence(dnaSeq, protein);
+ SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein);
if (peptide != null)
{
- Mapping map = mapping.getMappingBetween(dnaSeq, peptide);
+ int peptideLength = peptide.getLength();
+ 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 * CODON_LENGTH + CODON_LENGTH)
+ || (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()));
+ }
+
/*
- * traverse peptide adding gaps or codons to new cds sequence
+ * pre-fill the aligned cds sequence with gaps
*/
- MapList mapList = map.getMap();
+ char[] alignedCds = new char[peptideLength * CODON_LENGTH
+ + (addStopCodon ? CODON_LENGTH : 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 += CODON_LENGTH;
+ }
+ else
+ {
+ proteinPos++;
+ int[] codon = mapList.locateInTo(proteinPos, proteinPos);
+ if (codon == null)
+ {
+ // e.g. incomplete start codon, X in peptide
+ cdsCol += CODON_LENGTH;
+ }
+ 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 - CODON_LENGTH)
+ {
+ 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 - CODON_LENGTH; i < nucleotides.length; i++)
+ {
+ alignedCds[cdsCol++] = nucleotides[i];
+ }
+ }
+ cdsSeq.setSequence(new String(alignedCds));
+ return true;
}
}
}
int mappedFromLength = MappingUtils.getLength(aMapping.getMap()
.getFromRanges());
int dnaLength = seqDss.getLength();
- if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - 3)
+ if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - CODON_LENGTH)
{
return seqDss;
}
for (SequenceToSequenceMapping map : acf.getMappings())
{
Mapping mapping = map.getMapping();
- if (mapping != aMapping && mapping.getMap().getFromRatio() == 3
+ if (mapping != aMapping && mapping.getMap().getFromRatio() == CODON_LENGTH
&& proteinProduct == mapping.getTo()
&& seqDss != map.getFromSeq())
{
/*
* dna length should map to protein (or protein plus stop codon)
*/
- int codesForResidues = mappedDnaLength / 3;
+ int codesForResidues = mappedDnaLength / CODON_LENGTH;
if (codesForResidues == (proteinLength + 1))
{
// assuming extra codon is for STOP and not in peptide
if (codesForResidues == proteinLength)
{
proteinRange.add(new int[] { proteinStart, proteinEnd });
- return new MapList(ranges, proteinRange, 3, 1);
+ return new MapList(ranges, proteinRange, CODON_LENGTH, 1);
}
return null;
}
* are currently ignored here
*/
String trans = codon.contains("-") ? "-"
- : (codon.length() > 3 ? null : ResidueProperties
+ : (codon.length() > CODON_LENGTH ? null : ResidueProperties
.codonTranslate(codon));
if (trans != null && !trans.equals(residue))
{
// set score to 0f so 'graduated colour' option is offered! JAL-2060
SequenceFeature sf = new SequenceFeature(
SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
- peptidePos, 0f, "Jalview");
+ peptidePos, 0f, var.getSource());
StringBuilder attributes = new StringBuilder(32);
String id = (String) var.variant.getValue(ID);
if (id != null)
}
sf.setValue(ID, id);
attributes.append(ID).append("=").append(id);
- // TODO handle other species variants
+ // TODO handle other species variants JAL-2064
StringBuilder link = new StringBuilder(32);
try
{
* @param dnaToProtein
* @return
*/
+ @SuppressWarnings("unchecked")
static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
SequenceI dnaSeq, MapList dnaToProtein)
{
List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
if (codonVariants == null)
{
- codonVariants = new ArrayList[3];
+ codonVariants = new ArrayList[CODON_LENGTH];
codonVariants[0] = new ArrayList<DnaVariant>();
codonVariants[1] = new ArrayList<DnaVariant>();
codonVariants[2] = new ArrayList<DnaVariant>();
/*
* save nucleotide (and any variant) for each codon position
*/
- for (int codonPos = 0; codonPos < 3; codonPos++)
+ for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++)
{
String nucleotide = String.valueOf(
dnaSeq.getCharAt(codon[codonPos] - dnaStart))
*/
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);
int width = columnMap.size();
char gap = unaligned.getGapCharacter();
int realignedCount = 0;
+ // TODO: verify this loop scales sensibly for very wide/high alignments
for (SequenceI seq : unaligned.getSequences())
{
if (!unmapped.contains(seq))
{
char[] newSeq = new char[width];
- Arrays.fill(newSeq, gap);
+ Arrays.fill(newSeq, gap); // JBPComment - doubt this is faster than the
+ // Integer iteration below
int newCol = 0;
int lastCol = 0;
System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
newSeq = tmp;
}
+ // TODO: optimise SequenceI to avoid char[]->String->char[]
seq.setSequence(String.valueOf(newSeq));
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
+ * - sequences to be aligned based on aligned
+ * @param aligned
+ * - 'guide' alignment containing sequences derived from same dataset
+ * as unaligned
+ * @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 from dataset sequence to alignment sequence(s)
+ Map<SequenceI, List<SequenceI>> alignedDatasets = new HashMap<SequenceI, List<SequenceI>>();
+ for (SequenceI seq : aligned.getSequences())
+ {
+ SequenceI ds = seq.getDatasetSequence();
+ if (alignedDatasets.get(ds) == null)
+ {
+ alignedDatasets.put(ds, new ArrayList<SequenceI>());
+ }
+ alignedDatasets.get(ds).add(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;
+ * heuristic rule: pair off sequences in order for the case where
+ * more than one shares the same dataset sequence
+ */
+ for (SequenceI seq : unaligned.getSequences())
+ {
+ List<SequenceI> alignedSequences = alignedDatasets.get(seq
+ .getDatasetSequence());
+ // TODO: getSequenceAsString() will be deprecated in the future
+ // TODO: need to leave to SequenceI implementor to update gaps
+ seq.setSequence(alignedSequences.get(0).getSequenceAsString());
+ if (alignedSequences.size() > 0)
+ {
+ // pop off aligned sequences (except the last one)
+ alignedSequences.remove(0);
+ }
+ }
+
+ 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.
*
git://source.jalview.org / jalview.git / blobdiff