Mercurial > repos > rlegendre > ribo_tools
comparison get_codon_frequency.py @ 16:fcfdb2607cb8 draft
Uploaded
| author | rlegendre |
|---|---|
| date | Wed, 13 May 2015 04:34:59 -0400 |
| parents | 702e60e819c2 |
| children | c87c40e642af |
comparison
equal
deleted
inserted
replaced
| 15:702e60e819c2 | 16:fcfdb2607cb8 |
|---|---|
| 323 # # standard deviation = absolute value of diffence between replicats of cond1 | 323 # # standard deviation = absolute value of diffence between replicats of cond1 |
| 324 std_cond1.append(std(array([(cond1_1[i] * 100 / sum11), (cond1_2[i] * 100 / sum12)]))) | 324 std_cond1.append(std(array([(cond1_1[i] * 100 / sum11), (cond1_2[i] * 100 / sum12)]))) |
| 325 # # cond2 = moyenne of replicats cond1divided by max | 325 # # cond2 = moyenne of replicats cond1divided by max |
| 326 cond2_val[i] = ((cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2) | 326 cond2_val[i] = ((cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2) |
| 327 cond2[i] = ((cond2_1[i] + cond2_2[i]) / 2) | 327 cond2[i] = ((cond2_1[i] + cond2_2[i]) / 2) |
| 328 # # standard deviation = absolute value of diffence between replicats of cond2 | 328 # # standard deviation = absolute value of difference between replicats of cond2 |
| 329 std_cond2.append(std(array([((cond2_1[i]) * 100 / sum21), ((cond2_2[i]) * 100 / sum22)]))) | 329 std_cond2.append(std(array([((cond2_1[i]) * 100 / sum21), ((cond2_2[i]) * 100 / sum22)]))) |
| 330 # # max value for each codon | 330 # # max value for each codon |
| 331 max_val.append(max((cond1_1[i] / sum11 + cond1_2[i] / sum12) / 2, (cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2)) | 331 max_val.append(max((cond1_1[i] / sum11 + cond1_2[i] / sum12) / 2, (cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2)) |
| 332 | 332 |
| 333 # for graph design | 333 # for graph design |
| 377 expected = array(cond2_count) | 377 expected = array(cond2_count) |
| 378 observed = array(list(cond2.itervalues())) | 378 observed = array(list(cond2.itervalues())) |
| 379 | 379 |
| 380 # write result | 380 # write result |
| 381 with open(outfile, 'w') as out : | 381 with open(outfile, 'w') as out : |
| 382 out.write('Codon\tRaw_' + c1 + '\tRaw_' + c2 + '\tNorm_' + c1 + '\tNorm_' + c2 + '\tFC\tFC_' + c1 + '\tFC_' + c2 + '\n') | 382 out.write('Codon\tRaw_' + c1 + '\tRaw_' + c2 + '\tNorm_' + c1 + '\tNorm_' + c2 + '\tFC(Mut/WT)\n') |
| 383 for i in codon_sorted: | 383 for i in codon_sorted: |
| 384 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\t' + str((cond2_1[i] / sum21) / (cond1_1[i] / sum11)) + '\t' + str((cond2_2[i] / sum22) / (cond1_1[i] / sum11)) + '\n') | 384 ## if global foldchange is equal to zero |
| 385 if cond1_norm[i] == 0 and cond2_norm[i] == 0: | |
| 386 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t1.0\n') | |
| 387 elif cond1_norm[i] == 0 : | |
| 388 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t0.0\n') | |
| 389 else: | |
| 390 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\n') | |
| 385 chi = stats.chisquare(observed, expected) | 391 chi = stats.chisquare(observed, expected) |
| 386 out.write('Khi2 test\n') | 392 out.write('Khi2 test\n') |
| 387 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n') | 393 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n') |
| 388 | 394 |
| 389 | 395 |
| 482 | 488 |
| 483 # write result | 489 # write result |
| 484 with open(outfile, 'w') as out : | 490 with open(outfile, 'w') as out : |
| 485 out.write('Codon\tRaw_' + c1 + '\tRaw_' + c2 + '\tNorm_' + c1 + '\tNorm_' + c2 + '\tFC(Mut/WT)\n') | 491 out.write('Codon\tRaw_' + c1 + '\tRaw_' + c2 + '\tNorm_' + c1 + '\tNorm_' + c2 + '\tFC(Mut/WT)\n') |
| 486 for i in codon_sorted: | 492 for i in codon_sorted: |
| 487 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\n') | 493 if cond1_norm[i] == 0 and cond2_norm[i] == 0: |
| 494 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t1.0\n') | |
| 495 elif cond1_norm[i] == 0 : | |
| 496 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t0.0\n') | |
| 497 else: | |
| 498 out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\n') | |
| 488 out.write('Khi2 test\n') | 499 out.write('Khi2 test\n') |
| 489 chi = stats.chisquare(observed, expected) | 500 chi = stats.chisquare(observed, expected) |
| 490 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n') | 501 out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n') |
| 491 | 502 |
| 492 # # get max value for each codon for histogram | 503 # # get max value for each codon for histogram |
| 604 pl.axhline(y=1,color='r') | 615 pl.axhline(y=1,color='r') |
| 605 pl.xticks(index, label, rotation=90) | 616 pl.xticks(index, label, rotation=90) |
| 606 pl.ylabel('Foldchange of codon occupancy') | 617 pl.ylabel('Foldchange of codon occupancy') |
| 607 ax.yaxis.set_ticks_position('left') | 618 ax.yaxis.set_ticks_position('left') |
| 608 ax.xaxis.set_ticks_position('bottom') | 619 ax.xaxis.set_ticks_position('bottom') |
| 620 pl.ylim(-1,3) | |
| 609 pl.title(site+" site") | 621 pl.title(site+" site") |
| 610 pl.savefig(dirout + '/fc_codons.png', format="png", dpi=340) | 622 pl.savefig(dirout + '/fc_codons.png', format="png", dpi=340) |
| 611 | 623 |
| 612 | 624 |
| 613 def __main__(): | 625 def __main__(): |