X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=34fe221a98d9099fe4e072f403a7bf12b61d580e;hb=37de9310bec3501cbc6381e0c3dcb282fcaad812;hp=d1cd5a31fbacfd0681375ce886f015483e2a59db;hpb=6b5f772a654104d856da1651c3d411981be463df;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index d1cd5a3..34fe221 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -22,7 +22,6 @@ 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; @@ -74,7 +73,9 @@ public class AlignmentUtils { private static final int CODON_LENGTH = 3; + private static final String SEQUENCE_VARIANT = "sequence_variant:"; + private static final String ID = "ID"; /** @@ -447,8 +448,8 @@ public class AlignmentUtils */ if (cdnaLength != mappedLength && cdnaLength > 2) { - String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - CODON_LENGTH, CODON_LENGTH) - .toUpperCase(); + String lastCodon = String.valueOf(cdnaSeqChars, + cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase(); for (String stop : ResidueProperties.STOP) { if (lastCodon.equals(stop)) @@ -510,8 +511,7 @@ 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); @@ -937,7 +937,7 @@ public class AlignmentUtils .println("alignCdsSequenceAsProtein needs aligned sequence!"); return false; } - + List dnaMappings = MappingUtils .findMappingsForSequence(cdsSeq, mappings); for (AlignedCodonFrame mapping : dnaMappings) @@ -959,7 +959,8 @@ public class AlignmentUtils .getFromRanges()); int mappedToLength = MappingUtils .getLength(mapList.getToRanges()); - boolean addStopCodon = (cdsLength == mappedFromLength * CODON_LENGTH + CODON_LENGTH) + boolean addStopCodon = (cdsLength == mappedFromLength + * CODON_LENGTH + CODON_LENGTH) || (peptide.getDatasetSequence().getLength() == mappedFromLength - 1); if (cdsLength != mappedToLength && !addStopCodon) { @@ -1098,7 +1099,7 @@ public class AlignmentUtils // TODO resolve JAL-2022 so this fudge can be removed int mappedSequenceCount = protein.getHeight() - unmappedProtein.size(); addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount); - + return alignedCodons; } @@ -1689,11 +1690,20 @@ public class AlignmentUtils * 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(); + cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping, + dataset).deriveSequence(); + // cdsSeq has a name constructed as CDS| + // will be either the accession for the coding sequence, + // marked in the /via/ dbref to the protein product accession + // or it will be the original nucleotide accession. + SequenceI cdsSeqDss = cdsSeq.getDatasetSequence(); + cdsSeqs.add(cdsSeq); + if (!dataset.getSequences().contains(cdsSeqDss)) { + // check if this sequence is a newly created one + // so needs adding to the dataset dataset.addSequence(cdsSeqDss); } @@ -1702,10 +1712,12 @@ public class AlignmentUtils */ List cdsRange = Collections.singletonList(new int[] { 1, cdsSeq.getLength() }); - MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(), - mapList.getFromRatio(), mapList.getToRatio()); + MapList cdsToProteinMap = new MapList(cdsRange, + mapList.getToRanges(), mapList.getFromRatio(), + mapList.getToRatio()); AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame(); - cdsToProteinMapping.addMap(cdsSeq, proteinProduct, cdsToProteinMap); + cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct, + cdsToProteinMap); /* * guard against duplicating the mapping if repeating this action @@ -1715,31 +1727,15 @@ public class AlignmentUtils 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); - } - } - + propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(), + proteinProduct, aMapping); /* * add another mapping from original 'from' range to CDS */ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), - cdsRange, 1, - 1); - dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeq, + cdsRange, 1, 1); + dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, dnaToCdsMap); if (!mappings.contains(dnaToCdsMapping)) { @@ -1753,12 +1749,37 @@ public class AlignmentUtils * same source and accession, so need a different accession for * the CDS from the dna sequence */ - DBRefEntryI dnaRef = dnaDss.getSourceDBRef(); - if (dnaRef != null) + + // specific use case: + // Genomic contig ENSCHR:1, contains coding regions for ENSG01, + // ENSG02, ENSG03, with transcripts and products similarly named. + // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01 + + // JBPNote: ?? can't actually create an example that demonstrates we + // need to + // synthesize an xref. + + for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs()) { + // 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))); + cdsSeqDss.addDBRef(cdsCrossRef); + + // 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(dnaRef.getSource(), - dnaRef.getVersion(), cdsSeq.getName()); + + DBRefEntry proteinToCdsRef = new DBRefEntry( + primRef.getSource(), primRef.getVersion(), + cdsSeq.getName()); + // proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap .getInverse())); proteinProduct.addDBRef(proteinToCdsRef); @@ -1767,7 +1788,7 @@ public class AlignmentUtils /* * transfer any features on dna that overlap the CDS */ - transferFeatures(dnaSeq, cdsSeq, cdsToProteinMap, null, + transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null, SequenceOntologyI.CDS); } } @@ -1814,7 +1835,8 @@ public class AlignmentUtils int mappedFromLength = MappingUtils.getLength(aMapping.getMap() .getFromRanges()); int dnaLength = seqDss.getLength(); - if (mappedFromLength == dnaLength || mappedFromLength == dnaLength - CODON_LENGTH) + if (mappedFromLength == dnaLength + || mappedFromLength == dnaLength - CODON_LENGTH) { return seqDss; } @@ -1830,7 +1852,8 @@ public class AlignmentUtils for (SequenceToSequenceMapping map : acf.getMappings()) { Mapping mapping = map.getMapping(); - if (mapping != aMapping && mapping.getMap().getFromRatio() == CODON_LENGTH + if (mapping != aMapping + && mapping.getMap().getFromRatio() == CODON_LENGTH && proteinProduct == mapping.getTo() && seqDss != map.getFromSeq()) { @@ -1864,9 +1887,14 @@ public class AlignmentUtils * * @param seq * @param mapping + * @param dataset + * - existing dataset. We check for sequences that look like the CDS + * we are about to construct, if one exists already, then we will + * just return that one. * @return CDS sequence (as a dataset sequence) */ - static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping) + static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping, + AlignmentI dataset) { char[] seqChars = seq.getSequence(); List fromRanges = mapping.getMap().getFromRanges(); @@ -1901,12 +1929,124 @@ public class AlignmentUtils 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()); + if (matches != null) + { + boolean matched = false; + for (SequenceI mtch : matches) + { + if (mtch.getStart() != newSeq.getStart()) + { + continue; + } + if (mtch.getEnd() != newSeq.getEnd()) + { + continue; + } + if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence())) + { + continue; + } + if (!matched) + { + matched = true; + newSeq = mtch; + } + else + { + System.err + .println("JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):" + + mtch.toString()); + } + } + } + } // newSeq.setDescription(mapFromId); return newSeq; } /** + * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to + * the given mapping. + * + * @param cdsSeq + * @param contig + * @param mapping + * @return list of DBRefEntrys added. + */ + public 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) + { + for (DBRefEntry dbr : contig.getDBRefs()) + { + if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap()) + { + MapList map = dbr.getMap().getMap(); + // check if map is the CDS mapping + if (mapping.getMap().equals(map)) + { + direct.add(dbr); + directSources.add(dbr.getSource()); + } + } + } + } + DBRefEntry[] onSource = DBRefUtils.selectRefs( + proteinProduct.getDBRefs(), + directSources.toArray(new String[0])); + List propagated = new ArrayList(); + + // and generate appropriate mappings + for (DBRefEntry cdsref : direct) + { + // clone maplist and mapping + MapList cdsposmap = new MapList(Arrays.asList(new int[][] { new int[] + { cdsSeq.getStart(), cdsSeq.getEnd() } }), cdsref.getMap().getMap() + .getToRanges(), 3, 1); + Mapping cdsmap = new Mapping(cdsref.getMap().getTo(), cdsref.getMap() + .getMap()); + + // create dbref + DBRefEntry newref = new DBRefEntry(cdsref.getSource(), + cdsref.getVersion(), cdsref.getAccessionId(), new Mapping( + cdsmap.getTo(), cdsposmap)); + + // and see if we can map to the protein product for this mapping. + // onSource is the filtered set of accessions on protein that we are + // tranferring, so we assume accession is the same. + if (cdsmap.getTo() == null && onSource != null) + { + List sourceRefs = DBRefUtils.searchRefs(onSource, + cdsref.getAccessionId()); + if (sourceRefs != null) + { + for (DBRefEntry srcref : sourceRefs) + { + if (srcref.getSource().equalsIgnoreCase(cdsref.getSource())) + { + // we have found a complementary dbref on the protein product, so + // update mapping's getTo + newref.getMap().setTo(proteinProduct); + } + } + } + } + cdsSeq.addDBRef(newref); + propagated.add(newref); + } + return propagated; + } + + /** * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the * feature start/end ranges, optionally omitting specified feature types. * Returns the number of features copied. @@ -2347,7 +2487,9 @@ public class AlignmentUtils StringBuilder link = new StringBuilder(32); try { - link.append(desc).append(" ").append(id) + 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()); @@ -2356,8 +2498,7 @@ public class AlignmentUtils // as if } } - String clinSig = (String) var.variant - .getValue(CLINICAL_SIGNIFICANCE); + String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE); if (clinSig != null) { sf.setValue(CLINICAL_SIGNIFICANCE, clinSig); @@ -2513,7 +2654,7 @@ public class AlignmentUtils { AlignmentI copy = new Alignment(new Alignment(seqs)); copy.setDataset(dataset); - + boolean isProtein = !copy.isNucleotide(); SequenceIdMatcher matcher = new SequenceIdMatcher(seqs); if (xrefs != null) { @@ -2524,7 +2665,8 @@ public class AlignmentUtils { for (DBRefEntry dbref : dbrefs) { - if (dbref.getMap() == null || dbref.getMap().getTo() == null) + if (dbref.getMap() == null || dbref.getMap().getTo() == null + || dbref.getMap().getTo().isProtein() != isProtein) { continue; } @@ -2605,7 +2747,7 @@ public class AlignmentUtils } newCol++; } - + /* * trim trailing gaps */