f0041f697d3eb98aeefb7bbd9e8c1906b51a49ee
[jalview.git] / src / jalview / analysis / Dna.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.ArrayList;
21 import java.util.Hashtable;
22 import java.util.Vector;
23
24 import jalview.datamodel.AlignedCodonFrame;
25 import jalview.datamodel.Alignment;
26 import jalview.datamodel.AlignmentAnnotation;
27 import jalview.datamodel.AlignmentI;
28 import jalview.datamodel.Annotation;
29 import jalview.datamodel.DBRefEntry;
30 import jalview.datamodel.FeatureProperties;
31 import jalview.datamodel.Mapping;
32 import jalview.datamodel.Sequence;
33 import jalview.datamodel.SequenceFeature;
34 import jalview.datamodel.SequenceI;
35 import jalview.schemes.ResidueProperties;
36 import jalview.util.MapList;
37 import jalview.util.ShiftList;
38
39 public class Dna
40 {
41   /**
42    * 
43    * @param cdp1
44    * @param cdp2
45    * @return -1 if cdp1 aligns before cdp2, 0 if in the same column or cdp2 is
46    *         null, +1 if after cdp2
47    */
48   private static int compare_codonpos(int[] cdp1, int[] cdp2)
49   {
50     if (cdp2 == null
51             || (cdp1[0] == cdp2[0] && cdp1[1] == cdp2[1] && cdp1[2] == cdp2[2]))
52       return 0;
53     if (cdp1[0] < cdp2[0] || cdp1[1] < cdp2[1] || cdp1[2] < cdp2[2])
54       return -1; // one base in cdp1 precedes the corresponding base in the
55     // other codon
56     return 1; // one base in cdp1 appears after the corresponding base in the
57     // other codon.
58   }
59
60   /**
61    * DNA->mapped protein sequence alignment translation given set of sequences
62    * 1. id distinct coding regions within selected region for each sequence 2.
63    * generate peptides based on inframe (or given) translation or (optionally
64    * and where specified) out of frame translations (annotated appropriately) 3.
65    * align peptides based on codon alignment
66    */
67   /**
68    * id potential products from dna 1. search for distinct products within
69    * selected region for each selected sequence 2. group by associated DB type.
70    * 3. return as form for input into above function
71    */
72   /**
73    * 
74    */
75   /**
76    * create a new alignment of protein sequences by an inframe translation of
77    * the provided NA sequences
78    * 
79    * @param selection
80    * @param seqstring
81    * @param viscontigs
82    * @param gapCharacter
83    * @param annotations
84    * @param aWidth
85    * @param dataset
86    *          destination dataset for translated sequences and mappings
87    * @return
88    */
89   public static AlignmentI CdnaTranslate(SequenceI[] selection,
90           String[] seqstring, int viscontigs[], char gapCharacter,
91           AlignmentAnnotation[] annotations, int aWidth, Alignment dataset)
92   {
93     return CdnaTranslate(selection, seqstring, null, viscontigs,
94             gapCharacter, annotations, aWidth, dataset);
95   }
96
97   /**
98    * 
99    * @param selection
100    * @param seqstring
101    * @param product
102    *          - array of DbRefEntry objects from which exon map in seqstring is
103    *          derived
104    * @param viscontigs
105    * @param gapCharacter
106    * @param annotations
107    * @param aWidth
108    * @param dataset
109    * @return
110    */
111   public static AlignmentI CdnaTranslate(SequenceI[] selection,
112           String[] seqstring, DBRefEntry[] product, int viscontigs[],
113           char gapCharacter, AlignmentAnnotation[] annotations, int aWidth,
114           Alignment dataset)
115   {
116     AlignedCodonFrame codons = new AlignedCodonFrame(aWidth); // stores hash of
117     // subsequent
118     // positions for
119     // each codon
120     // start position
121     // in alignment
122     int s, sSize = selection.length;
123     Vector pepseqs = new Vector();
124     for (s = 0; s < sSize; s++)
125     {
126       SequenceI newseq = translateCodingRegion(selection[s], seqstring[s],
127               viscontigs, codons, gapCharacter,
128               (product != null) ? product[s] : null); // possibly anonymous
129       // product
130       if (newseq != null)
131       {
132         pepseqs.addElement(newseq);
133         SequenceI ds = newseq;
134         while (ds.getDatasetSequence() != null)
135         {
136           ds = ds.getDatasetSequence();
137         }
138         dataset.addSequence(ds);
139       }
140     }
141     if (codons.aaWidth == 0)
142       return null;
143     SequenceI[] newseqs = new SequenceI[pepseqs.size()];
144     pepseqs.copyInto(newseqs);
145     AlignmentI al = new Alignment(newseqs);
146     al.padGaps(); // ensure we look aligned.
147     al.setDataset(dataset);
148     translateAlignedAnnotations(annotations, al, codons);
149     al.addCodonFrame(codons);
150     return al;
151   }
152
153   /**
154    * fake the collection of DbRefs with associated exon mappings to identify if
155    * a translation would generate distinct product in the currently selected
156    * region.
157    * 
158    * @param selection
159    * @param viscontigs
160    * @return
161    */
162   public static boolean canTranslate(SequenceI[] selection,
163           int viscontigs[])
164   {
165     for (int gd = 0; gd < selection.length; gd++)
166     {
167       SequenceI dna = selection[gd];
168       jalview.datamodel.DBRefEntry[] dnarefs = jalview.util.DBRefUtils
169               .selectRefs(dna.getDBRef(),
170                       jalview.datamodel.DBRefSource.DNACODINGDBS);
171       if (dnarefs != null)
172       {
173         // intersect with pep
174         // intersect with pep
175         Vector mappedrefs = new Vector();
176         DBRefEntry[] refs = dna.getDBRef();
177         for (int d = 0; d < refs.length; d++)
178         {
179           if (refs[d].getMap() != null && refs[d].getMap().getMap() != null
180                   && refs[d].getMap().getMap().getFromRatio() == 3
181                   && refs[d].getMap().getMap().getToRatio() == 1)
182           {
183             mappedrefs.addElement(refs[d]); // add translated protein maps
184           }
185         }
186         dnarefs = new DBRefEntry[mappedrefs.size()];
187         mappedrefs.copyInto(dnarefs);
188         for (int d = 0; d < dnarefs.length; d++)
189         {
190           Mapping mp = dnarefs[d].getMap();
191           if (mp != null)
192           {
193             for (int vc = 0; vc < viscontigs.length; vc += 2)
194             {
195               int[] mpr = mp.locateMappedRange(viscontigs[vc],
196                       viscontigs[vc + 1]);
197               if (mpr != null)
198               {
199                 return true;
200               }
201             }
202           }
203         }
204       }
205     }
206     return false;
207   }
208
209   /**
210    * generate a set of translated protein products from annotated sequenceI
211    * 
212    * @param selection
213    * @param viscontigs
214    * @param gapCharacter
215    * @param dataset
216    *          destination dataset for translated sequences
217    * @param annotations
218    * @param aWidth
219    * @return
220    */
221   public static AlignmentI CdnaTranslate(SequenceI[] selection,
222           int viscontigs[], char gapCharacter, Alignment dataset)
223   {
224     int alwidth = 0;
225     Vector cdnasqs = new Vector();
226     Vector cdnasqi = new Vector();
227     Vector cdnaprod = new Vector();
228     for (int gd = 0; gd < selection.length; gd++)
229     {
230       SequenceI dna = selection[gd];
231       jalview.datamodel.DBRefEntry[] dnarefs = jalview.util.DBRefUtils
232               .selectRefs(dna.getDBRef(),
233                       jalview.datamodel.DBRefSource.DNACODINGDBS);
234       if (dnarefs != null)
235       {
236         // intersect with pep
237         Vector mappedrefs = new Vector();
238         DBRefEntry[] refs = dna.getDBRef();
239         for (int d = 0; d < refs.length; d++)
240         {
241           if (refs[d].getMap() != null && refs[d].getMap().getMap() != null
242                   && refs[d].getMap().getMap().getFromRatio() == 3
243                   && refs[d].getMap().getMap().getToRatio() == 1)
244           {
245             mappedrefs.addElement(refs[d]); // add translated protein maps
246           }
247         }
248         dnarefs = new DBRefEntry[mappedrefs.size()];
249         mappedrefs.copyInto(dnarefs);
250         for (int d = 0; d < dnarefs.length; d++)
251         {
252           Mapping mp = dnarefs[d].getMap();
253           StringBuffer sqstr = new StringBuffer();
254           if (mp != null)
255           {
256             Mapping intersect = mp.intersectVisContigs(viscontigs);
257             // generate seqstring for this sequence based on mapping
258
259             if (sqstr.length() > alwidth)
260               alwidth = sqstr.length();
261             cdnasqs.addElement(sqstr.toString());
262             cdnasqi.addElement(dna);
263             cdnaprod.addElement(intersect);
264           }
265         }
266       }
267       SequenceI[] cdna = new SequenceI[cdnasqs.size()];
268       DBRefEntry[] prods = new DBRefEntry[cdnaprod.size()];
269       String[] xons = new String[cdnasqs.size()];
270       cdnasqs.copyInto(xons);
271       cdnaprod.copyInto(prods);
272       cdnasqi.copyInto(cdna);
273       return CdnaTranslate(cdna, xons, prods, viscontigs, gapCharacter,
274               null, alwidth, dataset);
275     }
276     return null;
277   }
278
279   /**
280    * translate na alignment annotations onto translated amino acid alignment al
281    * using codon mapping codons
282    * 
283    * @param annotations
284    * @param al
285    * @param codons
286    */
287   public static void translateAlignedAnnotations(
288           AlignmentAnnotation[] annotations, AlignmentI al,
289           AlignedCodonFrame codons)
290   {
291     // //////////////////////////////
292     // Copy annotations across
293     //
294     // Can only do this for columns with consecutive codons, or where
295     // annotation is sequence associated.
296
297     int pos, a, aSize;
298     if (annotations != null)
299     {
300       for (int i = 0; i < annotations.length; i++)
301       {
302         // Skip any autogenerated annotation
303         if (annotations[i].autoCalculated)
304         {
305           continue;
306         }
307
308         aSize = codons.getaaWidth(); // aa alignment width.
309         jalview.datamodel.Annotation[] anots = (annotations[i].annotations == null) ? null
310                 : new jalview.datamodel.Annotation[aSize];
311         if (anots != null)
312         {
313           for (a = 0; a < aSize; a++)
314           {
315             // process through codon map.
316             if (codons.codons[a] != null
317                     && codons.codons[a][0] == (codons.codons[a][2] - 2))
318             {
319               anots[a] = getCodonAnnotation(codons.codons[a],
320                       annotations[i].annotations);
321             }
322           }
323         }
324
325         jalview.datamodel.AlignmentAnnotation aa = new jalview.datamodel.AlignmentAnnotation(
326                 annotations[i].label, annotations[i].description, anots);
327         aa.graph = annotations[i].graph;
328         aa.graphGroup = annotations[i].graphGroup;
329         aa.graphHeight = annotations[i].graphHeight;
330         if (annotations[i].getThreshold() != null)
331         {
332           aa.setThreshold(new jalview.datamodel.GraphLine(annotations[i]
333                   .getThreshold()));
334         }
335         if (annotations[i].hasScore)
336         {
337           aa.setScore(annotations[i].getScore());
338         }
339         if (annotations[i].sequenceRef != null)
340         {
341           SequenceI aaSeq = codons
342                   .getAaForDnaSeq(annotations[i].sequenceRef);
343           if (aaSeq != null)
344           {
345             // aa.compactAnnotationArray(); // throw away alignment annotation
346             // positioning
347             aa.setSequenceRef(aaSeq);
348             aa.createSequenceMapping(aaSeq, aaSeq.getStart(), true); // rebuild
349             // mapping
350             aa.adjustForAlignment();
351             aaSeq.addAlignmentAnnotation(aa);
352           }
353
354         }
355         al.addAnnotation(aa);
356       }
357     }
358   }
359
360   private static Annotation getCodonAnnotation(int[] is,
361           Annotation[] annotations)
362   {
363     // Have a look at all the codon positions for annotation and put the first
364     // one found into the translated annotation pos.
365     int contrib = 0;
366     Annotation annot = null;
367     for (int p = 0; p < 3; p++)
368     {
369       if (annotations[is[p]] != null)
370       {
371         if (annot == null)
372         {
373           annot = new Annotation(annotations[is[p]]);
374           contrib = 1;
375         }
376         else
377         {
378           // merge with last
379           Annotation cpy = new Annotation(annotations[is[p]]);
380           if (annot.colour == null)
381           {
382             annot.colour = cpy.colour;
383           }
384           if (annot.description == null || annot.description.length() == 0)
385           {
386             annot.description = cpy.description;
387           }
388           if (annot.displayCharacter == null)
389           {
390             annot.displayCharacter = cpy.displayCharacter;
391           }
392           if (annot.secondaryStructure == 0)
393           {
394             annot.secondaryStructure = cpy.secondaryStructure;
395           }
396           annot.value += cpy.value;
397           contrib++;
398         }
399       }
400     }
401     if (contrib > 1)
402     {
403       annot.value /= (float) contrib;
404     }
405     return annot;
406   }
407
408   /**
409    * Translate a na sequence
410    * 
411    * @param selection
412    *          sequence displayed under viscontigs visible columns
413    * @param seqstring
414    *          ORF read in some global alignment reference frame
415    * @param viscontigs
416    *          mapping from global reference frame to visible seqstring ORF read
417    * @param codons
418    *          Definition of global ORF alignment reference frame
419    * @param gapCharacter
420    * @param newSeq
421    * @return sequence ready to be added to alignment.
422    */
423   public static SequenceI translateCodingRegion(SequenceI selection,
424           String seqstring, int[] viscontigs, AlignedCodonFrame codons,
425           char gapCharacter, DBRefEntry product)
426   {
427     java.util.List skip = new ArrayList();
428     int skipint[] = null;
429     ShiftList vismapping = new ShiftList(); // map from viscontigs to seqstring
430     // intervals
431     int vc, scontigs[] = new int[viscontigs.length];
432     int npos = 0;
433     for (vc = 0; vc < viscontigs.length; vc += 2)
434     {
435       if (vc == 0)
436       {
437         vismapping.addShift(npos, viscontigs[vc]);
438       }
439       else
440       {
441         // hidden region
442         vismapping.addShift(npos, viscontigs[vc] - viscontigs[vc - 1] + 1);
443       }
444       scontigs[vc] = viscontigs[vc];
445       scontigs[vc + 1] = viscontigs[vc + 1];
446     }
447
448     StringBuffer protein = new StringBuffer();
449     String seq = seqstring.replace('U', 'T');
450     char codon[] = new char[3];
451     int cdp[] = new int[3], rf = 0, lastnpos = 0, nend;
452     int aspos = 0;
453     int resSize = 0;
454     for (npos = 0, nend = seq.length(); npos < nend; npos++)
455     {
456       if (!jalview.util.Comparison.isGap(seq.charAt(npos)))
457       {
458         cdp[rf] = npos; // store position
459         codon[rf++] = seq.charAt(npos); // store base
460       }
461       // filled an RF yet ?
462       if (rf == 3)
463       {
464         String aa = ResidueProperties.codonTranslate(new String(codon));
465         rf = 0;
466         if (aa == null)
467         {
468           aa = String.valueOf(gapCharacter);
469           if (skipint == null)
470           {
471             skipint = new int[]
472             { cdp[0], cdp[2] };
473           }
474           skipint[1] = cdp[2];
475         }
476         else
477         {
478           if (skipint != null)
479           {
480             // edit scontigs
481             skipint[0] = vismapping.shift(skipint[0]);
482             skipint[1] = vismapping.shift(skipint[1]);
483             for (vc = 0; vc < scontigs.length; vc += 2)
484             {
485               if (scontigs[vc + 1] < skipint[0])
486               {
487                 continue;
488               }
489               if (scontigs[vc] <= skipint[0])
490               {
491                 if (skipint[0] == scontigs[vc])
492                 {
493
494                 }
495                 else
496                 {
497                   int[] t = new int[scontigs.length + 2];
498                   System.arraycopy(scontigs, 0, t, 0, vc - 1);
499                   // scontigs[vc]; //
500                 }
501               }
502             }
503             skip.add(skipint);
504             skipint = null;
505           }
506           if (aa.equals("STOP"))
507           {
508             aa = "X";
509           }
510           resSize++;
511         }
512         // insert/delete gaps prior to this codon - if necessary
513         boolean findpos = true;
514         while (findpos)
515         {
516           // first ensure that the codons array is long enough.
517           codons.checkCodonFrameWidth(aspos);
518           // now check to see if we place the aa at the current aspos in the
519           // protein alignment
520           switch (Dna.compare_codonpos(cdp, codons.codons[aspos]))
521           {
522           case -1:
523             codons.insertAAGap(aspos, gapCharacter);
524             findpos = false;
525             break;
526           case +1:
527             // this aa appears after the aligned codons at aspos, so prefix it
528             // with a gap
529             aa = "" + gapCharacter + aa;
530             aspos++;
531             // if (aspos >= codons.aaWidth)
532             // codons.aaWidth = aspos + 1;
533             break; // check the next position for alignment
534           case 0:
535             // codon aligns at aspos position.
536             findpos = false;
537           }
538         }
539         // codon aligns with all other sequence residues found at aspos
540         protein.append(aa);
541         lastnpos = npos;
542         if (codons.codons[aspos] == null)
543         {
544           // mark this column as aligning to this aligned reading frame
545           codons.codons[aspos] = new int[]
546           { cdp[0], cdp[1], cdp[2] };
547         }
548         if (aspos >= codons.aaWidth)
549         {
550           // update maximum alignment width
551           // (we can do this without calling checkCodonFrameWidth because it was
552           // already done above)
553           codons.setAaWidth(aspos);
554         }
555         // ready for next translated reading frame alignment position (if any)
556         aspos++;
557       }
558     }
559     if (resSize > 0)
560     {
561       SequenceI newseq = new Sequence(selection.getName(),
562               protein.toString());
563       if (rf != 0)
564       {
565         jalview.bin.Cache.log
566                 .debug("trimming contigs for incomplete terminal codon.");
567         // map and trim contigs to ORF region
568         vc = scontigs.length - 1;
569         lastnpos = vismapping.shift(lastnpos); // place npos in context of
570         // whole dna alignment (rather
571         // than visible contigs)
572         // incomplete ORF could be broken over one or two visible contig
573         // intervals.
574         while (vc >= 0 && scontigs[vc] > lastnpos)
575         {
576           if (vc > 0 && scontigs[vc - 1] > lastnpos)
577           {
578             vc -= 2;
579           }
580           else
581           {
582             // correct last interval in list.
583             scontigs[vc] = lastnpos;
584           }
585         }
586
587         if (vc > 0 && (vc + 1) < scontigs.length)
588         {
589           // truncate map list to just vc elements
590           int t[] = new int[vc + 1];
591           System.arraycopy(scontigs, 0, t, 0, vc + 1);
592           scontigs = t;
593         }
594         if (vc <= 0)
595           scontigs = null;
596       }
597       if (scontigs != null)
598       {
599         npos = 0;
600         // map scontigs to actual sequence positions on selection
601         for (vc = 0; vc < scontigs.length; vc += 2)
602         {
603           scontigs[vc] = selection.findPosition(scontigs[vc]); // not from 1!
604           scontigs[vc + 1] = selection.findPosition(scontigs[vc + 1]); // exclusive
605           if (scontigs[vc + 1] == selection.getEnd())
606             break;
607         }
608         // trim trailing empty intervals.
609         if ((vc + 2) < scontigs.length)
610         {
611           int t[] = new int[vc + 2];
612           System.arraycopy(scontigs, 0, t, 0, vc + 2);
613           scontigs = t;
614         }
615         /*
616          * delete intervals in scontigs which are not translated. 1. map skip
617          * into sequence position intervals 2. truncate existing ranges and add
618          * new ranges to exclude untranslated regions. if (skip.size()>0) {
619          * Vector narange = new Vector(); for (vc=0; vc<scontigs.length; vc++) {
620          * narange.addElement(new int[] {scontigs[vc]}); } int sint=0,iv[]; vc =
621          * 0; while (sint<skip.size()) { skipint = (int[]) skip.elementAt(sint);
622          * do { iv = (int[]) narange.elementAt(vc); if (iv[0]>=skipint[0] &&
623          * iv[0]<=skipint[1]) { if (iv[0]==skipint[0]) { // delete beginning of
624          * range } else { // truncate range and create new one if necessary iv =
625          * (int[]) narange.elementAt(vc+1); if (iv[0]<=skipint[1]) { // truncate
626          * range iv[0] = skipint[1]; } else { } } } else if (iv[0]<skipint[0]) {
627          * iv = (int[]) narange.elementAt(vc+1); } } while (iv[0]) } }
628          */
629         MapList map = new MapList(scontigs, new int[]
630         { 1, resSize }, 3, 1);
631
632         // update newseq as if it was generated as mapping from product
633
634         if (product != null)
635         {
636           newseq.setName(product.getSource() + "|"
637                   + product.getAccessionId());
638           if (product.getMap() != null)
639           {
640             // Mapping mp = product.getMap();
641             // newseq.setStart(mp.getPosition(scontigs[0]));
642             // newseq.setEnd(mp
643             // .getPosition(scontigs[scontigs.length - 1]));
644           }
645         }
646         transferCodedFeatures(selection, newseq, map, null, null);
647         SequenceI rseq = newseq.deriveSequence(); // construct a dataset
648         // sequence for our new
649         // peptide, regardless.
650         // store a mapping (this actually stores a mapping between the dataset
651         // sequences for the two sequences
652         codons.addMap(selection, rseq, map);
653         return rseq;
654       }
655     }
656     // register the mapping somehow
657     //
658     return null;
659   }
660
661   /**
662    * Given a peptide newly translated from a dna sequence, copy over and set any
663    * features on the peptide from the DNA. If featureTypes is null, all features
664    * on the dna sequence are searched (rather than just the displayed ones), and
665    * similarly for featureGroups.
666    * 
667    * @param dna
668    * @param pep
669    * @param map
670    * @param featureTypes
671    *          hash who's keys are the displayed feature type strings
672    * @param featureGroups
673    *          hash where keys are feature groups and values are Boolean objects
674    *          indicating if they are displayed.
675    */
676   private static void transferCodedFeatures(SequenceI dna, SequenceI pep,
677           MapList map, Hashtable featureTypes, Hashtable featureGroups)
678   {
679     SequenceFeature[] sf = dna.getDatasetSequence().getSequenceFeatures();
680     Boolean fgstate;
681     jalview.datamodel.DBRefEntry[] dnarefs = jalview.util.DBRefUtils
682             .selectRefs(dna.getDBRef(),
683                     jalview.datamodel.DBRefSource.DNACODINGDBS);
684     if (dnarefs != null)
685     {
686       // intersect with pep
687       for (int d = 0; d < dnarefs.length; d++)
688       {
689         Mapping mp = dnarefs[d].getMap();
690         if (mp != null)
691         {
692         }
693       }
694     }
695     if (sf != null)
696     {
697       for (int f = 0; f < sf.length; f++)
698       {
699         fgstate = (featureGroups == null) ? null : ((Boolean) featureGroups
700                 .get(sf[f].featureGroup));
701         if ((featureTypes == null || featureTypes.containsKey(sf[f]
702                 .getType())) && (fgstate == null || fgstate.booleanValue()))
703         {
704           if (FeatureProperties.isCodingFeature(null, sf[f].getType()))
705           {
706             // if (map.intersectsFrom(sf[f].begin, sf[f].end))
707             {
708
709             }
710           }
711         }
712       }
713     }
714   }
715 }