X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FCrossRef.java;h=71b0aa0b9adb5cd92d39b5f7908fa3dcc349e4ed;hb=655b78299307682a4c7a6e5af0ed4618cbc9c924;hp=8fd07067962c515b62ac728c59a1578950312fb4;hpb=438af5436afff7bc1ce5ab380ec5813a7beb11b9;p=jalview.git diff --git a/src/jalview/analysis/CrossRef.java b/src/jalview/analysis/CrossRef.java index 8fd0706..71b0aa0 100644 --- a/src/jalview/analysis/CrossRef.java +++ b/src/jalview/analysis/CrossRef.java @@ -28,18 +28,17 @@ import jalview.datamodel.Mapping; import jalview.datamodel.Sequence; import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceI; -import jalview.util.Comparison; import jalview.util.DBRefUtils; import jalview.util.MapList; import jalview.ws.SequenceFetcherFactory; import jalview.ws.seqfetcher.ASequenceFetcher; import java.util.ArrayList; +import java.util.Iterator; import java.util.List; /** - * Functions for cross-referencing sequence databases. user must first specify - * if cross-referencing from protein or dna (set dna==true) + * Functions for cross-referencing sequence databases. * * @author JimP * @@ -47,28 +46,43 @@ import java.util.List; public class CrossRef { /* - * A sub-class that ignores Parent attribute when comparing sequence - * features. This avoids 'duplicate' CDS features that only - * differ in their parent Transcript ids. + * the dataset of the alignment for which we are searching for + * cross-references; in some cases we may resolve xrefs by + * searching in the dataset */ - class MySequenceFeature extends SequenceFeature - { - private SequenceFeature feat; + private AlignmentI dataset; - MySequenceFeature(SequenceFeature sf) - { - this.feat = sf; - } + /* + * the sequences for which we are seeking cross-references + */ + private SequenceI[] fromSeqs; - @Override - public boolean equals(Object o) - { - return feat.equals(o, true); - } - } + /** + * matcher built from dataset + */ + SequenceIdMatcher matcher; /** + * sequences found by cross-ref searches to fromSeqs + */ + List rseqs; + /** + * Constructor + * + * @param seqs + * the sequences for which we are seeking cross-references + * @param ds + * the containing alignment dataset (may be searched to resolve + * cross-references) + */ + public CrossRef(SequenceI[] seqs, AlignmentI ds) + { + fromSeqs = seqs; + dataset = ds.getDataset() == null ? ds : ds.getDataset(); + } + + /** * Returns a list of distinct database sources for which sequences have either * * * @param dna - * true if seqs are nucleotide - * @param seqs - * sequences whose xrefs we are seeking - * @param dataset - * an alignment to search for indirect references + * - when true, cross-references *from* dna returned. When false, + * cross-references *from* protein are returned * @return */ - public static List findXrefSourcesForSequences(boolean dna, - SequenceI[] seqs, AlignmentI dataset) + public List findXrefSourcesForSequences(boolean dna) { List sources = new ArrayList(); - for (SequenceI seq : seqs) + for (SequenceI seq : fromSeqs) { if (seq != null) { - findXrefSourcesForSequence(seq, dna, dataset, sources); + findXrefSourcesForSequence(seq, dna, sources); } } return sources; @@ -110,41 +120,41 @@ public class CrossRef * * @param seq * the sequence whose dbrefs we are searching against - * @param dna - * true if the sequence is nucleotide - * @param dataset - * an alignment to search for indirect references + * @param fromDna + * when true, context is DNA - so sources identifying protein + * products will be returned. * @param sources * a list of sources to add matches to */ - static void findXrefSourcesForSequence(SequenceI seq, boolean dna, - AlignmentI dataset, List sources) + void findXrefSourcesForSequence(SequenceI seq, boolean fromDna, + List sources) { /* * first find seq's xrefs (dna-to-peptide or peptide-to-dna) */ - DBRefEntry[] rfs = DBRefUtils.selectDbRefs(!dna, seq.getDBRefs()); + DBRefEntry[] rfs = DBRefUtils.selectDbRefs(!fromDna, seq.getDBRefs()); addXrefsToSources(rfs, sources); if (dataset != null) { /* * find sequence's direct (dna-to-dna, peptide-to-peptide) xrefs */ - DBRefEntry[] lrfs = DBRefUtils.selectDbRefs(dna, seq.getDBRefs()); - List rseqs = new ArrayList(); + DBRefEntry[] lrfs = DBRefUtils.selectDbRefs(fromDna, seq.getDBRefs()); + List foundSeqs = new ArrayList(); /* * find sequences in the alignment which xref one of these DBRefs * i.e. is xref-ed to a common sequence identifier */ - CrossRef.searchDatasetXrefs(seq, !dna, lrfs, dataset, rseqs, null); + searchDatasetXrefs(fromDna, seq, lrfs, foundSeqs, null); /* * add those sequences' (dna-to-peptide or peptide-to-dna) dbref sources */ - for (SequenceI rs : rseqs) + for (SequenceI rs : foundSeqs) { - DBRefEntry[] xrs = DBRefUtils.selectDbRefs(!dna, rs.getDBRefs()); + DBRefEntry[] xrs = DBRefUtils + .selectDbRefs(!fromDna, rs.getDBRefs()); addXrefsToSources(xrs, sources); } } @@ -157,13 +167,16 @@ public class CrossRef * @param xrefs * @param sources */ - static void addXrefsToSources(DBRefEntry[] xrefs, List sources) + void addXrefsToSources(DBRefEntry[] xrefs, List sources) { if (xrefs != null) { for (DBRefEntry ref : xrefs) { - String source = ref.getSource(); + /* + * avoid duplication e.g. ENSEMBL and Ensembl + */ + String source = DBRefUtils.getCanonicalName(ref.getSource()); if (!sources.contains(source)) { sources.add(source); @@ -173,24 +186,33 @@ public class CrossRef } /** + * Attempts to find cross-references from the sequences provided in the + * constructor to the given source database. Cross-references may be found + *
    + *
  • in dbrefs on the sequence which hold a mapping to a sequence + *
      + *
    • provided with a fetched sequence (e.g. ENA translation), or
      • + *
    • populated previously after getting cross-references
      • + *
    + *
  • as other sequences in the alignment which share a dbref identifier with + * the sequence
    • + *
  • by fetching from the remote database
    • + *
+ * The cross-referenced sequences, and mappings to them, are added to the + * alignment dataset. * - * @param seqs - * sequences whose xrefs are being retrieved - * @param dna - * true if sequences are nucleotide * @param source - * @param al - * alignment to search for cross-referenced sequences (and possibly - * add to) - * @return products (as dataset sequences) + * @return cross-referenced sequences (as dataset sequences) */ - public static Alignment findXrefSequences(SequenceI[] seqs, - final boolean dna, final String source, AlignmentI al) + public Alignment findXrefSequences(String source, boolean fromDna) { - AlignmentI dataset = al.getDataset() == null ? al : al.getDataset(); - List rseqs = new ArrayList(); + + rseqs = new ArrayList(); AlignedCodonFrame cf = new AlignedCodonFrame(); - for (SequenceI seq : seqs) + matcher = new SequenceIdMatcher( + dataset.getSequences()); + + for (SequenceI seq : fromSeqs) { SequenceI dss = seq; while (dss.getDatasetSequence() != null) @@ -198,245 +220,570 @@ public class CrossRef dss = dss.getDatasetSequence(); } boolean found = false; - DBRefEntry[] xrfs = DBRefUtils.selectDbRefs(!dna, dss.getDBRefs()); + DBRefEntry[] xrfs = DBRefUtils + .selectDbRefs(!fromDna, dss.getDBRefs()); + // ENST & ENSP comes in to both Protein and nucleotide, so we need to + // filter them + // out later. if ((xrfs == null || xrfs.length == 0) && dataset != null) { /* * found no suitable dbrefs on sequence - look for sequences in the * alignment which share a dbref with this one */ - DBRefEntry[] lrfs = DBRefUtils.selectDbRefs(dna, seq.getDBRefs()); + DBRefEntry[] lrfs = DBRefUtils.selectDbRefs(fromDna, + seq.getDBRefs()); /* * find sequences (except this one!), of complementary type, * which have a dbref to an accession id for this sequence, * and add them to the results */ - found = CrossRef.searchDatasetXrefs(dss, !dna, lrfs, dataset, - rseqs, cf); + found = searchDatasetXrefs(fromDna, dss, lrfs, rseqs, cf); } - for (int r = 0; xrfs != null && r < xrfs.length; r++) + if (xrfs == null && !found) { - DBRefEntry xref = xrfs[r]; - if (source != null && !source.equals(xref.getSource())) - { - continue; - } - if (xref.hasMap()) + /* + * no dbref to source on this sequence or matched + * complementary sequence in the dataset + */ + continue; + } + List sourceRefs = DBRefUtils.searchRefsForSource(xrfs, + source); + Iterator refIterator = sourceRefs.iterator(); + // At this point, if we are retrieving Ensembl, we still don't filter out + // ENST when looking for protein crossrefs. + while (refIterator.hasNext()) + { + DBRefEntry xref = refIterator.next(); + found = false; + // we're only interested in coding cross-references, not + // locus->transcript + if (xref.hasMap() && xref.getMap().getMap().isTripletMap()) { - if (xref.getMap().getTo() != null) + SequenceI mappedTo = xref.getMap().getTo(); + if (mappedTo != null) { + /* + * dbref contains the sequence it maps to; add it to the + * results unless we have done so already (could happen if + * fetching xrefs for sequences which have xrefs in common) + * for example: UNIPROT {P0CE19, P0CE20} -> EMBL {J03321, X06707} + */ found = true; - SequenceI rsq = new Sequence(xref.getMap().getTo()); + /* + * problem: matcher.findIdMatch() is lenient - returns a sequence + * with a dbref to the search arg e.g. ENST for ENSP - wrong + * but findInDataset() matches ENSP when looking for Uniprot... + */ + SequenceI matchInDataset = findInDataset(xref); + if (matchInDataset != null && xref.getMap().getTo() != null + && matchInDataset != xref.getMap().getTo()) + { + System.err + .println("Implementation problem (reopen JAL-2154): CrossRef.findInDataset seems to have recovered a different sequence than the one explicitly mapped for xref." + + "Found:" + + matchInDataset + + "\nExpected:" + + xref.getMap().getTo() + + "\nFor xref:" + + xref); + } + /*matcher.findIdMatch(mappedTo);*/ + if (matchInDataset != null) + { + if (!rseqs.contains(matchInDataset)) + { + rseqs.add(matchInDataset); + } + // even if rseqs contained matchInDataset - check mappings between + // these seqs are added + // need to try harder to only add unique mappings + if (xref.getMap().getMap().isTripletMap() + && dataset.getMapping(seq, matchInDataset) == null + && cf.getMappingBetween(seq, matchInDataset) == null) + { + // materialise a mapping for highlighting between these + // sequences + if (fromDna) + { + cf.addMap(dss, matchInDataset, xref.getMap().getMap(), + xref.getMap().getMappedFromId()); + } + else + { + cf.addMap(matchInDataset, dss, xref.getMap().getMap() + .getInverse(), xref.getMap().getMappedFromId()); + } + } + + refIterator.remove(); + continue; + } + // TODO: need to determine if this should be a deriveSequence + SequenceI rsq = new Sequence(mappedTo); rseqs.add(rsq); - if (xref.getMap().getMap().getFromRatio() != xref - .getMap().getMap().getToRatio()) + if (xref.getMap().getMap().isTripletMap()) { // get sense of map correct for adding to product alignment. - if (dna) + if (fromDna) { // map is from dna seq to a protein product - cf.addMap(dss, rsq, xref.getMap().getMap()); + cf.addMap(dss, rsq, xref.getMap().getMap(), xref.getMap() + .getMappedFromId()); } else { // map should be from protein seq to its coding dna - cf.addMap(rsq, dss, xref.getMap().getMap().getInverse()); + cf.addMap(rsq, dss, xref.getMap().getMap().getInverse(), + xref.getMap().getMappedFromId()); } } } } + if (!found) { - // do a bit more work - search for sequences with references matching - // xrefs on this sequence. - if (dataset != null) + SequenceI matchedSeq = matcher.findIdMatch(xref.getSource() + "|" + + xref.getAccessionId()); + // if there was a match, check it's at least the right type of + // molecule! + if (matchedSeq != null && matchedSeq.isProtein() == fromDna) { - found = searchDataset(dss, xref, dataset, rseqs, cf, false,/*true?*/ - !dna); - if (found) + if (constructMapping(seq, matchedSeq, xref, cf, fromDna)) { - xrfs[r] = null; // we've recovered seqs for this one. + found = true; } } } + + if (!found) + { + // do a bit more work - search for sequences with references matching + // xrefs on this sequence. + found = searchDataset(fromDna, dss, xref, rseqs, cf, false); + } + if (found) + { + refIterator.remove(); + } } - if (!found) + + /* + * fetch from source database any dbrefs we haven't resolved up to here + */ + if (!sourceRefs.isEmpty()) { - if (xrfs != null && xrfs.length > 0) + retrieveCrossRef(sourceRefs, seq, xrfs, fromDna, cf); + } + } + + Alignment ral = null; + if (rseqs.size() > 0) + { + ral = new Alignment(rseqs.toArray(new SequenceI[rseqs.size()])); + if (!cf.isEmpty()) + { + dataset.addCodonFrame(cf); + } + } + return ral; + } + + private void retrieveCrossRef(List sourceRefs, SequenceI seq, + DBRefEntry[] xrfs, boolean fromDna, AlignedCodonFrame cf) + { + ASequenceFetcher sftch = SequenceFetcherFactory.getSequenceFetcher(); + SequenceI[] retrieved = null; + SequenceI dss = seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence(); + // first filter in case we are retrieving crossrefs that have already been + // retrieved. this happens for cases where a database record doesn't yield + // protein products for CDS + removeAlreadyRetrievedSeqs(sourceRefs, fromDna); + if (sourceRefs.size() == 0) + { + // no more work to do! We already had all requested sequence records in + // the dataset. + return; + } + try + { + retrieved = sftch.getSequences(sourceRefs, !fromDna); + } catch (Exception e) + { + System.err + .println("Problem whilst retrieving cross references for Sequence : " + + seq.getName()); + e.printStackTrace(); + } + + if (retrieved != null) + { + boolean addedXref = false; + for (SequenceI retrievedSequence : retrieved) + { + // dataset gets contaminated ccwith non-ds sequences. why ??! + // try: Ensembl -> Nuc->Ensembl, Nuc->Uniprot-->Protein->EMBL-> + SequenceI retrievedDss = retrievedSequence.getDatasetSequence() == null ? retrievedSequence + : retrievedSequence.getDatasetSequence(); + addedXref |= importCrossRefSeq(cf, dss, retrievedDss); + } + if (!addedXref) + { + // try again, after looking for matching IDs + // shouldn't need to do this unless the dbref mechanism has broken. + updateDbrefMappings(seq, xrfs, retrieved, cf, fromDna); + for (SequenceI retrievedSequence : retrieved) { - ASequenceFetcher sftch = SequenceFetcherFactory - .getSequenceFetcher(); - SequenceI[] retrieved = null; - int l = xrfs.length; - for (int r = 0; r < xrfs.length; r++) + // dataset gets contaminated ccwith non-ds sequences. why ??! + // try: Ensembl -> Nuc->Ensembl, Nuc->Uniprot-->Protein->EMBL-> + SequenceI retrievedDss = retrievedSequence.getDatasetSequence() == null ? retrievedSequence + : retrievedSequence.getDatasetSequence(); + addedXref |= importCrossRefSeq(cf, dss, retrievedDss); + } + } + } + } + + /** + * Search dataset for sequences with a primary reference contained in + * sourceRefs. + * + * @param sourceRefs + * - list of references to filter. + * @param fromDna + * - type of sequence to search for matching primary reference. + */ + private void removeAlreadyRetrievedSeqs(List sourceRefs, + boolean fromDna) + { + DBRefEntry[] dbrSourceSet = sourceRefs.toArray(new DBRefEntry[0]); + for (SequenceI sq : dataset.getSequences()) + { + boolean dupeFound = false; + // !fromDna means we are looking only for nucleotide sequences, not + // protein + if (sq.isProtein() == fromDna) + { + for (DBRefEntry dbr : sq.getPrimaryDBRefs()) + { + for (DBRefEntry found : DBRefUtils.searchRefs(dbrSourceSet, dbr)) { - // filter out any irrelevant or irretrievable references - if (xrfs[r] == null - || ((source != null && !source.equals(xrfs[r] - .getSource())) || !sftch.isFetchable(xrfs[r] - .getSource()))) - { - l--; - xrfs[r] = null; - } + sourceRefs.remove(found); + dupeFound = true; } - if (l > 0) + } + } + if (dupeFound) + { + // rebuild the search array from the filtered sourceRefs list + dbrSourceSet = sourceRefs.toArray(new DBRefEntry[0]); + } + } + } + + /** + * process sequence retrieved via a dbref on source sequence to resolve and + * transfer data + * + * @param cf + * @param sourceSequence + * @param retrievedSequence + * @return true if retrieveSequence was imported + */ + private boolean importCrossRefSeq(AlignedCodonFrame cf, + SequenceI sourceSequence, SequenceI retrievedSequence) + { + /** + * set when retrievedSequence has been verified as a crossreference for + * sourceSequence + */ + boolean imported = false; + DBRefEntry[] dbr = retrievedSequence.getDBRefs(); + List newDsSeqs = new ArrayList(); + if (dbr != null) + { + for (DBRefEntry dbref : dbr) + { + SequenceI matched = findInDataset(dbref); + if (matched == sourceSequence) + { + // verified retrieved and source sequence cross-reference each other + imported = true; + } + // find any entry where we should put in the sequence being + // cross-referenced into the map + Mapping map = dbref.getMap(); + if (map != null) + { + if (map.getTo() != null && map.getMap() != null) { - // System.out - // .println("Attempting to retrieve cross referenced sequences."); - DBRefEntry[] t = new DBRefEntry[l]; - l = 0; - for (int r = 0; r < xrfs.length; r++) + if (map.getTo() == sourceSequence) { - if (xrfs[r] != null) - { - t[l++] = xrfs[r]; - } + // already called to import once, and most likely this sequence + // already imported ! + continue; } - xrfs = t; - try + if (matched == null) { - retrieved = sftch.getSequences(xrfs, !dna); - // problem here is we don't know which of xrfs resulted in which - // retrieved element - } catch (Exception e) - { - System.err - .println("Problem whilst retrieving cross references for Sequence : " - + seq.getName()); - e.printStackTrace(); + /* + * sequence is new to dataset, so save a reference so it can be added. + */ + newDsSeqs.add(map.getTo()); + continue; } - if (retrieved != null) - { - updateDbrefMappings(dna, seq, xrfs, retrieved, cf); + /* + * there was a matching sequence in dataset, so now, check to see if we can update the map.getTo() sequence to the existing one. + */ - SequenceIdMatcher matcher = new SequenceIdMatcher( - dataset.getSequences()); - List copiedFeatures = new ArrayList(); - CrossRef me = new CrossRef(); - for (int rs = 0; rs < retrieved.length; rs++) + try + { + // compare ms with dss and replace with dss in mapping + // if map is congruent + SequenceI ms = map.getTo(); + // TODO findInDataset requires exact sequence match but + // 'congruent' test is only for the mapped part + // maybe not a problem in practice since only ENA provide a + // mapping and it is to the full protein translation of CDS + // matcher.findIdMatch(map.getTo()); + // TODO addendum: if matched is shorter than getTo, this will fail + // - when it should really succeed. + int sf = map.getMap().getToLowest(); + int st = map.getMap().getToHighest(); + SequenceI mappedrg = ms.getSubSequence(sf, st); + if (mappedrg.getLength() > 0 + && ms.getSequenceAsString().equals( + matched.getSequenceAsString())) { - // TODO: examine each sequence for 'redundancy' - DBRefEntry[] dbr = retrieved[rs].getDBRefs(); - if (dbr != null && dbr.length > 0) + /* + * sequences were a match, + */ + String msg = "Mapping updated from " + ms.getName() + + " to retrieved crossreference " + + matched.getName(); + System.out.println(msg); + + DBRefEntry[] toRefs = map.getTo().getDBRefs(); + if (toRefs != null) + { + /* + * transfer database refs + */ + for (DBRefEntry ref : toRefs) + { + matched.addDBRef(ref); // add or update mapping + } + } + map.setTo(matched); + + /* + * give the reverse reference the inverse mapping + * (if it doesn't have one already) + */ + setReverseMapping(matched, dbref, cf); + + /* + * copy sequence features as well, avoiding + * duplication (e.g. same variation from two + * transcripts) + */ + SequenceFeature[] sfs = ms.getSequenceFeatures(); + if (sfs != null) { - for (int di = 0; di < dbr.length; di++) + for (SequenceFeature feat : sfs) { - // find any entry where we should put in the sequence being - // cross-referenced into the map - Mapping map = dbr[di].getMap(); - if (map != null) + /* + * make a flyweight feature object which ignores Parent + * attribute in equality test; this avoids creating many + * otherwise duplicate exon features on genomic sequence + */ + SequenceFeature newFeature = new SequenceFeature( + feat) { - if (map.getTo() != null && map.getMap() != null) + @Override + public boolean equals(Object o) { - SequenceI matched = matcher - .findIdMatch(map.getTo()); - if (matched != null) - { - /* - * already got an xref to this sequence; update this - * map to point to the same sequence, and add - * any new dbrefs to it - */ - for (DBRefEntry ref : map.getTo().getDBRefs()) - { - matched.addDBRef(ref); // add or update mapping - } - map.setTo(matched); - } - else - { - matcher.add(map.getTo()); - } - try - { - // compare ms with dss and replace with dss in mapping - // if map is congruent - SequenceI ms = map.getTo(); - int sf = map.getMap().getToLowest(); - int st = map.getMap().getToHighest(); - SequenceI mappedrg = ms.getSubSequence(sf, st); - // SequenceI loc = dss.getSubSequence(sf, st); - if (mappedrg.getLength() > 0 - && ms.getSequenceAsString().equals( - dss.getSequenceAsString())) - // && mappedrg.getSequenceAsString().equals( - // loc.getSequenceAsString())) - { - String msg = "Mapping updated from " - + ms.getName() - + " to retrieved crossreference " - + dss.getName(); - System.out.println(msg); - // method to update all refs of existing To on - // retrieved sequence with dss and merge any props - // on To onto dss. - map.setTo(dss); - /* - * copy sequence features as well, avoiding - * duplication (e.g. same variation from 2 - * transcripts) - */ - SequenceFeature[] sfs = ms - .getSequenceFeatures(); - if (sfs != null) - { - for (SequenceFeature feat : sfs) - { - /* - * we override SequenceFeature.equals here (but - * not elsewhere) to ignore Parent attribute - * TODO not quite working yet! - */ - if (!copiedFeatures - .contains(me.new MySequenceFeature( - feat))) - { - dss.addSequenceFeature(feat); - copiedFeatures.add(feat); - } - } - } - cf.addMap(retrieved[rs].getDatasetSequence(), - dss, map.getMap()); - } - // TODO remove this 'else' and the cf.addMap above? - else - { - cf.addMap(retrieved[rs].getDatasetSequence(), - map.getTo(), map.getMap()); - } - } catch (Exception e) - { - System.err - .println("Exception when consolidating Mapped sequence set..."); - e.printStackTrace(System.err); - } + return super.equals(o, true); } - } + }; + matched.addSequenceFeature(newFeature); } } - retrieved[rs].updatePDBIds(); - rseqs.add(retrieved[rs]); + } + cf.addMap(retrievedSequence, map.getTo(), map.getMap()); + } catch (Exception e) + { + System.err + .println("Exception when consolidating Mapped sequence set..."); + e.printStackTrace(System.err); } } } } } + if (imported) + { + retrievedSequence.updatePDBIds(); + rseqs.add(retrievedSequence); + if (dataset.findIndex(retrievedSequence) == -1) + { + dataset.addSequence(retrievedSequence); + matcher.add(retrievedSequence); + } + for (SequenceI newToSeq : newDsSeqs) + { - Alignment ral = null; - if (rseqs.size() > 0) + if (dataset.findIndex(newToSeq) == -1) + { + dataset.addSequence(newToSeq); + matcher.add(newToSeq); + } + } + } + return imported; + } + /** + * Sets the inverse sequence mapping in the corresponding dbref of the mapped + * to sequence (if any). This is used after fetching a cross-referenced + * sequence, if the fetched sequence has a mapping to the original sequence, + * to set the mapping in the original sequence's dbref. + * + * @param mapFrom + * the sequence mapped from + * @param dbref + * @param mappings + */ + void setReverseMapping(SequenceI mapFrom, DBRefEntry dbref, + AlignedCodonFrame mappings) + { + SequenceI mapTo = dbref.getMap().getTo(); + if (mapTo == null) { - ral = new Alignment(rseqs.toArray(new SequenceI[rseqs.size()])); - if (cf != null && !cf.isEmpty()) + return; + } + DBRefEntry[] dbrefs = mapTo.getDBRefs(); + if (dbrefs == null) + { + return; + } + for (DBRefEntry toRef : dbrefs) + { + if (toRef.hasMap() && mapFrom == toRef.getMap().getTo()) { - ral.addCodonFrame(cf); + /* + * found the reverse dbref; update its mapping if null + */ + if (toRef.getMap().getMap() == null) + { + MapList inverse = dbref.getMap().getMap().getInverse(); + toRef.getMap().setMap(inverse); + mappings.addMap(mapTo, mapFrom, inverse); + } } } - return ral; + } + + /** + * Returns null or the first sequence in the dataset which is identical to + * xref.mapTo, and has a) a primary dbref matching xref, or if none found, the + * first one with an ID source|xrefacc + * + * @param xref + * with map and mapped-to sequence + * @return + */ + SequenceI findInDataset(DBRefEntry xref) + { + if (xref == null || !xref.hasMap() || xref.getMap().getTo() == null) + { + return null; + } + SequenceI mapsTo = xref.getMap().getTo(); + String name = xref.getAccessionId(); + String name2 = xref.getSource() + "|" + name; + SequenceI dss = mapsTo.getDatasetSequence() == null ? mapsTo : mapsTo + .getDatasetSequence(); + // first check ds if ds is directly referenced + if (dataset.findIndex(dss) > -1) + { + return dss; + } + DBRefEntry template = new DBRefEntry(xref.getSource(), null, + xref.getAccessionId()); + /** + * remember the first ID match - in case we don't find a match to template + */ + SequenceI firstIdMatch = null; + for (SequenceI seq : dataset.getSequences()) + { + // first check primary refs. + List match = DBRefUtils.searchRefs(seq.getPrimaryDBRefs() + .toArray(new DBRefEntry[0]), template); + if (match != null && match.size() == 1 && sameSequence(seq, dss)) + { + return seq; + } + /* + * clumsy alternative to using SequenceIdMatcher which currently + * returns sequences with a dbref to the matched accession id + * which we don't want + */ + if (firstIdMatch == null + && (name.equals(seq.getName()) || seq.getName().startsWith( + name2))) + { + if (sameSequence(seq, dss)) + { + firstIdMatch = seq; + } + } + } + return firstIdMatch; + } + + /** + * Answers true if seq1 and seq2 contain exactly the same characters (ignoring + * case), else false. This method compares the lengths, then each character in + * turn, in order to 'fail fast'. For case-sensitive comparison, it would be + * possible to use Arrays.equals(seq1.getSequence(), seq2.getSequence()). + * + * @param seq1 + * @param seq2 + * @return + */ + // TODO move to Sequence / SequenceI + static boolean sameSequence(SequenceI seq1, SequenceI seq2) + { + if (seq1 == seq2) + { + return true; + } + if (seq1 == null || seq2 == null) + { + return false; + } + char[] c1 = seq1.getSequence(); + char[] c2 = seq2.getSequence(); + if (c1.length != c2.length) + { + return false; + } + for (int i = 0; i < c1.length; i++) + { + int diff = c1[i] - c2[i]; + /* + * same char or differ in case only ('a'-'A' == 32) + */ + if (diff != 0 && diff != 32 && diff != -32) + { + return false; + } + } + return true; } /** @@ -444,62 +791,129 @@ public class CrossRef * retrieved sequence if found, and adds any new mappings to the * AlignedCodonFrame * - * @param dna * @param mapFrom * @param xrefs * @param retrieved * @param acf */ - static void updateDbrefMappings(boolean dna, SequenceI mapFrom, - DBRefEntry[] xrefs, SequenceI[] retrieved, AlignedCodonFrame acf) + void updateDbrefMappings(SequenceI mapFrom, DBRefEntry[] xrefs, + SequenceI[] retrieved, AlignedCodonFrame acf, boolean fromDna) { - SequenceIdMatcher matcher = new SequenceIdMatcher(retrieved); + SequenceIdMatcher idMatcher = new SequenceIdMatcher(retrieved); for (DBRefEntry xref : xrefs) { if (!xref.hasMap()) { String targetSeqName = xref.getSource() + "|" + xref.getAccessionId(); - SequenceI[] matches = matcher.findAllIdMatches(targetSeqName); + SequenceI[] matches = idMatcher.findAllIdMatches(targetSeqName); if (matches == null) { return; } for (SequenceI seq : matches) { - MapList mapping = null; - if (dna) - { - mapping = AlignmentUtils.mapCdnaToProtein(seq, mapFrom); - } - else - { - mapping = AlignmentUtils.mapCdnaToProtein(mapFrom, seq); - if (mapping != null) - { - mapping = mapping.getInverse(); - } - } - if (mapping != null) - { - xref.setMap(new Mapping(seq, mapping)); - if (dna) - { - AlignmentUtils.computeProteinFeatures(mapFrom, seq, mapping); - } - if (dna) - { - acf.addMap(mapFrom, seq, mapping); - } - else - { - acf.addMap(seq, mapFrom, mapping.getInverse()); - } - continue; - } + constructMapping(mapFrom, seq, xref, acf, fromDna); + } + } + } + } + + /** + * Tries to make a mapping between sequences. If successful, adds the mapping + * to the dbref and the mappings collection and answers true, otherwise + * answers false. The following methods of making are mapping are tried in + * turn: + *
    + *
  • if 'mapTo' holds a mapping to 'mapFrom', take the inverse; this is, for + * example, the case after fetching EMBL cross-references for a Uniprot + * sequence
    • + *
  • else check if the dna translates exactly to the protein (give or take + * start and stop codons>
    • + *
  • else try to map based on CDS features on the dna sequence
    • + *
+ * + * @param mapFrom + * @param mapTo + * @param xref + * @param mappings + * @return + */ + boolean constructMapping(SequenceI mapFrom, SequenceI mapTo, + DBRefEntry xref, AlignedCodonFrame mappings, boolean fromDna) + { + MapList mapping = null; + SequenceI dsmapFrom = mapFrom.getDatasetSequence() == null ? mapFrom + : mapFrom.getDatasetSequence(); + SequenceI dsmapTo = mapTo.getDatasetSequence() == null ? mapTo + : mapTo.getDatasetSequence(); + /* + * look for a reverse mapping, if found make its inverse. + * Note - we do this on dataset sequences only. + */ + if (dsmapTo.getDBRefs() != null) + { + for (DBRefEntry dbref : dsmapTo.getDBRefs()) + { + String name = dbref.getSource() + "|" + dbref.getAccessionId(); + if (dbref.hasMap() && dsmapFrom.getName().startsWith(name)) + { + /* + * looks like we've found a map from 'mapTo' to 'mapFrom' + * - invert it to make the mapping the other way + */ + MapList reverse = dbref.getMap().getMap().getInverse(); + xref.setMap(new Mapping(dsmapTo, reverse)); + mappings.addMap(mapFrom, dsmapTo, reverse); + return true; } } } + + if (fromDna) + { + mapping = AlignmentUtils.mapCdnaToProtein(mapTo, mapFrom); + } + else + { + mapping = AlignmentUtils.mapCdnaToProtein(mapFrom, mapTo); + if (mapping != null) + { + mapping = mapping.getInverse(); + } + } + if (mapping == null) + { + return false; + } + xref.setMap(new Mapping(mapTo, mapping)); + + /* + * and add a reverse DbRef with the inverse mapping + */ + if (mapFrom.getDatasetSequence() != null && false) + // && mapFrom.getDatasetSequence().getSourceDBRef() != null) + { + // possible need to search primary references... except, why doesn't xref + // == getSourceDBRef ?? + // DBRefEntry dbref = new DBRefEntry(mapFrom.getDatasetSequence() + // .getSourceDBRef()); + // dbref.setMap(new Mapping(mapFrom.getDatasetSequence(), mapping + // .getInverse())); + // mapTo.addDBRef(dbref); + } + + if (fromDna) + { + AlignmentUtils.computeProteinFeatures(mapFrom, mapTo, mapping); + mappings.addMap(mapFrom, mapTo, mapping); + } + else + { + mappings.addMap(mapTo, mapFrom, mapping.getInverse()); + } + + return true; } /** @@ -507,15 +921,16 @@ public class CrossRef * dataset (that is not equal to sequenceI) Identifies matching DBRefEntry * based on source and accession string only - Map and Version are nulled. * + * @param fromDna + * - true if context was searching from Dna sequences, false if + * context was searching from Protein sequences * @param sequenceI * @param lrfs - * @param dataset - * @param rseqs + * @param foundSeqs * @return true if matches were found. */ - private static boolean searchDatasetXrefs(SequenceI sequenceI, - boolean dna, DBRefEntry[] lrfs, AlignmentI dataset, - List rseqs, AlignedCodonFrame cf) + private boolean searchDatasetXrefs(boolean fromDna, SequenceI sequenceI, + DBRefEntry[] lrfs, List foundSeqs, AlignedCodonFrame cf) { boolean found = false; if (lrfs == null) @@ -528,8 +943,7 @@ public class CrossRef // add in wildcards xref.setVersion(null); xref.setMap(null); - found |= searchDataset(sequenceI, xref, dataset, rseqs, cf, false, - dna); + found |= searchDataset(fromDna, sequenceI, xref, foundSeqs, cf, false); } return found; } @@ -538,25 +952,33 @@ public class CrossRef * Searches dataset for DBRefEntrys matching the given one (xrf) and adds the * associated sequence to rseqs * - * @param sequenceI + * @param fromDna + * true if context was searching for refs *from* dna sequence, false + * if context was searching for refs *from* protein sequence + * @param fromSeq * a sequence to ignore (start point of search) * @param xrf * a cross-reference to try to match - * @param dataset - * sequences to search in - * @param rseqs + * @param foundSeqs * result list to add to - * @param cf + * @param mappings * a set of sequence mappings to add to * @param direct - * - search all references or only subset - * @param dna - * search dna or protein xrefs (if direct=false) + * - indicates the type of relationship between returned sequences, + * xrf, and sequenceI that is required. + *
    + *
  • direct implies xrf is a primary reference for sequenceI AND + * the sequences to be located (eg a uniprot ID for a protein + * sequence, and a uniprot ref on a transcript sequence).
    • + *
  • indirect means xrf is a cross reference with respect to + * sequenceI or all the returned sequences (eg a genomic reference + * associated with a locus and one or more transcripts)
    • + *
* @return true if relationship found and sequence added. */ - public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf, - AlignmentI dataset, List rseqs, AlignedCodonFrame cf, - boolean direct, boolean dna) + boolean searchDataset(boolean fromDna, SequenceI fromSeq, DBRefEntry xrf, + List foundSeqs, AlignedCodonFrame mappings, + boolean direct) { boolean found = false; if (dataset == null) @@ -578,22 +1000,23 @@ public class CrossRef if (nxt.getDatasetSequence() != null) { System.err - .println("Implementation warning: getProducts passed a dataset alignment without dataset sequences in it!"); + .println("Implementation warning: CrossRef initialised with a dataset alignment with non-dataset sequences in it! (" + + nxt.getDisplayId(true) + + " has ds reference " + + nxt.getDatasetSequence().getDisplayId(true) + + ")"); } - if (nxt == sequenceI || nxt == sequenceI.getDatasetSequence()) + if (nxt == fromSeq || nxt == fromSeq.getDatasetSequence()) { continue; } - // check if this is the correct sequence type + /* + * only look at same molecule type if 'direct', or + * complementary type if !direct + */ { - // TODO 'direct' is always set to false - remove? - // or should it be 'true' from findXrefSequences? - // also its Javadoc conflicts with its use: - // test below implies 'direct' means find complementary sequences, - // !direct means select same molecule type - boolean isDna = Comparison - .isNucleotide(new SequenceI[] { nxt }); - if ((direct && isDna == dna) || (!direct && isDna != dna)) + boolean isDna = !nxt.isProtein(); + if (direct ? (isDna != fromDna) : (isDna == fromDna)) { // skip this sequence because it is wrong molecule type continue; @@ -602,55 +1025,55 @@ public class CrossRef // look for direct or indirect references in common DBRefEntry[] poss = nxt.getDBRefs(); - DBRefEntry[] cands = null; - /* - * TODO does this make any sense? - * if 'direct', search the dbrefs for xrf - * else, filter the dbrefs by type and then search for xrf - * - the result is the same isn't it? - */ - if (direct) - { - cands = DBRefUtils.searchRefs(poss, xrf); - } - else - { - poss = DBRefUtils.selectDbRefs(!dna, poss); - cands = DBRefUtils.searchRefs(poss, xrf); - } - if (cands != null) + List cands = null; + + // todo: indirect specifies we select either direct references to nxt + // that match xrf which is indirect to sequenceI, or indirect + // references to nxt that match xrf which is direct to sequenceI + cands = DBRefUtils.searchRefs(poss, xrf); + // else + // { + // poss = DBRefUtils.selectDbRefs(nxt.isProtein()!fromDna, poss); + // cands = DBRefUtils.searchRefs(poss, xrf); + // } + if (!cands.isEmpty()) { - if (!rseqs.contains(nxt)) + if (foundSeqs.contains(nxt)) { - found = true; - rseqs.add(nxt); - boolean foundmap = cf != null; - // don't search if we aren't given a codon map object - for (int r = 0; foundmap && r < cands.length; r++) + continue; + } + found = true; + foundSeqs.add(nxt); + if (mappings != null && !direct) + { + /* + * if the matched sequence has mapped dbrefs to + * protein product / cdna, add equivalent mappings to + * our source sequence + */ + for (DBRefEntry candidate : cands) { - if (cands[r].hasMap()) + Mapping mapping = candidate.getMap(); + if (mapping != null) { - Mapping mapping = cands[r].getMap(); MapList map = mapping.getMap(); if (mapping.getTo() != null && map.getFromRatio() != map.getToRatio()) { - // get sense of map correct for adding to product - // alignment. - if (dna) + /* + * add a mapping, as from dna to peptide sequence + */ + if (map.getFromRatio() == 3) { - // map is from dna seq to a protein product - cf.addMap(sequenceI, nxt, map); + mappings.addMap(nxt, fromSeq, map); } else { - // map should be from protein seq to its coding dna - cf.addMap(nxt, sequenceI, map.getInverse()); + mappings.addMap(nxt, fromSeq, map.getInverse()); } } } } - // TODO: add mapping between sequences if necessary } } }