JAL-2113 tidy up canonicalisation of Uniprot xrefs in EMBL data
[jalview.git] / src / jalview / datamodel / xdb / embl / EmblEntry.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.datamodel.xdb.embl;
22
23 import jalview.analysis.SequenceIdMatcher;
24 import jalview.bin.Cache;
25 import jalview.datamodel.DBRefEntry;
26 import jalview.datamodel.DBRefSource;
27 import jalview.datamodel.FeatureProperties;
28 import jalview.datamodel.Mapping;
29 import jalview.datamodel.Sequence;
30 import jalview.datamodel.SequenceFeature;
31 import jalview.datamodel.SequenceI;
32 import jalview.util.DBRefUtils;
33 import jalview.util.DnaUtils;
34 import jalview.util.MapList;
35 import jalview.util.MappingUtils;
36 import jalview.util.StringUtils;
37
38 import java.text.ParseException;
39 import java.util.Arrays;
40 import java.util.Hashtable;
41 import java.util.List;
42 import java.util.Map;
43 import java.util.Map.Entry;
44 import java.util.Vector;
45 import java.util.regex.Pattern;
46
47 /**
48  * Data model for one entry returned from an EMBL query, as marshalled by a
49  * Castor binding file
50  * 
51  * For example:
52  * http://www.ebi.ac.uk/Tools/dbfetch/dbfetch?db=ena_sequence&id=J03321
53  * &format=emblxml
54  * 
55  * @see embl_mapping.xml
56  */
57 public class EmblEntry
58 {
59   private static final Pattern SPACE_PATTERN = Pattern.compile(" ");
60
61   String accession;
62
63   String entryVersion;
64
65   String sequenceVersion;
66
67   String dataClass;
68
69   String moleculeType;
70
71   String topology;
72
73   String sequenceLength;
74
75   String taxonomicDivision;
76
77   String description;
78
79   String firstPublicDate;
80
81   String firstPublicRelease;
82
83   String lastUpdatedDate;
84
85   String lastUpdatedRelease;
86
87   Vector<String> keywords;
88
89   Vector<DBRefEntry> dbRefs;
90
91   Vector<EmblFeature> features;
92
93   EmblSequence sequence;
94
95   /**
96    * @return the accession
97    */
98   public String getAccession()
99   {
100     return accession;
101   }
102
103   /**
104    * @param accession
105    *          the accession to set
106    */
107   public void setAccession(String accession)
108   {
109     this.accession = accession;
110   }
111
112   /**
113    * @return the dbRefs
114    */
115   public Vector<DBRefEntry> getDbRefs()
116   {
117     return dbRefs;
118   }
119
120   /**
121    * @param dbRefs
122    *          the dbRefs to set
123    */
124   public void setDbRefs(Vector<DBRefEntry> dbRefs)
125   {
126     this.dbRefs = dbRefs;
127   }
128
129   /**
130    * @return the features
131    */
132   public Vector<EmblFeature> getFeatures()
133   {
134     return features;
135   }
136
137   /**
138    * @param features
139    *          the features to set
140    */
141   public void setFeatures(Vector<EmblFeature> features)
142   {
143     this.features = features;
144   }
145
146   /**
147    * @return the keywords
148    */
149   public Vector<String> getKeywords()
150   {
151     return keywords;
152   }
153
154   /**
155    * @param keywords
156    *          the keywords to set
157    */
158   public void setKeywords(Vector<String> keywords)
159   {
160     this.keywords = keywords;
161   }
162
163   /**
164    * @return the sequence
165    */
166   public EmblSequence getSequence()
167   {
168     return sequence;
169   }
170
171   /**
172    * @param sequence
173    *          the sequence to set
174    */
175   public void setSequence(EmblSequence sequence)
176   {
177     this.sequence = sequence;
178   }
179
180   /**
181    * Recover annotated sequences from EMBL file
182    * 
183    * @param sourceDb
184    * @param peptides
185    *          a list of protein products found so far (to add to)
186    * @return dna dataset sequence with DBRefs and features
187    */
188   public SequenceI getSequence(String sourceDb, List<SequenceI> peptides)
189   {
190     SequenceI dna = new Sequence(sourceDb + "|" + accession,
191             sequence.getSequence());
192     dna.setDescription(description);
193     DBRefEntry retrievedref = new DBRefEntry(sourceDb,
194             getSequenceVersion(), accession);
195     dna.addDBRef(retrievedref);
196     // add map to indicate the sequence is a valid coordinate frame for the
197     // dbref
198     retrievedref.setMap(new Mapping(null, new int[] { 1, dna.getLength() },
199             new int[] { 1, dna.getLength() }, 1, 1));
200
201     /*
202      * transform EMBL Database refs to canonical form
203      */
204     if (dbRefs != null)
205     {
206       for (DBRefEntry dbref : dbRefs)
207       {
208         dbref.setSource(DBRefUtils.getCanonicalName(dbref.getSource()));
209         dna.addDBRef(dbref);
210       }
211     }
212
213     SequenceIdMatcher matcher = new SequenceIdMatcher(peptides);
214     try
215     {
216       for (EmblFeature feature : features)
217       {
218         if (FeatureProperties.isCodingFeature(sourceDb, feature.getName()))
219         {
220           parseCodingFeature(feature, sourceDb, dna, peptides, matcher);
221         }
222       }
223     } catch (Exception e)
224     {
225       System.err.println("EMBL Record Features parsing error!");
226       System.err
227               .println("Please report the following to help@jalview.org :");
228       System.err.println("EMBL Record " + accession);
229       System.err.println("Resulted in exception: " + e.getMessage());
230       e.printStackTrace(System.err);
231     }
232
233     return dna;
234   }
235
236   /**
237    * Extracts coding region and product from a CDS feature and properly decorate
238    * it with annotations.
239    * 
240    * @param feature
241    *          coding feature
242    * @param sourceDb
243    *          source database for the EMBLXML
244    * @param dna
245    *          parent dna sequence for this record
246    * @param peptides
247    *          list of protein product sequences for Embl entry
248    * @param matcher
249    *          helper to match xrefs in already retrieved sequences
250    */
251   void parseCodingFeature(EmblFeature feature, String sourceDb,
252           SequenceI dna, List<SequenceI> peptides, SequenceIdMatcher matcher)
253   {
254     boolean isEmblCdna = sourceDb.equals(DBRefSource.EMBLCDS);
255
256     int[] exon = getCdsRanges(feature);
257
258     String prseq = null;
259     String prname = "";
260     String prid = null;
261     Map<String, String> vals = new Hashtable<String, String>();
262
263     /*
264      * codon_start 1/2/3 in EMBL corresponds to phase 0/1/2 in CDS
265      * (phase is required for CDS features in GFF3 format)
266      */
267     int codonStart = 1;
268
269     /*
270      * parse qualifiers, saving protein translation, protein id,
271      * codon start position, product (name), and 'other values'
272      */
273     if (feature.getQualifiers() != null)
274     {
275       for (Qualifier q : feature.getQualifiers())
276       {
277         String qname = q.getName();
278         if (qname.equals("translation"))
279         {
280           // remove all spaces (precompiled String.replaceAll(" ", ""))
281           prseq = SPACE_PATTERN.matcher(q.getValues()[0]).replaceAll("");
282         }
283         else if (qname.equals("protein_id"))
284         {
285           prid = q.getValues()[0];
286         }
287         else if (qname.equals("codon_start"))
288         {
289           try
290           {
291             codonStart = Integer.parseInt(q.getValues()[0]);
292           } catch (NumberFormatException e)
293           {
294             System.err.println("Invalid codon_start in XML for "
295                     + accession + ": " + e.getMessage());
296           }
297         }
298         else if (qname.equals("product"))
299         {
300           // sometimes name is returned e.g. for V00488
301           prname = q.getValues()[0];
302         }
303         else
304         {
305           // throw anything else into the additional properties hash
306           String[] qvals = q.getValues();
307           if (qvals != null)
308           {
309             String commaSeparated = StringUtils.arrayToSeparatorList(qvals,
310                     ",");
311             vals.put(qname, commaSeparated);
312           }
313         }
314       }
315     }
316
317     DBRefEntry protEMBLCDS = null;
318     exon = MappingUtils.removeStartPositions(codonStart - 1, exon);
319     boolean noProteinDbref = true;
320
321     SequenceI product = null;
322     Mapping map = null;
323     if (prseq != null && prname != null && prid != null)
324     {
325       /*
326        * look for product in peptides list, if not found, add it
327        */
328       product = matcher.findIdMatch(prid);
329       if (product == null)
330       {
331         product = new Sequence(prid, prseq, 1, prseq.length());
332         product.setDescription(((prname.length() == 0) ? "Protein Product from "
333                 + sourceDb
334                 : prname));
335         peptides.add(product);
336         matcher.add(product);
337       }
338
339       // we have everything - create the mapping and perhaps the protein
340       // sequence
341       if (exon == null || exon.length == 0)
342       {
343         System.err
344                 .println("Implementation Notice: EMBLCDS records not properly supported yet - Making up the CDNA region of this sequence... may be incorrect ("
345                         + sourceDb + ":" + getAccession() + ")");
346         if (prseq.length() * 3 == (1 - codonStart + dna.getSequence().length))
347         {
348           System.err
349                   .println("Not allowing for additional stop codon at end of cDNA fragment... !");
350           // this might occur for CDS sequences where no features are
351           // marked.
352           exon = new int[] { dna.getStart() + (codonStart - 1),
353               dna.getEnd() };
354           map = new Mapping(product, exon, new int[] { 1, prseq.length() },
355                   3, 1);
356         }
357         if ((prseq.length() + 1) * 3 == (1 - codonStart + dna.getSequence().length))
358         {
359           System.err
360                   .println("Allowing for additional stop codon at end of cDNA fragment... will probably cause an error in VAMSAs!");
361           exon = new int[] { dna.getStart() + (codonStart - 1),
362               dna.getEnd() - 3 };
363           map = new Mapping(product, exon, new int[] { 1, prseq.length() },
364                   3, 1);
365         }
366       }
367       else
368       {
369         // Trim the exon mapping if necessary - the given product may only be a
370         // fragment of a larger protein. (EMBL:AY043181 is an example)
371
372         if (isEmblCdna)
373         {
374           // TODO: Add a DbRef back to the parent EMBL sequence with the exon
375           // map
376           // if given a dataset reference, search dataset for parent EMBL
377           // sequence if it exists and set its map
378           // make a new feature annotating the coding contig
379         }
380         else
381         {
382           // final product length truncation check
383           // TODO should from range include stop codon even if not in protein
384           // in order to include stop codon in CDS sequence (as done for
385           // Ensembl)?
386           int[] cdsRanges = adjustForProteinLength(prseq.length(), exon);
387           map = new Mapping(product, cdsRanges, new int[] { 1,
388               prseq.length() }, 3, 1);
389           // reconstruct the EMBLCDS entry
390           // TODO: this is only necessary when there codon annotation is
391           // complete (I think JBPNote)
392           DBRefEntry pcdnaref = new DBRefEntry();
393           pcdnaref.setAccessionId(prid);
394           pcdnaref.setSource(DBRefSource.EMBLCDS);
395           pcdnaref.setVersion(getSequenceVersion()); // same as parent EMBL
396                                                      // version.
397           MapList mp = new MapList(new int[] { 1, prseq.length() },
398                   new int[] { 1 + (codonStart - 1),
399                       (codonStart - 1) + 3 * prseq.length() }, 1, 3);
400           pcdnaref.setMap(new Mapping(mp));
401           if (product != null)
402           {
403             product.addDBRef(pcdnaref);
404             protEMBLCDS = new DBRefEntry(pcdnaref);
405             protEMBLCDS.setSource(DBRefSource.EMBLCDSProduct);
406             product.addDBRef(protEMBLCDS);
407           }
408         }
409       }
410       // add cds feature to dna seq - this may include the stop codon
411       for (int xint = 0; exon != null && xint < exon.length; xint += 2)
412       {
413         SequenceFeature sf = makeCdsFeature(exon, xint, prname, prid, vals,
414                 codonStart);
415         sf.setType(feature.getName()); // "CDS"
416         sf.setEnaLocation(feature.getLocation());
417         sf.setFeatureGroup(sourceDb);
418         dna.addSequenceFeature(sf);
419       }
420     }
421
422     /*
423      * add dbRefs to sequence, and mappings for Uniprot xrefs
424      */
425     if (feature.dbRefs != null)
426     {
427       boolean mappingUsed = false;
428       for (DBRefEntry ref : feature.dbRefs)
429       {
430         /*
431          * ensure UniProtKB/Swiss-Prot converted to UNIPROT
432          */
433         ref.setSource(DBRefUtils.getCanonicalName(ref.getSource()));
434         if (ref.getSource().equals(DBRefSource.UNIPROT))
435         {
436           String proteinSeqName = DBRefSource.UNIPROT + "|"
437                   + ref.getAccessionId();
438           if (map != null && map.getTo() != null)
439           {
440             if (mappingUsed)
441             {
442               /*
443                * two or more Uniprot xrefs for the same CDS - 
444                * each needs a distinct Mapping (as to a different sequence)
445                */
446               map = new Mapping(map);
447             }
448             mappingUsed = true;
449
450             /*
451              * try to locate the protein mapped to (possibly by a 
452              * previous CDS feature)
453              */
454             SequenceI proteinSeq = matcher.findIdMatch(proteinSeqName);
455             if (proteinSeq == null)
456             {
457               proteinSeq = new Sequence(proteinSeqName,
458                       product.getSequenceAsString());
459               matcher.add(proteinSeq);
460               peptides.add(proteinSeq);
461             }
462             map.setTo(proteinSeq);
463             map.getTo().addDBRef(
464                     new DBRefEntry(ref.getSource(), ref.getVersion(), ref
465                             .getAccessionId()));
466             ref.setMap(map);
467           }
468           noProteinDbref = false;
469         }
470         if (product != null)
471         {
472           DBRefEntry pref = new DBRefEntry(ref.getSource(),
473                   ref.getVersion(), ref.getAccessionId());
474           pref.setMap(null); // reference is direct
475           product.addDBRef(pref);
476           // Add converse mapping reference
477           if (map != null)
478           {
479             Mapping pmap = new Mapping(dna, map.getMap().getInverse());
480             pref = new DBRefEntry(sourceDb, getSequenceVersion(),
481                     this.getAccession());
482             pref.setMap(pmap);
483             if (map.getTo() != null)
484             {
485               map.getTo().addDBRef(pref);
486             }
487           }
488         }
489         dna.addDBRef(ref);
490       }
491       if (noProteinDbref && product != null)
492       {
493         // add protein coding reference to dna sequence so xref matches
494         if (protEMBLCDS == null)
495         {
496           protEMBLCDS = new DBRefEntry();
497           protEMBLCDS.setAccessionId(prid);
498           protEMBLCDS.setSource(DBRefSource.EMBLCDSProduct);
499           protEMBLCDS.setVersion(getSequenceVersion());
500           protEMBLCDS
501                   .setMap(new Mapping(product, map.getMap().getInverse()));
502         }
503         product.addDBRef(protEMBLCDS);
504
505         // Add converse mapping reference
506         if (map != null)
507         {
508           Mapping pmap = new Mapping(product, protEMBLCDS.getMap().getMap()
509                   .getInverse());
510           DBRefEntry ncMap = new DBRefEntry(protEMBLCDS);
511           ncMap.setMap(pmap);
512           if (map.getTo() != null)
513           {
514             dna.addDBRef(ncMap);
515           }
516         }
517       }
518     }
519   }
520
521   /**
522    * Helper method to construct a SequenceFeature for one cds range
523    * 
524    * @param exons
525    *          array of cds [start, end, ...] positions
526    * @param exonStartIndex
527    *          offset into the exons array
528    * @param proteinName
529    * @param proteinAccessionId
530    * @param vals
531    *          map of 'miscellaneous values' for feature
532    * @param codonStart
533    *          codon start position for CDS (1/2/3, normally 1)
534    * @return
535    */
536   protected SequenceFeature makeCdsFeature(int[] exons, int exonStartIndex,
537           String proteinName, String proteinAccessionId,
538           Map<String, String> vals, int codonStart)
539   {
540     int exonNumber = exonStartIndex / 2 + 1;
541     SequenceFeature sf = new SequenceFeature();
542     sf.setBegin(Math.min(exons[exonStartIndex], exons[exonStartIndex + 1]));
543     sf.setEnd(Math.max(exons[exonStartIndex], exons[exonStartIndex + 1]));
544     sf.setDescription(String.format("Exon %d for protein '%s' EMBLCDS:%s",
545             exonNumber, proteinName, proteinAccessionId));
546     sf.setPhase(String.valueOf(codonStart - 1));
547     sf.setStrand(exons[exonStartIndex] <= exons[exonStartIndex + 1] ? "+"
548             : "-");
549     sf.setValue(FeatureProperties.EXONPOS, exonNumber);
550     sf.setValue(FeatureProperties.EXONPRODUCT, proteinName);
551     if (!vals.isEmpty())
552     {
553       StringBuilder sb = new StringBuilder();
554       boolean first = true;
555       for (Entry<String, String> val : vals.entrySet())
556       {
557         if (!first)
558         {
559           sb.append(";");
560         }
561         sb.append(val.getKey()).append("=").append(val.getValue());
562         first = false;
563         sf.setValue(val.getKey(), val.getValue());
564       }
565       sf.setAttributes(sb.toString());
566     }
567     return sf;
568   }
569
570   /**
571    * Returns the CDS positions as a single array of [start, end, start, end...]
572    * positions. If on the reverse strand, these will be in descending order.
573    * 
574    * @param feature
575    * @return
576    */
577   protected int[] getCdsRanges(EmblFeature feature)
578   {
579     if (feature.location == null)
580     {
581       return new int[] {};
582     }
583
584     try
585     {
586       List<int[]> ranges = DnaUtils.parseLocation(feature.location);
587       return listToArray(ranges);
588     } catch (ParseException e)
589     {
590       Cache.log.warn(String.format(
591               "Not parsing inexact CDS location %s in ENA %s",
592               feature.location, this.accession));
593       return new int[] {};
594     }
595   }
596
597   /**
598    * Converts a list of [start, end] ranges to a single array of [start, end,
599    * start, end ...]
600    * 
601    * @param ranges
602    * @return
603    */
604   int[] listToArray(List<int[]> ranges)
605   {
606     int[] result = new int[ranges.size() * 2];
607     int i = 0;
608     for (int[] range : ranges)
609     {
610       result[i++] = range[0];
611       result[i++] = range[1];
612     }
613     return result;
614   }
615
616   /**
617    * truncate the last exon interval to the prlength'th codon
618    * 
619    * @param prlength
620    * @param exon
621    * @return new exon
622    */
623   static int[] adjustForProteinLength(int prlength, int[] exon)
624   {
625     if (prlength <= 0 || exon == null)
626     {
627       return exon;
628     }
629     int desiredCdsLength = prlength * 3;
630     int exonLength = MappingUtils.getLength(Arrays.asList(exon));
631
632     /*
633      * assuming here exon might include stop codon in addition to protein codons
634      */
635     if (desiredCdsLength == exonLength
636             || desiredCdsLength == exonLength - 3)
637     {
638       return exon;
639     }
640
641     int origxon[];
642     int sxpos = -1;
643     int endxon = 0;
644     origxon = new int[exon.length];
645     System.arraycopy(exon, 0, origxon, 0, exon.length);
646     int cdspos = 0;
647     for (int x = 0; x < exon.length; x += 2)
648     {
649       cdspos += Math.abs(exon[x + 1] - exon[x]) + 1;
650       if (desiredCdsLength <= cdspos)
651       {
652         // advanced beyond last codon.
653         sxpos = x;
654         if (desiredCdsLength != cdspos)
655         {
656           // System.err
657           // .println("Truncating final exon interval on region by "
658           // + (cdspos - cdslength));
659         }
660
661         /*
662          * shrink the final exon - reduce end position if forward
663          * strand, increase it if reverse
664          */
665         if (exon[x + 1] >= exon[x])
666         {
667           endxon = exon[x + 1] - cdspos + desiredCdsLength;
668         }
669         else
670         {
671           endxon = exon[x + 1] + cdspos - desiredCdsLength;
672         }
673         break;
674       }
675     }
676
677     if (sxpos != -1)
678     {
679       // and trim the exon interval set if necessary
680       int[] nxon = new int[sxpos + 2];
681       System.arraycopy(exon, 0, nxon, 0, sxpos + 2);
682       nxon[sxpos + 1] = endxon; // update the end boundary for the new exon
683                                 // set
684       exon = nxon;
685     }
686     return exon;
687   }
688
689   public String getSequenceVersion()
690   {
691     return sequenceVersion;
692   }
693
694   public void setSequenceVersion(String sequenceVersion)
695   {
696     this.sequenceVersion = sequenceVersion;
697   }
698
699   public String getSequenceLength()
700   {
701     return sequenceLength;
702   }
703
704   public void setSequenceLength(String sequenceLength)
705   {
706     this.sequenceLength = sequenceLength;
707   }
708
709   public String getEntryVersion()
710   {
711     return entryVersion;
712   }
713
714   public void setEntryVersion(String entryVersion)
715   {
716     this.entryVersion = entryVersion;
717   }
718
719   public String getMoleculeType()
720   {
721     return moleculeType;
722   }
723
724   public void setMoleculeType(String moleculeType)
725   {
726     this.moleculeType = moleculeType;
727   }
728
729   public String getTopology()
730   {
731     return topology;
732   }
733
734   public void setTopology(String topology)
735   {
736     this.topology = topology;
737   }
738
739   public String getTaxonomicDivision()
740   {
741     return taxonomicDivision;
742   }
743
744   public void setTaxonomicDivision(String taxonomicDivision)
745   {
746     this.taxonomicDivision = taxonomicDivision;
747   }
748
749   public String getDescription()
750   {
751     return description;
752   }
753
754   public void setDescription(String description)
755   {
756     this.description = description;
757   }
758
759   public String getFirstPublicDate()
760   {
761     return firstPublicDate;
762   }
763
764   public void setFirstPublicDate(String firstPublicDate)
765   {
766     this.firstPublicDate = firstPublicDate;
767   }
768
769   public String getFirstPublicRelease()
770   {
771     return firstPublicRelease;
772   }
773
774   public void setFirstPublicRelease(String firstPublicRelease)
775   {
776     this.firstPublicRelease = firstPublicRelease;
777   }
778
779   public String getLastUpdatedDate()
780   {
781     return lastUpdatedDate;
782   }
783
784   public void setLastUpdatedDate(String lastUpdatedDate)
785   {
786     this.lastUpdatedDate = lastUpdatedDate;
787   }
788
789   public String getLastUpdatedRelease()
790   {
791     return lastUpdatedRelease;
792   }
793
794   public void setLastUpdatedRelease(String lastUpdatedRelease)
795   {
796     this.lastUpdatedRelease = lastUpdatedRelease;
797   }
798
799   public String getDataClass()
800   {
801     return dataClass;
802   }
803
804   public void setDataClass(String dataClass)
805   {
806     this.dataClass = dataClass;
807   }
808 }