X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=d04d810e3c102a12aea5c1e7c30c925bf7b1e91c;hb=refs%2Fheads%2Freleases%2FRelease_2_11_1_4_debianpatches;hp=77cc4d6218d8a94b5d768299431e94c7111f469a;hpb=b22d623e48f01e7c0ab274e30b1b05e289afc34f;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 77cc4d6..d04d810 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 @@ -965,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(); @@ -982,10 +979,10 @@ public class AlignmentUtils .getFromRanges()); int mappedToLength = MappingUtils .getLength(mapList.getToRanges()); - boolean addStopCodon = (cdsLength == mappedFromLength + boolean addStopCodon = peptide.getDatasetSequence().getEnd()==peptide.getEnd() && ((cdsLength == mappedFromLength * CODON_LENGTH + CODON_LENGTH) || (peptide.getDatasetSequence() - .getLength() == mappedFromLength - 1); + .getLength() == mappedFromLength - 1)); if (cdsLength != mappedToLength && !addStopCodon) { System.err.println(String.format( @@ -1736,15 +1733,8 @@ 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[] @@ -1753,16 +1743,26 @@ public class AlignmentUtils 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(), @@ -1996,45 +1996,31 @@ public class AlignmentUtils SequenceI newSeq = null; - final MapList maplist = mapping.getMap(); - if (maplist.isContiguous() && maplist.isFromForwardStrand()) - { - /* - * 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]; + /* + * construct CDS sequence by splicing mapped from ranges + */ + char[] seqChars = seq.getSequence(); + List fromRanges = mapping.getMap().getFromRanges(); + int cdsWidth = MappingUtils.getLength(fromRanges); + char[] newSeqChars = new char[cdsWidth]; - int newPos = 0; - for (int[] range : fromRanges) + int newPos = 0; + for (int[] range : fromRanges) + { + if (range[0] <= range[1]) { - 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 + // 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 + { + // reverse strand mapping - copy and complement one by one + for (int i = range[0]; i >= range[1]; i--) { - // 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]); - } + newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]); } } @@ -2068,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()); } } } @@ -2388,394 +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 - peptide.getStart())); - 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(",")) - { - if (!base1.equalsIgnoreCase(base)) - { - String codon = base.toUpperCase() + base2.toLowerCase() - + base3.toLowerCase(); - String canonical = base1.toUpperCase() + base2.toLowerCase() - + base3.toLowerCase(); - if (addPeptideVariant(peptide, peptidePos, residue, var, - codon, canonical)) - { - 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(",")) - { - if (!base2.equalsIgnoreCase(base)) - { - String codon = base1.toLowerCase() + base.toUpperCase() - + base3.toLowerCase(); - String canonical = base1.toLowerCase() + base2.toUpperCase() - + base3.toLowerCase(); - if (addPeptideVariant(peptide, peptidePos, residue, var, - codon, canonical)) - { - 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(",")) - { - if (!base3.equalsIgnoreCase(base)) - { - String codon = base1.toLowerCase() + base2.toLowerCase() - + base.toUpperCase(); - String canonical = base1.toLowerCase() + base2.toLowerCase() - + base3.toUpperCase(); - if (addPeptideVariant(peptide, peptidePos, residue, var, - codon, canonical)) - { - count++; - } - } - } - } - } - } - - return count; - } - - /** - * 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 canonical) - { - /* - * 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("-") ? 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)); - 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)) - { - id = id.substring(SEQUENCE_VARIANT.length()); - } - 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) - { - // 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; - } - - /** - * 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. @@ -2837,6 +2434,17 @@ public class AlignmentUtils */ public static int alignAs(AlignmentI unaligned, AlignmentI aligned) { + /* + * easy case - aligning a copy of aligned sequences + */ + if (alignAsSameSequences(unaligned, aligned)) + { + return unaligned.getHeight(); + } + + /* + * fancy case - aligning via mappings between sequences + */ List unmapped = new ArrayList<>(); Map> columnMap = buildMappedColumnsMap( unaligned, aligned, unmapped); @@ -2902,9 +2510,7 @@ public class AlignmentUtils * - 'guide' alignment containing sequences derived from same * dataset as unaligned * @return - * @deprecated probably obsolete and incomplete */ - @Deprecated static boolean alignAsSameSequences(AlignmentI unaligned, AlignmentI aligned) { @@ -2954,6 +2560,13 @@ public class AlignmentUtils { List alignedSequences = alignedDatasets .get(seq.getDatasetSequence()); + if (alignedSequences.isEmpty()) + { + /* + * defensive check - shouldn't happen! (JAL-3536) + */ + continue; + } SequenceI alignedSeq = alignedSequences.get(0); /* @@ -2975,6 +2588,12 @@ public class AlignmentUtils } } + /* + * finally remove gapped columns (e.g. introns) + */ + new RemoveGapColCommand("", unaligned.getSequencesArray(), 0, + unaligned.getWidth() - 1, unaligned); + return true; }