src/jalview/analysis/CrossRefs.java

   1 package jalview.analysis;
   2
   3 import jalview.analysis.CrossRef.MySequenceFeature;
   4 import jalview.datamodel.AlignedCodonFrame;
   5 import jalview.datamodel.Alignment;
   6 import jalview.datamodel.AlignmentI;
   7 import jalview.datamodel.DBRefEntry;
   8 import jalview.datamodel.Mapping;
   9 import jalview.datamodel.Sequence;
  10 import jalview.datamodel.SequenceFeature;
  11 import jalview.datamodel.SequenceI;
  12 import jalview.util.Comparison;
  13 import jalview.util.DBRefUtils;
  14 import jalview.util.MapList;
  15 import jalview.ws.SequenceFetcherFactory;
  16 import jalview.ws.seqfetcher.ASequenceFetcher;
  17
  18 import java.util.ArrayList;
  19 import java.util.Iterator;
  20 import java.util.List;
  21
  22 public class CrossRefs
  23 {
  24   /*
  25    * A sub-class that ignores Parent attribute when comparing sequence
  26    * features. This avoids 'duplicate' CDS features that only
  27    * differ in their parent Transcript ids.
  28    */
  29   class MySequenceFeature extends SequenceFeature
  30   {
  31     private SequenceFeature feat;
  32
  33     MySequenceFeature(SequenceFeature sf)
  34     {
  35       this.feat = sf;
  36     }
  37
  38     @Override
  39     public boolean equals(Object o)
  40     {
  41       return feat.equals(o, true);
  42     }
  43   }
  44
  45   /**
  46    * Finds cross-references for sequences from a specified source database.
  47    * These may be found in four ways:
  48    * <ul>
  49    * <li>as a DBRefEntry on the known sequence, which has a mapped-to sequence</li>
  50    * <li>a sequence of complementary type in the alignment dataset, which has a
  51    * DBRefEntry to one of the known sequence's 'direct' DBRefs</li>
  52    * <li>a sequence of complementary type in the alignment, which has a
  53    * DBRefEntry to one of the known sequence's 'cross-ref' DBRefs</li>
  54    * <li>by fetching the accession from the remote database</li>
  55    * </ul>
  56    *
  57    * @param seqs
  58    *          the sequences whose cross-references we are searching for
  59    * @param dna
  60    *          true if the sequences are from a nucleotide alignment, else false
  61    * @param source
  62    *          the database source we want cross-references to
  63    * @param dataset
  64    *          the alignment dataset the sequences belong to
  65    * @return an alignment containing cross-reference sequences, or null if none
  66    *         found
  67    */
  68   public static AlignmentI findXrefSequences(SequenceI[] seqs, boolean dna,
  69           String source, AlignmentI dataset)
  70   {
  71     /*
  72      * filter to only those sequences of the right type (nucleotide/protein)
  73      */
  74     List<SequenceI> fromSeqs = new ArrayList<SequenceI>();
  75     for (SequenceI seq : seqs)
  76     {
  77       if (dna == Comparison.isNucleotide(seq))
  78       {
  79         fromSeqs.add(seq);
  80       }
  81     }
  82     return findXrefSequences(fromSeqs, dna, source, dataset);
  83   }
  84
  85   /**
  86    * Finds cross-references for sequences from a specified source database.
  87    * These may be found in four ways:
  88    * <ul>
  89    * <li>as a DBRefEntry on the known sequence, which has a mapped-to sequence</li>
  90    * <li>a sequence of complementary type in the alignment dataset, which has a
  91    * DBRefEntry to one of the known sequence's 'direct' DBRefs</li>
  92    * <li>a sequence of complementary type in the alignment, which has a
  93    * DBRefEntry to one of the known sequence's 'cross-ref' DBRefs</li>
  94    * <li>by fetching the accession from the remote database</li>
  95    * </ul>
  96    *
  97    * @param seqs
  98    *          the sequences whose cross-references we are searching for,
  99    *          filtered to only those which are of the type denoted by 'dna'
 100    * @param dna
 101    *          true if the sequences are from a nucleotide alignment, else false
 102    * @param source
 103    *          the database source we want cross-references to
 104    * @param dataset
 105    *          the alignment dataset the sequences belong to
 106    * @return an alignment containing cross-reference sequences, or null if none
 107    *         found
 108    */
 109   static AlignmentI findXrefSequences(List<SequenceI> fromSeqs,
 110           boolean dna, String source, AlignmentI dataset)
 111   {
 112     List<SequenceI> foundSeqs = new ArrayList<SequenceI>();
 113     AlignedCodonFrame mappings = new AlignedCodonFrame();
 114
 115     List<DBRefEntry> unresolvedRefs = new ArrayList<DBRefEntry>();
 116
 117     /*
 118      * first extract any mapped sequences from sourceRefs
 119      * if successful, sequence is removed from fromSeqs
 120      * if unsuccessful, dbrefs are added to unresolvedRefs
 121      */
 122     findMappedDbrefs(fromSeqs, source, foundSeqs,
 123             unresolvedRefs, mappings);
 124
 125     /*
 126      * then search the alignment dataset for dbref resolutions
 127      */
 128     findIndirectCrossReferences(fromSeqs, source, dataset, foundSeqs,
 129             unresolvedRefs, mappings);
 130
 131     /*
 132      * fetch any remaining sourceRefs from the source database
 133      */
 134     fetchCrossReferences(fromSeqs, unresolvedRefs, foundSeqs, mappings,
 135             dna, dataset);
 136
 137     if (foundSeqs.isEmpty())
 138     {
 139       return null;
 140     }
 141     AlignmentI crossRefs = new Alignment(
 142             foundSeqs.toArray(new SequenceI[foundSeqs.size()]));
 143     crossRefs.addCodonFrame(mappings);
 144     return crossRefs;
 145   }
 146
 147   /**
 148    * Looks for DBRefEntrys to 'source' which have a mapping to a sequence. If
 149    * found, adds the sequence to foundSeqs and removes the dbref from the list.
 150    * DBRefs with no mapping are added to the 'unresolvedRefs' list (setting
 151    * version number to 0 i.e. use source and accession only).
 152    *
 153    * @param fromSeqs
 154    *          the dataset sequences we are searching from
 155    * @param source
 156    *          the database source we are searching dbrefs for
 157    * @param foundSeqs
 158    *          a list of found sequences to add to
 159    * @param unresolvedRefs
 160    *          a list of unresolved cross-references to add to
 161    * @param mappings
 162    *          a set of sequence mappings to add to
 163    * @return
 164    */
 165   static void findMappedDbrefs(List<SequenceI> fromSeqs, String source,
 166           List<SequenceI> foundSeqs, List<DBRefEntry> unresolvedRefs,
 167           AlignedCodonFrame mappings)
 168   {
 169     Iterator<SequenceI> it = fromSeqs.iterator();
 170     while (it.hasNext())
 171     {
 172       SequenceI seq = it.next();
 173       SequenceI dss = seq.getDatasetSequence();
 174       dss = dss == null ? seq : dss;
 175
 176       DBRefEntry[] dbRefs = seq.getDBRefs();
 177       if (dbRefs == null)
 178       {
 179         continue;
 180       }
 181       boolean resolved = false;
 182       for (DBRefEntry dbref : dbRefs)
 183       {
 184         if (!source.equals(dbref.getSource()))
 185         {
 186           continue;
 187         }
 188         DBRefEntry todo = new DBRefEntry(dbref.getSource(), "0",
 189                 dbref.getAccessionId());
 190         Mapping map = dbref.getMap();
 191         if (map != null)
 192         {
 193           unresolvedRefs.remove(todo);
 194           resolved = true;
 195           SequenceI mappedTo = map.getTo();
 196           if (mappedTo != null)
 197           {
 198             foundSeqs.add(new Sequence(mappedTo));
 199
 200             /*
 201              * check mapping is not 'direct' (it shouldn't be if we reach here)
 202              * and add mapping (dna-to-peptide or vice versa) to the set
 203              */
 204             MapList mapList = map.getMap();
 205             int fromRatio = mapList.getFromRatio();
 206             int toRatio = mapList.getToRatio();
 207             if (fromRatio != toRatio)
 208             {
 209               if (fromRatio == 3)
 210               {
 211                 mappings.addMap(dss, mappedTo, mapList);
 212               }
 213               else
 214               {
 215                 mappings.addMap(mappedTo, dss, mapList.getInverse());
 216               }
 217             }
 218           }
 219         }
 220         else
 221         {
 222           /*
 223            * no mapping to resolve dbref - add source+accession to list to resolve
 224            */
 225           if (!unresolvedRefs.contains(todo))
 226           {
 227             unresolvedRefs.add(todo);
 228           }
 229         }
 230       }
 231       if (resolved)
 232       {
 233         it.remove();
 234       }
 235     }
 236   }
 237
 238   /**
 239    * Tries to fetch seq's database references to 'source' database, and add them
 240    * to the foundSeqs list. If found, tries to make a mapping between seq and
 241    * the retrieved sequence and insert it into the database reference.
 242    *
 243    * @param fromSeqs
 244    * @param sourceRefs
 245    * @param foundSeqs
 246    * @param mappings
 247    * @param dna
 248    */
 249   static void fetchCrossReferences(List<SequenceI> fromSeqs,
 250           List<DBRefEntry> sourceRefs, List<SequenceI> foundSeqs,
 251           AlignedCodonFrame mappings, boolean dna, AlignmentI dataset)
 252   {
 253     ASequenceFetcher sftch = SequenceFetcherFactory.getSequenceFetcher();
 254     SequenceI[] retrieved;
 255     try
 256     {
 257       retrieved = sftch.getSequences(sourceRefs, !dna);
 258     } catch (Exception e)
 259     {
 260       System.err.println("Problem whilst retrieving cross references: "
 261               + e.getMessage());
 262       e.printStackTrace();
 263       return;
 264     }
 265
 266     if (retrieved == null)
 267     {
 268       return;
 269     }
 270     updateDbrefMappings(dna, fromSeqs, sourceRefs, retrieved, mappings);
 271
 272     SequenceIdMatcher matcher = new SequenceIdMatcher(
 273             dataset.getSequences());
 274     List<SequenceFeature> copiedFeatures = new ArrayList<SequenceFeature>();
 275     CrossRefs me = new CrossRefs();
 276     for (int rs = 0; rs < retrieved.length; rs++)
 277     {
 278       // TODO: examine each sequence for 'redundancy'
 279       DBRefEntry[] dbr = retrieved[rs].getDBRefs();
 280       if (dbr != null && dbr.length > 0)
 281       {
 282         for (int di = 0; di < dbr.length; di++)
 283         {
 284           // find any entry where we should put in the sequence being
 285           // cross-referenced into the map
 286           Mapping map = dbr[di].getMap();
 287           if (map != null)
 288           {
 289             if (map.getTo() != null && map.getMap() != null)
 290             {
 291               SequenceI matched = matcher.findIdMatch(map.getTo());
 292               if (matched != null)
 293               {
 294                 /*
 295                  * already got an xref to this sequence; update this
 296                  * map to point to the same sequence, and add
 297                  * any new dbrefs to it
 298                  */
 299                 for (DBRefEntry ref : map.getTo().getDBRefs())
 300                 {
 301                   matched.addDBRef(ref); // add or update mapping
 302                 }
 303                 map.setTo(matched);
 304               }
 305               else
 306               {
 307                 matcher.add(map.getTo());
 308               }
 309               try
 310               {
 311                 // compare ms with dss and replace with dss in mapping
 312                 // if map is congruent
 313                 SequenceI ms = map.getTo();
 314                 int sf = map.getMap().getToLowest();
 315                 int st = map.getMap().getToHighest();
 316                 SequenceI mappedrg = ms.getSubSequence(sf, st);
 317                 // SequenceI loc = dss.getSubSequence(sf, st);
 318                 if (mappedrg.getLength() > 0
 319                         && ms.getSequenceAsString().equals(
 320                                 fromSeqs.getSequenceAsString()))
 321                 // && mappedrg.getSequenceAsString().equals(
 322                 // loc.getSequenceAsString()))
 323                 {
 324                   String msg = "Mapping updated from " + ms.getName()
 325                           + " to retrieved crossreference "
 326                           + fromSeqs.getName();
 327                   System.out.println(msg);
 328                   // method to update all refs of existing To on
 329                   // retrieved sequence with dss and merge any props
 330                   // on To onto dss.
 331                   map.setTo(fromSeqs);
 332                   /*
 333                    * copy sequence features as well, avoiding
 334                    * duplication (e.g. same variation from 2
 335                    * transcripts)
 336                    */
 337                   SequenceFeature[] sfs = ms.getSequenceFeatures();
 338                   if (sfs != null)
 339                   {
 340                     for (SequenceFeature feat : sfs)
 341                     {
 342                       /*
 343                        * we override SequenceFeature.equals here (but
 344                        * not elsewhere) to ignore Parent attribute
 345                        * TODO not quite working yet!
 346                        */
 347                       if (!copiedFeatures
 348                               .contains(me.new MySequenceFeature(feat)))
 349                       {
 350                         fromSeqs.addSequenceFeature(feat);
 351                         copiedFeatures.add(feat);
 352                       }
 353                     }
 354                   }
 355                 }
 356                 mappings.addMap(retrieved[rs].getDatasetSequence(),
 357                         map.getTo(), map.getMap());
 358               } catch (Exception e)
 359               {
 360                 System.err
 361                         .println("Exception when consolidating Mapped sequence set...");
 362                 e.printStackTrace(System.err);
 363               }
 364             }
 365           }
 366         }
 367       }
 368       retrieved[rs].updatePDBIds();
 369       foundSeqs.add(retrieved[rs]);
 370     }
 371   }
 372
 373   /**
 374    * Searches the alignment for a sequence of complementary type to 'seq' which
 375    * shares a DBRefEntry with it. If found, adds the sequence to foundSeqs and
 376    * removes the resolved sourceRef from the search list.
 377    *
 378    * @param fromSeqs
 379    * @param source
 380    * @param unresolvedRefs
 381    * @param foundSeqs
 382    * @param unresolvedRefs
 383    * @param mappings
 384    * @return
 385    */
 386   static void findIndirectCrossReferences(List<SequenceI> fromSeqs,
 387           String source, AlignmentI dataset,
 388           List<SequenceI> foundSeqs, List<DBRefEntry> unresolvedRefs,
 389           AlignedCodonFrame mappings)
 390   {
 391     Iterator<DBRefEntry> refs = unresolvedRefs.iterator();
 392     while (refs.hasNext())
 393     {
 394       DBRefEntry dbref = refs.next();
 395       boolean found = false;
 396       // boolean found = searchDatasetForCrossReference(fromSeqs, dbref,
 397       // foundSeqs,
 398       // unresolvedRefs, mappings);
 399       if (found)
 400       {
 401         refs.remove();
 402       }
 403     }
 404   }
 405
 406   /**
 407    * Searches the dataset for a sequence of opposite type to 'excluding', which
 408    * has a cross-reference matching dbref. If found, adds the sequence to
 409    * foundSeqs and removes dbref from the search list.
 410    *
 411    * @param excluding
 412    *          a sequence to ignore (start point of search)
 413    * @param dbref
 414    *          a cross-reference to try to match
 415    * @param dataset
 416    *          sequences to search in
 417    * @param foundSeqs
 418    *          result list to add to
 419    * @param mappings
 420    *          a set of sequence mappings to add to
 421    * @return true if relationship found and sequence added
 422    */
 423   static boolean searchDatasetForCrossReference(SequenceI excluding,
 424           DBRefEntry dbref, AlignmentI dataset, List<SequenceI> foundSeqs,
 425           AlignedCodonFrame mappings)
 426   {
 427     boolean fromNucleotide = Comparison.isNucleotide(excluding);
 428     boolean found = false;
 429     if (dataset == null)
 430     {
 431       return false;
 432     }
 433     if (dataset.getSequences() == null)
 434     {
 435       return false;
 436     }
 437     List<SequenceI> ds;
 438     synchronized (ds = dataset.getSequences())
 439     {
 440       for (SequenceI nxt : ds)
 441       {
 442         if (nxt != null)
 443         {
 444           if (nxt.getDatasetSequence() != null)
 445           {
 446             System.err
 447                     .println("Implementation warning: getProducts passed a dataset alignment without dataset sequences in it!");
 448           }
 449           if (nxt == excluding || nxt == excluding.getDatasetSequence())
 450           {
 451             continue;
 452           }
 453           if (foundSeqs.contains(nxt))
 454           {
 455             /*
 456              * already added this sequence to cross-refs
 457              */
 458             continue;
 459           }
 460           boolean isDna = Comparison.isNucleotide(nxt);
 461           if (isDna == fromNucleotide)
 462           {
 463             /*
 464              * skip this sequence - wrong molecule type
 465              */
 466             continue;
 467           }
 468
 469           /*
 470            * check if this sequence has any dbref matching source and accession
 471            * (version and mapping may differ)
 472            */
 473           List<DBRefEntry> candidates = DBRefUtils.searchRefs(
 474                   nxt.getDBRefs(), dbref);
 475
 476           if (candidates.isEmpty())
 477           {
 478             continue;
 479           }
 480           found = true;
 481           foundSeqs.add(nxt);
 482           if (mappings != null)
 483           {
 484             // don't search if we aren't given a codon map object
 485             for (DBRefEntry candidate : candidates)
 486             {
 487               if (candidate.hasMap())
 488               {
 489                 Mapping mapping = candidate.getMap();
 490                 MapList map = mapping.getMap();
 491                 if (mapping.getTo() != null
 492                         && map.getFromRatio() != map.getToRatio())
 493                 {
 494                   if (fromNucleotide)
 495                   {
 496                     // map is from dna seq to a protein product
 497                     mappings.addMap(excluding, nxt, map);
 498                   }
 499                   else
 500                   {
 501                     // map is from protein seq to its coding dna
 502                     mappings.addMap(nxt, excluding, map.getInverse());
 503                   }
 504                 }
 505               }
 506             }
 507           }
 508         }
 509       }
 510     }
 511     return found;
 512   }
 513
 514   /**
 515    * Updates any empty mappings in the cross-references with one to a compatible
 516    * retrieved sequence if found, and adds any new mappings to the
 517    * AlignedCodonFrame
 518    *
 519    * @param dna
 520    * @param fromSeqs
 521    * @param xrefs
 522    * @param retrieved
 523    * @param mappings
 524    */
 525   static void updateDbrefMappings(boolean dna, List<SequenceI> fromSeqs,
 526           List<DBRefEntry> xrefs, SequenceI[] retrieved,
 527           AlignedCodonFrame mappings)
 528   {
 529     SequenceIdMatcher matcher = new SequenceIdMatcher(retrieved);
 530     for (DBRefEntry xref : xrefs)
 531     {
 532       if (!xref.hasMap())
 533       {
 534         String targetSeqName = xref.getSource() + "|"
 535                 + xref.getAccessionId();
 536         SequenceI[] matches = matcher.findAllIdMatches(targetSeqName);
 537         if (matches == null)
 538         {
 539           return;
 540         }
 541         for (SequenceI seq : matches)
 542         {
 543           MapList mapping = null;
 544           if (dna)
 545           {
 546             mapping = AlignmentUtils.mapCdnaToProtein(seq, fromSeqs);
 547           }
 548           else
 549           {
 550             mapping = AlignmentUtils.mapCdnaToProtein(fromSeqs, seq);
 551             if (mapping != null)
 552             {
 553               mapping = mapping.getInverse();
 554             }
 555           }
 556           if (mapping != null)
 557           {
 558             xref.setMap(new Mapping(seq, mapping));
 559             if (dna)
 560             {
 561               AlignmentUtils.computeProteinFeatures(fromSeqs, seq, mapping);
 562             }
 563             if (dna)
 564             {
 565               mappings.addMap(fromSeqs, seq, mapping);
 566             }
 567             else
 568             {
 569               mappings.addMap(seq, fromSeqs, mapping.getInverse());
 570             }
 571             continue;
 572           }
 573         }
 574       }
 575     }
 576   }
 577 }