|
| 1 | +function calculateCodonBias() { |
| 2 | + |
| 3 | + // Script to calculate amino acid composition in reference sequences |
| 4 | + |
| 5 | + var codingFeatures = {}; |
| 6 | + var resultMap = glue.command(["list", "feature","-w", "featureMetatags.name = 'CODES_AMINO_ACIDS'"]); |
| 7 | + var featureList = resultMap["listResult"]; |
| 8 | + var codingFeatureList = featureList["row"]; |
| 9 | + _.each(codingFeatureList,function(featureObj){ |
| 10 | + |
| 11 | + //glue.log("INFO", "RESULT WAS ", featureObj); |
| 12 | + |
| 13 | + var valueArray = featureObj["value"]; |
| 14 | + var codingFeatureName = valueArray[0]; |
| 15 | + //glue.log("INFO", "NAME WAS ", featureName) |
| 16 | + codingFeatures[codingFeatureName] = featureObj; |
| 17 | + |
| 18 | + }); |
| 19 | + |
| 20 | + //glue.log("INFO", "RESULT WAS ", codingFeatures); |
| 21 | + |
| 22 | + // get list of reference sequences from GLUE |
| 23 | + var referencesResult = glue.command(["list","reference"]); |
| 24 | + //glue.log("INFO", "RESULT WAS ", referencesResult); |
| 25 | + |
| 26 | + var listResult = referencesResult["listResult"]; |
| 27 | + var referencesList = listResult["row"]; |
| 28 | + //glue.log("INFO", "RESULT WAS ", referencesList); |
| 29 | + |
| 30 | + |
| 31 | + // iterate through reference list and get AA composition of each coding feature |
| 32 | + codonCompositionResults = {} |
| 33 | + _.each(referencesList, function(refObj) { |
| 34 | + |
| 35 | + //glue.log("INFO", "RESULT WAS ", refObj); |
| 36 | + |
| 37 | + var refseqResults = {}; |
| 38 | + var referenceProperties = refObj["value"]; |
| 39 | + var referenceName = referenceProperties[0]; |
| 40 | + var sequenceID = referenceProperties[2]; |
| 41 | + |
| 42 | + //glue.log("INFO", "Reference name result was:", referenceName); |
| 43 | + |
| 44 | + // list all features annotated in this reference |
| 45 | + // GLUE COMMAND: reference [referenceName] list feature-location |
| 46 | + glue.inMode("/reference/"+referenceName, function() { |
| 47 | + |
| 48 | + var featuresResult = glue.tableToObjects(glue.command(["list", "feature-location"])); |
| 49 | + //glue.log("INFO", "RESULT WAS ", featuresResult); |
| 50 | + |
| 51 | + // iterate through features |
| 52 | + _.each(featuresResult, function(featureObj) { |
| 53 | + |
| 54 | + //glue.log("INFO", "RESULT WAS ", featureObj); |
| 55 | + |
| 56 | + var featureResults = {}; |
| 57 | + var aaCounts = {}; |
| 58 | + |
| 59 | + // get amino acid sequence |
| 60 | + var featureName = featureObj["feature.name"]; |
| 61 | + //glue.log("INFO", "Feature name result was:", featureName); |
| 62 | + |
| 63 | + // get amino acid sequence of the features |
| 64 | + if (codingFeatures[featureName]) { |
| 65 | + |
| 66 | + glue.inMode("/feature-location/"+featureName, function() { |
| 67 | + |
| 68 | + var codonSequenceResult = glue.tableToObjects(glue.command(["amino-acid"])); |
| 69 | + |
| 70 | + //glue.log("INFO", "Table result was:", codonSequenceResult); |
| 71 | + |
| 72 | + // iterate through and get codon composition of feature |
| 73 | + var length = 0; |
| 74 | + _.each(codonSequenceResult, function(codonObj) { |
| 75 | + |
| 76 | + var codon = codonObj["codonNts"]; |
| 77 | + var aa = codonObj["aminoAcid"]; |
| 78 | + |
| 79 | + // increment amino acid counts |
| 80 | + if (aaCounts[aa]) { |
| 81 | + |
| 82 | + aaCounts[aa] += 1; |
| 83 | + |
| 84 | + } else { |
| 85 | + aaCounts[aa] = 1; |
| 86 | + } |
| 87 | + |
| 88 | + length += 1; |
| 89 | + |
| 90 | + //glue.log("INFO", "Amino acid result:", aa); |
| 91 | + //glue.log("INFO", " Codon result:", codon); |
| 92 | + |
| 93 | + if (featureResults[aa]) { |
| 94 | + |
| 95 | + // get the map with codon counts for this aa |
| 96 | + aaCodonMap = featureResults[aa]; |
| 97 | + aaCodonMap[codon] += 1; |
| 98 | + |
| 99 | + } |
| 100 | + else { |
| 101 | + |
| 102 | + // create a map with all codons for this aa as keys |
| 103 | + aaCodonMap = create_codon_map(aa); |
| 104 | + // increment the map for the current codon |
| 105 | + |
| 106 | + aaCodonMap[codon] = 1; |
| 107 | + |
| 108 | + // associate this data structure with the genome feature map |
| 109 | + // by using the current aa as a key |
| 110 | + featureResults[aa] = aaCodonMap; |
| 111 | + |
| 112 | + } |
| 113 | + |
| 114 | + }); |
| 115 | + |
| 116 | + // store the counts for each amino acid |
| 117 | + featureResults["aa_counts"] = aaCounts; |
| 118 | + |
| 119 | + // store the total sequence length |
| 120 | + featureResults["length"] = length; |
| 121 | + |
| 122 | + // Store the feature name |
| 123 | + featureResults["feature_name"] = featureName; |
| 124 | + |
| 125 | + }); |
| 126 | + |
| 127 | + //glue.log("INFO", "FEATURE RESULT WAS ", featureResults); die; |
| 128 | + |
| 129 | + // store feature result |
| 130 | + refseqResults[featureName] = featureResults; |
| 131 | + refseqResults["reference_name"] = referenceName; |
| 132 | + refseqResults["sequence_id"] = sequenceID; |
| 133 | + |
| 134 | + } |
| 135 | + |
| 136 | + }); |
| 137 | + |
| 138 | + // store reference result |
| 139 | + codonCompositionResults[referenceName] = refseqResults; |
| 140 | + |
| 141 | + }); |
| 142 | + |
| 143 | + }); |
| 144 | + |
| 145 | + // glue.log("INFO", "FINAL RESULT WAS ", codonCompositionResults); |
| 146 | + var outputArray = []; |
| 147 | + |
| 148 | + // get each reference |
| 149 | + |
| 150 | + _.each(codonCompositionResults, function(refseqResults) { |
| 151 | + |
| 152 | + //glue.log("INFO", "RESULT WAS ", refseqResults); |
| 153 | + |
| 154 | + var referenceName = refseqResults["reference_name"]; |
| 155 | + var sequenceID = refseqResults["sequence_id"]; |
| 156 | + |
| 157 | + // get each feature |
| 158 | + _.each(refseqResults, function(featureResults) { |
| 159 | + |
| 160 | + var featureName = featureResults["feature_name"]; |
| 161 | + |
| 162 | + // get each amino |
| 163 | + aaCounts = featureResults["aa_counts"]; |
| 164 | + var referenceName = refseqResults["reference_name"]; |
| 165 | + _.each(_.keys(aaCounts), function(aa) { |
| 166 | + |
| 167 | + var aaCount = aaCounts[aa]; |
| 168 | + // glue.log("INFO", "Got count '"+count+"' for aa character'"+aa+"'"); |
| 169 | + |
| 170 | + var codonCounts = featureResults[aa]; |
| 171 | + //glue.log("INFO", "RESULT WAS ", codonCounts); die; |
| 172 | + _.each(_.keys(codonCounts), function(codon) { |
| 173 | + |
| 174 | + var codonCount = codonCounts[codon]; |
| 175 | + |
| 176 | + // write line: reference - feature - aa - aa count - codon - codon count |
| 177 | + outputArray.push({ |
| 178 | + |
| 179 | + referenceName: referenceName, |
| 180 | + sequenceID: sequenceID, |
| 181 | + featureName: featureName, |
| 182 | + aa: aa, |
| 183 | + aaCount: aaCount, |
| 184 | + codon: codon, |
| 185 | + codonCount: codonCount |
| 186 | + |
| 187 | + }); |
| 188 | + |
| 189 | + |
| 190 | + }); |
| 191 | + |
| 192 | + }); |
| 193 | + |
| 194 | + }); |
| 195 | + |
| 196 | + }); |
| 197 | + |
| 198 | + return outputArray; |
| 199 | +} |
| 200 | + |
| 201 | + |
| 202 | +function create_codon_map(aa){ |
| 203 | + |
| 204 | + var codonMap = {}; |
| 205 | + |
| 206 | + if (aa == 'A') { |
| 207 | + |
| 208 | + codonMap['GCT'] = 0; |
| 209 | + codonMap['GCC'] = 0; |
| 210 | + codonMap['GCA'] = 0; |
| 211 | + codonMap['GCG'] = 0; |
| 212 | + |
| 213 | + } else if (aa == 'C') { |
| 214 | + |
| 215 | + codonMap['TGT'] = 0; |
| 216 | + codonMap['TGC'] = 0; |
| 217 | + |
| 218 | + |
| 219 | + } else if (aa == 'D') { |
| 220 | + |
| 221 | + codonMap['GAT'] = 0; |
| 222 | + codonMap['GAC'] = 0; |
| 223 | + |
| 224 | + } else if (aa == 'E') { |
| 225 | + |
| 226 | + codonMap['GAA'] = 0; |
| 227 | + codonMap['GAG'] = 0; |
| 228 | + |
| 229 | + } else if (aa == 'F') { |
| 230 | + |
| 231 | + codonMap['TTT'] = 0; |
| 232 | + codonMap['TTC'] = 0; |
| 233 | + |
| 234 | + } else if (aa == 'G') { |
| 235 | + |
| 236 | + codonMap['GGT'] = 0; |
| 237 | + codonMap['GGC'] = 0; |
| 238 | + codonMap['GGA'] = 0; |
| 239 | + codonMap['GGG'] = 0; |
| 240 | + |
| 241 | + } else if (aa == 'H') { |
| 242 | + |
| 243 | + codonMap['CAT'] = 0; |
| 244 | + codonMap['CAC'] = 0; |
| 245 | + |
| 246 | + } else if (aa == 'I') { |
| 247 | + |
| 248 | + codonMap['ATT'] = 0; |
| 249 | + codonMap['ATC'] = 0; |
| 250 | + codonMap['ATA'] = 0; |
| 251 | + |
| 252 | + } else if (aa == 'K') { |
| 253 | + |
| 254 | + codonMap['AAA'] = 0; |
| 255 | + codonMap['AAG'] = 0; |
| 256 | + |
| 257 | + } else if (aa == 'L') { |
| 258 | + |
| 259 | + codonMap['CTT'] = 0; |
| 260 | + codonMap['CTC'] = 0; |
| 261 | + codonMap['CTA'] = 0; |
| 262 | + codonMap['CTG'] = 0; |
| 263 | + codonMap['TTA'] = 0; |
| 264 | + codonMap['TTG'] = 0; |
| 265 | + |
| 266 | + |
| 267 | + } else if (aa == 'M') { |
| 268 | + |
| 269 | + codonMap['ATG'] = 0; |
| 270 | + |
| 271 | + } else if (aa == 'N') { |
| 272 | + |
| 273 | + codonMap['AAT'] = 0; |
| 274 | + codonMap['AAC'] = 0; |
| 275 | + |
| 276 | + } else if (aa == 'P') { |
| 277 | + |
| 278 | + codonMap['CCT'] = 0; |
| 279 | + codonMap['CCC'] = 0; |
| 280 | + codonMap['CCA'] = 0; |
| 281 | + codonMap['CCG'] = 0; |
| 282 | + |
| 283 | + } else if (aa == 'Q') { |
| 284 | + |
| 285 | + codonMap['CAA'] = 0; |
| 286 | + codonMap['CAG'] = 0; |
| 287 | + |
| 288 | + } else if (aa == 'R') { |
| 289 | + |
| 290 | + codonMap['CGT'] = 0; |
| 291 | + codonMap['CGC'] = 0; |
| 292 | + codonMap['CGA'] = 0; |
| 293 | + codonMap['CGG'] = 0; |
| 294 | + codonMap['AGA'] = 0; |
| 295 | + codonMap['AGG'] = 0; |
| 296 | + |
| 297 | + } else if (aa == 'S') { |
| 298 | + |
| 299 | + codonMap['TCT'] = 0; |
| 300 | + codonMap['TCC'] = 0; |
| 301 | + codonMap['TCA'] = 0; |
| 302 | + codonMap['TCG'] = 0; |
| 303 | + codonMap['AGT'] = 0; |
| 304 | + codonMap['AGC'] = 0; |
| 305 | + |
| 306 | + } else if (aa == 'T') { |
| 307 | + |
| 308 | + codonMap['ACT'] = 0; |
| 309 | + codonMap['ACC'] = 0; |
| 310 | + codonMap['ACA'] = 0; |
| 311 | + codonMap['ACG'] = 0; |
| 312 | + |
| 313 | + } else if (aa == 'V') { |
| 314 | + |
| 315 | + codonMap['GTT'] = 0; |
| 316 | + codonMap['GTC'] = 0; |
| 317 | + codonMap['GTA'] = 0; |
| 318 | + codonMap['GTG'] = 0; |
| 319 | + |
| 320 | + } else if (aa == 'W') { |
| 321 | + |
| 322 | + codonMap['TGG'] = 0; |
| 323 | + |
| 324 | + } else if (aa == 'Y') { |
| 325 | + |
| 326 | + codonMap['TAT'] = 0; |
| 327 | + codonMap['TAC'] = 0; |
| 328 | + |
| 329 | + } else if (aa == '*') { |
| 330 | + |
| 331 | + codonMap['TAA'] = 0; |
| 332 | + codonMap['TAG'] = 0; |
| 333 | + codonMap['TGA'] = 0; |
| 334 | + |
| 335 | + } else if (aa == 'X') { |
| 336 | + |
| 337 | + codonMap['NNN'] = 0; |
| 338 | + |
| 339 | + } else { |
| 340 | + |
| 341 | + glue.log("INFO", "AMINO RESULT WAS ", aa); |
| 342 | + die; |
| 343 | + |
| 344 | + |
| 345 | + } |
| 346 | + |
| 347 | + return codonMap; |
| 348 | + |
| 349 | +} |
| 350 | + |
| 351 | + |
| 352 | + |
0 commit comments