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