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