X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=20fcf25b9b13c433e8b83698a84963b419fa4d3b;hb=2b237dfb6fc6d232af79673700b2f57ae8d5a10f;hp=eb1ee4bb04bbd0c9eef128d0cc67c57b98146ab2;hpb=6ed535f7ef953468f8827255ec6ebcd5a6e54d8d;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index eb1ee4b..20fcf25 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -20,8 +20,12 @@ */ package jalview.analysis; +import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE; + +import jalview.api.DBRefEntryI; import jalview.datamodel.AlignedCodon; import jalview.datamodel.AlignedCodonFrame; +import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping; import jalview.datamodel.Alignment; import jalview.datamodel.AlignmentAnnotation; import jalview.datamodel.AlignmentI; @@ -36,10 +40,13 @@ import jalview.io.gff.SequenceOntologyFactory; import jalview.io.gff.SequenceOntologyI; import jalview.schemes.ResidueProperties; import jalview.util.Comparison; +import jalview.util.DBRefUtils; 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; @@ -66,6 +73,31 @@ import java.util.TreeMap; public class AlignmentUtils { + private static final String SEQUENCE_VARIANT = "sequence_variant:"; + private static final String ID = "ID"; + + /** + * A data model to hold the 'normal' base value at a position, and an optional + * sequence variant feature + */ + static class DnaVariant + { + String base; + + SequenceFeature variant; + + DnaVariant(String nuc) + { + base = nuc; + } + + DnaVariant(String nuc, SequenceFeature var) + { + base = nuc; + variant = var; + } + } + /** * given an existing alignment, create a new alignment including all, or up to * flankSize additional symbols from each sequence's dataset sequence @@ -820,6 +852,11 @@ public class AlignmentUtils */ public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna) { + if (protein.isNucleotide() || !dna.isNucleotide()) + { + System.err.println("Wrong alignment type in alignProteinAsDna"); + return 0; + } List unmappedProtein = new ArrayList(); Map> alignedCodons = buildCodonColumnsMap( protein, dna, unmappedProtein); @@ -827,6 +864,162 @@ public class AlignmentUtils } /** + * Realigns the given dna to match the alignment of the protein, using codon + * mappings to translate aligned peptide positions to codons. + * + * @param dna + * the alignment whose sequences are realigned by this method + * @param protein + * the protein alignment whose alignment we are 'copying' + * @return the number of sequences that were realigned + */ + public static int alignCdsAsProtein(AlignmentI dna, AlignmentI protein) + { + if (protein.isNucleotide() || !dna.isNucleotide()) + { + System.err.println("Wrong alignment type in alignProteinAsDna"); + return 0; + } + // todo: implement this + List mappings = protein.getCodonFrames(); + int alignedCount = 0; + for (SequenceI dnaSeq : dna.getSequences()) + { + if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings, + dna.getGapCharacter())) + { + alignedCount++; + } + } + return alignedCount; + } + + /** + * Helper method to align (if possible) the dna sequence to match the + * alignment of a mapped protein sequence. This is currently limited to + * handling coding sequence only. + * + * @param cdsSeq + * @param protein + * @param mappings + * @param gapChar + * @return + */ + static boolean alignCdsSequenceAsProtein(SequenceI cdsSeq, + AlignmentI protein, List mappings, char gapChar) + { + SequenceI cdsDss = cdsSeq.getDatasetSequence(); + if (cdsDss == null) + { + System.err + .println("alignCdsSequenceAsProtein needs aligned sequence!"); + return false; + } + + List dnaMappings = MappingUtils + .findMappingsForSequence(cdsSeq, mappings); + for (AlignedCodonFrame mapping : dnaMappings) + { + SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein); + int peptideLength = peptide.getLength(); + if (peptide != null) + { + Mapping map = mapping.getMappingBetween(cdsSeq, peptide); + if (map != null) + { + MapList mapList = map.getMap(); + if (map.getTo() == peptide.getDatasetSequence()) + { + mapList = mapList.getInverse(); + } + int cdsLength = cdsDss.getLength(); + int mappedFromLength = MappingUtils.getLength(mapList + .getFromRanges()); + int mappedToLength = MappingUtils + .getLength(mapList.getToRanges()); + boolean addStopCodon = (cdsLength == mappedFromLength * 3 + 3) + || (peptide.getDatasetSequence().getLength() == mappedFromLength - 1); + if (cdsLength != mappedToLength && !addStopCodon) + { + System.err + .println(String + .format("Can't align cds as protein (length mismatch %d/%d): %s", + cdsLength, mappedToLength, + cdsSeq.getName())); + } + + /* + * pre-fill the aligned cds sequence with gaps + */ + char[] alignedCds = new char[peptideLength * 3 + + (addStopCodon ? 3 : 0)]; + Arrays.fill(alignedCds, gapChar); + + /* + * walk over the aligned peptide sequence and insert mapped + * codons for residues in the aligned cds sequence + */ + char[] alignedPeptide = peptide.getSequence(); + char[] nucleotides = cdsDss.getSequence(); + int copiedBases = 0; + int cdsStart = cdsDss.getStart(); + int proteinPos = peptide.getStart() - 1; + int cdsCol = 0; + for (char residue : alignedPeptide) + { + if (Comparison.isGap(residue)) + { + cdsCol += 3; + } + else + { + proteinPos++; + int[] codon = mapList.locateInTo(proteinPos, proteinPos); + if (codon == null) + { + // e.g. incomplete start codon, X in peptide + cdsCol += 3; + } + else + { + for (int j = codon[0]; j <= codon[1]; j++) + { + char mappedBase = nucleotides[j - cdsStart]; + alignedCds[cdsCol++] = mappedBase; + copiedBases++; + } + } + } + } + + /* + * append stop codon if not mapped from protein, + * closing it up to the end of the mapped sequence + */ + if (copiedBases == nucleotides.length - 3) + { + for (int i = alignedCds.length - 1; i >= 0; i--) + { + if (!Comparison.isGap(alignedCds[i])) + { + cdsCol = i + 1; // gap just after end of sequence + break; + } + } + for (int i = nucleotides.length - 3; i < nucleotides.length; i++) + { + alignedCds[cdsCol++] = nucleotides[i]; + } + } + cdsSeq.setSequence(new String(alignedCds)); + return true; + } + } + } + return false; + } + + /** * Builds a map whose key is an aligned codon position (3 alignment column * numbers base 0), and whose value is a map from protein sequence to each * protein's peptide residue for that codon. The map generates an ordering of @@ -1303,15 +1496,19 @@ public class AlignmentUtils Collection types, List forSequences, boolean anyType, boolean doShow) { - for (AlignmentAnnotation aa : al.getAlignmentAnnotation()) + AlignmentAnnotation[] anns = al.getAlignmentAnnotation(); + if (anns != null) { - if (anyType || types.contains(aa.label)) + for (AlignmentAnnotation aa : anns) { - if ((aa.sequenceRef != null) - && (forSequences == null || forSequences - .contains(aa.sequenceRef))) + if (anyType || types.contains(aa.label)) { - aa.visible = doShow; + if ((aa.sequenceRef != null) + && (forSequences == null || forSequences + .contains(aa.sequenceRef))) + { + aa.visible = doShow; + } } } } @@ -1366,78 +1563,276 @@ public class AlignmentUtils * Constructs an alignment consisting of the mapped (CDS) regions in the given * nucleotide sequences, and updates mappings to match. The CDS sequences are * added to the original alignment's dataset, which is shared by the new - * alignment. + * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are + * added to the alignment dataset. * * @param dna - * aligned dna sequences - * @param mappings - * from dna to protein; these are replaced with new mappings - * @param al + * aligned nucleotide (dna or cds) sequences + * @param dataset + * the alignment dataset the sequences belong to + * @param products + * (optional) to restrict results to CDS that map to specified + * protein products * @return an alignment whose sequences are the cds-only parts of the dna * sequences (or null if no mappings are found) */ public static AlignmentI makeCdsAlignment(SequenceI[] dna, - List mappings, AlignmentI al) + AlignmentI dataset, SequenceI[] products) { + if (dataset == null || dataset.getDataset() != null) + { + throw new IllegalArgumentException( + "IMPLEMENTATION ERROR: dataset.getDataset() must be null!"); + } + List foundSeqs = new ArrayList(); List cdsSeqs = new ArrayList(); - - for (SequenceI seq : dna) + List mappings = dataset.getCodonFrames(); + HashSet productSeqs = null; + if (products != null) + { + productSeqs = new HashSet(); + for (SequenceI seq : products) + { + productSeqs.add(seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence()); + } + } + + /* + * Construct CDS sequences from mappings on the alignment dataset. + * The logic is: + * - find the protein product(s) mapped to from each dna sequence + * - if the mapping covers the whole dna sequence (give or take start/stop + * codon), take the dna as the CDS sequence + * - else search dataset mappings for a suitable dna sequence, i.e. one + * whose whole sequence is mapped to the protein + * - if no sequence found, construct one from the dna sequence and mapping + * (and add it to dataset so it is found if this is repeated) + */ + for (SequenceI dnaSeq : dna) { - AlignedCodonFrame cdsMappings = new AlignedCodonFrame(); + SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + List seqMappings = MappingUtils - .findMappingsForSequence(seq, mappings); - List alignmentMappings = al.getCodonFrames(); + .findMappingsForSequence(dnaSeq, mappings); for (AlignedCodonFrame mapping : seqMappings) { - for (Mapping aMapping : mapping.getMappingsFromSequence(seq)) + List mappingsFromSequence = mapping + .getMappingsFromSequence(dnaSeq); + + for (Mapping aMapping : mappingsFromSequence) { - SequenceI cdsSeq = makeCdsSequence(seq.getDatasetSequence(), - aMapping); + MapList mapList = aMapping.getMap(); + if (mapList.getFromRatio() == 1) + { + /* + * not a dna-to-protein mapping (likely dna-to-cds) + */ + continue; + } + + /* + * skip if mapping is not to one of the target set of proteins + */ + SequenceI proteinProduct = aMapping.getTo(); + if (productSeqs != null && !productSeqs.contains(proteinProduct)) + { + continue; + } + + /* + * try to locate the CDS from the dataset mappings; + * guard against duplicate results (for the case that protein has + * dbrefs to both dna and cds sequences) + */ + SequenceI cdsSeq = findCdsForProtein(mappings, dnaSeq, + seqMappings, aMapping); + if (cdsSeq != null) + { + if (!foundSeqs.contains(cdsSeq)) + { + foundSeqs.add(cdsSeq); + SequenceI derivedSequence = cdsSeq.deriveSequence(); + cdsSeqs.add(derivedSequence); + if (!dataset.getSequences().contains(cdsSeq)) + { + dataset.addSequence(cdsSeq); + } + } + continue; + } + + /* + * didn't find mapped CDS sequence - construct it and add + * its dataset sequence to the dataset + */ + cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping); + SequenceI cdsSeqDss = cdsSeq.createDatasetSequence(); cdsSeqs.add(cdsSeq); - + if (!dataset.getSequences().contains(cdsSeqDss)) + { + dataset.addSequence(cdsSeqDss); + } + /* * add a mapping from CDS to the (unchanged) mapped to range */ List cdsRange = Collections.singletonList(new int[] { 1, cdsSeq.getLength() }); - MapList map = new MapList(cdsRange, aMapping.getMap() - .getToRanges(), aMapping.getMap().getFromRatio(), - aMapping.getMap().getToRatio()); - cdsMappings.addMap(cdsSeq, aMapping.getTo(), map); + MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(), + mapList.getFromRatio(), mapList.getToRatio()); + AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); + cdsToProteinMapping.addMap(cdsSeq, proteinProduct, cdsToProteinMap); + + /* + * guard against duplicating the mapping if repeating this action + */ + if (!mappings.contains(cdsToProteinMapping)) + { + mappings.add(cdsToProteinMapping); + } + + /* + * copy protein's dbrefs to CDS sequence + * this enables Get Cross-References from CDS alignment + */ + DBRefEntry[] proteinRefs = DBRefUtils.selectDbRefs(false, + proteinProduct.getDBRefs()); + if (proteinRefs != null) + { + for (DBRefEntry ref : proteinRefs) + { + DBRefEntry cdsToProteinRef = new DBRefEntry(ref); + cdsToProteinRef.setMap(new Mapping(proteinProduct, + cdsToProteinMap)); + cdsSeqDss.addDBRef(cdsToProteinRef); + } + } /* * add another mapping from original 'from' range to CDS */ - map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1, + AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); + MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), + cdsRange, 1, 1); - cdsMappings.addMap(seq.getDatasetSequence(), cdsSeq, map); + dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq, + dnaToCdsMap); + if (!mappings.contains(dnaToCdsMapping)) + { + mappings.add(dnaToCdsMapping); + } - alignmentMappings.add(cdsMappings); + /* + * add DBRef with mapping from protein to CDS + * (this enables Get Cross-References from protein alignment) + * This is tricky because we can't have two DBRefs with the + * same source and accession, so need a different accession for + * the CDS from the dna sequence + */ + DBRefEntryI dnaRef = dnaDss.getSourceDBRef(); + if (dnaRef != null) + { + // assuming cds version same as dna ?!? + DBRefEntry proteinToCdsRef = new DBRefEntry(dnaRef.getSource(), + dnaRef.getVersion(), cdsSeq.getName()); + proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap + .getInverse())); + proteinProduct.addDBRef(proteinToCdsRef); + } /* * transfer any features on dna that overlap the CDS */ - transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS); + transferFeatures(dnaSeq, cdsSeq, cdsToProteinMap, null, + SequenceOntologyI.CDS); } } } + AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs + .size()])); + cds.setDataset(dataset); + + return cds; + } + + /** + * A helper method that finds a CDS sequence in the alignment dataset that is + * mapped to the given protein sequence, and either is, or has a mapping from, + * the given dna sequence. + * + * @param mappings + * set of all mappings on the dataset + * @param dnaSeq + * a dna (or cds) sequence we are searching from + * @param seqMappings + * the set of mappings involving dnaSeq + * @param aMapping + * an initial candidate from seqMappings + * @return + */ + static SequenceI findCdsForProtein(List mappings, + SequenceI dnaSeq, List seqMappings, + Mapping aMapping) + { + /* + * TODO a better dna-cds-protein mapping data representation to allow easy + * navigation; until then this clunky looping around lists of mappings + */ + SequenceI seqDss = dnaSeq.getDatasetSequence() == null ? dnaSeq + : dnaSeq.getDatasetSequence(); + SequenceI proteinProduct = aMapping.getTo(); + /* - * add CDS seqs to shared dataset + * is this mapping from the whole dna sequence (i.e. CDS)? + * allowing for possible stop codon on dna but not peptide */ - Alignment dataset = al.getDataset(); - for (SequenceI seq : cdsSeqs) + int mappedFromLength = MappingUtils.getLength(aMapping.getMap() + .getFromRanges()); + int dnaLength = seqDss.getLength(); + if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - 3) { - if (!dataset.getSequences().contains(seq.getDatasetSequence())) + return seqDss; + } + + /* + * looks like we found the dna-to-protein mapping; search for the + * corresponding cds-to-protein mapping + */ + List mappingsToPeptide = MappingUtils + .findMappingsForSequence(proteinProduct, mappings); + for (AlignedCodonFrame acf : mappingsToPeptide) + { + for (SequenceToSequenceMapping map : acf.getMappings()) { - dataset.addSequence(seq.getDatasetSequence()); + Mapping mapping = map.getMapping(); + if (mapping != aMapping && mapping.getMap().getFromRatio() == 3 + && proteinProduct == mapping.getTo() + && seqDss != map.getFromSeq()) + { + 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 + */ + SequenceI cdsSeq = map.getFromSeq(); + List dnaToCdsMaps = MappingUtils + .findMappingsForSequence(cdsSeq, seqMappings); + if (!dnaToCdsMaps.isEmpty()) + { + return cdsSeq; + } + } + } } } - AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs - .size()])); - cds.setDataset(dataset); - - return cds; + return null; } /** @@ -1447,7 +1842,7 @@ public class AlignmentUtils * * @param seq * @param mapping - * @return + * @return CDS sequence (as a dataset sequence) */ static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) { @@ -1477,9 +1872,15 @@ public class AlignmentUtils } } - SequenceI newSeq = new Sequence(seq.getName() + "|" - + mapping.getTo().getName(), newSeqChars, 1, newPos); - newSeq.createDatasetSequence(); + /* + * 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); + // newSeq.setDescription(mapFromId); + return newSeq; } @@ -1743,66 +2144,230 @@ public class AlignmentUtils * /ENSP00000288602?feature=transcript_variation;content-type=text/xml * which would be a bit slower but possibly more reliable */ - LinkedHashMap variants = buildDnaVariantsMap( + + /* + * 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()) + for (Entry[]> variant : variants.entrySet()) { int peptidePos = variant.getKey(); - String[][] codonVariants = variant.getValue(); - String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); // 0-based - List peptideVariants = computePeptideVariants(codonVariants, - residue); - if (!peptideVariants.isEmpty()) - { - String desc = residue + "," // include canonical residue in description - + StringUtils.listToDelimitedString(peptideVariants, ", "); - SequenceFeature sf = new SequenceFeature( - SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos, - peptidePos, 0f, null); - peptide.addSequenceFeature(sf); - count++; + 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++; + } + } + } } } /* - * ugly sort to get sequence features in start position order - * - would be better to store in Sequence as a TreeSet instead? + * variants in second codon base */ - Arrays.sort(peptide.getSequenceFeatures(), - new Comparator() + 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)) { - @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; - } - }); + 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; } /** - * Builds a map whose key is position in the protein sequence, and value is an - * array of all variants for the coding codon positions + * 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() > 3 ? 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, "Jalview"); + StringBuilder attributes = new StringBuilder(32); + String id = (String) var.variant.getValue(ID); + if (id != null) + { + if (id.startsWith(SEQUENCE_VARIANT)) + { + id = id.substring(SEQUENCE_VARIANT.length()); + } + sf.setValue(ID, id); + attributes.append(ID).append("=").append(id); + // TODO handle other species variants + StringBuilder link = new StringBuilder(32); + try + { + link.append(desc).append(" ").append(id) + .append("|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=") + .append(URLEncoder.encode(id, "UTF-8")); + sf.addLink(link.toString()); + } catch (UnsupportedEncodingException e) + { + // as if + } + } + String clinSig = (String) var.variant + .getValue(CLINICAL_SIGNIFICANCE); + if (clinSig != null) + { + 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 */ - static LinkedHashMap buildDnaVariantsMap( + static LinkedHashMap[]> buildDnaVariantsMap( SequenceI dnaSeq, MapList dnaToProtein) { /* - * map from peptide position to all variant features of the codon for it - * LinkedHashMap ensures we add the peptide features in sequence order + * 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(); + LinkedHashMap[]> variants = new LinkedHashMap[]>(); SequenceOntologyI so = SequenceOntologyFactory.getInstance(); SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures(); @@ -1835,10 +2400,13 @@ public class AlignmentUtils continue; } int peptidePosition = mapsTo[0]; - String[][] codonVariants = variants.get(peptidePosition); + List[] codonVariants = variants.get(peptidePosition); if (codonVariants == null) { - codonVariants = new String[3][]; + codonVariants = new ArrayList[3]; + codonVariants[0] = new ArrayList(); + codonVariants[1] = new ArrayList(); + codonVariants[2] = new ArrayList(); variants.put(peptidePosition, codonVariants); } @@ -1867,31 +2435,38 @@ public class AlignmentUtils lastCodon = codon; /* - * save nucleotide (and this variant) for each codon position + * save nucleotide (and any variant) for each codon position */ for (int codonPos = 0; codonPos < 3; codonPos++) { String nucleotide = String.valueOf( dnaSeq.getCharAt(codon[codonPos] - dnaStart)) .toUpperCase(); - if (codonVariants[codonPos] == null) + List codonVariant = codonVariants[codonPos]; + if (codon[codonPos] == dnaCol) { - /* - * record current dna base - */ - codonVariants[codonPos] = new String[] { nucleotide }; + 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)); + } } - if (codon[codonPos] == dnaCol) + else if (codonVariant.isEmpty()) { /* - * add alleles to dna base (and any previously found alleles) + * record (possibly non-varying) base value */ - String[] known = codonVariants[codonPos]; - String[] dnaVariants = new String[alleles.length + known.length]; - System.arraycopy(known, 0, dnaVariants, 0, known.length); - System.arraycopy(alleles, 0, dnaVariants, known.length, - alleles.length); - codonVariants[codonPos] = dnaVariants; + codonVariant.add(new DnaVariant(nucleotide)); } } } @@ -1900,83 +2475,21 @@ public class AlignmentUtils } /** - * Returns a sorted, non-redundant list of all peptide translations generated - * by the given dna variants, excluding the current residue value - * - * @param codonVariants - * an array of base values (acgtACGT) for codon positions 1, 2, 3 - * @param residue - * the current residue translation - * @return - */ - static List computePeptideVariants(String[][] codonVariants, - String residue) - { - List result = new ArrayList(); - for (String base1 : codonVariants[0]) - { - for (String base2 : codonVariants[1]) - { - for (String base3 : codonVariants[2]) - { - String codon = base1 + base2 + base3; - /* - * 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 ignored here - */ - String peptide = codon.contains("-") ? "-" - : (codon.length() > 3 ? null : ResidueProperties - .codonTranslate(codon)); - if (peptide != null && !result.contains(peptide) - && !peptide.equalsIgnoreCase(residue)) - { - result.add(peptide); - } - } - } - } - - /* - * sort alphabetically with STOP at the end - */ - Collections.sort(result, new Comparator() - { - - @Override - public int compare(String o1, String o2) - { - if ("STOP".equals(o1)) - { - return 1; - } - else if ("STOP".equals(o2)) - { - return -1; - } - else - { - return o1.compareTo(o2); - } - } - }); - return result; - } - - /** * 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. * * @param seqs * @param xrefs + * @param dataset + * the alignment dataset shared by the new copy * @return */ public static AlignmentI makeCopyAlignment(SequenceI[] seqs, - SequenceI[] xrefs) + SequenceI[] xrefs, AlignmentI dataset) { AlignmentI copy = new Alignment(new Alignment(seqs)); + copy.setDataset(dataset); SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) @@ -2087,13 +2600,13 @@ public class AlignmentUtils { /* * Map will hold, for each aligned column position, a map of - * {unalignedSequence, sequenceCharacter} at that position. + * {unalignedSequence, characterPerSequence} at that position. * TreeMap keeps the entries in ascending column order. */ Map> map = new TreeMap>(); /* - * report any sequences that have no mapping so can't be realigned + * record any sequences that have no mapping so can't be realigned */ unmapped.addAll(unaligned.getSequences()); @@ -2106,7 +2619,7 @@ public class AlignmentUtils SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned); if (fromSeq != null) { - Mapping seqMap = mapping.getMappingForSequence(seq); + Mapping seqMap = mapping.getMappingBetween(fromSeq, seq); if (addMappedPositions(seq, fromSeq, seqMap, map)) { unmapped.remove(seq); @@ -2137,6 +2650,20 @@ public class AlignmentUtils static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq, Mapping seqMap, Map> map) { + if (seqMap == null) + { + return false; + } + + /* + * invert mapping if it is from unaligned to aligned sequence + */ + if (seqMap.getTo() == fromSeq.getDatasetSequence()) + { + seqMap = new Mapping(seq.getDatasetSequence(), seqMap.getMap() + .getInverse()); + } + char[] fromChars = fromSeq.getSequence(); int toStart = seq.getStart(); char[] toChars = seq.getSequence(); @@ -2170,7 +2697,8 @@ public class AlignmentUtils * of the next character of the mapped-to sequence; stop when all * the characters of the range have been counted */ - while (mappedCharPos <= range[1]) + while (mappedCharPos <= range[1] && fromCol <= fromChars.length + && fromCol >= 0) { if (!Comparison.isGap(fromChars[fromCol - 1])) { @@ -2193,4 +2721,19 @@ public class AlignmentUtils } return true; } + + // strictly temporary hack until proper criteria for aligning protein to cds + // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot + public static boolean looksLikeEnsembl(AlignmentI alignment) + { + for (SequenceI seq : alignment.getSequences()) + { + String name = seq.getName(); + if (!name.startsWith("ENSG") && !name.startsWith("ENST")) + { + return false; + } + } + return true; + } }