src/jalview/analysis/Dna.java

   1 package jalview.analysis;\r
   2 \r
   3 import java.util.Hashtable;\r
   4 import java.util.Vector;\r
   5 \r
   6 import jalview.datamodel.Alignment;\r
   7 import jalview.datamodel.AlignmentAnnotation;\r
   8 import jalview.datamodel.AlignmentI;\r
   9 import jalview.datamodel.Annotation;\r
  10 import jalview.datamodel.ColumnSelection;\r
  11 import jalview.datamodel.Sequence;\r
  12 import jalview.datamodel.SequenceFeature;\r
  13 import jalview.datamodel.SequenceI;\r
  14 import jalview.schemes.ResidueProperties;\r
  15 import jalview.util.MapList;\r
  16 \r
  17 public class Dna {\r
  18   /**\r
  19    * \r
  20    * @param cdp1\r
  21    * @param cdp2\r
  22    * @return -1 if cdp1 aligns before cdp2, 0 if in the same column or cdp2 is null, +1 if after cdp2\r
  23    */\r
  24   private static int compare_codonpos(int[] cdp1, int[] cdp2) {\r
  25     if (cdp2==null || (cdp1[0]==cdp2[0] && cdp1[1] == cdp2[1] && cdp1[2] == cdp2[2]))\r
  26       return 0;\r
  27     if (cdp1[0]<cdp2[0] || cdp1[1]<cdp2[1] || cdp1[2]<cdp2[2])\r
  28       return -1; // one base in cdp1 precedes the corresponding base in the other codon\r
  29     return 1; // one base in cdp1 appears after the corresponding base in the other codon.\r
  30   }\r
  31   /**\r
  32    * create a new alignment of protein sequences\r
  33    * by an inframe translation of the provided NA sequences\r
  34    * @param selection\r
  35    * @param seqstring\r
  36    * @param viscontigs\r
  37    * @param gapCharacter\r
  38    * @param annotations\r
  39    * @param aWidth\r
  40    * @return\r
  41    */\r
  42   public static AlignmentI CdnaTranslate(SequenceI[] selection, String[] seqstring, int viscontigs[], char gapCharacter, \r
  43       AlignmentAnnotation[] annotations, int aWidth) {\r
  44     int s, sSize = selection.length;\r
  45     SequenceI [] newSeq = new SequenceI[sSize];\r
  46     int res, resSize;\r
  47     StringBuffer protein;\r
  48     String seq;\r
  49 \r
  50     int[][] codons = new int[aWidth][]; // stores hash of subsequent positions for each codon start position in alignment\r
  51 \r
  52     for (res=0;res<aWidth;res++)\r
  53       codons[res]=null;\r
  54     int aslen=0; // final width of aligned translated aa sequences\r
  55     for(s=0; s<sSize; s++)\r
  56     {\r
  57       int vc,scontigs[]=new int[viscontigs.length];\r
  58 \r
  59       for (vc=0;vc<scontigs.length; vc+=2)\r
  60       {\r
  61         scontigs[vc]=selection[s].findPosition(viscontigs[vc]); // not from 1!\r
  62         scontigs[vc+1]=selection[s].findPosition(viscontigs[vc+1]-1); // exclusive\r
  63         if (scontigs[vc+1]==selection[s].getEnd())\r
  64           break;\r
  65       }\r
  66       if ((vc+2)<scontigs.length) {\r
  67         int t[] = new int[vc+2];\r
  68         System.arraycopy(scontigs, 0, t, 0, vc+2);\r
  69         scontigs = t;\r
  70       }\r
  71       protein = new StringBuffer();\r
  72       seq = seqstring[s].replace('U', 'T');\r
  73       char codon[]=new char[3];\r
  74       int cdp[]=new int[3],rf=0,gf=0,nend,npos;\r
  75       int aspos=0;\r
  76       resSize=0;\r
  77       for (npos=0,nend=seq.length(); npos<nend; npos++) {\r
  78         if (!jalview.util.Comparison.isGap(seq.charAt(npos))) { \r
  79           cdp[rf] = npos; // store position\r
  80           codon[rf++]=seq.charAt(npos); // store base\r
  81         }\r
  82         // filled an RF yet ?\r
  83         if (rf==3) {\r
  84           String aa = ResidueProperties.codonTranslate(new String(codon));\r
  85           rf=0;\r
  86           if(aa==null)\r
  87             aa=String.valueOf(gapCharacter);\r
  88           else {\r
  89             if(aa.equals("STOP"))\r
  90             {\r
  91               aa="X";\r
  92             }\r
  93             resSize++;\r
  94           }\r
  95           // insert/delete gaps prior to this codon - if necessary\r
  96           boolean findpos=true;\r
  97           while (findpos) \r
  98           {\r
  99             // first ensure that the codons array is long enough.\r
 100             if (codons.length<=aslen+1) {\r
 101               // probably never have to do this ?\r
 102               int[][] c = new int[codons.length+10][];\r
 103               for (int i=0; i<codons.length; i++) {\r
 104                 c[i] = codons[i];\r
 105                 codons[i]=null;\r
 106               }\r
 107               codons = c;\r
 108             }\r
 109             // now check to see if we place the aa at the current aspos in the protein alignment\r
 110             switch (Dna.compare_codonpos(cdp, codons[aspos])) \r
 111             {\r
 112             case -1:\r
 113               // this aa appears before the aligned codons at aspos - so shift them.\r
 114               aslen++;\r
 115               for (int sq=0;sq<s; sq++) {\r
 116                 newSeq[sq].insertCharAt(aspos, gapCharacter);\r
 117               }\r
 118               System.arraycopy(codons, aspos, codons, aspos+1, aslen-aspos);\r
 119               codons[aspos]=null; // clear so new codon position can be marked.\r
 120               findpos=false;\r
 121               break;\r
 122             case +1:\r
 123               // this aa appears after the aligned codons at aspos, so prefix it with a gap\r
 124               aa = ""+gapCharacter+aa;\r
 125               aspos++;\r
 126               if (aspos>=aslen)\r
 127                 aslen=aspos+1;\r
 128               break; // check the next position for alignment\r
 129             case 0:\r
 130               // codon aligns at aspos position.\r
 131               findpos = false;\r
 132             }\r
 133           }\r
 134           // codon aligns with all other sequence residues found at aspos\r
 135           protein.append(aa);\r
 136           if (codons[aspos]==null) \r
 137           {\r
 138             // mark this column as aligning to this aligned reading frame \r
 139             codons[aspos] = new int[] { cdp[0], cdp[1], cdp[2] };\r
 140           }\r
 141           aspos++;\r
 142           if (aspos>=aslen)\r
 143             aslen=aspos+1;\r
 144         }\r
 145       }\r
 146       if (resSize>0) \r
 147       {\r
 148         newSeq[s] = new Sequence(selection[s].getName(),\r
 149             protein.toString());\r
 150         if (rf!=0) \r
 151         {\r
 152           jalview.bin.Cache.log.debug("trimming contigs for incomplete terminal codon.");\r
 153           // trim contigs\r
 154           vc=scontigs.length-1;\r
 155           nend-=rf;\r
 156           // incomplete ORF could be broken over one or two visible contig intervals.\r
 157           while (vc>0 && scontigs[vc]>nend)\r
 158           {\r
 159             if (scontigs[vc-1]>nend) \r
 160             {\r
 161               vc-=2;\r
 162             } else {\r
 163               // correct last interval in list.\r
 164               scontigs[vc]=nend;\r
 165             }\r
 166           }\r
 167           if ((vc+2)<scontigs.length) {\r
 168             // truncate map list\r
 169             int t[] = new int[vc+1];\r
 170             System.arraycopy(scontigs,0,t,0,vc+1);\r
 171             scontigs=t;\r
 172           }\r
 173         }\r
 174         MapList map = new MapList(scontigs, new int[] { 1, resSize },3,1); // TODO: store mapping on newSeq for linked DNA/Protein viewing.\r
 175       }\r
 176       // register the mapping somehow\r
 177       // \r
 178     }\r
 179     if (aslen==0)\r
 180       return null;\r
 181     AlignmentI al = new Alignment(newSeq);\r
 182     al.padGaps();  // ensure we look aligned.\r
 183     al.setDataset(null);\r
 184 \r
 185 \r
 186     ////////////////////////////////\r
 187     // Copy annotations across\r
 188     //\r
 189     // Can only do this for columns with consecutive codons, or where\r
 190     // annotation is sequence associated.\r
 191     \r
 192     int pos,a,aSize;\r
 193     if(annotations!=null)\r
 194     {\r
 195       for (int i = 0; i < annotations.length; i++)\r
 196       {\r
 197         // Skip any autogenerated annotation\r
 198         if (annotations[i].autoCalculated) {\r
 199           continue;\r
 200         }\r
 201         \r
 202         aSize = aslen; // aa alignment width.\r
 203         jalview.datamodel.Annotation[] anots = \r
 204           (annotations[i].annotations==null) \r
 205           ? null :\r
 206             new jalview.datamodel.Annotation[aSize];\r
 207         if (anots!=null)\r
 208         {\r
 209           for (a = 0; a < aSize; a++)\r
 210           {\r
 211             // process through codon map.\r
 212             if (codons[a]!=null && codons[a][0]==(codons[a][2]-2))\r
 213             {\r
 214               pos = codons[a][0];\r
 215               if (annotations[i].annotations[pos] == null\r
 216                       || annotations[i].annotations[pos] == null)\r
 217                 continue;\r
 218             \r
 219               anots[a] = new Annotation(annotations[i].annotations[pos]);\r
 220             }\r
 221           }\r
 222         }\r
 223 \r
 224         jalview.datamodel.AlignmentAnnotation aa\r
 225         = new jalview.datamodel.AlignmentAnnotation(annotations[i].label,\r
 226             annotations[i].description, anots);\r
 227         if (annotations[i].hasScore)\r
 228         {\r
 229           aa.setScore(annotations[i].getScore());\r
 230         }\r
 231         al.addAnnotation(aa);\r
 232       }\r
 233     }\r
 234     return al;\r
 235   }\r
 236 }\r