src/jalview/analysis/CrossRef.java

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