X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=d8c1cdff513f4bd63ec26f58c1d18a1e134941c6;hb=68ab92fb78438a63d70043b1b7f5740380668d7a;hp=228446bc7cdb8162d58d4e40b349141ebaf7c6ed;hpb=e464da0ea20b0322846ce623598a16f2a1cba7ae;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 228446b..d8c1cdf 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,8 +20,24 @@ */ package jalview.analysis; -import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collection; +import java.util.Collections; +import java.util.HashMap; +import java.util.HashSet; +import java.util.Iterator; +import java.util.LinkedHashMap; +import java.util.List; +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; +import jalview.bin.Cache; +import jalview.commands.RemoveGapColCommand; import jalview.datamodel.AlignedCodon; import jalview.datamodel.AlignedCodonFrame; import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping; @@ -37,7 +53,6 @@ import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceGroup; import jalview.datamodel.SequenceI; import jalview.datamodel.features.SequenceFeatures; -import jalview.io.gff.Gff3Helper; import jalview.io.gff.SequenceOntologyI; import jalview.schemes.ResidueProperties; import jalview.util.Comparison; @@ -45,25 +60,6 @@ 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.HashMap; -import java.util.HashSet; -import java.util.Iterator; -import java.util.LinkedHashMap; -import java.util.List; -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; /** * grab bag of useful alignment manipulation operations Expect these to be @@ -74,12 +70,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 @@ -128,7 +127,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()) { @@ -258,7 +257,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(); @@ -267,7 +266,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); @@ -294,8 +293,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 @@ -465,7 +464,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)) { @@ -536,7 +535,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; } @@ -568,7 +568,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; } @@ -881,7 +882,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); @@ -960,11 +961,12 @@ public class AlignmentUtils .findMappingsForSequence(cdsSeq, mappings); for (AlignedCodonFrame mapping : dnaMappings) { - SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); + List foundMap=new ArrayList<>(); + SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein,foundMap); if (peptide != null) { final int peptideLength = peptide.getLength(); - Mapping map = mapping.getMappingBetween(cdsSeq, peptide); + Mapping map = foundMap.get(0).getMapping(); if (map != null) { MapList mapList = map.getMap(); @@ -1092,7 +1094,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()) @@ -1138,9 +1140,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()) @@ -1319,7 +1321,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); @@ -1456,7 +1458,7 @@ public class AlignmentUtils { continue; } - final List result = new ArrayList(); + final List result = new ArrayList<>(); for (AlignmentAnnotation dsann : datasetAnnotations) { /* @@ -1638,13 +1640,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 @@ -1731,31 +1733,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(), @@ -1865,7 +1872,7 @@ public class AlignmentUtils return; } - MapList newMap = targetToFrom.traverse(fromLoci.getMap()); + MapList newMap = targetToFrom.traverse(fromLoci.getMapping()); if (newMap != null) { @@ -1886,7 +1893,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, @@ -1911,7 +1918,15 @@ public class AlignmentUtils 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; + } } /* @@ -1936,10 +1951,12 @@ public class AlignmentUtils { /* * 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()) @@ -1969,6 +1986,19 @@ public class AlignmentUtils static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, AlignmentI dataset) { + /* + * 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; + + /* + * construct CDS sequence by splicing mapped from ranges + */ char[] seqChars = seq.getSequence(); List fromRanges = mapping.getMap().getFromRanges(); int cdsWidth = MappingUtils.getLength(fromRanges); @@ -1993,15 +2023,10 @@ public class AlignmentUtils 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()); @@ -2029,9 +2054,8 @@ public class AlignmentUtils } else { - System.err.println( - "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" - + mtch.toString()); + Cache.log.error( + "JAL-2154 regression: warning - found (and ignored) a duplicate CDS sequence:" + mtch.toString()); } } } @@ -2054,9 +2078,11 @@ public class AlignmentUtils protected static List propagateDBRefsToCDS(SequenceI cdsSeq, SequenceI contig, SequenceI proteinProduct, Mapping mapping) { + // gather direct refs from contig congruent with mapping - List direct = new ArrayList(); - HashSet directSources = new HashSet(); + List direct = new ArrayList<>(); + HashSet directSources = new HashSet<>(); + if (contig.getDBRefs() != null) { for (DBRefEntry dbr : contig.getDBRefs()) @@ -2076,7 +2102,7 @@ public class AlignmentUtils DBRefEntry[] 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) @@ -2142,6 +2168,10 @@ 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 @@ -2217,7 +2247,10 @@ public class AlignmentUtils /** * 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 @@ -2229,6 +2262,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(); @@ -2243,7 +2286,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) @@ -2252,8 +2295,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 }); @@ -2264,7 +2311,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. @@ -2274,7 +2321,7 @@ public class AlignmentUtils */ protected static List findCdsPositions(SequenceI dnaSeq) { - List result = new ArrayList(); + List result = new ArrayList<>(); List sfs = dnaSeq.getFeatures().getFeaturesByOntology( SequenceOntologyI.CDS); @@ -2283,7 +2330,6 @@ public class AlignmentUtils return result; } SequenceFeatures.sortFeatures(sfs, true); - int startPhase = 0; for (SequenceFeature sf : sfs) { @@ -2301,7 +2347,7 @@ public class AlignmentUtils */ int begin = sf.getBegin(); int end = sf.getEnd(); - if (result.isEmpty()) + if (result.isEmpty() && phase > 0) { begin += phase; if (begin > end) @@ -2316,16 +2362,6 @@ public class AlignmentUtils } /* - * 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 @@ -2337,356 +2373,6 @@ public class AlignmentUtils } /** - * 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); - } - - 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(Gff3Helper.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(Gff3Helper.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(Gff3Helper.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; - SequenceFeature sf = new SequenceFeature( - SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, - peptidePos, 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. - *

- * This depends on dna variants being held as a comma-separated list as - * property "alleles" on variant features. - * - * @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[]>(); - - List dnaFeatures = dnaSeq.getFeatures() - .getFeaturesByOntology(SequenceOntologyI.SEQUENCE_VARIANT); - if (dnaFeatures.isEmpty()) - { - 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; - } - - /* - * ignore variant if not a SNP - */ - String alls = (String) sf.getValue(Gff3Helper.ALLELES); - if (alls == null) - { - continue; // non-SNP VCF variant perhaps - can't process this - } - - String[] alleles = alls.toUpperCase().split(","); - boolean isSnp = true; - for (String allele : alleles) - { - if (allele.trim().length() > 1) - { - isSnp = false; - } - } - if (!isSnp) - { - continue; - } - - 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); - } - - /* - * 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. @@ -2759,7 +2445,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(); @@ -2819,10 +2505,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, @@ -2834,7 +2520,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(); @@ -2846,15 +2532,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); } /* @@ -2862,13 +2555,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()); + 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) @@ -2876,6 +2588,12 @@ public class AlignmentUtils } } + /* + * finally remove gapped columns (e.g. introns) + */ + new RemoveGapColCommand("", unaligned.getSequencesArray(), 0, + unaligned.getWidth() - 1, unaligned); + return true; } @@ -2897,7 +2615,7 @@ public class AlignmentUtils * {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 @@ -3002,7 +2720,7 @@ public class AlignmentUtils Map seqsMap = map.get(fromCol); if (seqsMap == null) { - seqsMap = new HashMap(); + seqsMap = new HashMap<>(); map.put(fromCol, seqsMap); } seqsMap.put(seq, seq.getCharAt(mappedCharPos - toStart));