<jupyter_start><jupyter_text>Loading the big dataframe:<jupyter_code>data = pd.read_csv('../../Data/New_data_table_Xist.txt',\
                    sep="\t",\
                    header=0,\
                    index_col=0)<jupyter_output><empty_output><jupyter_text>Removing data with no coverage:<jupyter_code>positionstomask = data[(data["SHAPE_reactivity_ex_vivo_1"] < -900) | 
                       (data["SHAPE_reactivity_ex_vivo_2"] < -900) |
                       (data["SHAPE_reactivity_in_cell_1"] < -900) | 
                       (data["SHAPE_reactivity_in_cell_2"] < -900)].index
print(len(positionstomask))
data.loc[positionstomask,data.columns[11:]] = np.nan<jupyter_output>2517<jupyter_text>Checking what nucleotides were found to be modified in the ex_vivo and in_vivo data:<jupyter_code>samples  = ['Xist_sequence','diffBUM_HMM_ex_vivo_norm','diffBUM_HMM_ex_vivo','diffBUM_HMM_in_vivo_norm','diffBUM_HMM_in_vivo',
            'diffBUM_HMM_ex_vivo_norm_unique','diffBUM_HMM_ex_vivo_unique',
            'diffBUM_HMM_in_vivo_norm_unique','diffBUM_HMM_in_vivo_unique',
            'deltaSHAPE_ex_vivo_1','deltaSHAPE_ex_vivo_2','deltaSHAPE_in_vivo_1','deltaSHAPE_in_vivo_2']
nucleotides = ["A","T","C","G"]
modinfo = pd.DataFrame(0,index=nucleotides,columns=samples)
modinfo.columns<jupyter_output><empty_output><jupyter_text>Xist sequence:<jupyter_code>for i in data.nucleotide:
    modinfo.loc[i,'Xist_sequence'] += 1<jupyter_output><empty_output><jupyter_text>Setting a threshold for calling DRNs in the diffBUM-HMM data:<jupyter_code>threshold = 0.95<jupyter_output><empty_output><jupyter_text>First I need to figure out how long the stretches of DRNs are in the diffBUM-HMM and deltaSHAPE data:This will produce the data for Figure 6A.<jupyter_code>intervallengthcounter = defaultdict(lambda: defaultdict(int))<jupyter_output><empty_output><jupyter_text>diffBUM-HMM:<jupyter_code>### How many DRNs are there in the data?:
ex_vivo_pos = data[data.ex_vivo >= threshold].index

### contigousarray2Intervals groups the DRNs that are next to each other in intervals.
intervals = contigousarray2Intervals(ex_vivo_pos)

### This prints the number of intervals that were detected and how many DRNs there were in the data.
print(len(intervals))
print(len(ex_vivo_pos))

for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["diffBUM_HMM_ex_vivo"][length] += 1
        
in_vivo_pos = data[data.in_vivo >= threshold].index
intervals = contigousarray2Intervals(in_vivo_pos)

print(len(intervals))
print(len(in_vivo_pos))

for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["diffBUM_HMM_in_vivo"][length] += 1<jupyter_output>849
1164
127
142<jupyter_text>deltaSHAPE rep 1:<jupyter_code>ex_vivo_pos = data[data.deltaSHAPE_rep1 > 0].index
intervals = contigousarray2Intervals(ex_vivo_pos)

print(len(intervals))
print(len(ex_vivo_pos))


for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["deltaSHAPE_rep1_ex_vivo"][length] += 1
    
in_vivo_pos = data[data.deltaSHAPE_rep1 < 0].index
intervals = contigousarray2Intervals(in_vivo_pos)

print(len(intervals))
print(len(in_vivo_pos))

count = 0

for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["deltaSHAPE_rep1_in_vivo"][length] += 1<jupyter_output>162
460
120
337<jupyter_text>deltaSHAPE rep 2:<jupyter_code>ex_vivo_pos = data[data.deltaSHAPE_rep2 > 0].index
intervals = contigousarray2Intervals(ex_vivo_pos)

print(len(intervals))
print(len(ex_vivo_pos))


for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["deltaSHAPE_rep2_ex_vivo"][length] += 1
        
in_vivo_pos = data[data.deltaSHAPE_rep2 < 0].index
intervals = contigousarray2Intervals(in_vivo_pos)

print(len(intervals))
print(len(in_vivo_pos))

count = 0

for (i,j) in intervals:
    length = (j-i)+1
    intervallengthcounter["deltaSHAPE_rep2_in_vivo"][length] += 1<jupyter_output>186
510
166
428<jupyter_text>Now making the plots for Figure 6A:<jupyter_code>numberofplots = len(intervallengthcounter.keys())
samples = list(intervallengthcounter.keys())
fig,ax = plt.subplots(numberofplots,figsize=[3,10],sharex=True)
for i in range(numberofplots):
    sample = samples[i]
    counts = intervallengthcounter[sample]
    x = list(counts.keys())
    y = list(counts.values())
    sns.barplot(x,y,ax=ax[i],color="blue")
    ax[i].set_title(sample)
    ax[i].set_ylabel('Counts')
plt.tight_layout()
fig.savefig("Figure_6A.pdf",dpi=400)<jupyter_output><empty_output><jupyter_text>Now calculating the nucleotide composition of the DRNs and how frequently they are located in single-stranded regions:  diffBUM_HMM:<jupyter_code>in_vivo = data[data.in_vivo >= threshold]["nucleotide"].values
ex_vivo = data[data.ex_vivo >= threshold]["nucleotide"].values
in_vivo_scaled = data[data.scaled_in_vivo >= threshold]["nucleotide"].values
ex_vivo_scaled = data[data.scaled_ex_vivo >= threshold]["nucleotide"].values
for nuc in in_vivo:
    modinfo.loc[nuc,"diffBUM_HMM_in_vivo"] += 1
for nuc in in_vivo_scaled:
    modinfo.loc[nuc,"diffBUM_HMM_in_vivo_norm"] += 1
for nuc in ex_vivo:
    modinfo.loc[nuc,"diffBUM_HMM_ex_vivo"] += 1
for nuc in ex_vivo_scaled:
    modinfo.loc[nuc,"diffBUM_HMM_ex_vivo_norm"] += 1<jupyter_output><empty_output><jupyter_text>diffBUM_HMM unique<jupyter_code>in_vivo_unique = data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.in_vivo >= threshold)]["nucleotide"].values
ex_vivo_unique = data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.ex_vivo >= threshold)]["nucleotide"].values
in_vivo_scaled_unique = data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                             (data.scaled_in_vivo >= threshold)]["nucleotide"].values
ex_vivo_scaled_unique = data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                             (data.scaled_ex_vivo >= threshold)]["nucleotide"].values
for nuc in in_vivo_unique:
    modinfo.loc[nuc,"diffBUM_HMM_in_vivo_unique"] += 1
for nuc in in_vivo_scaled_unique:
    modinfo.loc[nuc,"diffBUM_HMM_in_vivo_norm_unique"] += 1
for nuc in ex_vivo_unique:
    modinfo.loc[nuc,"diffBUM_HMM_ex_vivo_unique"] += 1
for nuc in ex_vivo_scaled_unique:
    modinfo.loc[nuc,"diffBUM_HMM_ex_vivo_norm_unique"] += 1<jupyter_output><empty_output><jupyter_text>Now deltaSHAPE data:<jupyter_code>in_vivo_1 = data[data.deltaSHAPE_rep1 > 0]["nucleotide"].values
ex_vivo_1 = data[data.deltaSHAPE_rep1 < 0]["nucleotide"].values
in_vivo_2 = data[data.deltaSHAPE_rep2 > 0]["nucleotide"].values
ex_vivo_2 = data[data.deltaSHAPE_rep2 < 0]["nucleotide"].values
for nuc in in_vivo_1:
    modinfo.loc[nuc,"deltaSHAPE_in_vivo_1"] += 1
for nuc in in_vivo_2:
    modinfo.loc[nuc,"deltaSHAPE_in_vivo_2"] += 1
for nuc in ex_vivo_1:
    modinfo.loc[nuc,"deltaSHAPE_ex_vivo_1"] += 1
for nuc in ex_vivo_2:
    modinfo.loc[nuc,"deltaSHAPE_ex_vivo_2"] += 1
modinfo<jupyter_output><empty_output><jupyter_text>Now correct for total number of modifications:<jupyter_code>modinfonorm = modinfo.div(modinfo.sum(axis=0), axis=1)<jupyter_output><empty_output><jupyter_text>Now check if there is a statistical significant difference between the nucleotide frequencies in each dataset using a chi-square test:<jupyter_code>modinfonorm.head()
for i in modinfonorm.columns[1:]:
    f_exp = modinfonorm['Xist_sequence'].values
    f_obs = modinfonorm[i].values
    #print(f_exp,f_obs)
    p_value = chisquare(f_obs,f_exp)[1]
    print("%s:\t%s" % (i,p_value))<jupyter_output>diffBUM_HMM_ex_vivo_norm:	0.9921921939576517
diffBUM_HMM_ex_vivo:	0.9861833528301358
diffBUM_HMM_in_vivo_norm:	0.8774512930475329
diffBUM_HMM_in_vivo:	0.8584003893945844
diffBUM_HMM_ex_vivo_norm_unique:	0.9916034877939462
diffBUM_HMM_ex_vivo_unique:	0.9857397075770127
diffBUM_HMM_in_vivo_norm_unique:	0.8734742568696059
diffBUM_HMM_in_vivo_unique:	0.8557288380621534
deltaSHAPE_ex_vivo_1:	0.99726258943411
deltaSHAPE_ex_vivo_2:	0.9854985448954785
deltaSHAPE_in_vivo_1:	0.9910509186704817
deltaSHAPE_in_vivo_2:	0.9942966910237633<jupyter_text>So the difference is not statistically significant but the sample size for this test is too small. Should be at least 5. There must be a better approach.  Now making the heatmap:<jupyter_code>modinfonorm.head()<jupyter_output><empty_output><jupyter_text>Now adding data for how many of the DRNs are located in regions predicted to be single-stranded in Xist:<jupyter_code>totalss = len(data[data.base_paired_to == 0].index)
deltaSHAPE_ex_vivo_rep1_total = len(data[data.deltaSHAPE_rep1 > 0].index)
deltaSHAPE_ex_vivo_rep2_total = len(data[data.deltaSHAPE_rep2 > 0].index)
deltaSHAPE_in_vivo_rep1_total = len(data[data.deltaSHAPE_rep1 < 0].index)
deltaSHAPE_in_vivo_rep2_total = len(data[data.deltaSHAPE_rep2 < 0].index)
in_vivo_total = len(data[(data.in_vivo >= threshold)].index)
ex_vivo_total = len(data[(data.ex_vivo >= threshold)].index)
in_vivo_unique_total = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.in_vivo >= threshold)].index)
ex_vivo_unique_total = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.ex_vivo >= threshold)].index)
in_vivo_scaled_total = len(data[(data.scaled_in_vivo >= threshold)].index)
ex_vivo_scaled_total = len(data[(data.scaled_ex_vivo >= threshold)].index)
in_vivo_scaled_unique_total = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                             (data.scaled_in_vivo >= threshold)].index)
ex_vivo_scaled_unique_total = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                             (data.scaled_ex_vivo >= threshold)].index)
deltaSHAPE_ex_vivo_rep1 = len(data[(data.deltaSHAPE_rep1 > 0) & (data.base_paired_to == 0)].index)
deltaSHAPE_ex_vivo_rep2 = len(data[(data.deltaSHAPE_rep2 > 0) & (data.base_paired_to == 0)].index)
deltaSHAPE_in_vivo_rep1 = len(data[(data.deltaSHAPE_rep1 < 0) & (data.base_paired_to == 0)].index)
deltaSHAPE_in_vivo_rep2 = len(data[(data.deltaSHAPE_rep2 < 0) & (data.base_paired_to == 0)].index)
in_vivo = len(data[(data.in_vivo >= threshold) & (data.base_paired_to == 0)].index)
ex_vivo = len(data[(data.ex_vivo >= threshold) & (data.base_paired_to == 0)].index)
in_vivo_unique = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.in_vivo >= threshold) & (data.base_paired_to == 0)].index)
ex_vivo_unique = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                      (data.ex_vivo >= threshold) & (data.base_paired_to == 0)].index)
in_vivo_scaled = len(data[(data.scaled_in_vivo >= threshold) & (data.base_paired_to == 0)].index)
ex_vivo_scaled = len(data[(data.scaled_ex_vivo >= threshold) & (data.base_paired_to == 0)].index)
in_vivo_scaled_unique = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                                 (data.scaled_in_vivo >= threshold) & (data.base_paired_to == 0)].index)
ex_vivo_scaled_unique = len(data[(data.deltaSHAPE_rep1 == 0) & (data.deltaSHAPE_rep2 == 0) & \
                                 (data.scaled_ex_vivo >= threshold) & (data.base_paired_to == 0)].index)
samples = modinfonorm.columns
samples
frequencies = [ totalss/len(data.index),
                ex_vivo_scaled/ex_vivo_scaled_total,
                ex_vivo/ex_vivo_total,
                in_vivo_scaled/in_vivo_scaled_total,
                in_vivo/in_vivo_total,
                ex_vivo_scaled_unique/ex_vivo_scaled_unique_total,
                ex_vivo_unique/ex_vivo_unique_total,
                in_vivo_scaled_unique/in_vivo_scaled_unique_total,
                in_vivo_unique/in_vivo_unique_total,
                deltaSHAPE_ex_vivo_rep1/deltaSHAPE_ex_vivo_rep1_total,
                deltaSHAPE_ex_vivo_rep2/deltaSHAPE_ex_vivo_rep2_total,
                deltaSHAPE_in_vivo_rep1/deltaSHAPE_in_vivo_rep1_total,
                deltaSHAPE_in_vivo_rep2/deltaSHAPE_in_vivo_rep2_total]
frequencies
len(samples)
df = pd.DataFrame(frequencies,columns=["frequencies"],index=samples)
df
fig,ax = plt.subplots(figsize=[9,3.5])
ax = sns.barplot(df.index,df.frequencies,color="blue")
ax.set_yticks(np.arange(0,1.2,0.2))
ax.set_yticklabels(np.arange(0,1.2,0.2))
ax.set_xticklabels(df.index,rotation=45,ha='right')
fig.savefig("Nucleotide_ssOverlap_counts.pdf",dpi=400)
modinfonorm = modinfonorm.append(df.T)
fig,ax = plt.subplots(figsize=[9,3.5])
ax = sns.heatmap(modinfonorm,\
            cmap='Blues',\
            linewidths=.5,\
            annot=True)
ax.set_yticklabels(modinfonorm.index,rotation=0,fontsize=14)

fig.savefig("Figure_6B.pdf",dpi=400)<jupyter_output><empty_output><jupyter_text>There seem to be a lot of nucleotides in the diffBUM-HMM in vivo data that are A-rich but according to the secondary structure data do not seem to be part of single-stranded regions? Could this have biological significance?<jupyter_code>A_and_double_stranded = data[(data.nucleotide == "A") & \
                             (data.in_vivo >= threshold) & \
                             (data.base_paired_to > 0)]
print(len(A_and_double_stranded.index))<jupyter_output>36<jupyter_text>How many DRNs were in the in vivo data?<jupyter_code>totaldrnsinvivo = data[data.in_vivo >= threshold]
print(len(totaldrnsinvivo.index))
totaldrnsinvivo.to_csv("diffBUM-HMM_DRNs_in_vivo_locations.txt",sep="\t")<jupyter_output><empty_output><jupyter_text>That is quite a reasonable fraction.Where are they located?<jupyter_code>A_and_double_stranded.index<jupyter_output><empty_output><jupyter_text>They don't seem to cluster very well...Perhaps add something about this to the discussion.<jupyter_code>A_and_double_stranded.to_csv("diffBUM-HMM_DRNs_in_As_in_vivo_locations.txt",sep="\t")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12181528662420382, 'fraction_numerical': 0.03488708743485813, 'mean_word_length': 5.105612019443217, 'pattern_counts': {'":': 0, '<': 91, '<?xml version=': 0, '>': 116, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, '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': '442472', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5731, 'n_tokens_neox': 5341, 'n_words': 761, 'path': 'Jupyter_notebooks/Figure_6/Nucleotide_analyses.ipynb', 'repo_name': 'marangiop/my_binder', 'size': '13816'}