<jupyter_start><jupyter_text>RNAi recoding   Reading in the Smed transcriptomeWe read in the Smed_v6 transcriptome orfs that were extracted using orfipy. We then join them all into one string and obtain the codon frequencies.<jupyter_code>fname = 'data/dd_Smed_v6_transcripts_orfs.fa'
descriptors, seqs = read_many_fasta(fname)
# join all ORFS into one large transcriptome 
transcriptome = ''.join(seqs)
# get aminoacidweights and codon weights 

codon_frequencies_dic_sw = get_codon_frequencies(transcriptome)<jupyter_output><empty_output><jupyter_text>I also found a published version of amino acid frequencies:<jupyter_code>df = pd.read_csv('data/codon_usage_smed.csv')


AAs = df['codon'].values
freqs = df['frequency'].values/1000.

codon_frequencies_dic = {}
for i in range(len(AAs)):
    codon_frequencies_dic[AAs[i]] = freqs[i]
print(sum(freqs))<jupyter_output>1.00000000001<jupyter_text>We use our codon frequencies dictionary to compute CAI weights (based on the weight definition for the CAI) for all codons $$w_i = \frac{f_i}{\max (f_j)} i,j \in [ \text{synonymouse codons for amino acid} ]$$Where $f_i$ is the frequency of codon $i$. We obtain two dictionaries: aminoacidweights: keys are amino acids, values are arrays of $w_i$ for all synonymous codons. The order of the codons is the as those used in aminoacidcode.     gencodeweights: keys are codons, values are $w_i$ for each codon<jupyter_code>aminoacidweights, gencodeweights = get_codon_weights(codon_frequencies_dic)<jupyter_output><empty_output><jupyter_text>We pickle dump everything so we do not have to repeat the above line later.<jupyter_code>pickle.dump( aminoacidweights,
            open( "data/Smed_transcriptome_aminoacidweights.p", "wb" ) )
pickle.dump( gencodeweights, 
            open( "data/Smed_transcriptome_gencodeweights.p", "wb" ) )
pickle.dump( aminoacidcode,
            open( "data/aminoacidcode.p", "wb" ))<jupyter_output><empty_output><jupyter_text>We reload everything with pickle because why not.<jupyter_code>aminoacidweights = pickle.load( open( "data/Smed_transcriptome_aminoacidweights.p",
                                     "rb" ) )
gencodeweights = pickle.load( open( "data/Smed_transcriptome_gencodeweights.p", 
                                   "rb" ) )
aminoacidcode = pickle.load(open("data/aminoacidcode.p", 'rb'))<jupyter_output><empty_output><jupyter_text>Recoding the H2B ORF for RNAi Here we define the ORF for the histone 2B gene.<jupyter_code>H2B_ORF = 'atggcaattaaaggtaagatcgctgctaagtctgttaagaagatttcaaaggaagttgctcctaaaacagacaaaaagaaaaggatacataaacgcaaagaaagttatggtatttacatctacaaagtgttgagacaagttcatccagatactggaatttctggcaaagcaatgtctattatgaatagctttgtcaacgatgttttcgaaagaatcgcttcggaggctagcaagttggcaacttataataagaaatcaaccataaccagcagagagattcaaactgcagtcaggttaattttaccaggagaattggctaaacacgcagttagtgaaggaaccaaagctgttacaaaatacacaggatccaaa'<jupyter_output><empty_output><jupyter_text>We add on the stop codon.<jupyter_code>H2B_ORF = H2B_ORF.upper() + 'TAA'<jupyter_output><empty_output><jupyter_text>Here we define the protein for the histone 2B gene.<jupyter_code>H2B_protein = 'MAIKGKIAAKSVKKISKEVAPKTDKKKRIHKRKESYGIYIYKVLRQVHPDTGISGKAMSIMNSFVNDVFERIASEASKLATYNKKSTITSREIQTAVRLILPGELAKHAVSEGTKAVTKYTGSK*'<jupyter_output><empty_output><jupyter_text>Now we can use the function get_RNAi_seq to randomly sample different recoded H2B proteins. The function get_RNAi_seq requires the ORF, protein sequence, an aminoacidweights and gencodeweights dictionary. We run 1000 random samples and do not enforce that every codon be different. It returns the list of tested sequences (seqs), scores ($CAI + D$/2) for each sequence, codon adaptation indices (CAIs), and Hamming distances (dists = $D$).<jupyter_code>seqs, scores, cais, dists = get_RNAi_seq(H2B_ORF, H2B_protein, aminoacidweights,
                            gencodeweights, trials = 5000,  enforce_different_codons = False)<jupyter_output><empty_output><jupyter_text>We redo the process but enforce that every codon must be different.<jupyter_code>seqs_diff, scores_diff, cais_diff, dists_diff = get_RNAi_seq(H2B_ORF, H2B_protein, aminoacidweights,
                            gencodeweights, trials = 5000,  enforce_different_codons = True)<jupyter_output><empty_output><jupyter_text>We define a function to compute ECDFs<jupyter_code>def ecdf_vals(data):
    """Return x and y values for an ECDF."""
    return np.sort(data), np.arange(1, len(data)+1) / len(data)<jupyter_output><empty_output><jupyter_text>We plot ECDFs of the CAIs.<jupyter_code># Make figure
p = bokeh.plotting.figure(
    frame_width=400,
    frame_height=300,
    x_axis_label='CAI',
    y_axis_label='ECDF',
)
cais, ecdf_cais = ecdf_vals(cais)
p.circle(cais, ecdf_cais, legend_label = 'Not all different ')

cais_diff, ecdf_cais_diff = ecdf_vals(cais_diff)
p.circle(cais_diff, ecdf_cais_diff, legend_label = 'all different', color = 'orange')
p.legend.location = 'bottom_right'
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>We plot ECDFs of the hamming distances<jupyter_code># Make figure
p = bokeh.plotting.figure(
    frame_width=400,
    frame_height=300,
    x_axis_label='Hamming Distance',
    y_axis_label='ECDF',
    
)
dists, ecdf_dists = ecdf_vals(dists)
p.circle(dists, ecdf_dists, legend_label = 'Not all different ')

dists_diff, ecdf_dists_diff = ecdf_vals(dists_diff)
p.circle(dists_diff, ecdf_dists_diff, legend_label = 'all different', color = 'orange')
p.legend.location = 'bottom_right'
p.x_range = bokeh.models.Range1d(.1, .6)
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>Now we recode the luc ORFS!!!!<jupyter_code>SmedNluc2 = 'ATGGTGTTTACTTTGGAAGATTTTGTTGGAGATTGGAGACAAACTGCTGGTTACAATCTGGATCAGGTACTGGAACAAGGCGGTGTTAGTTCATTATTCCAAAACCTGGGTGTGAGTGTAACTCCGATTCAGCGAATAGTGTTGTCTGGAGAAAATGGGCTGAAGATTGATATACACGTCATAATTCCATACGAAGGCTTAAGCGGTGATCAAATGGGACAAATTGAAAAAATTTTTAAAGTAGTTTACCCAGTTGACGACCATCATTTTAAAGTTATCCTTCATTACGGTACACTGGTTATAGATGGTGTAACTCCAAATATGATCGATTATTTCGGAAGACCTTACGAAGGCATAGCCGTTTTTGATGGAAAAAAGATTACAGTAACAGGTACATTGTGGAACGGAAATAAGATTATTGACGAACGTTTAATTAACCCAGATGGAAGTTTGCTCTTTAGAGTTACAATTAATGGTGTGACAGGATGGAGATTATGCGAACGGATACTCGCGTAA'
smednluc2_protein = 'MVFTLEDFVGDWRQTAGYNLDQVLEQGGVSSLFQNLGVSVTPIQRIVLSGENGLKIDIHVIIPYEGLSGDQMGQIEKIFKVVYPVDDHHFKVILHYGTLVIDGVTPNMIDYFGRPYEGIAVFDGKKITVTGTLWNGNKIIDERLINPDGSLLFRVTINGVTGWRLCERILA*'
hluc = 'ATGGTCTTCACACTCGAAGATTTCGTTGGGGACTGGCGACAGACAGCCGGCTACAACCTGGACCAAGTCCTTGAACAGGGAGGTGTGTCCAGTTTGTTTCAGAATCTCGGGGTGTCCGTAACTCCGATCCAAAGGATTGTCCTGAGCGGTGAAAATGGGCTGAAGATCGACATCCATGTCATCATCCCGTATGAAGGTCTGAGCGGCGACCAAATGGGCCAGATCGAAAAAATTTTTAAGGTGGTGTACCCTGTGGATGATCATCACTTTAAGGTGATCCTGCACTATGGCACACTGGTAATCGACGGGGTTACGCCGAACATGATCGACTATTTCGGACGGCCGTATGAAGGCATCGCCGTGTTCGACGGCAAAAAGATCACTGTAACAGGGACCCTGTGGAACGGCAACAAAATTATCGACGAGCGCCTGATCAACCCCGACGGCTCCCTGCTGTTCCGAGTAACCATCAACGGAGTGACCGGCTGGCGGCTGTGCGAACGCATTCTGGCGTAA'<jupyter_output><empty_output><jupyter_text>I wonder what the CAI for each ORF is?<jupyter_code>print('CAI for SMed Nuc:', get_CAI(SmedNluc2, gencodeweights))
print('CAI for Human Nuc:', get_CAI(hluc, gencodeweights))
print('Hamming Distance vs Smed vs Human Nuc', get_hamming_dist(SmedNluc2, hluc))<jupyter_output>CAI for SMed Nuc: 0.708949323659455
CAI for Human Nuc: 0.49758650479795064
Hamming Distance vs Smed vs Human Nuc 0.25775193798449614<jupyter_text>Seems ok, now lets recoede the Smed ORF<jupyter_code>seqs, scores, cais, dists = get_RNAi_seq(SmedNluc2, smednluc2_protein, aminoacidweights,
                            gencodeweights, trials = 5000,  enforce_different_codons = False)<jupyter_output><empty_output><jupyter_text>We redo the process but enforce that every codon must be different.<jupyter_code>seqs_diff, scores_diff, cais_diff, dists_diff = get_RNAi_seq(SmedNluc2, smednluc2_protein, aminoacidweights,
                            gencodeweights, trials = 5000,  enforce_different_codons = True)

print('Max CAI for any sequence: ', max(cais))
print('Max CAI for all codons different sequence: ', max(cais_diff))
print('Max Hamming Distance for any sequence: ', max(dists))
print('Max Hamming Distance for all codons different sequence: ', max(dists_diff))<jupyter_output>Max Hamming Distance for any sequence:  0.2596899224806202
Max Hamming Distance for all codons different sequence:  0.374031007751938<jupyter_text>We plot ECDFs of the CAIs.<jupyter_code># Make figure
p = bokeh.plotting.figure(
    frame_width=400,
    frame_height=300,
    x_axis_label='CAI',
    y_axis_label='ECDF',
)
cais_sort, ecdf_cais = ecdf_vals(cais)
p.circle(cais_sort, ecdf_cais, legend_label = 'Not all different ')

cais_diff_sort, ecdf_cais_diff = ecdf_vals(cais_diff)
p.circle(cais_diff_sort, ecdf_cais_diff, legend_label = 'all different', color = 'orange')
p.legend.location = 'bottom_right'
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>We plot ECDFs of the hamming distances<jupyter_code># Make figure
p = bokeh.plotting.figure(
    frame_width=400,
    frame_height=300,
    x_axis_label='Hamming Distance',
    y_axis_label='ECDF',
    
)
dists_sort, ecdf_dists = ecdf_vals(dists)
p.circle(dists_sort, ecdf_dists, legend_label = 'Not all different ')

dists_diff_sort, ecdf_dists_diff = ecdf_vals(dists_diff)
p.circle(dists_diff_sort, ecdf_dists_diff, legend_label = 'all different', color = 'orange')
p.legend.location = 'bottom_right'
p.x_range = bokeh.models.Range1d(.1, .6)
bokeh.io.show(p)<jupyter_output><empty_output><jupyter_text>It looks like if they are all different the hamming distance is pretty similar between all sequences, which suggests we should mainly choose by CAI for thehse cases.<jupyter_code>sorted_inds = np.argsort(cais_diff)
seqs_sort = np.array(seqs)[sorted_inds]<jupyter_output><empty_output><jupyter_text>For the sequences in which not every codon is enforced to be the same distance matters.<jupyter_code># Make figure
p = bokeh.plotting.figure(
    frame_width=400,
    frame_height=300,
    x_axis_label='CAIs',
    y_axis_label='Dist',
    
)
p.circle(cais_diff[::10], dists_diff[::10])

p.circle(cais[::10], dists[::10], color = 'orange')
bokeh.io.show(p)
ge<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08465189873417721, 'fraction_numerical': 0.018789556962025316, 'mean_word_length': 7.007125890736342, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '407425', 'license': 'MIT', 'n_tokens_mistral': 4206, 'n_tokens_neox': 3854, 'n_words': 701, 'path': 'cloning-scripts-testing.ipynb', 'repo_name': 'sophiejwalton/random-useful-scripts', 'size': '10112'}