X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=55efaa51c71d83467b20609ef0284ffe96463d34;hb=9b85e1552fa57f02cf6cd312cfbd7efdfd079ea3;hp=74066d778593d437ce44a31d947b716d0681b08d;hpb=d2299844ae932a515a5007f30caf766a2c83ad97;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 74066d7..55efaa5 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -22,6 +22,7 @@ package jalview.analysis; import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE; +import jalview.commands.RemoveGapColCommand; import jalview.datamodel.AlignedCodon; import jalview.datamodel.AlignedCodonFrame; import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping; @@ -29,16 +30,20 @@ import jalview.datamodel.Alignment; import jalview.datamodel.AlignmentAnnotation; import jalview.datamodel.AlignmentI; import jalview.datamodel.DBRefEntry; +import jalview.datamodel.GeneLociI; import jalview.datamodel.IncompleteCodonException; import jalview.datamodel.Mapping; import jalview.datamodel.Sequence; import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceGroup; import jalview.datamodel.SequenceI; -import jalview.io.gff.SequenceOntologyFactory; +import jalview.datamodel.features.SequenceFeatures; +import jalview.io.gff.Gff3Helper; import jalview.io.gff.SequenceOntologyI; import jalview.schemes.ResidueProperties; import jalview.util.Comparison; +import jalview.util.DBRefUtils; +import jalview.util.IntRangeComparator; import jalview.util.MapList; import jalview.util.MappingUtils; import jalview.util.StringUtils; @@ -49,7 +54,6 @@ import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; import java.util.Collections; -import java.util.Comparator; import java.util.HashMap; import java.util.HashSet; import java.util.Iterator; @@ -59,6 +63,7 @@ import java.util.Map; import java.util.Map.Entry; import java.util.NoSuchElementException; import java.util.Set; +import java.util.SortedMap; import java.util.TreeMap; /** @@ -70,23 +75,30 @@ import java.util.TreeMap; */ public class AlignmentUtils { + private static final int CODON_LENGTH = 3; private static final String SEQUENCE_VARIANT = "sequence_variant:"; - private static final String ID = "ID"; + + /* + * the 'id' attribute is provided for variant features fetched from + * Ensembl using its REST service with JSON format + */ + public static final String VARIANT_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) @@ -94,6 +106,20 @@ public class AlignmentUtils base = nuc; variant = var; } + + public String getSource() + { + return variant == null ? null : variant.getFeatureGroup(); + } + + /** + * toString for aid in the debugger only + */ + @Override + public String toString() + { + return base + ":" + (variant == null ? "" : variant.getDescription()); + } } /** @@ -106,7 +132,7 @@ public class AlignmentUtils */ public static AlignmentI expandContext(AlignmentI core, int flankSize) { - List sq = new ArrayList(); + List sq = new ArrayList<>(); int maxoffset = 0; for (SequenceI s : core.getSequences()) { @@ -159,10 +185,12 @@ public class AlignmentUtils } } // TODO use Character.toLowerCase to avoid creating String objects? - char[] upstream = new String(ds.getSequence(s.getStart() - 1 - - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray(); - char[] downstream = new String(ds.getSequence(s_end - 1, s_end - + dstream_ds)).toLowerCase().toCharArray(); + char[] upstream = new String(ds + .getSequence(s.getStart() - 1 - ustream_ds, s.getStart() - 1)) + .toLowerCase().toCharArray(); + char[] downstream = new String( + ds.getSequence(s_end - 1, s_end + dstream_ds)).toLowerCase() + .toCharArray(); char[] coreseq = s.getSequence(); char[] nseq = new char[offset + upstream.length + downstream.length + coreseq.length]; @@ -177,8 +205,8 @@ public class AlignmentUtils System.arraycopy(upstream, 0, nseq, p, upstream.length); System.arraycopy(coreseq, 0, nseq, p + upstream.length, coreseq.length); - System.arraycopy(downstream, 0, nseq, p + coreseq.length - + upstream.length, downstream.length); + System.arraycopy(downstream, 0, nseq, + p + coreseq.length + upstream.length, downstream.length); s.setSequence(new String(nseq)); s.setStart(s.getStart() - ustream_ds); s.setEnd(s_end + downstream.length); @@ -234,7 +262,7 @@ public class AlignmentUtils public static Map> getSequencesByName( AlignmentI al) { - Map> theMap = new LinkedHashMap>(); + Map> theMap = new LinkedHashMap<>(); for (SequenceI seq : al.getSequences()) { String name = seq.getName(); @@ -243,7 +271,7 @@ public class AlignmentUtils List seqs = theMap.get(name); if (seqs == null) { - seqs = new ArrayList(); + seqs = new ArrayList<>(); theMap.put(name, seqs); } seqs.add(seq); @@ -270,8 +298,8 @@ public class AlignmentUtils return false; } - Set mappedDna = new HashSet(); - Set mappedProtein = new HashSet(); + Set mappedDna = new HashSet<>(); + Set mappedProtein = new HashSet<>(); /* * First pass - map sequences where cross-references exist. This include @@ -305,9 +333,9 @@ public class AlignmentUtils * @return */ protected static boolean mapProteinToCdna( - final AlignmentI proteinAlignment, - final AlignmentI cdnaAlignment, Set mappedDna, - Set mappedProtein, boolean xrefsOnly) + final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment, + Set mappedDna, Set mappedProtein, + boolean xrefsOnly) { boolean mappingExistsOrAdded = false; List thisSeqs = proteinAlignment.getSequences(); @@ -336,9 +364,8 @@ public class AlignmentUtils * Don't map non-xrefd sequences more than once each. This heuristic * allows us to pair up similar sequences in ordered alignments. */ - if (!xrefsOnly - && (mappedProtein.contains(aaSeq) || mappedDna - .contains(cdnaSeq))) + if (!xrefsOnly && (mappedProtein.contains(aaSeq) + || mappedDna.contains(cdnaSeq))) { continue; } @@ -372,7 +399,7 @@ public class AlignmentUtils * Answers true if the mappings include one between the given (dataset) * sequences. */ - public static boolean mappingExists(List mappings, + protected static boolean mappingExists(List mappings, SequenceI aaSeq, SequenceI cdnaSeq) { if (mappings != null) @@ -391,7 +418,8 @@ public class AlignmentUtils /** * Builds a mapping (if possible) of a cDNA to a protein sequence. *
    - *
  • first checks if the cdna translates exactly to the protein sequence
    • + *
  • first checks if the cdna translates exactly to the protein + * sequence
    • *
  • else checks for translation after removing a STOP codon
    • *
  • else checks for translation after removing a START codon
    • *
  • if that fails, inspect CDS features on the cDNA sequence
    • @@ -413,8 +441,9 @@ public class AlignmentUtils * String objects. */ final SequenceI proteinDataset = proteinSeq.getDatasetSequence(); - char[] aaSeqChars = proteinDataset != null ? proteinDataset - .getSequence() : proteinSeq.getSequence(); + char[] aaSeqChars = proteinDataset != null + ? proteinDataset.getSequence() + : proteinSeq.getSequence(); final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence(); char[] cdnaSeqChars = cdnaDataset != null ? cdnaDataset.getSequence() : cdnaSeq.getSequence(); @@ -426,7 +455,7 @@ public class AlignmentUtils /* * 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(); @@ -438,14 +467,14 @@ public class AlignmentUtils */ if (cdnaLength != mappedLength && cdnaLength > 2) { - String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3) - .toUpperCase(); - for (String stop : ResidueProperties.STOP) + String lastCodon = String.valueOf(cdnaSeqChars, + cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(); + for (String stop : ResidueProperties.STOP_CODONS) { if (lastCodon.equals(stop)) { - cdnaEnd -= 3; - cdnaLength -= 3; + cdnaEnd -= CODON_LENGTH; + cdnaLength -= CODON_LENGTH; break; } } @@ -455,14 +484,13 @@ public class AlignmentUtils * If lengths still don't match, try ignoring start codon. */ int startOffset = 0; - if (cdnaLength != mappedLength - && cdnaLength > 2 - && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase() + if (cdnaLength != mappedLength && cdnaLength > 2 + && 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)) @@ -470,8 +498,9 @@ public class AlignmentUtils /* * protein is translation of dna (+/- start/stop codons) */ - MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] - { proteinStart, proteinEnd }, 3, 1); + MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, + new int[] + { proteinStart, proteinEnd }, CODON_LENGTH, 1); return map; } @@ -502,16 +531,17 @@ public class AlignmentUtils 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); /* * allow * in protein to match untranslatable in dna */ final char aaRes = aaSeqChars[aaPos]; - if ((translated == null || "STOP".equals(translated)) && aaRes == '*') + if ((translated == null || ResidueProperties.STOP.equals(translated)) + && aaRes == '*') { continue; } @@ -540,10 +570,11 @@ public class AlignmentUtils { return true; } - if (dnaPos == cdnaSeqChars.length - 3) + if (dnaPos == cdnaSeqChars.length - CODON_LENGTH) { - String codon = String.valueOf(cdnaSeqChars, dnaPos, 3); - if ("STOP".equals(ResidueProperties.codonTranslate(codon))) + String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); + if (ResidueProperties.STOP + .equals(ResidueProperties.codonTranslate(codon))) { return true; } @@ -621,10 +652,9 @@ public class AlignmentUtils * @param preserveUnmappedGaps * @param preserveMappedGaps */ - public static void alignSequenceAs(SequenceI alignTo, - SequenceI alignFrom, AlignedCodonFrame mapping, String myGap, - char sourceGap, boolean preserveMappedGaps, - boolean preserveUnmappedGaps) + public static void alignSequenceAs(SequenceI alignTo, SequenceI alignFrom, + AlignedCodonFrame mapping, String myGap, char sourceGap, + boolean preserveMappedGaps, boolean preserveUnmappedGaps) { // TODO generalise to work for Protein-Protein, dna-dna, dna-protein @@ -640,15 +670,16 @@ public class AlignmentUtils int toOffset = alignTo.getStart() - 1; int sourceGapMappedLength = 0; boolean inExon = false; - final char[] thisSeq = alignTo.getSequence(); - final char[] thatAligned = alignFrom.getSequence(); - StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length); + final int toLength = alignTo.getLength(); + final int fromLength = alignFrom.getLength(); + StringBuilder thisAligned = new StringBuilder(2 * toLength); /* * Traverse the 'model' aligned sequence */ - for (char sourceChar : thatAligned) + for (int i = 0; i < fromLength; i++) { + char sourceChar = alignFrom.getCharAt(i); if (sourceChar == sourceGap) { sourceGapMappedLength += ratio; @@ -688,9 +719,9 @@ public class AlignmentUtils */ int intronLength = 0; while (basesWritten + toOffset < mappedCodonEnd - && thisSeqPos < thisSeq.length) + && thisSeqPos < toLength) { - final char c = thisSeq[thisSeqPos++]; + final char c = alignTo.getCharAt(thisSeqPos++); if (c != myGapChar) { basesWritten++; @@ -716,7 +747,7 @@ public class AlignmentUtils int gapsToAdd = calculateGapsToInsert(preserveMappedGaps, preserveUnmappedGaps, sourceGapMappedLength, inExon, trailingCopiedGap.length(), intronLength, startOfCodon); - for (int i = 0; i < gapsToAdd; i++) + for (int k = 0; k < gapsToAdd; k++) { thisAligned.append(myGapChar); } @@ -744,9 +775,9 @@ public class AlignmentUtils * At end of model aligned sequence. Copy any remaining target sequence, optionally * including (intron) gaps. */ - while (thisSeqPos < thisSeq.length) + while (thisSeqPos < toLength) { - final char c = thisSeq[thisSeqPos++]; + final char c = alignTo.getCharAt(thisSeqPos++); if (c != myGapChar || preserveUnmappedGaps) { thisAligned.append(c); @@ -819,8 +850,9 @@ public class AlignmentUtils } else { - gapsToAdd = Math.min(intronLength + trailingGapLength - - sourceGapMappedLength, trailingGapLength); + gapsToAdd = Math.min( + intronLength + trailingGapLength - sourceGapMappedLength, + trailingGapLength); } } } @@ -850,13 +882,191 @@ public class AlignmentUtils */ public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) { - List unmappedProtein = new ArrayList(); + if (protein.isNucleotide() || !dna.isNucleotide()) + { + System.err.println("Wrong alignment type in alignProteinAsDna"); + return 0; + } + List unmappedProtein = new ArrayList<>(); Map> alignedCodons = buildCodonColumnsMap( protein, dna, unmappedProtein); return alignProteinAs(protein, alignedCodons, unmappedProtein); } /** + * 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 + * 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 mappings = protein.getCodonFrames(); + int alignedCount = 0; + int width = 0; // alignment width for padding CDS + for (SequenceI dnaSeq : dna.getSequences()) + { + 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 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 mappings, + char gapChar) + { + SequenceI cdsDss = cdsSeq.getDatasetSequence(); + if (cdsDss == null) + { + System.err + .println("alignCdsSequenceAsProtein needs aligned sequence!"); + return false; + } + + List dnaMappings = MappingUtils + .findMappingsForSequence(cdsSeq, mappings); + for (AlignedCodonFrame mapping : dnaMappings) + { + SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); + if (peptide != null) + { + final 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(); + } + final 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())); + } + + /* + * pre-fill the aligned cds sequence with gaps + */ + 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 + */ + int copiedBases = 0; + int cdsStart = cdsDss.getStart(); + int proteinPos = peptide.getStart() - 1; + int cdsCol = 0; + + for (int col = 0; col < peptideLength; col++) + { + char residue = peptide.getCharAt(col); + + 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 = cdsDss.getCharAt(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 == cdsLength - 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 = cdsLength - CODON_LENGTH; i < cdsLength; i++) + { + alignedCds[cdsCol++] = cdsDss.getCharAt(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 @@ -888,7 +1098,7 @@ public class AlignmentUtils * {dnaSequence, {proteinSequence, codonProduct}} at that position. The * comparator keeps the codon positions ordered. */ - Map> alignedCodons = new TreeMap>( + Map> alignedCodons = new TreeMap<>( new CodonComparator()); for (SequenceI dnaSeq : dna.getSequences()) @@ -899,8 +1109,8 @@ public class AlignmentUtils if (prot != null) { Mapping seqMap = mapping.getMappingForSequence(dnaSeq); - addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), - seqMap, alignedCodons); + addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), seqMap, + alignedCodons); unmappedProtein.remove(prot); } } @@ -913,7 +1123,7 @@ public class AlignmentUtils // TODO resolve JAL-2022 so this fudge can be removed int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); - + return alignedCodons; } @@ -934,9 +1144,9 @@ public class AlignmentUtils // TODO delete this ugly hack once JAL-2022 is resolved // i.e. we can model startPhase > 0 (incomplete start codon) - List sequencesChecked = new ArrayList(); + List sequencesChecked = new ArrayList<>(); AlignedCodon lastCodon = null; - Map toAdd = new HashMap(); + Map toAdd = new HashMap<>(); for (Entry> entry : alignedCodons .entrySet()) @@ -953,8 +1163,8 @@ public class AlignmentUtils AlignedCodon codon = sequenceCodon.getValue(); if (codon.peptideCol > 1) { - System.err - .println("Problem mapping protein with >1 unmapped start positions: " + System.err.println( + "Problem mapping protein with >1 unmapped start positions: " + seq.getName()); } else if (codon.peptideCol == 1) @@ -965,8 +1175,8 @@ public class AlignmentUtils if (lastCodon != null) { AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1, - lastCodon.pos2, lastCodon.pos3, String.valueOf(seq - .getCharAt(0)), 0); + lastCodon.pos2, lastCodon.pos3, + String.valueOf(seq.getCharAt(0)), 0); toAdd.put(seq, firstPeptide); } else @@ -1015,21 +1225,26 @@ public class AlignmentUtils List unmappedProtein) { /* - * Prefill aligned sequences with gaps before inserting aligned protein - * residues. + * prefill peptide sequences with gaps */ int alignedWidth = alignedCodons.size(); char[] gaps = new char[alignedWidth]; Arrays.fill(gaps, protein.getGapCharacter()); - String allGaps = String.valueOf(gaps); + Map peptides = new HashMap<>(); for (SequenceI seq : protein.getSequences()) { if (!unmappedProtein.contains(seq)) { - seq.setSequence(allGaps); + peptides.put(seq, Arrays.copyOf(gaps, gaps.length)); } } + /* + * Traverse the codons left to right (as defined by CodonComparator) + * and insert peptides in each column where the sequence is mapped. + * This gives a peptide 'alignment' where residues are aligned if their + * corresponding codons occupy the same columns in the cdna alignment. + */ int column = 0; for (AlignedCodon codon : alignedCodons.keySet()) { @@ -1037,12 +1252,20 @@ public class AlignmentUtils .get(codon); for (Entry entry : columnResidues.entrySet()) { - // place translated codon at its column position in sequence - entry.getKey().getSequence()[column] = entry.getValue().product - .charAt(0); + char residue = entry.getValue().product.charAt(0); + peptides.get(entry.getKey())[column] = residue; } column++; } + + /* + * and finally set the constructed sequences + */ + for (Entry entry : peptides.entrySet()) + { + entry.getKey().setSequence(new String(entry.getValue())); + } + return 0; } @@ -1102,7 +1325,7 @@ public class AlignmentUtils Map seqProduct = alignedCodons.get(codon); if (seqProduct == null) { - seqProduct = new HashMap(); + seqProduct = new HashMap<>(); alignedCodons.put(codon, seqProduct); } seqProduct.put(protein, codon); @@ -1115,7 +1338,8 @@ public class AlignmentUtils *
      *
    • One alignment must be nucleotide, and the other protein
      • *
    • At least one pair of sequences must be already mapped, or mappable
      • - *
    • Mappable means the nucleotide translation matches the protein sequence
      • + *
    • Mappable means the nucleotide translation matches the protein + * sequence
      • *
    • The translation may ignore start and stop codons if present in the * nucleotide
      • *
    @@ -1171,9 +1395,10 @@ public class AlignmentUtils return false; } - SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq - .getDatasetSequence(); - SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq + SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + SequenceI proteinDs = proteinSeq.getDatasetSequence() == null + ? proteinSeq : proteinSeq.getDatasetSequence(); for (AlignedCodonFrame mapping : mappings) @@ -1210,8 +1435,7 @@ public class AlignmentUtils * the alignment to check for presence of annotations */ public static void findAddableReferenceAnnotations( - List sequenceScope, - Map labelForCalcId, + List sequenceScope, Map labelForCalcId, final Map> candidates, AlignmentI al) { @@ -1238,7 +1462,7 @@ public class AlignmentUtils { continue; } - final List result = new ArrayList(); + final List result = new ArrayList<>(); for (AlignmentAnnotation dsann : datasetAnnotations) { /* @@ -1315,8 +1539,8 @@ public class AlignmentUtils /** * Set visibility of alignment annotations of specified types (labels), for - * specified sequences. This supports controls like - * "Show all secondary structure", "Hide all Temp factor", etc. + * specified sequences. This supports controls like "Show all secondary + * structure", "Hide all Temp factor", etc. * * @al the alignment to scan for annotations * @param types @@ -1340,9 +1564,8 @@ public class AlignmentUtils { if (anyType || types.contains(aa.label)) { - if ((aa.sequenceRef != null) - && (forSequences == null || forSequences - .contains(aa.sequenceRef))) + if ((aa.sequenceRef != null) && (forSequences == null + || forSequences.contains(aa.sequenceRef))) { aa.visible = doShow; } @@ -1404,43 +1627,64 @@ public class AlignmentUtils * added to the alignment dataset. * * @param dna - * aligned dna sequences + * aligned nucleotide (dna or cds) sequences * @param dataset - * - throws error if not given a dataset + * the alignment dataset the sequences belong to + * @param products + * (optional) to restrict results to CDS that map to specified + * protein products * @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, - AlignmentI dataset) + AlignmentI dataset, SequenceI[] products) { - if (dataset.getDataset() != null) + if (dataset == null || dataset.getDataset() != null) { - throw new Error( + throw new IllegalArgumentException( "IMPLEMENTATION ERROR: dataset.getDataset() must be null!"); } - List cdsSeqs = new ArrayList(); + List foundSeqs = new ArrayList<>(); + List cdsSeqs = new ArrayList<>(); List mappings = dataset.getCodonFrames(); - + HashSet productSeqs = null; + if (products != null) + { + productSeqs = new HashSet<>(); + for (SequenceI seq : products) + { + productSeqs.add(seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence()); + } + } /* - * construct CDS sequences from the (cds-to-protein) mappings made earlier; - * this makes it possible to model multiple products from dna (e.g. EMBL); - * however it does mean we don't have the EMBL protein_id (a property on - * the CDS features) in order to make the CDS sequence name :-( + * Construct CDS sequences from mappings on the alignment dataset. + * The logic is: + * - find the protein product(s) mapped to from each dna sequence + * - if the mapping covers the whole dna sequence (give or take start/stop + * codon), take the dna as the CDS sequence + * - else search dataset mappings for a suitable dna sequence, i.e. one + * whose whole sequence is mapped to the protein + * - if no sequence found, construct one from the dna sequence and mapping + * (and add it to dataset so it is found if this is repeated) */ - for (SequenceI seq : dna) + for (SequenceI dnaSeq : dna) { - SequenceI seqDss = seq.getDatasetSequence() == null ? seq : seq - .getDatasetSequence(); + SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + List seqMappings = MappingUtils - .findMappingsForSequence(seq, mappings); + .findMappingsForSequence(dnaSeq, mappings); for (AlignedCodonFrame mapping : seqMappings) { - List mappingsFromSequence = mapping.getMappingsFromSequence(seq); + List mappingsFromSequence = mapping + .getMappingsFromSequence(dnaSeq); for (Mapping aMapping : mappingsFromSequence) { - if (aMapping.getMap().getFromRatio() == 1) + MapList mapList = aMapping.getMap(); + if (mapList.getFromRatio() == 1) { /* * not a dna-to-protein mapping (likely dna-to-cds) @@ -1449,43 +1693,32 @@ public class AlignmentUtils } /* - * check for an existing CDS sequence i.e. a 3:1 mapping to - * the dna mapping's product + * skip if mapping is not to one of the target set of proteins */ - SequenceI cdsSeq = null; - // TODO better mappings collection data model so we can do - // a table lookup instead of double loops to find mappings SequenceI proteinProduct = aMapping.getTo(); - for (AlignedCodonFrame acf : MappingUtils - .findMappingsForSequence(proteinProduct, mappings)) + if (productSeqs != null && !productSeqs.contains(proteinProduct)) { - for (SequenceToSequenceMapping map : acf.getMappings()) - { - if (map.getMapping().getMap().getFromRatio() == 3 - && proteinProduct == map.getMapping().getTo() - && seqDss != map.getFromSeq()) - { - /* - * found a 3:1 mapping to the protein product which is not - * from the dna sequence...assume it is from the CDS sequence - * TODO mappings data model that brings together related - * dna-cds-protein mappings in one object - */ - cdsSeq = map.getFromSeq(); - } - } + 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) { - /* - * mappings are always to dataset sequences so create an aligned - * sequence to own it; add the dataset sequence to the dataset - */ - SequenceI derivedSequence = cdsSeq.deriveSequence(); - cdsSeqs.add(derivedSequence); - if (!dataset.getSequences().contains(cdsSeq)) + if (!foundSeqs.contains(cdsSeq)) { - dataset.addSequence(cdsSeq); + foundSeqs.add(cdsSeq); + SequenceI derivedSequence = cdsSeq.deriveSequence(); + cdsSeqs.add(derivedSequence); + if (!dataset.getSequences().contains(cdsSeq)) + { + dataset.addSequence(cdsSeq); + } } continue; } @@ -1494,24 +1727,36 @@ public class AlignmentUtils * didn't find mapped CDS sequence - construct it and add * its dataset sequence to the dataset */ - cdsSeq = makeCdsSequence(seq.getDatasetSequence(), aMapping); - SequenceI cdsSeqDss = cdsSeq.createDatasetSequence(); + cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping, + dataset).deriveSequence(); + // cdsSeq has a name constructed as CDS| + // will be either the accession for the coding sequence, + // marked in the /via/ dbref to the protein product accession + // or it will be the original nucleotide accession. + SequenceI cdsSeqDss = cdsSeq.getDatasetSequence(); + cdsSeqs.add(cdsSeq); + if (!dataset.getSequences().contains(cdsSeqDss)) { + // check if this sequence is a newly created one + // so needs adding to the dataset dataset.addSequence(cdsSeqDss); } /* * add a mapping from CDS to the (unchanged) mapped to range */ - List cdsRange = Collections.singletonList(new int[] { 1, - cdsSeq.getLength() }); - MapList map = new MapList(cdsRange, aMapping.getMap() - .getToRanges(), aMapping.getMap().getFromRatio(), - aMapping.getMap().getToRatio()); + List cdsRange = Collections + .singletonList(new int[] + { cdsSeq.getStart(), + cdsSeq.getLength() + cdsSeq.getStart() - 1 }); + MapList cdsToProteinMap = new MapList(cdsRange, + mapList.getToRanges(), mapList.getFromRatio(), + mapList.getToRatio()); AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); - cdsToProteinMapping.addMap(cdsSeq, proteinProduct, map); + cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct, + cdsToProteinMap); /* * guard against duplicating the mapping if repeating this action @@ -1521,482 +1766,645 @@ public class AlignmentUtils mappings.add(cdsToProteinMapping); } + propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(), + proteinProduct, aMapping); /* * add another mapping from original 'from' range to CDS */ - AlignedCodonFrame dnaToProteinMapping = new AlignedCodonFrame(); - map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1, - 1); - dnaToProteinMapping.addMap(seq.getDatasetSequence(), cdsSeq, map); - if (!mappings.contains(dnaToProteinMapping)) + AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); + final MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), + cdsRange, 1, 1); + dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, + dnaToCdsMap); + if (!mappings.contains(dnaToCdsMapping)) { - mappings.add(dnaToProteinMapping); + mappings.add(dnaToCdsMapping); } + /* + * transfer dna chromosomal loci (if known) to the CDS + * sequence (via the mapping) + */ + final MapList cdsToDnaMap = dnaToCdsMap.getInverse(); + transferGeneLoci(dnaSeq, cdsToDnaMap, cdsSeq); + + /* + * 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 + */ + + // specific use case: + // Genomic contig ENSCHR:1, contains coding regions for ENSG01, + // ENSG02, ENSG03, with transcripts and products similarly named. + // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01 + + // JBPNote: ?? can't actually create an example that demonstrates we + // need to + // synthesize an xref. + + for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs()) + { + /* + * create a cross-reference from CDS to the source sequence's + * primary reference and vice versa + */ + String source = primRef.getSource(); + String version = primRef.getVersion(); + DBRefEntry cdsCrossRef = new DBRefEntry(source, source + ":" + + version, primRef.getAccessionId()); + cdsCrossRef.setMap(new Mapping(dnaDss, new MapList(cdsToDnaMap))); + cdsSeqDss.addDBRef(cdsCrossRef); + + dnaSeq.addDBRef(new DBRefEntry(source, version, cdsSeq + .getName(), new Mapping(cdsSeqDss, dnaToCdsMap))); + + // problem here is that the cross-reference is synthesized - + // cdsSeq.getName() may be like 'CDS|dnaaccession' or + // 'CDS|emblcdsacc' + // assuming cds version same as dna ?!? + + DBRefEntry proteinToCdsRef = new DBRefEntry(source, version, + 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, dnaToCdsMap, null, + SequenceOntologyI.CDS); } } } AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs .size()])); - cds.setDataset((Alignment) dataset); + cds.setDataset(dataset); return cds; } /** - * 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). + * Tries to transfer gene loci (dbref to chromosome positions) from fromSeq to + * toSeq, mediated by the given mapping between the sequences * - * @param seq - * @param mapping - * @return CDS sequence (as a dataset sequence) + * @param fromSeq + * @param targetToFrom + * Map + * @param targetSeq + */ + protected static void transferGeneLoci(SequenceI fromSeq, + MapList targetToFrom, SequenceI targetSeq) + { + if (targetSeq.getGeneLoci() != null) + { + // already have - don't override + return; + } + GeneLociI fromLoci = fromSeq.getGeneLoci(); + if (fromLoci == null) + { + return; + } + + MapList newMap = targetToFrom.traverse(fromLoci.getMapping()); + + if (newMap != null) + { + targetSeq.setGeneLoci(fromLoci.getSpeciesId(), + fromLoci.getAssemblyId(), fromLoci.getChromosomeId(), newMap); + } + } + + /** + * 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 + * a transcript-to-peptide mapping + * @return */ - static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) + static SequenceI findCdsForProtein(List mappings, + SequenceI dnaSeq, List seqMappings, + Mapping aMapping) { - char[] seqChars = seq.getSequence(); - List fromRanges = mapping.getMap().getFromRanges(); - int cdsWidth = MappingUtils.getLength(fromRanges); - char[] newSeqChars = new char[cdsWidth]; + /* + * 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(); - int newPos = 0; - for (int[] range : fromRanges) + /* + * 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 - CODON_LENGTH) { - if (range[0] <= range[1]) + /* + * if sequence has CDS features, this is a transcript with no UTR + * - do not take this as the CDS sequence! (JAL-2789) + */ + if (seqDss.getFeatures().getFeaturesByOntology(SequenceOntologyI.CDS) + .isEmpty()) { - // 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; + return seqDss; } - else + } + + /* + * 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()) { - // reverse strand mapping - copy and complement one by one - for (int i = range[0]; i >= range[1]; i--) + Mapping mapping = map.getMapping(); + if (mapping != aMapping + && mapping.getMap().getFromRatio() == CODON_LENGTH + && proteinProduct == mapping.getTo() + && seqDss != map.getFromSeq()) { - newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); + 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 the CDS + * is mapped from the given dna start sequence + */ + SequenceI cdsSeq = map.getFromSeq(); + // todo this test is weak if seqMappings contains multiple mappings; + // we get away with it if transcript:cds relationship is 1:1 + List dnaToCdsMaps = MappingUtils + .findMappingsForSequence(cdsSeq, seqMappings); + if (!dnaToCdsMaps.isEmpty()) + { + return cdsSeq; + } + } } } } - - SequenceI newSeq = new Sequence(seq.getName() + "|" - + mapping.getTo().getName(), newSeqChars, 1, newPos); - return newSeq; + return null; } /** - * 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. + * 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 fromSeq - * @param toSeq - * @param select - * if not null, only features of this type are copied (including - * subtypes in the Sequence Ontology) + * @param seq * @param mapping - * the mapping from 'fromSeq' to 'toSeq' - * @param omitting + * @param dataset + * - existing dataset. We check for sequences that look like the CDS + * we are about to construct, if one exists already, then we will + * just return that one. + * @return CDS sequence (as a dataset sequence) */ - public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, - MapList mapping, String select, String... omitting) + static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, + AlignmentI dataset) { - SequenceI copyTo = toSeq; - while (copyTo.getDatasetSequence() != null) + /* + * construct CDS sequence name as "CDS|" with 'from id' held in the mapping + * if set (e.g. EMBL protein_id), else sequence name appended + */ + String mapFromId = mapping.getMappedFromId(); + final String seqId = "CDS|" + + (mapFromId != null ? mapFromId : seq.getName()); + + SequenceI newSeq = null; + + final MapList maplist = mapping.getMap(); + if (maplist.isContiguous() && maplist.isFromForwardStrand()) { - copyTo = copyTo.getDatasetSequence(); + /* + * just a subsequence, keep same dataset sequence + */ + int start = maplist.getFromLowest(); + int end = maplist.getFromHighest(); + newSeq = seq.getSubSequence(start - 1, end); + newSeq.setName(seqId); } - - SequenceOntologyI so = SequenceOntologyFactory.getInstance(); - int count = 0; - SequenceFeature[] sfs = fromSeq.getSequenceFeatures(); - if (sfs != null) + else { - for (SequenceFeature sf : sfs) + /* + * construct by splicing mapped from ranges + */ + char[] seqChars = seq.getSequence(); + List fromRanges = maplist.getFromRanges(); + int cdsWidth = MappingUtils.getLength(fromRanges); + char[] newSeqChars = new char[cdsWidth]; + + int newPos = 0; + for (int[] range : fromRanges) { - String type = sf.getType(); - if (select != null && !so.isA(type, select)) + if (range[0] <= range[1]) { - continue; + // 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; } - boolean omit = false; - for (String toOmit : omitting) + else { - if (type.equals(toOmit)) + // reverse strand mapping - copy and complement one by one + for (int i = range[0]; i >= range[1]; i--) { - omit = true; + newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); } } - if (omit) - { - continue; - } + } - /* - * locate the mapped range - null if either start or end is - * not mapped (no partial overlaps are calculated) - */ - 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) + newSeq = new Sequence(seqId, newSeqChars, 1, newPos); + } + + if (dataset != null) + { + SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName()); + if (matches != null) + { + boolean matched = false; + for (SequenceI mtch : matches) { - mappedTo = mapping.locateInTo(end, end); - if (mappedTo != null) + if (mtch.getStart() != newSeq.getStart()) { - /* - * end of exon is in CDS range - 5' overlap - * to a range from the start of the peptide - */ - mappedTo[0] = 1; + continue; } - } - if (mappedTo == null) - { - mappedTo = mapping.locateInTo(start, start); - if (mappedTo != null) + if (mtch.getEnd() != newSeq.getEnd()) { - /* - * start of exon is in CDS range - 3' overlap - * to a range up to the end of the peptide - */ - mappedTo[1] = toSeq.getLength(); + continue; + } + if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence())) + { + continue; + } + if (!matched) + { + matched = true; + newSeq = mtch; + } + else + { + System.err.println( + "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + + mtch.toString()); } - } - 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; + // newSeq.setDescription(mapFromId); + + return newSeq; } /** - * Returns a mapping from dna to protein by inspecting sequence features of - * type "CDS" on the dna. + * Adds any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to + * the given mapping. * - * @param dnaSeq - * @param proteinSeq - * @return + * @param cdsSeq + * @param contig + * @param proteinProduct + * @param mapping + * @return list of DBRefEntrys added */ - public static MapList mapCdsToProtein(SequenceI dnaSeq, - SequenceI proteinSeq) + protected static List propagateDBRefsToCDS(SequenceI cdsSeq, + SequenceI contig, SequenceI proteinProduct, Mapping mapping) { - List ranges = findCdsPositions(dnaSeq); - int mappedDnaLength = MappingUtils.getLength(ranges); - int proteinLength = proteinSeq.getLength(); - int proteinStart = proteinSeq.getStart(); - int proteinEnd = proteinSeq.getEnd(); + // gather direct refs from contig congruent with mapping + List direct = new ArrayList<>(); + HashSet directSources = new HashSet<>(); - /* - * incomplete start codon may mean X at start of peptide - * we ignore both for mapping purposes - */ - if (proteinSeq.getCharAt(0) == 'X') + if (contig.getDBRefs() != null) { - // todo JAL-2022 support startPhase > 0 - proteinStart++; - proteinLength--; - } - List proteinRange = new ArrayList(); - - /* - * dna length should map to protein (or protein plus stop codon) - */ - int codesForResidues = mappedDnaLength / 3; - if (codesForResidues == (proteinLength + 1)) - { - // assuming extra codon is for STOP and not in peptide - codesForResidues--; + for (DBRefEntry dbr : contig.getDBRefs()) + { + if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap()) + { + MapList map = dbr.getMap().getMap(); + // check if map is the CDS mapping + if (mapping.getMap().equals(map)) + { + direct.add(dbr); + directSources.add(dbr.getSource()); + } + } + } } - if (codesForResidues == proteinLength) - { - proteinRange.add(new int[] { proteinStart, proteinEnd }); - return new MapList(ranges, proteinRange, 3, 1); + DBRefEntry[] onSource = DBRefUtils.selectRefs( + proteinProduct.getDBRefs(), + directSources.toArray(new String[0])); + List propagated = new ArrayList<>(); + + // and generate appropriate mappings + for (DBRefEntry cdsref : direct) + { + // clone maplist and mapping + MapList cdsposmap = new MapList( + Arrays.asList(new int[][] + { new int[] { cdsSeq.getStart(), cdsSeq.getEnd() } }), + cdsref.getMap().getMap().getToRanges(), 3, 1); + Mapping cdsmap = new Mapping(cdsref.getMap().getTo(), + cdsref.getMap().getMap()); + + // create dbref + DBRefEntry newref = new DBRefEntry(cdsref.getSource(), + cdsref.getVersion(), cdsref.getAccessionId(), + new Mapping(cdsmap.getTo(), cdsposmap)); + + // and see if we can map to the protein product for this mapping. + // onSource is the filtered set of accessions on protein that we are + // tranferring, so we assume accession is the same. + if (cdsmap.getTo() == null && onSource != null) + { + List sourceRefs = DBRefUtils.searchRefs(onSource, + cdsref.getAccessionId()); + if (sourceRefs != null) + { + for (DBRefEntry srcref : sourceRefs) + { + if (srcref.getSource().equalsIgnoreCase(cdsref.getSource())) + { + // we have found a complementary dbref on the protein product, so + // update mapping's getTo + newref.getMap().setTo(proteinProduct); + } + } + } + } + cdsSeq.addDBRef(newref); + propagated.add(newref); } - return null; + return propagated; } /** - * 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. + * 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 dnaSeq - * @return + * @param fromSeq + * @param toSeq + * @param mapping + * the mapping from 'fromSeq' to 'toSeq' + * @param select + * if not null, only features of this type are copied (including + * subtypes in the Sequence Ontology) + * @param omitting */ - public static List findCdsPositions(SequenceI dnaSeq) + protected static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, + MapList mapping, String select, String... omitting) { - List result = new ArrayList(); - SequenceFeature[] sfs = dnaSeq.getSequenceFeatures(); - if (sfs == null) + SequenceI copyTo = toSeq; + while (copyTo.getDatasetSequence() != null) { - return result; + copyTo = copyTo.getDatasetSequence(); + } + if (fromSeq == copyTo || fromSeq.getDatasetSequence() == copyTo) + { + return 0; // shared dataset sequence } - SequenceOntologyI so = SequenceOntologyFactory.getInstance(); - int startPhase = 0; + /* + * get features, optionally restricted by an ontology term + */ + List sfs = select == null ? fromSeq.getFeatures() + .getPositionalFeatures() : fromSeq.getFeatures() + .getFeaturesByOntology(select); + int count = 0; for (SequenceFeature sf : sfs) { + String type = sf.getType(); + boolean omit = false; + for (String toOmit : omitting) + { + if (type.equals(toOmit)) + { + omit = true; + } + } + if (omit) + { + continue; + } + /* - * process a CDS feature (or a sub-type of CDS) + * locate the mapped range - null if either start or end is + * not mapped (no partial overlaps are calculated) */ - if (so.isA(sf.getType(), SequenceOntologyI.CDS)) + 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) { - int phase = 0; - try - { - phase = Integer.parseInt(sf.getPhase()); - } catch (NumberFormatException e) + mappedTo = mapping.locateInTo(end, end); + if (mappedTo != null) { - // ignore + /* + * end of exon is in CDS range - 5' overlap + * to a range from the start of the peptide + */ + mappedTo[0] = 1; } - /* - * 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()) + } + if (mappedTo == null) + { + mappedTo = mapping.locateInTo(start, start); + if (mappedTo != null) { - begin += phase; - if (begin > end) - { - // shouldn't happen! - System.err - .println("Error: start phase extends beyond start CDS in " - + dnaSeq.getName()); - } + /* + * start of exon is in CDS range - 3' overlap + * to a range up to the end of the peptide + */ + mappedTo[1] = toSeq.getLength(); } - 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) + if (mappedTo != null) { - return Integer.compare(o1[0], o2[0]); + int newBegin = Math.min(mappedTo[0], mappedTo[1]); + int newEnd = Math.max(mappedTo[0], mappedTo[1]); + SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd, + sf.getFeatureGroup(), sf.getScore()); + copyTo.addSequenceFeature(copy); + count++; } - }); - return result; + } + return count; } /** - * 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. + * Returns a mapping from dna to protein by inspecting sequence features of + * type "CDS" on the dna. A mapping is constructed if the total CDS feature + * length is 3 times the peptide length (optionally after dropping a trailing + * stop codon). This method does not check whether the CDS nucleotide sequence + * translates to the peptide sequence. * * @param dnaSeq - * @param peptide - * @param dnaToProtein + * @param proteinSeq + * @return */ - public static int computeProteinFeatures(SequenceI dnaSeq, - SequenceI peptide, MapList dnaToProtein) + public static MapList mapCdsToProtein(SequenceI dnaSeq, + SequenceI proteinSeq) { - while (dnaSeq.getDatasetSequence() != null) - { - dnaSeq = dnaSeq.getDatasetSequence(); - } - while (peptide.getDatasetSequence() != null) - { - peptide = peptide.getDatasetSequence(); - } - - transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON); + List ranges = findCdsPositions(dnaSeq); + int mappedDnaLength = MappingUtils.getLength(ranges); /* - * 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 + * if not a whole number of codons, truncate mapping */ + int codonRemainder = mappedDnaLength % CODON_LENGTH; + if (codonRemainder > 0) + { + mappedDnaLength -= codonRemainder; + MappingUtils.removeEndPositions(codonRemainder, ranges); + } - /* - * build a map with codon variations for each potentially varying peptide - */ - LinkedHashMap[]> variants = buildDnaVariantsMap( - dnaSeq, dnaToProtein); + int proteinLength = proteinSeq.getLength(); + int proteinStart = proteinSeq.getStart(); + int proteinEnd = proteinSeq.getEnd(); /* - * scan codon variations, compute peptide variants and add to peptide sequence + * incomplete start codon may mean X at start of peptide + * we ignore both for mapping purposes */ - int count = 0; - for (Entry[]> variant : variants.entrySet()) + if (proteinSeq.getCharAt(0) == 'X') { - int peptidePos = variant.getKey(); - List[] codonVariants = variant.getValue(); - count += computePeptideVariants(peptide, peptidePos, codonVariants); + // todo JAL-2022 support startPhase > 0 + proteinStart++; + proteinLength--; } + List proteinRange = new ArrayList<>(); /* - * sort to get sequence features in start position order - * - would be better to store in Sequence as a TreeSet or NCList? + * dna length should map to protein (or protein plus stop codon) */ - if (peptide.getSequenceFeatures() != null) + int codesForResidues = mappedDnaLength / CODON_LENGTH; + if (codesForResidues == (proteinLength + 1)) { - 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; - } - }); + // assuming extra codon is for STOP and not in peptide + // todo: check trailing codon is indeed a STOP codon + codesForResidues--; + mappedDnaLength -= CODON_LENGTH; + MappingUtils.removeEndPositions(CODON_LENGTH, ranges); } - return count; + + if (codesForResidues == proteinLength) + { + proteinRange.add(new int[] { proteinStart, proteinEnd }); + return new MapList(ranges, proteinRange, CODON_LENGTH, 1); + } + return null; } /** - * 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. + * 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 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 + * @param dnaSeq + * @return */ - static int computePeptideVariants(SequenceI peptide, int peptidePos, - List[] codonVariants) + protected static List findCdsPositions(SequenceI dnaSeq) { - 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; + List result = new ArrayList<>(); - /* - * variants in first codon base - */ - for (DnaVariant var : codonVariants[0]) + List sfs = dnaSeq.getFeatures().getFeaturesByOntology( + SequenceOntologyI.CDS); + if (sfs.isEmpty()) { - 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++; - } - } - } - } + return result; } + SequenceFeatures.sortFeatures(sfs, true); - /* - * variants in second codon base - */ - for (DnaVariant var : codonVariants[1]) + for (SequenceFeature sf : sfs) { - if (var.variant != null) + int phase = 0; + try { - 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++; - } - } - } + phase = Integer.parseInt(sf.getPhase()); + } catch (NumberFormatException e) + { + // ignore } - } - - /* - * variants in third codon base - */ - for (DnaVariant var : codonVariants[2]) - { - if (var.variant != null) + /* + * 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() && phase > 0) { - String alleles = (String) var.variant.getValue("alleles"); - if (alleles != null) + begin += phase; + if (begin > end) { - for (String base : alleles.split(",")) - { - String codon = base1 + base2 + base; - if (addPeptideVariant(peptide, peptidePos, residue, var, codon)) - { - count++; - } - } + // shouldn't happen! + System.err + .println("Error: start phase extends beyond start CDS in " + + dnaSeq.getName()); } } + result.add(new int[] { begin, end }); } - return count; + /* + * 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, IntRangeComparator.ASCENDING); + return result; } /** - * 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. + * Helper method that adds a peptide variant feature. 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 + * the variant codon e.g. aCg + * @param canonical + * the 'normal' codon e.g. aTg * @return true if a feature was added, else false */ static boolean addPeptideVariant(SequenceI peptide, int peptidePos, - String residue, DnaVariant var, String codon) + String residue, DnaVariant var, String codon, String canonical) { /* * get peptide translation of codon e.g. GAT -> D @@ -2004,181 +2412,71 @@ public class AlignmentUtils * 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 trans = codon.contains("-") ? null + : (codon.length() > CODON_LENGTH ? null + : ResidueProperties.codonTranslate(codon)); + if (trans == null) + { + return false; + } + String desc = canonical + "/" + codon; + String featureType = ""; + if (trans.equals(residue)) + { + featureType = SequenceOntologyI.SYNONYMOUS_VARIANT; + } + else if (ResidueProperties.STOP.equals(trans)) + { + featureType = SequenceOntologyI.STOP_GAINED; + } + else { 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) + desc = "p." + residue3Char + peptidePos + trans3Char; + featureType = SequenceOntologyI.NONSYNONYMOUS_VARIANT; + } + SequenceFeature sf = new SequenceFeature(featureType, desc, peptidePos, + peptidePos, var.getSource()); + + StringBuilder attributes = new StringBuilder(32); + String id = (String) var.variant.getValue(VARIANT_ID); + if (id != null) + { + if (id.startsWith(SEQUENCE_VARIANT)) { - sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); - attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=") - .append(clinSig); + id = id.substring(SEQUENCE_VARIANT.length()); } - peptide.addSequenceFeature(sf); - if (attributes.length() > 0) + sf.setValue(VARIANT_ID, id); + attributes.append(VARIANT_ID).append("=").append(id); + // TODO handle other species variants JAL-2064 + 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) { - sf.setAttributes(attributes.toString()); + // as if } - 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[]> 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) + String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE); + if (clinSig != null) { - return variants; + sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); + attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=") + .append(clinSig); } - - int dnaStart = dnaSeq.getStart(); - int[] lastCodon = null; - int lastPeptidePostion = 0; - - /* - * build a map of codon variations for peptides - */ - for (SequenceFeature sf : dnaFeatures) + peptide.addSequenceFeature(sf); + if (attributes.length() > 0) { - 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)); - } - } - } + sf.setAttributes(attributes.toString()); } - return variants; + return true; } /** @@ -2188,33 +2486,16 @@ public class AlignmentUtils * * @param seqs * @param xrefs + * @param dataset + * the alignment dataset shared by the new copy * @return */ public static AlignmentI makeCopyAlignment(SequenceI[] seqs, - SequenceI[] xrefs) + SequenceI[] xrefs, AlignmentI dataset) { AlignmentI copy = new Alignment(new Alignment(seqs)); - - /* - * add mappings between sequences to the new alignment - */ - AlignedCodonFrame mappings = new AlignedCodonFrame(); - copy.addCodonFrame(mappings); - for (int i = 0; i < copy.getHeight(); i++) - { - SequenceI from = seqs[i]; - SequenceI to = copy.getSequenceAt(i); - if (to.getDatasetSequence() != null) - { - to = to.getDatasetSequence(); - } - int start = from.getStart(); - int end = from.getEnd(); - MapList map = new MapList(new int[] { start, end }, new int[] { - start, end }, 1, 1); - mappings.addMap(to, from, map); - } - + copy.setDataset(dataset); + boolean isProtein = !copy.isNucleotide(); SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) { @@ -2225,7 +2506,8 @@ public class AlignmentUtils { for (DBRefEntry dbref : dbrefs) { - if (dbref.getMap() == null || dbref.getMap().getTo() == null) + if (dbref.getMap() == null || dbref.getMap().getTo() == null + || dbref.getMap().getTo().isProtein() != isProtein) { continue; } @@ -2259,19 +2541,32 @@ public class AlignmentUtils */ public static int alignAs(AlignmentI unaligned, AlignmentI aligned) { - List unmapped = new ArrayList(); + /* + * 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; + // 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; @@ -2293,7 +2588,7 @@ public class AlignmentUtils } newCol++; } - + /* * trim trailing gaps */ @@ -2303,6 +2598,7 @@ public class AlignmentUtils System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1); newSeq = tmp; } + // TODO: optimise SequenceI to avoid char[]->String->char[] seq.setSequence(String.valueOf(newSeq)); realignedCount++; } @@ -2311,6 +2607,97 @@ public class AlignmentUtils } /** + * 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> alignedDatasets = new HashMap<>(); + for (SequenceI seq : aligned.getSequences()) + { + SequenceI ds = seq.getDatasetSequence(); + if (alignedDatasets.get(ds) == null) + { + alignedDatasets.put(ds, new ArrayList()); + } + alignedDatasets.get(ds).add(seq); + } + + /* + * first pass - check whether all sequences to be aligned share a + * dataset sequence with an aligned sequence; also note the leftmost + * ungapped column from which to copy + */ + int leftmost = Integer.MAX_VALUE; + for (SequenceI seq : unaligned.getSequences()) + { + final SequenceI ds = seq.getDatasetSequence(); + if (!alignedDatasets.containsKey(ds)) + { + return false; + } + SequenceI alignedSeq = alignedDatasets.get(ds) + .get(0); + int startCol = alignedSeq.findIndex(seq.getStart()); // 1.. + leftmost = Math.min(leftmost, startCol); + } + + /* + * second pass - copy aligned sequences; + * heuristic rule: pair off sequences in order for the case where + * more than one shares the same dataset sequence + */ + final char gapCharacter = aligned.getGapCharacter(); + for (SequenceI seq : unaligned.getSequences()) + { + List alignedSequences = alignedDatasets + .get(seq.getDatasetSequence()); + SequenceI alignedSeq = alignedSequences.get(0); + + /* + * gap fill for leading (5') UTR if any + */ + // TODO this copies intron columns - wrong! + int startCol = alignedSeq.findIndex(seq.getStart()); // 1.. + int endCol = alignedSeq.findIndex(seq.getEnd()); + char[] seqchars = new char[endCol - leftmost + 1]; + Arrays.fill(seqchars, gapCharacter); + char[] toCopy = alignedSeq.getSequence(startCol - 1, endCol); + System.arraycopy(toCopy, 0, seqchars, startCol - leftmost, + toCopy.length); + seq.setSequence(String.valueOf(seqchars)); + if (alignedSequences.size() > 0) + { + // pop off aligned sequences (except the last one) + alignedSequences.remove(0); + } + } + + /* + * finally remove gapped columns (e.g. introns) + */ + new RemoveGapColCommand("", unaligned.getSequencesArray(), 0, + unaligned.getWidth() - 1, unaligned); + + 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. * @@ -2319,18 +2706,19 @@ public class AlignmentUtils * @param unmapped * @return */ - static Map> buildMappedColumnsMap( - AlignmentI unaligned, AlignmentI aligned, List unmapped) + static SortedMap> buildMappedColumnsMap( + AlignmentI unaligned, AlignmentI aligned, + List unmapped) { /* * 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> map = new TreeMap>(); + SortedMap> map = new TreeMap<>(); /* - * 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()); @@ -2355,7 +2743,8 @@ public class AlignmentUtils } /** - * Helper method that adds to a map the mapped column positions of a sequence.
    + * 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. @@ -2379,9 +2768,16 @@ public class AlignmentUtils return false; } - char[] fromChars = fromSeq.getSequence(); + /* + * invert mapping if it is from unaligned to aligned sequence + */ + if (seqMap.getTo() == fromSeq.getDatasetSequence()) + { + seqMap = new Mapping(seq.getDatasetSequence(), + seqMap.getMap().getInverse()); + } + int toStart = seq.getStart(); - char[] toChars = seq.getSequence(); /* * traverse [start, end, start, end...] ranges in fromSeq @@ -2412,9 +2808,10 @@ public class AlignmentUtils * 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 <= fromSeq.getLength() + && fromCol >= 0) { - if (!Comparison.isGap(fromChars[fromCol - 1])) + if (!Comparison.isGap(fromSeq.getCharAt(fromCol - 1))) { /* * mapped from sequence has a character in this column @@ -2423,10 +2820,10 @@ public class AlignmentUtils Map seqsMap = map.get(fromCol); if (seqsMap == null) { - seqsMap = new HashMap(); + seqsMap = new HashMap<>(); map.put(fromCol, seqsMap); } - seqsMap.put(seq, toChars[mappedCharPos - toStart]); + seqsMap.put(seq, seq.getCharAt(mappedCharPos - toStart)); mappedCharPos++; } fromCol += (forward ? 1 : -1);