src/jalview/analysis/AlignmentUtil.java

   1 package jalview.analysis;\r
   2 \r
   3 import jalview.datamodel.*;\r
   4 import jalview.io.*;\r
   5 import jalview.util.*;\r
   6 \r
   7 import java.util.*;\r
   8 \r
   9 public class AlignmentUtil {\r
  10 \r
  11   private AlignmentUtil() {\r
  12   }\r
  13 \r
  14   public static int[][] percentIdentity2(AlignmentI align) {\r
  15     return percentIdentity2(align,0,align.getWidth()-1);\r
  16   }\r
  17 \r
  18   public static int[][] percentIdentity2(AlignmentI align, int start, int end) {\r
  19     int [][] cons2 = new int[align.getWidth()][24];\r
  20     // Initialize the array\r
  21     for (int j=0;j<24;j++) {\r
  22       for (int i=0; i < align.getWidth();i++) {\r
  23         cons2[i][j] = 0;\r
  24       }\r
  25     }\r
  26 \r
  27     return cons2;\r
  28   }\r
  29   public static int getPixelHeight(int i, int j,int charHeight) {\r
  30     int h=0;\r
  31     while (i < j) {\r
  32       h += charHeight;\r
  33       i++;\r
  34     }\r
  35     return h;\r
  36   }\r
  37 \r
  38 \r
  39     public static Vector substitution_rates (AlignmentI align, int start, int end) {\r
  40 \r
  41         Vector rates = new Vector();\r
  42 \r
  43         int len = (end-start+1);\r
  44 \r
  45         // Turn seqs into char arrays\r
  46 \r
  47         int[][] seqint = new int[align.getHeight()][len];\r
  48 \r
  49 \r
  50         for (int j = 0; j < align.getHeight(); j++) {\r
  51 \r
  52             SequenceI seq = align.getSequenceAt(j);\r
  53 \r
  54             for (int i = 0 ; i < len; i++) {\r
  55                 char c = seq.getCharAt(start + i - 1);\r
  56 \r
  57                 if (c == 'A') {\r
  58                     seqint[j][i] = 0;\r
  59                 } else if (c == 'C') {\r
  60                     seqint[j][i] = 1;\r
  61                 } else if (c == 'T') {\r
  62                     seqint[j][i] = 2;\r
  63                 } else if (c == 'G') {\r
  64                     seqint[j][i] = 3;\r
  65                 } else {\r
  66                     seqint[j][i] = -1;\r
  67                 }\r
  68 \r
  69             }\r
  70 \r
  71         }\r
  72 \r
  73 \r
  74         //      print_matrix(seqint,2,len);     for (int j = 0; j < align.getHeight(); j++) {\r
  75 \r
  76             for (int j = 0; j < align.getHeight(); j++) {\r
  77 \r
  78                 for (int k = 0; k < align.getHeight(); k++) {\r
  79 \r
  80                 int counts[][] = new int[4][4];\r
  81                 int tot  = 0;\r
  82                 int tots[] = new int[4];\r
  83                 int fulltot = 0;\r
  84                 int fulltots[] = new int[4];\r
  85 \r
  86                 for (int i = 0 ; i < len; i++) {\r
  87 \r
  88                     if (k != j) {\r
  89 \r
  90                         //                      System.out.println("Seq " + j + " " + k + " " + i + " " + seqint[j][i] + " " + seqint[k][i]);\r
  91                         if (seqint[j][i] >= 0 &&\r
  92                             seqint[k][i] >= 0) {\r
  93                             counts[seqint[k][i]][seqint[j][i]]++;\r
  94 \r
  95                             //  print_matrix(counts,4,4);\r
  96                             tots[seqint[j][i]]++;\r
  97                             tot++;\r
  98                         }\r
  99                         if (seqint[j][i] != -1) {\r
 100                             fulltots[seqint[j][i]]++;\r
 101                             fulltot++;\r
 102                         }\r
 103                     }\r
 104                 }\r
 105 \r
 106                 if (k != j) {\r
 107 \r
 108                     System.out.println();\r
 109 \r
 110                     System.out.println("Sequence " + align.getSequenceAt(j).getName() + " "  + align.getSequenceAt(k).getName());\r
 111 \r
 112                     System.out.println();\r
 113                     print_matrix(counts,4,4);\r
 114                     System.out.println();\r
 115 \r
 116 \r
 117 \r
 118                     double[][] out = new double[4][4];// = constant_multiply_matrix(counts,1.0/tot,4,4);\r
 119 \r
 120                     for (int i = 0; i < 4; i++) {\r
 121                         for (int jj = 0; jj < 4; jj++) {\r
 122                             out[i][jj] = (double)counts[i][jj]/tots[jj];\r
 123                         }\r
 124                     }\r
 125 \r
 126                     print_matrix(out,4,4);\r
 127                     System.out.println();\r
 128 \r
 129 \r
 130                     System.out.print("RATES\t");\r
 131                     System.out.print(align.getSequenceAt(j).getName() + "\t"  + align.getSequenceAt(k).getName() + "\t");\r
 132 \r
 133                     for (int i = 0; i < 4; i++) {\r
 134                         for (int jj = 0; jj < 4; jj++) {\r
 135                             Format.print(System.out,"%4.3f\t",out[i][jj]);\r
 136                         }\r
 137                     }\r
 138                     System.out.println();\r
 139 \r
 140                     for (int i = 0; i < 4; i++) {\r
 141                         Format.print(System.out,"%4.3f\t",(double)fulltots[i]*1.0/fulltot);\r
 142                     }\r
 143 \r
 144                     System.out.println();\r
 145                     System.out.print("\nGC ");\r
 146 \r
 147                     Format.print(System.out,"%4.3f\t",(double)(100*fulltots[1]+fulltots[3])/fulltot);\r
 148 \r
 149                     System.out.print((fulltots[1]+fulltots[3]) + "\t" + fulltot);\r
 150 \r
 151 \r
 152                     System.out.println();\r
 153 \r
 154                     rates.addElement(out);\r
 155                 }\r
 156 \r
 157             }\r
 158         }\r
 159         return rates;\r
 160     }\r
 161 \r
 162     public static double[][]  constant_multiply_matrix(int[][] matrix, double c,int n, int m) {\r
 163         double[][] out = new double[n][m];\r
 164 \r
 165         for (int i = 0; i < n; i++) {\r
 166             for (int j = 0; j < m; j++) {\r
 167 \r
 168                 out[i][j] = matrix[i][j]*c;\r
 169             }\r
 170         }\r
 171         return out;\r
 172     }\r
 173 \r
 174 \r
 175     public static void print_matrix(int[][] matrix, int n, int m) {\r
 176 \r
 177 \r
 178         for (int i = 0; i < n; i++) {\r
 179             for (int j = 0; j < m; j++) {\r
 180 \r
 181                 System.out.print(matrix[i][j] + "\t");\r
 182             }\r
 183             System.out.println();\r
 184         }\r
 185     }\r
 186     public static void print_matrix(double[][] matrix, int n, int m) {\r
 187 \r
 188 \r
 189         for (int i = 0; i < n; i++) {\r
 190             for (int j = 0; j < m; j++) {\r
 191 \r
 192                 Format.print(System.out,"%4.3f\t",matrix[i][j]);\r
 193 \r
 194             }\r
 195             System.out.println();\r
 196         }\r
 197     }\r
 198 \r
 199     public static Hashtable findKmers(SequenceI seq, int start, int end, Vector kmers) {\r
 200 \r
 201          Hashtable pos = new Hashtable();\r
 202 \r
 203          for (int j = 0; j < kmers.size(); j++) {\r
 204 \r
 205            String   kmer    = ((Sequence)kmers.elementAt(j)).getSequence();\r
 206            Sequence kmerseq = (Sequence)kmers.elementAt(j);\r
 207 \r
 208            if (end < seq.getLength()) {\r
 209 \r
 210                String str = seq.getSequence(start,end);\r
 211 \r
 212                int i = 0;\r
 213 \r
 214                while (str.indexOf(kmer,i) != -1) {\r
 215 \r
 216                //               System.out.println("STring " + str + " " + i);\r
 217                    int coord = str.indexOf(kmer,i);\r
 218 \r
 219                    pos.put(new Integer(coord),new Integer(kmerseq.getLength()));\r
 220 \r
 221                    i  = coord + 1;\r
 222 \r
 223                }\r
 224            }\r
 225          }\r
 226          return pos;\r
 227     }\r
 228 \r
 229     class FeatureThread extends Thread\r
 230     {\r
 231 \r
 232     }\r
 233 \r
 234     public static void addUniprotFeatures(AlignmentI align)\r
 235     {\r
 236       EBIFetchClient ebi = new EBIFetchClient();\r
 237       int seqIndex = 0;\r
 238       Vector sequences = align.getSequences();\r
 239       SequenceI sequence;\r
 240 \r
 241       while (seqIndex < sequences.size())\r
 242       {\r
 243         StringBuffer ids = new StringBuffer("uniprot:");\r
 244         for (int i=0; seqIndex<sequences.size() && i<50; seqIndex++, i++)\r
 245         {\r
 246           sequence = (SequenceI) sequences.get(seqIndex);\r
 247           ids.append(sequence.getName() + ";");\r
 248         }\r
 249 \r
 250         String[] result = ebi.fetchData(ids.toString(), "xml", null);\r
 251         sequence = null;\r
 252         Vector features = null;\r
 253         String type, description, status, start, end, pdb = null;\r
 254 \r
 255         for (int r = 0; r < result.length; r++)\r
 256         {\r
 257           if(sequence==null && result[r].indexOf("<name>")>-1)\r
 258           {\r
 259             sequence = align.findName( parseElement( result[r], "<name>" )) ;\r
 260             features = new Vector();\r
 261             type=""; start="0"; end="0"; description=""; status=""; pdb="";\r
 262           }\r
 263 \r
 264           if(sequence==null)\r
 265             continue;\r
 266 \r
 267            if( result[r].indexOf("<property type=\"pdb accession\"")>-1)\r
 268            {\r
 269              pdb = parseValue( result[r], "value=" );\r
 270              sequence.setPDBId(pdb);\r
 271            }\r
 272 \r
 273            if(result[r].indexOf("feature type")>-1)\r
 274            {\r
 275              type = parseValue( result[r], "type=" );\r
 276              description = parseValue( result[r], "description=" );\r
 277              status = parseValue ( result[r], "status=");\r
 278 \r
 279              while( result[r].indexOf("<location>")==-1)\r
 280              {\r
 281                  r++;  //<location>\r
 282              }\r
 283              r++;\r
 284              if(result[r].indexOf("begin")>-1)\r
 285              {\r
 286                start = parseValue( result[r], "position=" );\r
 287                end = parseValue( result[++r], "position=" );\r
 288              }\r
 289              else\r
 290              {\r
 291                start = parseValue( result[r], "position=" );\r
 292                end = parseValue(   result[r], "position=" );\r
 293              }\r
 294              int sstart = Integer.parseInt(start);\r
 295              int eend = Integer.parseInt(end);\r
 296 \r
 297              if(sstart>=sequence.getStart() && eend<=sequence.getEnd())\r
 298              {\r
 299                SequenceFeature sf = new SequenceFeature(type,\r
 300                    sstart,\r
 301                    eend,\r
 302                    description,\r
 303                    status);\r
 304                features.add(sf);\r
 305              }\r
 306            }\r
 307 \r
 308            if(result[r].indexOf("</entry>")>-1)\r
 309            {\r
 310              sequence.setSequenceFeatures( features );\r
 311              features = null;\r
 312              sequence = null;\r
 313            }\r
 314         }\r
 315       }\r
 316     }\r
 317 \r
 318     static String parseValue(String line, String tag)\r
 319     {\r
 320 \r
 321       int index = line.indexOf(tag)+tag.length()+1;\r
 322       if(index==tag.length())\r
 323         return "";\r
 324 \r
 325       return line.substring( index, line.indexOf("\"", index+1) );\r
 326     }\r
 327 \r
 328 \r
 329     static String parseElement(String line, String tag)\r
 330     {\r
 331       int index = line.indexOf(tag)+tag.length();\r
 332       return line.substring( index, line.indexOf("</") ) ;\r
 333     }\r
 334 \r
 335 \r
 336 \r
 337 \r
 338 \r
 339 }\r