src/jalview/analysis/CrossRef.java

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