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