Copy
import numpy as np
import pandas as pd
import re
from numba import jit
from scipy.stats import ks_2samp
import statsmodels.stats.multitest as multi
####################函数定义部分####################
@jit(nopython=False)
def get_gini(R_, nucleotide_, cut_AC=False, ratio_nan=0.9, ratio_nan_AC=0.8):
R_nonull = R_[np.isnan(R_)==False]
ratio = len(R_nonull)/len(R_)
if ratio <= ratio_nan:
gini = np.nan
else:
if cut_AC:
R_AC = R_[(nucleotide_==b'A')|(nucleotide_==b'C')]
R_AC_nonull = R_AC[np.isnan(R_AC)==False]
ratio_AC = len(R_AC_nonull) / len(R_AC)
if (ratio_AC <= ratio_nan_AC)|len(R_AC)<=1|len(R_AC_nonull)<=1:
gini = np.nan
else:
sorted = np.sort(R_AC_nonull)
height, area = 0, 0
for i in range(0,len(sorted)):
height += sorted[i]
area += height - sorted[i] / 2.
fair_area = height * len(sorted) / 2.
if fair_area ==0:
gini=np.nan
else:
gini = (fair_area - area) / fair_area
else:
sorted = np.sort(R_nonull)
height, area = 0, 0
for i in range(0, len(sorted)):
height += sorted[i]
area += height - sorted[i] / 2.
fair_area = height * len(sorted) / 2.
if fair_area == 0:
gini = np.nan
else:
gini = (fair_area - area) / fair_area
return gini
@jit(nopython=False)
def get_p(R_1,R_2,ratio_nan=0.9):
R_1_nonull = R_1[np.isnan(R_1) == False]
ratio_1 = len(R_1_nonull) / len(R_1)
R_2_nonull = R_2[np.isnan(R_2) == False]
ratio_2 = len(R_2_nonull) / len(R_2)
if (ratio_1 <= ratio_nan)|(ratio_2 <= ratio_nan):
p = np.nan
else:
s,p=ks_2samp(R_1, R_2)
return p
@jit(nopython=False)
def get_window(X,location,window=50,step=1,cut_AC=False,ratio_nan=0.9,ratio_nan_AC=0.8):
#########读取数据#############
name = X[0]
nucleotide = np.array(list(X[1]))
modified = np.array(X[2].split(',')).astype('int')
modified_depth = np.array(X[3].split(',')).astype('int')
unmodified = np.array(X[4].split(',')).astype('int')
unmodified_depth = np.array(X[5].split(',')).astype('int')
Mutrs = modified/modified_depth
Mutru = unmodified/unmodified_depth
# R = Mutrs - Mutru
R = np.array(X[6].split(',')).astype('float')
R[R == -999] = np.nan
n = len(nucleotide)
index = np.array(range(0, n))
location_index = np.array(np.array(location['location'])[0].split(',')).astype('int')
#####滑动窗口计算change##################
nucleotide_ = np.zeros([len(range(0, n-(window), step)), window],dtype=np.string_)
R_ = np.zeros([len(range(0, n-(window), step)), window])
gini = np.zeros([len(range(0, n-(window), step)),])
index_ = np.zeros([len(range(0, n-(window), step)), window])
location_index_ = np.zeros([len(range(0, n - (window), step)), window])
j =0
for i in range(0, n-(window), step):
index_[j,:]=index[i:i+window]
location_index_[j,:] = location_index[i:i+window]
nucleotide_[j,:] = nucleotide[i:i+window]
R_[j,:] = R[i:i+window]
gini[j] = get_gini(R[i:i+window],nucleotide[i:i+window],cut_AC=cut_AC,ratio_nan=ratio_nan,ratio_nan_AC=ratio_nan_AC)
j=j+1
return index_,location_index_,nucleotide_,R_,gini
@jit(nopython=False)
def get_window_p(X1,X2, window=50, step=1,ratio_nan=0.9):
#########读取数据#############
name = X1[0]
nucleotide = np.array(list(X1[1]))
R_z = np.array(X1[6].split(',')).astype('float')
R_z[R_z == -999] = np.nan
R_f = np.array(X2[6].split(',')).astype('float')
R_f[R_f == -999] = np.nan
n = len(nucleotide)
#####滑动窗口计算change##################
nucleotide_ = np.zeros([len(range(0, n - (window), step)), window], dtype=np.string_)
R_z_ = np.zeros([len(range(0, n - (window), step)), window])
R_f_ = np.zeros([len(range(0, n - (window), step)), window])
p = np.zeros([len(range(0, n - (window), step)), ])
j = 0
for i in range(0, n - (window), step):
nucleotide_[j, :] = nucleotide[i:i + window]
R_z_[j, :] = R_z[i:i + window]
R_f_[j, :] = R_f[i:i + window]
p[j] = get_p(R_z[i:i + window], R_f[i:i + window],ratio_nan=ratio_nan,)
j = j + 1
p_=p.copy()
if len(p[~np.isnan(p)])>0:
a,p_bh,b,c =multi.multipletests(p[~np.isnan(p)],method='fdr_bh')
p_[~np.isnan(p_)]=p_bh
return p_,p
@jit(nopython=False)
def calculate_delta_gini(R_1, R_2, gini_1, gini_2, ratio_nan=0.9):
'''Calculates Standard RMSD on two vectors of numbers of the same length'''
# Check to see the vectors are of equal length.
if len(R_1) != len(R_2):
return np.nan
else:
R = R_1 - R_2
if len(R[np.isnan(R) == False]) / len(R) <= ratio_nan:
return np.nan
else:
delta_gini = gini_1 - gini_2
return delta_gini
@jit(nopython=False)
def get_gff(X,chr,strand,data_gff_):
gff=[]
if len(data_gff_.loc[data_gff_['location'] == 'CDS', 'strat']) == 0:
for i in range(len(X)):
gff_='erro'
gff.append(gff_)
else:
CDS_strat = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'strat'])[0]
CDS_end = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'end'])[0]
for i in range(len(X)):
if strand == '+':
if X[i]<CDS_strat:
gff_='five_prime_UTR'
elif X[i]>CDS_end:
gff_ = 'three_prime_UTR'
else:
gff_='CDS'
else:
if X[i]<CDS_strat:
gff_='three_prime_UTR'
elif X[i]>CDS_end:
gff_ = 'five_prime_UTR'
else:
gff_='CDS'
gff.append(gff_)
return gff
@jit(nopython=False)
def merge_gini(X1,X2,strat,end):
#########读取数据#############
name = X1[0]
nucleotide = np.array(list(X1[1]))
R_z = np.array(X1[6].split(',')).astype('float')
R_z[R_z == -999] = np.nan
R_f = np.array(X2[6].split(',')).astype('float')
R_f[R_f == -999] = np.nan
n = len(nucleotide)
R_z_=R_z[strat:end+1]
R_f_=R_f[strat:end+1]
gini_z=get_gini(R_z_,nucleotide_=[],ratio_nan=0)
gini_f = get_gini(R_f_, nucleotide_=[],ratio_nan=0)
delta=gini_z-gini_f
p=get_p(R_z_,R_f_,ratio_nan=0)
R_f_mean = R_f_[R_f_>=0].mean()
R_z_mean = R_z_[R_z_ >= 0].mean()
return delta,p,R_f_mean,R_z_mean,gini_z,gini_f
@jit(nopython=False)
def get_statistics(data_z,data_f,data_location_,data_gff_,strand):
data_z = np.array(data_z).reshape([7,])
data_f = np.array(data_f).reshape([7,])
chr=list(data_location_['chr'])[0]
index_z,location_index_z,nucleotide_z, R_z, gini_z = get_window(data_z,data_location_)
index_f,location_index_f,nucleotide_f, R_f, gini_f = get_window(data_f,data_location_)
# print('ok')
if len(R_z)!=len(R_f):
print("error")
else:
location_exon = get_gff(location_index_z[:,0],chr,strand,data_gff_)
# print('ok')
location_exon = np.array(location_exon)
p_bh, p = get_window_p(data_z, data_f)
delta_gini = np.zeros([len(R_z), ])
index_list = []
nucleotide_list = []
location_index_list = []
for i in range(len(R_z)):
index_list.append( ','.join(list(index_z[i,:].astype('int').astype('str'))))
nucleotide_list.append(','.join(list(nucleotide_z[i, :].astype('str'))))
location_index_list.append( ','.join(list(location_index_z[i,:].astype('int').astype('str'))))
delta_gini[i] = calculate_delta_gini(R_z[i, :], R_f[i, :], gini_z[i], gini_f[i])
index_list = np.array(index_list)
location_index_list = np.array(location_index_list)
nucleotide_list = np.array(nucleotide_list)
return delta_gini,p_bh,p, gini_z, gini_f,index_list,location_index_list,nucleotide_list,location_exon
# @jit(nopython=False,error_model="numpy")
def merge_data(delta_gini,index_list,data_z,data_f,location_list,chr,strand,data_gff_,windows=50):
data_z = np.array(data_z).reshape([7,])
data_f = np.array(data_f).reshape([7,])
nucleotide = np.array(list(data_z[1]))
index = np.array([x.split(',')[0] for x in index_list]).astype('int')
delta_gini_ =delta_gini.copy()
delta_gini_[np.isnan(delta_gini_)]=0
i=0
windows_num=[]
n=index[abs(delta_gini_)>0.1][0]+windows-1
for j in range(len(index)-1):
if abs(delta_gini_[j])<0.1:
windows_num.append(np.nan)
else:
if index[j]>n+1:
i=i+1
n=index[j]+windows-1
windows_num.append(i)
if abs(delta_gini_[len(index)-1])<0.1:
windows_num.append(np.nan)
else:
if index[j] > n + 1:
i = i + 1
windows_num.append(i)
data_all= pd.DataFrame(np.zeros([len(set(windows_num) - set([np.nan])),9]))
for i in list(set(windows_num) - set([np.nan])):
windows_num_=np.array(windows_num)
strat = index[windows_num_==i].astype('int')[0]
end = index[windows_num_==i].astype('int')[-1]+windows-1
strat_chr=location_list[strat]
end_chr = location_list[end]
delta,p=merge_gini(data_z,data_f,strat,end)
location = get_gff(np.array([strat_chr,end_chr]),chr,strand,data_gff_)
nucleotide_ =nucleotide[strat:end+1]
nucleotide_str = ','.join(list(nucleotide_.astype('str')))
data_all.iloc[i,:]=[strat,end,strat_chr,end_chr,location[0],location[1],delta,p,nucleotide_str]
data_all.columns=['start','end','start_chr','end_chr','location_start','location_end','delta','p','nucleotide']
return data_all
def merge_data_2(delta_gini,index_list,data_z,data_f,location_list,chr,strand,data_gff_,windows=50):
data_z = np.array(data_z).reshape([7,])
data_f = np.array(data_f).reshape([7,])
nucleotide = np.array(list(data_z[1]))
index = np.array([x.split(',')[0] for x in index_list]).astype('int')
delta_gini_ =delta_gini.copy()
delta_gini_[np.isnan(delta_gini_)]=0
i=0
windows_num=[]
if abs(delta_gini_[0]) < 0.1:
windows_num.append(np.nan)
else:
i = i + 1
windows_num.append(i)
for j in range(1,len(index)):
if abs(delta_gini_[j])<0.1:
windows_num.append(np.nan)
else:
if abs(delta_gini_[j-1])>0.1:
windows_num.append(i)
else:
i=i+1
windows_num.append(i)
data_all= pd.DataFrame(np.zeros([len(set(windows_num) - set([np.nan])),13]))
for i in list(set(windows_num) - set([np.nan])):
windows_num_=np.array(windows_num)
strat = index[windows_num_==i].astype('int')[0]
end = index[windows_num_==i].astype('int')[-1]+windows-1
delta, p, R_f_mean, R_z_mean,gini_z,gini_f=merge_gini(data_z,data_f,strat,end)
if abs(delta)<=0.1:
num = len(index[windows_num_==i].astype('int'))
for n in range(1,num):
end_=end-n
delta, p, R_f_mean, R_z_mean,gini_z,gini_f = merge_gini(data_z, data_f, strat, end_)
if abs(delta)>0.1:
end = end_
break
strat_chr=location_list[strat]
end_chr = location_list[end]
location = get_gff(np.array([strat_chr,end_chr]),chr,strand,data_gff_)
nucleotide_ =nucleotide[strat:end+1]
nucleotide_str = ','.join(list(nucleotide_.astype('str')))
data_all.iloc[i-1,:]=[strat,end,strat_chr,end_chr,location[0],location[1],delta,p,R_f_mean, R_z_mean,gini_z,gini_f,nucleotide_str]
data_all.columns=['start','end','start_chr','end_chr','location_start','location_end','delta','p','R_f_mean','R_z_mean','gini_z','gini_f','nucleotide']
return data_all
def main_sum(data_uv_z,data_uv_f,transcript_id_list,data_location,data_gff,save_path):
i=0
data_all = pd.DataFrame(columns=['transcript_id','start','end','start_chr','end_chr','location_start','location_end','delta','p','R_f_mean','R_z_mean','gini_z','gini_f','nucleotide'])
for transcript in transcript_id_list:
data_z = data_uv_z.loc[data_uv_z['transcript_id']==transcript,:]
data_f = data_uv_f.loc[data_uv_f['transcript_id'] == transcript,:]
data_location_ = data_location.loc[data_location['transcript_id'] == transcript,:]
data_gff_ = data_gff.loc[data_gff['transcript_id']== transcript,:]
location_list=np.array(list(data_location_['location'])[0].split(',')).astype('int')
chr = list(data_location_['chr'])[0]
strand = list(data_location_['strand'])[0]
delta_gini, p_bh,p,gini_z, gini_f, index_list, location_index_list, nucleotide_list,location_exon = get_statistics(data_z,data_f,data_location_,data_gff_,strand)
# print(delta_gini)
data_new = pd.DataFrame(np.vstack([location_exon,index_list,location_index_list,nucleotide_list,gini_z, gini_f,delta_gini,p_bh,p])).T
data_new.columns =['location_exon','index','location_index','nucleotide','gini_z', 'gini_f','delta_gini','p_bh','p']
data_new['transcript']=transcript
# pd.DataFrame(data_new).to_csv(save_path+'/'+transcript+'_gini.csv',sep='\t',header=True,index=False)
data_all_=merge_data_2(delta_gini, index_list, data_z, data_f, location_list,chr,strand,data_gff_)
data_all_['transcript_id']=transcript
data_all = pd.concat([data_all,data_all_])
i = i+1
print(i)
return data_all
####################主函数部分####################
if __name__ == '__main__':
col_uv_f = '/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_nouv/'
col_uv_z = '/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_uv/'
###col###
result = '/data/TA_QUIZ_RNA_regulation/result/PartIII.SHAPE-seq_analysis/col'
uv_z = col_uv_z
uv_f = col_uv_f
data_uv_z = pd.read_csv(uv_z + '/final.modified_unmodified_new', sep='\t')
data_uv_z.columns = ['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']
data_uv_z = data_uv_z[['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']]
data_uv_z = data_uv_z.drop_duplicates()
data_uv_f = pd.read_csv(uv_f + '/final.modified_unmodified_new', sep='\t', header=None)
data_uv_f.columns = ['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']
data_uv_f = data_uv_f[['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']]
data_uv_f = data_uv_f.drop_duplicates()
statistics_sum=pd.read_csv(result + '/gini_summary_50_1.csv', sep='\t')
# statistics_sum_new = pd.read_csv(result + '/summary_result_merge_50_1_0.csv',sep='\t')
transcript_all = list(statistics_sum.loc[statistics_sum['num_0.1']>0,'transcript_id'])
# transcript_all =transcript_all[:1000]
# transcript_all=['AT5G26000.1','AT5G26000.2']
# transcript_all = ['AT5G45260.2']
print(len(transcript_all))
gff_path = '/data/TA_QUIZ_RNA_regulation/data/ATH/GFF/Ath_genes.gff'
data_gff = pd.read_csv(gff_path, sep='\t')
# data_gff.columns = ['exosome','name','location','strat','end','.','+/-','num','id']
data_location = pd.read_csv('/data/TA_QUIZ_RNA_regulation/data/ATH/GTF/shape_map/result/transcript_exon_location.csv', sep='\t')
data_all = main_sum(data_uv_z,data_uv_f,transcript_all,data_location,data_gff,result+'/transcript_gini_merge')
data_all.to_csv(result+'/summary_result_merge_50_1.csv',sep='\t',index=False)
hit_level_coil_uv_f = pd.read_csv(uv_f+'cutoff.hit.group',sep='\t')
#hit_level_coil_uv_f = hit_level_coil_uv_f.reset_index()
hit_level_coil_uv_f.columns =['group','transcript_id','modified.median','unmodified.median','modified.sum','unmodified.sum','hit_f']
hit_level_coil_uv_z = pd.read_csv(uv_z+'cutoff.hit.group',sep='\t')
#hit_level_coil_uv_z = hit_level_coil_uv_z.reset_index()
hit_level_coil_uv_z.columns =['group','transcript_id','modified.median','unmodified.median','modified.sum','unmodified.sum','hit_z']
shape_data_coil =pd.merge(pd.merge(data_all,hit_level_coil_uv_f[['transcript_id','hit_f']],on='transcript_id',how='left'),hit_level_coil_uv_z[['transcript_id','hit_z']],on='transcript_id',how='left')
shape_data_coil_2 =shape_data_coil.loc[(shape_data_coil['hit_f']>2)&(shape_data_coil['hit_z']>2),:]
shape_data_coil.to_csv(result+'/summary_result_merge_50_1_new.csv',sep='\t',index=False) 
