X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=0c40873154035f5cbb2b1a6d264709add75810df;hb=006890b02106eb31841e6e84d75f1027434823e0;hp=ea96b3b29efeca02a5e236062ffc93a10e9daad7;hpb=9390152fd31fd0d9c52833f1e5c63d222cb38d82;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index ea96b3b..0c40873 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,8 +20,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,28 +28,26 @@ 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.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; -import java.io.UnsupportedEncodingException; -import java.net.URLEncoder; 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; @@ -72,12 +69,15 @@ 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 @@ -105,6 +105,15 @@ public class AlignmentUtils { return variant == null ? null : variant.getFeatureGroup(); } + + /** + * toString for aid in the debugger only + */ + @Override + public String toString() + { + return base + ":" + (variant == null ? "" : variant.getDescription()); + } } /** @@ -117,7 +126,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()) { @@ -170,10 +179,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]; @@ -188,8 +199,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); @@ -245,7 +256,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(); @@ -254,7 +265,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); @@ -281,8 +292,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 @@ -316,9 +327,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(); @@ -347,9 +358,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; } @@ -383,7 +393,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) @@ -402,7 +412,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
    • @@ -424,8 +435,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(); @@ -451,7 +463,7 @@ public class AlignmentUtils { String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(); - for (String stop : ResidueProperties.STOP) + for (String stop : ResidueProperties.STOP_CODONS) { if (lastCodon.equals(stop)) { @@ -466,8 +478,7 @@ public class AlignmentUtils * If lengths still don't match, try ignoring start codon. */ int startOffset = 0; - if (cdnaLength != mappedLength - && cdnaLength > 2 + if (cdnaLength != mappedLength && cdnaLength > 2 && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase() .equals(ResidueProperties.START)) { @@ -481,8 +492,9 @@ public class AlignmentUtils /* * protein is translation of dna (+/- start/stop codons) */ - MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] - { proteinStart, proteinEnd }, CODON_LENGTH, 1); + MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, + new int[] + { proteinStart, proteinEnd }, CODON_LENGTH, 1); return map; } @@ -512,7 +524,8 @@ public class AlignmentUtils int aaPos = 0; int dnaPos = cdnaStart; - for (; dnaPos < cdnaSeqChars.length - 2 && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++) + for (; dnaPos < cdnaSeqChars.length - 2 + && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++) { String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); final String translated = ResidueProperties.codonTranslate(codon); @@ -521,7 +534,8 @@ public class AlignmentUtils * 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; } @@ -553,7 +567,8 @@ public class AlignmentUtils if (dnaPos == cdnaSeqChars.length - CODON_LENGTH) { String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH); - if ("STOP".equals(ResidueProperties.codonTranslate(codon))) + if (ResidueProperties.STOP + .equals(ResidueProperties.codonTranslate(codon))) { return true; } @@ -631,10 +646,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 @@ -650,15 +664,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; @@ -698,9 +713,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++; @@ -726,7 +741,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); } @@ -754,9 +769,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); @@ -829,8 +844,9 @@ public class AlignmentUtils } else { - gapsToAdd = Math.min(intronLength + trailingGapLength - - sourceGapMappedLength, trailingGapLength); + gapsToAdd = Math.min( + intronLength + trailingGapLength - sourceGapMappedLength, + trailingGapLength); } } } @@ -865,7 +881,7 @@ public class AlignmentUtils System.err.println("Wrong alignment type in alignProteinAsDna"); return 0; } - List unmappedProtein = new ArrayList(); + List unmappedProtein = new ArrayList<>(); Map> alignedCodons = buildCodonColumnsMap( protein, dna, unmappedProtein); return alignProteinAs(protein, alignedCodons, unmappedProtein); @@ -929,7 +945,8 @@ public class AlignmentUtils * @return */ static boolean alignCdsSequenceAsProtein(SequenceI cdsSeq, - AlignmentI protein, List mappings, char gapChar) + AlignmentI protein, List mappings, + char gapChar) { SequenceI cdsDss = cdsSeq.getDatasetSequence(); if (cdsDss == null) @@ -946,7 +963,7 @@ public class AlignmentUtils SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); if (peptide != null) { - int peptideLength = peptide.getLength(); + final int peptideLength = peptide.getLength(); Mapping map = mapping.getMappingBetween(cdsSeq, peptide); if (map != null) { @@ -955,21 +972,20 @@ public class AlignmentUtils { mapList = mapList.getInverse(); } - int cdsLength = cdsDss.getLength(); + 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); + || (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())); + System.err.println(String.format( + "Can't align cds as protein (length mismatch %d/%d): %s", + cdsLength, mappedToLength, cdsSeq.getName())); } /* @@ -983,14 +999,15 @@ public class AlignmentUtils * 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) + + for (int col = 0; col < peptideLength; col++) { + char residue = peptide.getCharAt(col); + if (Comparison.isGap(residue)) { cdsCol += CODON_LENGTH; @@ -1008,7 +1025,7 @@ public class AlignmentUtils { for (int j = codon[0]; j <= codon[1]; j++) { - char mappedBase = nucleotides[j - cdsStart]; + char mappedBase = cdsDss.getCharAt(j - cdsStart); alignedCds[cdsCol++] = mappedBase; copiedBases++; } @@ -1020,7 +1037,7 @@ public class AlignmentUtils * append stop codon if not mapped from protein, * closing it up to the end of the mapped sequence */ - if (copiedBases == nucleotides.length - CODON_LENGTH) + if (copiedBases == cdsLength - CODON_LENGTH) { for (int i = alignedCds.length - 1; i >= 0; i--) { @@ -1030,9 +1047,9 @@ public class AlignmentUtils break; } } - for (int i = nucleotides.length - CODON_LENGTH; i < nucleotides.length; i++) + for (int i = cdsLength - CODON_LENGTH; i < cdsLength; i++) { - alignedCds[cdsCol++] = nucleotides[i]; + alignedCds[cdsCol++] = cdsDss.getCharAt(i); } } cdsSeq.setSequence(new String(alignedCds)); @@ -1075,7 +1092,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()) @@ -1086,8 +1103,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); } } @@ -1121,9 +1138,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()) @@ -1140,8 +1157,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) @@ -1152,8 +1169,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 @@ -1202,21 +1219,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()) { @@ -1224,12 +1246,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; } @@ -1289,7 +1319,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); @@ -1302,7 +1332,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
      • *
    @@ -1358,9 +1389,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) @@ -1397,8 +1429,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) { @@ -1425,7 +1456,7 @@ public class AlignmentUtils { continue; } - final List result = new ArrayList(); + final List result = new ArrayList<>(); for (AlignmentAnnotation dsann : datasetAnnotations) { /* @@ -1502,8 +1533,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 @@ -1527,9 +1558,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; } @@ -1567,11 +1597,12 @@ public class AlignmentUtils return false; } String name = seq2.getName(); - final DBRefEntry[] xrefs = seq1.getDBRefs(); + final List xrefs = seq1.getDBRefs(); if (xrefs != null) { - for (DBRefEntry xref : xrefs) + for (int ix = 0, nx = xrefs.size(); ix < nx; ix++) { + DBRefEntry xref = xrefs.get(ix); String xrefName = xref.getSource() + "|" + xref.getAccessionId(); // case-insensitive test, consistent with DBRefEntry.equalRef() if (xrefName.equalsIgnoreCase(name)) @@ -1608,13 +1639,13 @@ public class AlignmentUtils throw new IllegalArgumentException( "IMPLEMENTATION ERROR: dataset.getDataset() must be null!"); } - List foundSeqs = new ArrayList(); - 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(); + productSeqs = new HashSet<>(); for (SequenceI seq : products) { productSeqs.add(seq.getDatasetSequence() == null ? seq : seq @@ -1701,31 +1732,36 @@ public class AlignmentUtils 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 + * build the mapping from CDS to protein */ - List cdsRange = Collections.singletonList(new int[] { 1, - cdsSeq.getLength() }); + 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(cdsSeqDss, proteinProduct, - cdsToProteinMap); - /* - * guard against duplicating the mapping if repeating this action - */ - if (!mappings.contains(cdsToProteinMapping)) + if (!dataset.getSequences().contains(cdsSeqDss)) { - mappings.add(cdsToProteinMapping); + /* + * if this sequence is a newly created one, add it to the dataset + * and made a CDS to protein mapping (if sequence already exists, + * CDS-to-protein mapping _is_ the transcript-to-protein mapping) + */ + dataset.addSequence(cdsSeqDss); + AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); + cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct, + cdsToProteinMap); + + /* + * guard against duplicating the mapping if repeating this action + */ + if (!mappings.contains(cdsToProteinMapping)) + { + mappings.add(cdsToProteinMapping); + } } propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(), @@ -1734,7 +1770,7 @@ public class AlignmentUtils * add another mapping from original 'from' range to CDS */ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); - MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), + final MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), cdsRange, 1, 1); dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, dnaToCdsMap); @@ -1744,6 +1780,13 @@ public class AlignmentUtils } /* + * 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 @@ -1760,32 +1803,35 @@ public class AlignmentUtils // need to // synthesize an xref. - for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs()) + List primrefs = dnaDss.getPrimaryDBRefs(); + for (int ip = 0, np = primrefs.size(); ip < np; ip++) { - // creates a complementary cross-reference to the source sequence's - // primary reference. - - DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(), - primRef.getSource() + ":" + primRef.getVersion(), - primRef.getAccessionId()); - cdsCrossRef - .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap))); + DBRefEntry primRef = primrefs.get(ip); + /* + * 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( - primRef.getSource(), primRef.getVersion(), + 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 */ @@ -1803,6 +1849,38 @@ public class AlignmentUtils } /** + * Tries to transfer gene loci (dbref to chromosome positions) from fromSeq to + * toSeq, mediated by the given mapping between the sequences + * + * @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. @@ -1814,7 +1892,7 @@ public class AlignmentUtils * @param seqMappings * the set of mappings involving dnaSeq * @param aMapping - * an initial candidate from seqMappings + * a transcript-to-peptide mapping * @return */ static SequenceI findCdsForProtein(List mappings, @@ -1833,13 +1911,21 @@ public class AlignmentUtils * 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 mappedFromLength = MappingUtils + .getLength(aMapping.getMap().getFromRanges()); int dnaLength = seqDss.getLength(); if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - CODON_LENGTH) { - return seqDss; + /* + * 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()) + { + return seqDss; + } } /* @@ -1858,16 +1944,18 @@ public class AlignmentUtils && proteinProduct == mapping.getTo() && seqDss != map.getFromSeq()) { - mappedFromLength = MappingUtils.getLength(mapping.getMap() - .getFromRanges()); + mappedFromLength = MappingUtils + .getLength(mapping.getMap().getFromRanges()); if (mappedFromLength == map.getFromSeq().getLength()) { /* * found a 3:1 mapping to the protein product which covers - * the whole dna sequence i.e. is from CDS; finally check it - * is from the dna start sequence + * 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()) @@ -1897,39 +1985,61 @@ public class AlignmentUtils static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, AlignmentI dataset) { - char[] seqChars = seq.getSequence(); - List fromRanges = mapping.getMap().getFromRanges(); - int cdsWidth = MappingUtils.getLength(fromRanges); - char[] newSeqChars = new char[cdsWidth]; + /* + * 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()); - int newPos = 0; - for (int[] range : fromRanges) + SequenceI newSeq = null; + + final MapList maplist = mapping.getMap(); + if (maplist.isContiguous() && maplist.isFromForwardStrand()) { - if (range[0] <= range[1]) - { - // forward strand mapping - just copy the range - int length = range[1] - range[0] + 1; - System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos, - length); - newPos += length; - } - else + /* + * just a subsequence, keep same dataset sequence + */ + int start = maplist.getFromLowest(); + int end = maplist.getFromHighest(); + newSeq = seq.getSubSequence(start - 1, end); + newSeq.setName(seqId); + } + else + { + /* + * 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) { - // reverse strand mapping - copy and complement one by one - for (int i = range[0]; i >= range[1]; i--) + if (range[0] <= range[1]) + { + // forward strand mapping - just copy the range + int length = range[1] - range[0] + 1; + System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos, + length); + newPos += length; + } + else { - newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); + // reverse strand mapping - copy and complement one by one + for (int i = range[0]; i >= range[1]; i--) + { + newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); + } } } + + newSeq = new Sequence(seqId, newSeqChars, 1, newPos); } - /* - * assign 'from id' held in the mapping if set (e.g. EMBL protein_id), - * else generate a sequence name - */ - String mapFromId = mapping.getMappedFromId(); - String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName()); - SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos); if (dataset != null) { SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName()); @@ -1957,8 +2067,8 @@ public class AlignmentUtils } else { - System.err - .println("JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + System.err.println( + "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + mtch.toString()); } } @@ -1970,28 +2080,32 @@ public class AlignmentUtils } /** - * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to + * Adds any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to * the given mapping. * * @param cdsSeq * @param contig + * @param proteinProduct * @param mapping - * @return list of DBRefEntrys added. + * @return list of DBRefEntrys added */ - public static List propagateDBRefsToCDS(SequenceI cdsSeq, + protected static List propagateDBRefsToCDS(SequenceI cdsSeq, SequenceI contig, SequenceI proteinProduct, Mapping mapping) { - // gather direct refs from contig congrent with mapping - List direct = new ArrayList(); - HashSet directSources = new HashSet(); - if (contig.getDBRefs() != null) + // gather direct refs from contig congruent with mapping + List direct = new ArrayList<>(); + HashSet directSources = new HashSet<>(); + + List refs = contig.getDBRefs(); + if (refs != null) { - for (DBRefEntry dbr : contig.getDBRefs()) + for (int ib = 0, nb = refs.size(); ib < nb; ib++) { - if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap()) + DBRefEntry dbr = refs.get(ib); + MapList map; + if (dbr.hasMap() && (map = dbr.getMap().getMap()).isTripletMap()) { - MapList map = dbr.getMap().getMap(); // check if map is the CDS mapping if (mapping.getMap().equals(map)) { @@ -2001,25 +2115,27 @@ public class AlignmentUtils } } } - DBRefEntry[] onSource = DBRefUtils.selectRefs( + List onSource = DBRefUtils.selectRefs( proteinProduct.getDBRefs(), directSources.toArray(new String[0])); - List propagated = new ArrayList(); + List propagated = new ArrayList<>(); // and generate appropriate mappings - for (DBRefEntry cdsref : direct) + for (int ic = 0, nc = direct.size(); ic < nc; ic++) { + DBRefEntry cdsref = direct.get(ic); + Mapping m = cdsref.getMap(); // 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()); + MapList cdsposmap = new MapList( + Arrays.asList(new int[][] + { new int[] { cdsSeq.getStart(), cdsSeq.getEnd() } }), + m.getMap().getToRanges(), 3, 1); + Mapping cdsmap = new Mapping(m.getTo(), m.getMap()); // create dbref DBRefEntry newref = new DBRefEntry(cdsref.getSource(), - cdsref.getVersion(), cdsref.getAccessionId(), new Mapping( - cdsmap.getTo(), cdsposmap)); + 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 @@ -2054,14 +2170,14 @@ public class AlignmentUtils * * @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 mapping - * the mapping from 'fromSeq' to 'toSeq' * @param omitting */ - public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, + protected static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, MapList mapping, String select, String... omitting) { SequenceI copyTo = toSeq; @@ -2069,83 +2185,89 @@ public class AlignmentUtils { copyTo = copyTo.getDatasetSequence(); } + if (fromSeq == copyTo || fromSeq.getDatasetSequence() == copyTo) + { + return 0; // shared dataset sequence + } + + /* + * get features, optionally restricted by an ontology term + */ + List sfs = select == null ? fromSeq.getFeatures() + .getPositionalFeatures() : fromSeq.getFeatures() + .getFeaturesByOntology(select); - SequenceOntologyI so = SequenceOntologyFactory.getInstance(); int count = 0; - SequenceFeature[] sfs = fromSeq.getSequenceFeatures(); - if (sfs != null) + for (SequenceFeature sf : sfs) { - for (SequenceFeature sf : sfs) + String type = sf.getType(); + boolean omit = false; + for (String toOmit : omitting) { - String type = sf.getType(); - if (select != null && !so.isA(type, select)) + if (type.equals(toOmit)) { - continue; - } - boolean omit = false; - for (String toOmit : omitting) - { - if (type.equals(toOmit)) - { - omit = true; - } - } - if (omit) - { - continue; + omit = true; } + } + 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) - { - mappedTo = mapping.locateInTo(end, end); - if (mappedTo != null) - { - /* - * end of exon is in CDS range - 5' overlap - * to a range from the start of the peptide - */ - mappedTo[0] = 1; - } - } - if (mappedTo == null) + /* + * 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) + { + mappedTo = mapping.locateInTo(end, end); + if (mappedTo != null) { - mappedTo = mapping.locateInTo(start, start); - if (mappedTo != null) - { - /* - * start of exon is in CDS range - 3' overlap - * to a range up to the end of the peptide - */ - mappedTo[1] = toSeq.getLength(); - } + /* + * end of exon is in CDS range - 5' overlap + * to a range from the start of the peptide + */ + mappedTo[0] = 1; } + } + if (mappedTo == null) + { + mappedTo = mapping.locateInTo(start, start); if (mappedTo != null) { - 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++; + /* + * start of exon is in CDS range - 3' overlap + * to a range up to the end of the peptide + */ + mappedTo[1] = toSeq.getLength(); } } + if (mappedTo != null) + { + 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 count; } /** * Returns a mapping from dna to protein by inspecting sequence features of - * type "CDS" on the dna. + * 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 proteinSeq @@ -2157,6 +2279,16 @@ public class AlignmentUtils List ranges = findCdsPositions(dnaSeq); int mappedDnaLength = MappingUtils.getLength(ranges); + /* + * if not a whole number of codons, truncate mapping + */ + int codonRemainder = mappedDnaLength % CODON_LENGTH; + if (codonRemainder > 0) + { + mappedDnaLength -= codonRemainder; + MappingUtils.removeEndPositions(codonRemainder, ranges); + } + int proteinLength = proteinSeq.getLength(); int proteinStart = proteinSeq.getStart(); int proteinEnd = proteinSeq.getEnd(); @@ -2171,7 +2303,7 @@ public class AlignmentUtils proteinStart++; proteinLength--; } - List proteinRange = new ArrayList(); + List proteinRange = new ArrayList<>(); /* * dna length should map to protein (or protein plus stop codon) @@ -2180,8 +2312,12 @@ public class AlignmentUtils if (codesForResidues == (proteinLength + 1)) { // 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); } + if (codesForResidues == proteinLength) { proteinRange.add(new int[] { proteinStart, proteinEnd }); @@ -2192,7 +2328,7 @@ public class AlignmentUtils /** * 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 + * [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. @@ -2200,62 +2336,50 @@ public class AlignmentUtils * @param dnaSeq * @return */ - public static List findCdsPositions(SequenceI dnaSeq) + protected static List findCdsPositions(SequenceI dnaSeq) { - List result = new ArrayList(); - SequenceFeature[] sfs = dnaSeq.getSequenceFeatures(); - if (sfs == null) + List result = new ArrayList<>(); + + List sfs = dnaSeq.getFeatures().getFeaturesByOntology( + SequenceOntologyI.CDS); + if (sfs.isEmpty()) { return result; } - - SequenceOntologyI so = SequenceOntologyFactory.getInstance(); - int startPhase = 0; + SequenceFeatures.sortFeatures(sfs, true); for (SequenceFeature sf : sfs) { + int phase = 0; + try + { + String s = sf.getPhase(); + if (s != null) + { + phase = Integer.parseInt(s); + } + } catch (NumberFormatException e) + { + // leave as zero + } /* - * process a CDS feature (or a sub-type of CDS) + * phase > 0 on first codon means 5' incomplete - skip to the start + * of the next codon; example ENST00000496384 */ - if (so.isA(sf.getType(), SequenceOntologyI.CDS)) + int begin = sf.getBegin(); + int end = sf.getEnd(); + if (result.isEmpty() && phase > 0) { - int phase = 0; - try - { - phase = Integer.parseInt(sf.getPhase()); - } catch (NumberFormatException e) + begin += phase; + if (begin > end) { - // ignore - } - /* - * phase > 0 on first codon means 5' incomplete - skip to the start - * of the next codon; example ENST00000496384 - */ - int begin = sf.getBegin(); - int end = sf.getEnd(); - if (result.isEmpty()) - { - begin += phase; - if (begin > end) - { - // shouldn't happen! - System.err - .println("Error: start phase extends beyond start CDS in " - + dnaSeq.getName()); - } + // 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; + result.add(new int[] { begin, end }); } /* @@ -2265,381 +2389,11 @@ public class AlignmentUtils * ranges are assembled in order. Other cases should not use this method, * but instead construct an explicit mapping for CDS (e.g. EMBL parsing). */ - Collections.sort(result, new Comparator() - { - @Override - public int compare(int[] o1, int[] o2) - { - return Integer.compare(o1[0], o2[0]); - } - }); + Collections.sort(result, IntRangeComparator.ASCENDING); return result; } /** - * Maps exon features from dna to protein, and computes variants in peptide - * product generated by variants in dna, and adds them as sequence_variant - * features on the protein sequence. Returns the number of variant features - * added. - * - * @param dnaSeq - * @param peptide - * @param dnaToProtein - */ - public static int computeProteinFeatures(SequenceI dnaSeq, - SequenceI peptide, MapList dnaToProtein) - { - while (dnaSeq.getDatasetSequence() != null) - { - dnaSeq = dnaSeq.getDatasetSequence(); - } - while (peptide.getDatasetSequence() != null) - { - peptide = peptide.getDatasetSequence(); - } - - transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON); - - /* - * compute protein variants from dna variants and codon mappings; - * NB - alternatively we could retrieve this using the REST service e.g. - * http://rest.ensembl.org/overlap/translation - * /ENSP00000288602?feature=transcript_variation;content-type=text/xml - * which would be a bit slower but possibly more reliable - */ - - /* - * build a map with codon variations for each potentially varying peptide - */ - LinkedHashMap[]> variants = buildDnaVariantsMap( - dnaSeq, dnaToProtein); - - /* - * scan codon variations, compute peptide variants and add to peptide sequence - */ - int count = 0; - for (Entry[]> variant : variants.entrySet()) - { - int peptidePos = variant.getKey(); - List[] codonVariants = variant.getValue(); - count += computePeptideVariants(peptide, peptidePos, codonVariants); - } - - /* - * sort to get sequence features in start position order - * - would be better to store in Sequence as a TreeSet or NCList? - */ - if (peptide.getSequenceFeatures() != null) - { - Arrays.sort(peptide.getSequenceFeatures(), - new Comparator() - { - @Override - public int compare(SequenceFeature o1, SequenceFeature o2) - { - int c = Integer.compare(o1.getBegin(), o2.getBegin()); - return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd()) - : c; - } - }); - } - return count; - } - - /** - * Computes non-synonymous peptide variants from codon variants and adds them - * as sequence_variant features on the protein sequence (one feature per - * allele variant). Selected attributes (variant id, clinical significance) - * are copied over to the new features. - * - * @param peptide - * the protein sequence - * @param peptidePos - * the position to compute peptide variants for - * @param codonVariants - * a list of dna variants per codon position - * @return the number of features added - */ - static int computePeptideVariants(SequenceI peptide, int peptidePos, - List[] codonVariants) - { - 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() > CODON_LENGTH ? null : ResidueProperties - .codonTranslate(codon)); - if (trans != null && !trans.equals(residue)) - { - 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, var.getSource()); - 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 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) - { - // 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 - */ - @SuppressWarnings("unchecked") - 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) - { - return variants; - } - - int dnaStart = dnaSeq.getStart(); - int[] lastCodon = null; - int lastPeptidePostion = 0; - - /* - * build a map of codon variations for peptides - */ - for (SequenceFeature sf : dnaFeatures) - { - int dnaCol = sf.getBegin(); - if (dnaCol != sf.getEnd()) - { - // not handling multi-locus variant features - continue; - } - if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT)) - { - int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol); - if (mapsTo == null) - { - // feature doesn't lie within coding region - continue; - } - int peptidePosition = mapsTo[0]; - List[] codonVariants = variants.get(peptidePosition); - if (codonVariants == null) - { - codonVariants = new ArrayList[CODON_LENGTH]; - 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 < CODON_LENGTH; codonPos++) - { - String nucleotide = String.valueOf( - dnaSeq.getCharAt(codon[codonPos] - dnaStart)) - .toUpperCase(); - List codonVariant = codonVariants[codonPos]; - if (codon[codonPos] == dnaCol) - { - if (!codonVariant.isEmpty() - && codonVariant.get(0).variant == null) - { - /* - * already recorded base value, add this variant - */ - codonVariant.get(0).variant = sf; - } - else - { - /* - * add variant with base value - */ - codonVariant.add(new DnaVariant(nucleotide, sf)); - } - } - else if (codonVariant.isEmpty()) - { - /* - * record (possibly non-varying) base value - */ - codonVariant.add(new DnaVariant(nucleotide)); - } - } - } - } - return variants; - } - - /** * Makes an alignment with a copy of the given sequences, adding in any * non-redundant sequences which are mapped to by the cross-referenced * sequences. @@ -2659,19 +2413,25 @@ public class AlignmentUtils SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) { - for (SequenceI xref : xrefs) + // BH 2019.01.25 recoded to remove iterators + + for (int ix = 0, nx = xrefs.length; ix < nx; ix++) { - DBRefEntry[] dbrefs = xref.getDBRefs(); + SequenceI xref = xrefs[ix]; + List dbrefs = xref.getDBRefs(); if (dbrefs != null) { - for (DBRefEntry dbref : dbrefs) + for (int ir = 0, nir = dbrefs.size(); ir < nir; ir++) { - if (dbref.getMap() == null || dbref.getMap().getTo() == null - || dbref.getMap().getTo().isProtein() != isProtein) + DBRefEntry dbref = dbrefs.get(ir); + Mapping map = dbref.getMap(); + SequenceI mto; + if (map == null || (mto = map.getTo()) == null + || mto.isProtein() != isProtein) { continue; } - SequenceI mappedTo = dbref.getMap().getTo(); + SequenceI mappedTo = mto; SequenceI match = matcher.findIdMatch(mappedTo); if (match == null) { @@ -2712,7 +2472,7 @@ public class AlignmentUtils /* * fancy case - aligning via mappings between sequences */ - List unmapped = new ArrayList(); + List unmapped = new ArrayList<>(); Map> columnMap = buildMappedColumnsMap( unaligned, aligned, unmapped); int width = columnMap.size(); @@ -2772,10 +2532,10 @@ public class AlignmentUtils * true; else returns false * * @param unaligned - * - sequences to be aligned based on aligned + * - sequences to be aligned based on aligned * @param aligned - * - 'guide' alignment containing sequences derived from same dataset - * as unaligned + * - 'guide' alignment containing sequences derived from same + * dataset as unaligned * @return */ static boolean alignAsSameSequences(AlignmentI unaligned, @@ -2787,7 +2547,7 @@ public class AlignmentUtils } // map from dataset sequence to alignment sequence(s) - Map> alignedDatasets = new HashMap>(); + Map> alignedDatasets = new HashMap<>(); for (SequenceI seq : aligned.getSequences()) { SequenceI ds = seq.getDatasetSequence(); @@ -2799,15 +2559,22 @@ public class AlignmentUtils } /* - * first pass - check whether all sequences to be aligned share a dataset - * sequence with an aligned sequence + * 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()) { - if (!alignedDatasets.containsKey(seq.getDatasetSequence())) + 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); } /* @@ -2815,13 +2582,32 @@ public class AlignmentUtils * 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()); - // TODO: getSequenceAsString() will be deprecated in the future - // TODO: need to leave to SequenceI implementor to update gaps - seq.setSequence(alignedSequences.get(0).getSequenceAsString()); + List alignedSequences = alignedDatasets + .get(seq.getDatasetSequence()); + if (alignedSequences.isEmpty()) + { + /* + * defensive check - shouldn't happen! (JAL-3536) + */ + continue; + } + 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) @@ -2829,6 +2615,12 @@ public class AlignmentUtils } } + /* + * finally remove gapped columns (e.g. introns) + */ + new RemoveGapColCommand("", unaligned.getSequencesArray(), 0, + unaligned.getWidth() - 1, unaligned); + return true; } @@ -2842,14 +2634,15 @@ public class AlignmentUtils * @return */ static SortedMap> buildMappedColumnsMap( - AlignmentI unaligned, AlignmentI aligned, List unmapped) + AlignmentI unaligned, AlignmentI aligned, + List unmapped) { /* * Map will hold, for each aligned column position, a map of * {unalignedSequence, characterPerSequence} at that position. * TreeMap keeps the entries in ascending column order. */ - SortedMap> map = new TreeMap>(); + SortedMap> map = new TreeMap<>(); /* * record any sequences that have no mapping so can't be realigned @@ -2877,7 +2670,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. @@ -2906,13 +2700,11 @@ public class AlignmentUtils */ if (seqMap.getTo() == fromSeq.getDatasetSequence()) { - seqMap = new Mapping(seq.getDatasetSequence(), seqMap.getMap() - .getInverse()); + seqMap = new Mapping(seq.getDatasetSequence(), + seqMap.getMap().getInverse()); } - char[] fromChars = fromSeq.getSequence(); int toStart = seq.getStart(); - char[] toChars = seq.getSequence(); /* * traverse [start, end, start, end...] ranges in fromSeq @@ -2943,10 +2735,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] && fromCol <= fromChars.length + 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 @@ -2955,10 +2747,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);