src/jalview/analysis/CrossRef.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import jalview.datamodel.AlignedCodonFrame;
  24 import jalview.datamodel.Alignment;
  25 import jalview.datamodel.AlignmentI;
  26 import jalview.datamodel.DBRefEntry;
  27 import jalview.datamodel.DBRefSource;
  28 import jalview.datamodel.Mapping;
  29 import jalview.datamodel.Sequence;
  30 import jalview.datamodel.SequenceI;
  31 import jalview.util.DBRefUtils;
  32 import jalview.ws.SequenceFetcher;
  33 import jalview.ws.seqfetcher.ASequenceFetcher;
  34
  35 import java.util.ArrayList;
  36 import java.util.List;
  37 import java.util.Vector;
  38
  39 /**
  40  * Functions for cross-referencing sequence databases. user must first specify
  41  * if cross-referencing from protein or dna (set dna==true)
  42  *
  43  * @author JimP
  44  *
  45  */
  46 public class CrossRef
  47 {
  48   /**
  49    * Select just the DNA or protein references for a protein or dna sequence
  50    *
  51    * @param fromDna
  52    *          if true, select references from DNA (i.e. Protein databases), else
  53    *          DNA database references
  54    * @param refs
  55    *          a set of references to select from
  56    * @return
  57    */
  58   public static DBRefEntry[] findXDbRefs(boolean fromDna, DBRefEntry[] refs)
  59   {
  60     return DBRefUtils.selectRefs(refs, fromDna ? DBRefSource.PROTEINDBS
  61             : DBRefSource.DNACODINGDBS);
  62     // could attempt to find other cross
  63     // refs here - ie PDB xrefs
  64     // (not dna, not protein seq)
  65   }
  66
  67   /**
  68    * @param dna
  69    *          true if seqs are DNA seqs
  70    * @param seqs
  71    * @return a list of sequence database cross reference source types
  72    */
  73   public static String[] findSequenceXrefTypes(boolean dna, SequenceI[] seqs)
  74   {
  75     return findSequenceXrefTypes(dna, seqs, null);
  76   }
  77
  78   /**
  79    * Indirect references are references from other sequences from the dataset to
  80    * any of the direct DBRefEntrys on the given sequences.
  81    *
  82    * @param dna
  83    *          true if seqs are DNA seqs
  84    * @param seqs
  85    * @return a list of sequence database cross reference source types
  86    */
  87   public static String[] findSequenceXrefTypes(boolean dna,
  88           SequenceI[] seqs, AlignmentI dataset)
  89   {
  90     String[] dbrefs = null;
  91     List<String> refs = new ArrayList<String>();
  92     for (SequenceI seq : seqs)
  93     {
  94       if (seq != null)
  95       {
  96         SequenceI dss = seq;
  97         while (dss.getDatasetSequence() != null)
  98         {
  99           dss = dss.getDatasetSequence();
 100         }
 101         DBRefEntry[] rfs = findXDbRefs(dna, dss.getDBRef());
 102         if (rfs != null)
 103         {
 104           for (DBRefEntry ref : rfs)
 105           {
 106             if (!refs.contains(ref.getSource()))
 107             {
 108               refs.add(ref.getSource());
 109             }
 110           }
 111         }
 112         if (dataset != null)
 113         {
 114           // search for references to this sequence's direct references.
 115           DBRefEntry[] lrfs = CrossRef.findXDbRefs(!dna, seq.getDBRef());
 116           List<SequenceI> rseqs = new ArrayList<SequenceI>();
 117           CrossRef.searchDatasetXrefs(seq, !dna, lrfs, dataset, rseqs,
 118                   null); // don't need to specify codon frame for mapping here
 119           for (SequenceI rs : rseqs)
 120           {
 121             DBRefEntry[] xrs = findXDbRefs(dna, rs.getDBRef());
 122             if (xrs != null)
 123             {
 124               for (DBRefEntry ref : xrs)
 125               {
 126                 if (!refs.contains(ref.getSource()))
 127                 {
 128                   refs.add(ref.getSource());
 129                 }
 130               }
 131             }
 132             // looks like copy and paste - change rfs to xrs?
 133             // for (int r = 0; rfs != null && r < rfs.length; r++)
 134             // {
 135             // if (!refs.contains(rfs[r].getSource()))
 136             // {
 137             // refs.add(rfs[r].getSource());
 138             // }
 139             // }
 140           }
 141         }
 142       }
 143     }
 144     if (refs.size() > 0)
 145     {
 146       dbrefs = new String[refs.size()];
 147       refs.toArray(dbrefs);
 148     }
 149     return dbrefs;
 150   }
 151
 152   public static boolean hasCdnaMap(SequenceI[] seqs)
 153   {
 154     // TODO unused - remove?
 155     String[] reftypes = findSequenceXrefTypes(false, seqs);
 156     for (int s = 0; s < reftypes.length; s++)
 157     {
 158       if (reftypes.equals(DBRefSource.EMBLCDS))
 159       {
 160         return true;
 161         // no map
 162       }
 163     }
 164     return false;
 165   }
 166
 167   public static SequenceI[] getCdnaMap(SequenceI[] seqs)
 168   {
 169     // TODO unused - remove?
 170     Vector cseqs = new Vector();
 171     for (int s = 0; s < seqs.length; s++)
 172     {
 173       DBRefEntry[] cdna = findXDbRefs(true, seqs[s].getDBRef());
 174       for (int c = 0; c < cdna.length; c++)
 175       {
 176         if (cdna[c].getSource().equals(DBRefSource.EMBLCDS))
 177         {
 178           System.err
 179                   .println("TODO: unimplemented sequence retrieval for coding region sequence.");
 180           // TODO: retrieve CDS dataset sequences
 181           // need global dataset sequence retriever/resolver to reuse refs
 182           // and construct Mapping entry.
 183           // insert gaps in CDS according to peptide gaps.
 184           // add gapped sequence to cseqs
 185         }
 186       }
 187     }
 188     if (cseqs.size() > 0)
 189     {
 190       SequenceI[] rsqs = new SequenceI[cseqs.size()];
 191       cseqs.copyInto(rsqs);
 192       return rsqs;
 193     }
 194     return null;
 195
 196   }
 197
 198   /**
 199    *
 200    * @param dna
 201    * @param seqs
 202    * @return
 203    */
 204   public static Alignment findXrefSequences(SequenceI[] seqs, boolean dna,
 205           String source)
 206   {
 207     return findXrefSequences(seqs, dna, source, null);
 208   }
 209
 210   /**
 211    *
 212    * @param seqs
 213    *          sequences whose xrefs are being retrieved
 214    * @param dna
 215    *          true if sequences are nucleotide
 216    * @param source
 217    * @param dataset
 218    *          alignment to search for product sequences.
 219    * @return products (as dataset sequences)
 220    */
 221   public static Alignment findXrefSequences(SequenceI[] seqs, boolean dna,
 222           String source, AlignmentI dataset)
 223   {
 224     List<SequenceI> rseqs = new ArrayList<SequenceI>();
 225     AlignedCodonFrame cf = new AlignedCodonFrame();
 226     for (SequenceI seq : seqs)
 227     {
 228       SequenceI dss = seq;
 229       while (dss.getDatasetSequence() != null)
 230       {
 231         dss = dss.getDatasetSequence();
 232       }
 233       boolean found = false;
 234       DBRefEntry[] xrfs = CrossRef.findXDbRefs(dna, dss.getDBRef());
 235       if ((xrfs == null || xrfs.length == 0) && dataset != null)
 236       {
 237         System.out.println("Attempting to find ds Xrefs refs.");
 238         // FIXME should be dss not seq here?
 239         DBRefEntry[] lrfs = CrossRef.findXDbRefs(!dna, seq.getDBRef());
 240         // less ambiguous would be a 'find primary dbRefEntry' method.
 241         // filter for desired source xref here
 242         found = CrossRef.searchDatasetXrefs(dss, !dna, lrfs, dataset,
 243                 rseqs, cf);
 244       }
 245       for (int r = 0; xrfs != null && r < xrfs.length; r++)
 246       {
 247         if (source != null && !source.equals(xrfs[r].getSource()))
 248         {
 249           continue;
 250         }
 251         if (xrfs[r].hasMap())
 252         {
 253           if (xrfs[r].getMap().getTo() != null)
 254           {
 255             SequenceI rsq = new Sequence(xrfs[r].getMap().getTo());
 256             rseqs.add(rsq);
 257             if (xrfs[r].getMap().getMap().getFromRatio() != xrfs[r]
 258                     .getMap().getMap().getToRatio())
 259             {
 260               // get sense of map correct for adding to product alignment.
 261               if (dna)
 262               {
 263                 // map is from dna seq to a protein product
 264                 cf.addMap(dss, rsq, xrfs[r].getMap().getMap());
 265               }
 266               else
 267               {
 268                 // map should be from protein seq to its coding dna
 269                 cf.addMap(rsq, dss, xrfs[r].getMap().getMap().getInverse());
 270               }
 271             }
 272             found = true;
 273           }
 274         }
 275         if (!found)
 276         {
 277           // do a bit more work - search for sequences with references matching
 278           // xrefs on this sequence.
 279           if (dataset != null)
 280           {
 281             found |= searchDataset(dss, xrfs[r], dataset, rseqs, cf); // ,false,!dna);
 282             if (found)
 283             {
 284               xrfs[r] = null; // we've recovered seqs for this one.
 285             }
 286           }
 287         }
 288       }
 289       if (!found)
 290       {
 291         if (xrfs != null && xrfs.length > 0)
 292         {
 293           // Try and get the sequence reference...
 294           /*
 295            * Ideal world - we ask for a sequence fetcher implementation here if
 296            * (jalview.io.RunTimeEnvironment.getSequenceFetcher()) (
 297            */
 298           ASequenceFetcher sftch = new SequenceFetcher();
 299           SequenceI[] retrieved = null;
 300           int l = xrfs.length;
 301           for (int r = 0; r < xrfs.length; r++)
 302           {
 303             // filter out any irrelevant or irretrievable references
 304             if (xrfs[r] == null
 305                     || ((source != null && !source.equals(xrfs[r]
 306                             .getSource())) || !sftch.isFetchable(xrfs[r]
 307                             .getSource())))
 308             {
 309               l--;
 310               xrfs[r] = null;
 311             }
 312           }
 313           if (l > 0)
 314           {
 315             System.out
 316                     .println("Attempting to retrieve cross referenced sequences.");
 317             DBRefEntry[] t = new DBRefEntry[l];
 318             l = 0;
 319             for (int r = 0; r < xrfs.length; r++)
 320             {
 321               if (xrfs[r] != null)
 322               {
 323                 t[l++] = xrfs[r];
 324               }
 325             }
 326             xrfs = t;
 327             try
 328             {
 329               retrieved = sftch.getSequences(xrfs); // problem here is we don't
 330               // know which of xrfs
 331               // resulted in which
 332               // retrieved element
 333             } catch (Exception e)
 334             {
 335               System.err
 336                       .println("Problem whilst retrieving cross references for Sequence : "
 337                               + seq.getName());
 338               e.printStackTrace();
 339             }
 340             if (retrieved != null)
 341             {
 342               for (int rs = 0; rs < retrieved.length; rs++)
 343               {
 344                 // TODO: examine each sequence for 'redundancy'
 345                 DBRefEntry[] dbr = retrieved[rs].getDBRef();
 346                 if (dbr != null && dbr.length > 0)
 347                 {
 348                   for (int di = 0; di < dbr.length; di++)
 349                   {
 350                     // find any entry where we should put in the sequence being
 351                     // cross-referenced into the map
 352                     Mapping map = dbr[di].getMap();
 353                     if (map != null)
 354                     {
 355                       if (map.getTo() != null && map.getMap() != null)
 356                       {
 357                         // should search the local dataset to find any existing
 358                         // candidates for To !
 359                         try
 360                         {
 361                           // compare ms with dss and replace with dss in mapping
 362                           // if map is congruent
 363                           SequenceI ms = map.getTo();
 364                           int sf = map.getMap().getToLowest();
 365                           int st = map.getMap().getToHighest();
 366                           SequenceI mappedrg = ms.getSubSequence(sf, st);
 367                           SequenceI loc = dss.getSubSequence(sf, st);
 368                           if (mappedrg.getLength() > 0
 369                                   && mappedrg.getSequenceAsString().equals(
 370                                           loc.getSequenceAsString()))
 371                           {
 372                             System.err
 373                                     .println("Mapping updated for retrieved crossreference");
 374                             // method to update all refs of existing To on
 375                             // retrieved sequence with dss and merge any props
 376                             // on To onto dss.
 377                             map.setTo(dss);
 378                           }
 379                         } catch (Exception e)
 380                         {
 381                           System.err
 382                                   .println("Exception when consolidating Mapped sequence set...");
 383                           e.printStackTrace(System.err);
 384                         }
 385                       }
 386                     }
 387                   }
 388                 }
 389                 retrieved[rs].updatePDBIds();
 390                 rseqs.add(retrieved[rs]);
 391               }
 392             }
 393           }
 394         }
 395       }
 396     }
 397
 398     Alignment ral = null;
 399     if (rseqs.size() > 0)
 400     {
 401       SequenceI[] rsqs = new SequenceI[rseqs.size()];
 402       rseqs.toArray(rsqs);
 403       ral = new Alignment(rsqs);
 404       if (cf != null && cf.getProtMappings() != null)
 405       {
 406         ral.addCodonFrame(cf);
 407       }
 408     }
 409     return ral;
 410   }
 411
 412   /**
 413    * find references to lrfs in the cross-reference set of each sequence in
 414    * dataset (that is not equal to sequenceI) Identifies matching DBRefEntry
 415    * based on source and accession string only - Map and Version are nulled.
 416    *
 417    * @param sequenceI
 418    * @param lrfs
 419    * @param dataset
 420    * @param rseqs
 421    * @return true if matches were found.
 422    */
 423   private static boolean searchDatasetXrefs(SequenceI sequenceI,
 424           boolean dna, DBRefEntry[] lrfs, AlignmentI dataset,
 425           List<SequenceI> rseqs, AlignedCodonFrame cf)
 426   {
 427     boolean found = false;
 428     if (lrfs == null)
 429     {
 430       return false;
 431     }
 432     for (int i = 0; i < lrfs.length; i++)
 433     {
 434       DBRefEntry xref = new DBRefEntry(lrfs[i]);
 435       // add in wildcards
 436       xref.setVersion(null);
 437       xref.setMap(null);
 438       found = searchDataset(sequenceI, xref, dataset, rseqs, cf, false, dna);
 439     }
 440     return found;
 441   }
 442
 443   /**
 444    * search a given sequence dataset for references matching cross-references to
 445    * the given sequence
 446    *
 447    * @param sequenceI
 448    * @param xrf
 449    * @param dataset
 450    * @param rseqs
 451    *          set of unique sequences
 452    * @param cf
 453    * @return true if one or more unique sequences were found and added
 454    */
 455   public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf,
 456           AlignmentI dataset, List<SequenceI> rseqs, AlignedCodonFrame cf)
 457   {
 458     return searchDataset(sequenceI, xrf, dataset, rseqs, cf, true, false);
 459   }
 460
 461   /**
 462    * TODO: generalise to different protein classifications Search dataset for
 463    * DBRefEntrys matching the given one (xrf) and add the associated sequence to
 464    * rseq.
 465    *
 466    * @param sequenceI
 467    * @param xrf
 468    * @param dataset
 469    * @param rseqs
 470    * @param direct
 471    *          - search all references or only subset
 472    * @param dna
 473    *          search dna or protein xrefs (if direct=false)
 474    * @return true if relationship found and sequence added.
 475    */
 476   public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf,
 477           AlignmentI dataset, List<SequenceI> rseqs, AlignedCodonFrame cf,
 478           boolean direct, boolean dna)
 479   {
 480     boolean found = false;
 481     SequenceI[] typer = new SequenceI[1];
 482     if (dataset == null)
 483     {
 484       return false;
 485     }
 486     if (dataset.getSequences() == null)
 487     {
 488       System.err.println("Empty dataset sequence set - NO VECTOR");
 489       return false;
 490     }
 491     List<SequenceI> ds;
 492     synchronized (ds = dataset.getSequences())
 493     {
 494       for (SequenceI nxt : ds)
 495       {
 496         if (nxt != null)
 497         {
 498           if (nxt.getDatasetSequence() != null)
 499           {
 500             System.err
 501                     .println("Implementation warning: getProducts passed a dataset alignment without dataset sequences in it!");
 502           }
 503           if (nxt != sequenceI && nxt != sequenceI.getDatasetSequence())
 504           {
 505             // check if this is the correct sequence type
 506             {
 507               typer[0] = nxt;
 508               boolean isDna = jalview.util.Comparison.isNucleotide(typer);
 509               if ((direct && isDna == dna) || (!direct && isDna != dna))
 510               {
 511                 // skip this sequence because it is same molecule type
 512                 continue;
 513               }
 514             }
 515
 516             // look for direct or indirect references in common
 517             DBRefEntry[] poss = nxt.getDBRef(), cands = null;
 518             if (direct)
 519             {
 520               cands = jalview.util.DBRefUtils.searchRefs(poss, xrf);
 521             }
 522             else
 523             {
 524               poss = CrossRef.findXDbRefs(dna, poss); //
 525               cands = jalview.util.DBRefUtils.searchRefs(poss, xrf);
 526             }
 527             if (cands != null)
 528             {
 529               if (!rseqs.contains(nxt))
 530               {
 531                 rseqs.add(nxt);
 532                 boolean foundmap = cf != null;
 533                 // don't search if we aren't given a codon map object
 534                 for (int r = 0; foundmap && r < cands.length; r++)
 535                 {
 536                   if (cands[r].hasMap())
 537                   {
 538                     if (cands[r].getMap().getTo() != null
 539                             && cands[r].getMap().getMap().getFromRatio() != cands[r]
 540                                     .getMap().getMap().getToRatio())
 541                     {
 542                       foundmap = true;
 543                       // get sense of map correct for adding to product
 544                       // alignment.
 545                       if (dna)
 546                       {
 547                         // map is from dna seq to a protein product
 548                         cf.addMap(sequenceI, nxt, cands[r].getMap()
 549                                 .getMap());
 550                       }
 551                       else
 552                       {
 553                         // map should be from protein seq to its coding dna
 554                         cf.addMap(nxt, sequenceI, cands[r].getMap()
 555                                 .getMap().getInverse());
 556                       }
 557                     }
 558                   }
 559                 }
 560                 // TODO: add mapping between sequences if necessary
 561                 found = true;
 562               }
 563             }
 564
 565           }
 566         }
 567       }
 568     }
 569     return found;
 570   }
 571
 572   /**
 573    * precalculate different products that can be found for seqs in dataset and
 574    * return them.
 575    *
 576    * @param dna
 577    * @param seqs
 578    * @param dataset
 579    * @param fake
 580    *          - don't actually build lists - just get types
 581    * @return public static Object[] buildXProductsList(boolean dna, SequenceI[]
 582    *         seqs, AlignmentI dataset, boolean fake) { String types[] =
 583    *         jalview.analysis.CrossRef.findSequenceXrefTypes( dna, seqs,
 584    *         dataset); if (types != null) { System.out.println("Xref Types for:
 585    *         "+(dna ? "dna" : "prot")); for (int t = 0; t < types.length; t++) {
 586    *         System.out.println("Type: " + types[t]); SequenceI[] prod =
 587    *         jalview.analysis.CrossRef.findXrefSequences(seqs, dna, types[t]);
 588    *         System.out.println("Found " + ((prod == null) ? "no" : "" +
 589    *         prod.length) + " products"); if (prod!=null) { for (int p=0;
 590    *         p<prod.length; p++) { System.out.println("Prod "+p+":
 591    *         "+prod[p].getDisplayId(true)); } } } } else {
 592    *         System.out.println("Trying getProducts for
 593    *         "+al.getSequenceAt(0).getDisplayId(true));
 594    *         System.out.println("Search DS Xref for: "+(dna ? "dna" : "prot"));
 595    *         // have a bash at finding the products amongst all the retrieved
 596    *         sequences. SequenceI[] prod =
 597    *         jalview.analysis.CrossRef.findXrefSequences(al
 598    *         .getSequencesArray(), dna, null, ds); System.out.println("Found " +
 599    *         ((prod == null) ? "no" : "" + prod.length) + " products"); if
 600    *         (prod!=null) { // select non-equivalent sequences from dataset list
 601    *         for (int p=0; p<prod.length; p++) { System.out.println("Prod "+p+":
 602    *         "+prod[p].getDisplayId(true)); } } } }
 603    */
 604 }