72382393811ab4b8a7e97b6bcae75726958dec0c
[jalview.git] / 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 java.util.ArrayList;
24 import java.util.List;
25 import java.util.Vector;
26
27 import jalview.datamodel.AlignedCodonFrame;
28 import jalview.datamodel.Alignment;
29 import jalview.datamodel.AlignmentI;
30 import jalview.datamodel.DBRefEntry;
31 import jalview.datamodel.DBRefSource;
32 import jalview.datamodel.Sequence;
33 import jalview.datamodel.SequenceI;
34 import jalview.util.DBRefUtils;
35 import jalview.ws.SequenceFetcher;
36 import jalview.ws.seqfetcher.ASequenceFetcher;
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,
412           AlignedCodonFrame cf)
413   {
414     boolean found = false;
415     if (lrfs == null)
416     {
417       return false;
418     }
419     for (int i = 0; i < lrfs.length; i++)
420     {
421       DBRefEntry xref = new DBRefEntry(lrfs[i]);
422       // add in wildcards
423       xref.setVersion(null);
424       xref.setMap(null);
425       found = searchDataset(sequenceI, xref, dataset, rseqs, cf, false, dna);
426     }
427     return found;
428   }
429
430   /**
431    * search a given sequence dataset for references matching cross-references to
432    * the given sequence
433    * 
434    * @param sequenceI
435    * @param xrf
436    * @param dataset
437    * @param rseqs
438    *          set of unique sequences
439    * @param cf
440    * @return true if one or more unique sequences were found and added
441    */
442   public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf,
443           AlignmentI dataset, List<SequenceI> rseqs, AlignedCodonFrame cf)
444   {
445     return searchDataset(sequenceI, xrf, dataset, rseqs, cf, true, false);
446   }
447
448   /**
449    * TODO: generalise to different protein classifications Search dataset for
450    * DBRefEntrys matching the given one (xrf) and add the associated sequence to
451    * rseq.
452    * 
453    * @param sequenceI
454    * @param xrf
455    * @param dataset
456    * @param rseqs
457    * @param direct
458    *          - search all references or only subset
459    * @param dna
460    *          search dna or protein xrefs (if direct=false)
461    * @return true if relationship found and sequence added.
462    */
463   public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf,
464           AlignmentI dataset, List<SequenceI> rseqs, AlignedCodonFrame cf,
465           boolean direct, boolean dna)
466   {
467     boolean found = false;
468     SequenceI[] typer = new SequenceI[1];
469     if (dataset == null)
470     {
471       return false;
472     }
473     if (dataset.getSequences() == null)
474     {
475       System.err.println("Empty dataset sequence set - NO VECTOR");
476       return false;
477     }
478     List<SequenceI> ds;
479     synchronized (ds = dataset.getSequences())
480     {
481       for (SequenceI nxt : ds)
482       {
483         if (nxt != null)
484         {
485           if (nxt.getDatasetSequence() != null)
486           {
487             System.err
488                     .println("Implementation warning: getProducts passed a dataset alignment without dataset sequences in it!");
489           }
490           if (nxt != sequenceI && nxt != sequenceI.getDatasetSequence())
491           {
492             // check if this is the correct sequence type
493             {
494               typer[0] = nxt;
495               boolean isDna = jalview.util.Comparison.isNucleotide(typer);
496               if ((direct && isDna == dna) || (!direct && isDna != dna))
497               {
498                 // skip this sequence because it is same molecule type
499                 continue;
500               }
501             }
502
503             // look for direct or indirect references in common
504             DBRefEntry[] poss = nxt.getDBRef(), cands = null;
505             if (direct)
506             {
507               cands = jalview.util.DBRefUtils.searchRefs(poss, xrf);
508             }
509             else
510             {
511               poss = CrossRef.findXDbRefs(dna, poss); //
512               cands = jalview.util.DBRefUtils.searchRefs(poss, xrf);
513             }
514             if (cands != null)
515             {
516               if (!rseqs.contains(nxt))
517               {
518                 rseqs.add(nxt);
519                 boolean foundmap = cf != null;
520                 // don't search if we aren't given a codon map object
521                 for (int r = 0; foundmap && r < cands.length; r++)
522                 {
523                   if (cands[r].hasMap())
524                   {
525                     if (cands[r].getMap().getTo() != null
526                             && cands[r].getMap().getMap().getFromRatio() != cands[r]
527                                     .getMap().getMap().getToRatio())
528                     {
529                       foundmap = true;
530                       // get sense of map correct for adding to product
531                       // alignment.
532                       if (dna)
533                       {
534                         // map is from dna seq to a protein product
535                         cf.addMap(sequenceI, nxt, cands[r].getMap()
536                                 .getMap());
537                       }
538                       else
539                       {
540                         // map should be from protein seq to its coding dna
541                         cf.addMap(nxt, sequenceI, cands[r].getMap()
542                                 .getMap().getInverse());
543                       }
544                     }
545                   }
546                 }
547                 // TODO: add mapping between sequences if necessary
548                 found = true;
549               }
550             }
551
552           }
553         }
554       }
555     }
556     return found;
557   }
558
559   /**
560    * precalculate different products that can be found for seqs in dataset and
561    * return them.
562    * 
563    * @param dna
564    * @param seqs
565    * @param dataset
566    * @param fake
567    *          - don't actually build lists - just get types
568    * @return public static Object[] buildXProductsList(boolean dna, SequenceI[]
569    *         seqs, AlignmentI dataset, boolean fake) { String types[] =
570    *         jalview.analysis.CrossRef.findSequenceXrefTypes( dna, seqs,
571    *         dataset); if (types != null) { System.out.println("Xref Types for:
572    *         "+(dna ? "dna" : "prot")); for (int t = 0; t < types.length; t++) {
573    *         System.out.println("Type: " + types[t]); SequenceI[] prod =
574    *         jalview.analysis.CrossRef.findXrefSequences(seqs, dna, types[t]);
575    *         System.out.println("Found " + ((prod == null) ? "no" : "" +
576    *         prod.length) + " products"); if (prod!=null) { for (int p=0;
577    *         p<prod.length; p++) { System.out.println("Prod "+p+":
578    *         "+prod[p].getDisplayId(true)); } } } } else {
579    *         System.out.println("Trying getProducts for
580    *         "+al.getSequenceAt(0).getDisplayId(true));
581    *         System.out.println("Search DS Xref for: "+(dna ? "dna" : "prot"));
582    *         // have a bash at finding the products amongst all the retrieved
583    *         sequences. SequenceI[] prod =
584    *         jalview.analysis.CrossRef.findXrefSequences(al
585    *         .getSequencesArray(), dna, null, ds); System.out.println("Found " +
586    *         ((prod == null) ? "no" : "" + prod.length) + " products"); if
587    *         (prod!=null) { // select non-equivalent sequences from dataset list
588    *         for (int p=0; p<prod.length; p++) { System.out.println("Prod "+p+":
589    *         "+prod[p].getDisplayId(true)); } } } }
590    */
591 }