src/jalview/analysis/CrossRef.java

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