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               findpos=false;\r
 120               break;\r
 121             case +1:\r
 122               // this aa appears after the aligned codons at aspos, so prefix it with a gap\r
 123               aa = ""+gapCharacter+aa;\r
 124               aspos++;\r
 125               if (aspos>=aslen)\r
 126                 aslen=aspos+1;\r
 127               break; // check the next position for alignment\r
 128             case 0:\r
 129               // codon aligns at aspos position.\r
 130               findpos = false;\r
 131             }\r
 132           }\r
 133           // codon aligns with all other sequence residues found at aspos\r
 134           protein.append(aa);\r
 135           if (codons[aspos]==null) \r
 136           {\r
 137             // mark this column as aligning to this aligned reading frame \r
 138             codons[aspos] = new int[] { cdp[0], cdp[1], cdp[2] };\r
 139           }\r
 140           aspos++;\r
 141           if (aspos>=aslen)\r
 142             aslen=aspos+1;\r
 143         }\r
 144       }\r
 145       if (resSize>0) \r
 146       {\r
 147         newSeq[s] = new Sequence(selection[s].getName(),\r
 148             protein.toString());\r
 149         if (rf!=0) \r
 150         {\r
 151           jalview.bin.Cache.log.debug("trimming contigs for incomplete terminal codon.");\r
 152           // trim contigs\r
 153           vc=scontigs.length-1;\r
 154           nend-=rf;\r
 155           // incomplete ORF could be broken over one or two visible contig intervals.\r
 156           while (vc>0 && scontigs[vc]>nend)\r
 157           {\r
 158             if (scontigs[vc-1]>nend) \r
 159             {\r
 160               vc-=2;\r
 161             } else {\r
 162               // correct last interval in list.\r
 163               scontigs[vc]=nend;\r
 164             }\r
 165           }\r
 166           if ((vc+2)<scontigs.length) {\r
 167             // truncate map list\r
 168             int t[] = new int[vc+1];\r
 169             System.arraycopy(scontigs,0,t,0,vc+1);\r
 170             scontigs=t;\r
 171           }\r
 172         }\r
 173         MapList map = new MapList(scontigs, new int[] { 1, resSize },3,1); // TODO: store mapping on newSeq for linked DNA/Protein viewing.\r
 174       }\r
 175       // register the mapping somehow\r
 176       // \r
 177     }\r
 178     if (aslen==0)\r
 179       return null;\r
 180     AlignmentI al = new Alignment(newSeq);\r
 181     al.padGaps();  // ensure we look aligned.\r
 182     al.setDataset(null);\r
 183 \r
 184 \r
 185     ////////////////////////////////\r
 186     // Copy annotations across\r
 187     //\r
 188     // Can only do this for columns with consecutive codons, or where\r
 189     // annotation is sequence associated.\r
 190     \r
 191     int pos,a,aSize;\r
 192     if(annotations!=null)\r
 193     {\r
 194       for (int i = 0; i < annotations.length; i++)\r
 195       {\r
 196         // Skip any autogenerated annotation\r
 197         if (annotations[i].autoCalculated) {\r
 198           continue;\r
 199         }\r
 200 \r
 201         aSize = aslen; // aa alignment width.\r
 202         jalview.datamodel.Annotation[] anots =\r
 203           new jalview.datamodel.Annotation[aSize];\r
 204         \r
 205         for (a = 0; a < aSize; a++)\r
 206         {\r
 207           // process through codon map.\r
 208           if (codons[a]!=null && codons[a][0]==(codons[a][2]-2))\r
 209           {\r
 210             pos = codons[a][0];\r
 211             if (annotations[i].annotations[pos] == null\r
 212                 || annotations[i].annotations[pos] == null)\r
 213               continue;\r
 214             \r
 215             anots[a] = new Annotation(annotations[i].annotations[pos]);\r
 216           }\r
 217         }\r
 218 \r
 219         jalview.datamodel.AlignmentAnnotation aa\r
 220         = new jalview.datamodel.AlignmentAnnotation(annotations[i].label,\r
 221             annotations[i].description, anots);\r
 222         al.addAnnotation(aa);\r
 223       }\r
 224     }\r
 225     return al;\r
 226   }\r
 227 }\r