这章主要内容是不同类型数据之间的整合分析 (data integration)：
1 转录本丰度变化和翻译效率(TE)变化程度之间关系
绘制Log2FoldChange(TE)和Log2FoldChange(RNA-seq)的关系图。
2 转录本结构改变程度和翻译效率TE变化程度之间关系
通过假设检验计算结构改变程度是否影响翻译效率TE变化程度。
3 翻译效率变化基因的motif分析
TE上/下调的3'UTR和5'UTR的motif富集。
1) 转录本丰度变化和翻译效率(TE)变化程度之间关系
我们通过散点图表示转录本丰度光照前后变化和翻译效率(TE)光照前后变化程度之间关系，横轴为log2FC(RNA-seq)，这里仅比较WT。如下图
1.1) 计算所需文件
2.1.RNA-seq Analysis 中(DESeq)计算得到WT的差异表达数据(使用自己计算得到文件)，参考文件/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/5.DESeq2/wt/wt_rawdata.csv
2.2.Ribo-seq Analysis 中计算得到的差异翻译矩阵(使用自己计算得到文件)，参考文件/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv
1.2) 计算过程
这里我们给出参考脚本。
import numpy as np
import pandas as pd
import re
import seaborn as sns
import matplotlib.pyplot as plt
RF_data=pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv',sep='\t')
RF_data=RF_data.reset_index()
RF_data=RF_data.rename(columns={'index':'gene'})
RS_data=pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/5.DESeq2/wt/wt_rawdata.csv',sep=',')
RS_data=RS_data.rename(columns={'Row.names':'gene'})
data=pd.merge(RS_data[['gene','pvalue','padj','log2FoldChange']],RF_data[['gene','pvalue_final','pvalue.adjust','log2FC_TE_final']],on='gene',how='right')
data.columns=['gene','pvalue(RNA-seq)','padj(RNA-seq)','log2FoldChange(RNA-seq)','pvalue(TE)','padj(TE)','log2FoldChange(TE)']
data['group'] = 'darkgray'
result_1=data.loc[(data['pvalue(TE)']<0.05)&(data['log2FoldChange(TE)']>0)&(data['padj(RNA-seq)']>0.05),:]
result_2=data.loc[(data['pvalue(TE)']<0.05)&(data['log2FoldChange(TE)']<0)&(data['padj(RNA-seq)']>0.05),:]
data_= data[["log2FoldChange(RNA-seq)","log2FoldChange(TE)"]]
data_corr = data_.corr().iloc[0,1]
#确定坐标轴显示范围（自己定义调整）
xmin=-3
xmax=6
ymin=-8
ymax=8
#绘制散点图
fig = plt.figure(figsize=plt.figaspect(5/6)) #确定fig比例（h/w）
ax = fig.add_subplot()
ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax), title='')
ax.scatter(data['log2FoldChange(RNA-seq)'], data['log2FoldChange(TE)'], s=15, c=data['group'])
ax.scatter(result_1['log2FoldChange(RNA-seq)'], result_1['log2FoldChange(TE)'], s=20, marker='.',c='#cc0000',label = str(len(result_1))+' TE up mRNAs')
ax.scatter(result_2['log2FoldChange(RNA-seq)'], result_2['log2FoldChange(TE)'], s=20, marker='.',c='steelblue',label = str(len(result_2))+' TE down mRNAs')
ax.spines['right'].set_visible(False) #去掉右边框
ax.spines['top'].set_visible(False) #去掉上边框
ax.spines['bottom'].set_linewidth(2)###设置底部坐标轴的粗细
ax.spines['left'].set_linewidth(2)####设置左边坐标轴的粗细
plt.tick_params(labelsize=15)
ax.set_xticks(range(xmin,xmax,1)) #设置x轴刻度起点和步长
ax.set_yticks(range(ymin,ymax,2)) #设置y轴刻度起点和步长
# font2 = {'family': 'Times New Roman','weight': 'normal','size': 15}
font2 = {'weight': 'normal','size': 18}
plt.xlabel('log2FoldChange(RNA-seq)',font2)
plt.ylabel('log2FoldChange(TE)',font2)
plt.legend(fontsize=10)
font3 = {'weight': 'normal','size': 16}
plt.text(3,4, 'r = '+str(round(data_corr,2)),font3)
plt.tight_layout()
# plt.show()
# 图片输出路径
plt.savefig('logFC_TE_corr.png')
plt.close()
Questions：
1)你得到的转录本丰度变化和翻译效率(TE)变化关系图是怎样的？
2)二者是否存在相关关系，若存在，试解释产生这种相关性的原因。
提示：可以尝试利用Riboseq count matrix计算Riboseq丰度的差异，然后画出了logFC(RNA-seq)和logFC(Ribo-seq)之间的关系。 提示：可思考翻译效率计算方式解释。
2) 结构改变程度是否影响翻译效率TE变化程度
这里我们利用假设检验来检验结构改变程度是否影响翻译效率TE变化程度。
Whether translation up-regulation is related to structural down-regulation
Whether translation down-regulation is related to structural up-regulation
我们需定义结构改变基因为hit level>1且Have structure changed region(|delta gini index|>0.1)。当delta gini index>0.1时为structural up-regulation，当delta gini index<-0.1时为structural down-regulation。
我们需定义翻译效率改变基因为P0.5。
这里我们仅计算WT数据。
2.1) 计算所需文件
2.2.Ribo-seq Analysis计算得到的差异翻译矩阵(使用自己计算得到文件)，参考文件/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv
2.3.SHAPE Data Analysis计算得到的结构改变区域数据(使用自己计算得到文件)，参考文件/data/TA_QUIZ_RNA_regulation/result/PartIII.SHAPE-seq_analysis/merge/merge_data_WT.csv
2.2) 计算过程
这里我们给出参考脚本。
import pandas as pd # Data analysis
import numpy as np # Scientific computing
import matplotlib.pyplot as plt # Plotting
import matplotlib.colors as colors # Coloring
from scipy.stats import chi2_contingency
from scipy.stats import fisher_exact
#####翻译改变###################
exp_data_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/5.DESeq2/wt/wt_rawdata.csv',sep=',')
exp_data_wt =exp_data_wt.rename(columns={'Row.names':'gene','log2FoldChange':'log2FoldChange(EXP)','pvalue':'pvalue(EXP)','padj':'FDR(EXP)'})
exp_data_wt = exp_data_wt[['gene','log2FoldChange(EXP)','pvalue(EXP)','FDR(EXP)']]
ribo_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv',sep='\t')
ribo_wt = ribo_wt.reset_index()
ribo_wt =ribo_wt.rename(columns={'index':'gene','log2FC_TE_final':'log2FoldChange(TE)','pvalue_final':'pvalue(TE)','pvalue.adjust':'FDR(TE)'})
ribo_wt = ribo_wt[['gene','log2FoldChange(TE)','pvalue(TE)','FDR(TE)']]
result = pd.merge(ribo_wt,exp_data_wt,on='gene',how='left')
gene_all_=set(result.loc[((result['log2FoldChange(TE)']>0.5)|(result['log2FoldChange(TE)']<-0.5))&(result['pvalue(TE)']<0.05)&(result['FDR(EXP)']>0.05),'gene'])
gene_TE_up_2=set(result.loc[(result['log2FoldChange(TE)']>0.5)&(result['pvalue(TE)']<0.05)&(result['FDR(EXP)']>0.05),'gene'])
gene_TE_no_up_2=gene_all_-gene_TE_up_2
gene_TE_down_2=set(result.loc[(result['log2FoldChange(TE)']<-0.5)&(result['pvalue(TE)']<0.05)&(result['FDR(EXP)']>0.05),'gene'])
gene_TE_no_down_2=gene_all_-gene_TE_down_2
#####结结构改变###################
#Have structure changed region(|delta gini index|>0.1)
shape_data_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartIII.SHAPE-seq_analysis/merge/merge_data_WT.csv',sep='\t')
shape_data_wt = shape_data_wt.loc[(shape_data_wt['hit_z']>1)|(shape_data_wt['hit_f']>1),:]
shape_data_wt['up_down']= shape_data_wt['delta'].map(lambda x: 'up' if x>0 else 'down')
shape_data_wt['gene']=shape_data_wt['transcript_id'].map(lambda x:x.split('.')[0])
shape_up_1=set(shape_data_wt.loc[(shape_data_wt['up_down']=='up'),'gene'])&gene_all_
shape_no_up_1=gene_all_-shape_up_1
shape_down_1=set(shape_data_wt.loc[(shape_data_wt['up_down']=='down'),'gene'])&gene_all_
shape_no_down_1=gene_all_-shape_down_1
########################
dtest1=np.array([[len(gene_TE_down_2&shape_up_1),len(gene_TE_down_2&shape_no_up_1)],
[len(gene_TE_no_down_2&shape_up_1),len(gene_TE_no_down_2&shape_no_up_1)]])
# k,p,f,expctd =chi2_contingency(dtest)
o,p=fisher_exact(dtest1,alternative='greater')
print('shape up,TE down')
print(p)
dtest3=np.array([[len(gene_TE_up_2&shape_down_1),len(gene_TE_up_2&shape_no_down_1)],
[len(gene_TE_no_up_2&shape_down_1),len(gene_TE_no_up_2&shape_no_down_1)]])
o,p=fisher_exact(dtest3,alternative='greater')
print('shape down,TE up')
print(p)
给出假设检验结果
shape up,TE down
5.387478374073877e-23
shape down,TE up
0.7150732345900125
由此可知，结构变多可导致翻译效率下降。我们接下来用图片表示富集程度。计算脚本参考如下。
import pandas as pd # Data analysis
import numpy as np # Scientific computing
import matplotlib.pyplot as plt # Plotting
import matplotlib.colors as colors # Coloring
exp_data_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/5.DESeq2/wt/wt_rawdata.csv',sep=',')
exp_data_wt =exp_data_wt.rename(columns={'Row.names':'gene','log2FoldChange':'log2FoldChange(EXP)','pvalue':'pvalue(EXP)','padj':'FDR(EXP)'})
exp_data_wt = exp_data_wt[['gene','log2FoldChange(EXP)','pvalue(EXP)','FDR(EXP)']]
exp_wt_list=list(exp_data_wt.loc[exp_data_wt['FDR(EXP)']>0.05,'gene'])
ribo_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv',sep='\t')
ribo_wt = ribo_wt.reset_index()
ribo_wt =ribo_wt.rename(columns={'index':'gene','log2FC_TE_final':'log2FoldChange(TE)','pvalue_final':'pvalue(TE)','pvalue.adjust':'FDR(TE)'})
ribo_wt = ribo_wt[['gene','log2FoldChange(TE)','pvalue(TE)','FDR(TE)']]
shape_data_wt = pd.read_csv('/data/TA_QUIZ_RNA_regulation/result/PartIII.SHAPE-seq_analysis/merge/merge_data_WT.csv',sep='\t')
shape_data_wt = shape_data_wt.loc[(shape_data_wt['hit_z']>1)|(shape_data_wt['hit_f']>1),:]
shape_data_wt['up_down']= shape_data_wt['delta'].map(lambda x: 'up' if x>0 else 'down')
shape_data_wt['gene']=shape_data_wt['transcript_id'].map(lambda x:x.split('.')[0])
shape_up_list=list(shape_data_wt.loc[(shape_data_wt['up_down']=='up'),'gene'])
shape_down_list=list(shape_data_wt.loc[(shape_data_wt['up_down']=='down'),'gene'])
result = pd.merge(ribo_wt,exp_data_wt,on='gene',how='left')
result=result.loc[result['FDR(EXP)']>0.05,:]
result['shape_up']=result['gene'].map(lambda x:1 if x in shape_up_list else 0)
result['shape_down']=result['gene'].map(lambda x:1 if x in shape_down_list else 0)
result['x'] = result['log2FoldChange(TE)']
result['y'] = -np.log10(result['pvalue(TE)'])
#设置显著性阈值
x_threshold=0.5
y_threshold=-np.log10(0.05)
result_2=result.loc[((result['x']<-0.5)|(result['x'] >0.5))&(result['y'] >y_threshold)&(result['shape_up']==1),:]
result_3=result.loc[((result['x']<-0.5)|(result['x'] >0.5))&(result['y'] >y_threshold)&(result['shape_down']==1),:]
#分组为up, normal, down
result['group'] = 'dimgrey'
# result.loc[(result.x > x_threshold)&(result.y > y_threshold),'group'] = 'tab:red' #x=-+x_threshold直接截断
# result.loc[(result.x < -x_threshold)&(result.y > y_threshold),'group'] = 'tab:blue' #x=-+x_threshold直接截断
# result.loc[result.y < y_threshold,'group'] = 'dimgrey' #阈值以下点为灰色
#确定坐标轴显示范围
xmin=-8
xmax=8
ymin=-1
ymax=8
#绘制散点图
fig = plt.figure(figsize=plt.figaspect(5/7)) #确定fig比例（h/w）
ax = fig.add_subplot()
ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax), title='')
ax.scatter(result['x'], result['y'], s=2, c=result['group'])
# ax.scatter(result_2['x'], result_2['y'], s=10, marker='o',c='',edgecolors='tab:purple',label = 'More Structure in UV+')
ax.scatter(result_2['x'], result_2['y'], s=8, marker='o',c='',edgecolors='fuchsia',label = 'More Structure in UV+')
ax.scatter(result_3['x'], result_3['y'], s=5, marker='o',c='',edgecolors='tab:orange',label = 'Less Structure in UV+')
ax.spines['right'].set_visible(False) #去掉右边框
ax.spines['top'].set_visible(False) #去掉上边框
#水平和竖直线
ax.vlines(-x_threshold, ymin, ymax, color='dimgrey',linestyle='dashed', linewidth=1) #画竖直线
ax.vlines(x_threshold, ymin, ymax, color='dimgrey',linestyle='dashed', linewidth=1) #画竖直线
ax.hlines(y_threshold, xmin, xmax, color='dimgrey',linestyle='dashed', linewidth=1) #画竖水平线
plt.tick_params(labelsize=13)
ax.set_xticks(range(xmin,xmax,2)) #设置x轴刻度起点和步长
ax.set_yticks(range(ymin,ymax,2)) #设置y轴刻度起点和步长
# font2 = {'family': 'Times New Roman','weight': 'normal','size': 15}
font2 = {'weight': 'normal','size': 18}
plt.xlabel('Log2FoldChange(TE)',font2)
plt.ylabel('Log10Pvalue(TE)',font2)
plt.legend(fontsize=13)
plt.tight_layout()
plt.savefig('/data/liuxiaofan/test/volcano.png')
plt.close()
结果如下图
3) 翻译效率变化基因的motif分析
本节我们将分别对TE上调和TE下调的所有基因的3‘UTR和5'UTR富集序列motif。
3.1) 数据说明
2.2.Ribo-seq Analysis中计算得到的差异翻译矩阵(使用自己计算得到文件)，参考文件/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/7.TE/wt.0-vs-1.TE_new.csv
所有基因3’UTR/5'UTR的fastq文件
5'UTR的fastq文件：/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/merge_5UTR.fasta
3'UTR的fastq文件：/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/merge_3UTR.fasta
3.2) 提取所有3’UTR/5'UTR的fastq文件(可跳过)
input：
/data/TA_QUIZ_RNA_regulation/data/ATH/GFF/Ath_genes.gff：参考基因组注释文件
/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_nouv/：包含转录本全长序列信息文件
/data/TA_QUIZ_RNA_regulation/data/ATH/GTF/shape_map/result/transcript_exon_location.csv：转录本位置信息文件
import numpy as np
import pandas as pd
col_uv_f = '/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_nouv/'
gff_path = '/data/TA_QUIZ_RNA_regulation/data/ATH/GFF/Ath_genes.gff'
data_gff = pd.read_csv(gff_path, sep='\t')
data_location = pd.read_csv( '/data/TA_QUIZ_RNA_regulation/data/ATH/GTF/shape_map/result/transcript_exon_location.csv',sep='\t')
data_gff=pd.merge(data_gff,data_location[['transcript_id','strand']],on='transcript_id',how='left')
gene=pd.read_csv('/data/TA_QUIZ_RNA_regulation/data/gene_list/wt/ribo_wt_gene_list.txt',sep='\t',header=None)
gene_list=list(gene[0])
data_1 = pd.read_csv(col_uv_f + '/final.modified_unmodified_new', sep='\t')
data_1.columns = ['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']
data_1 = data_1[['transcript_id', 'Nucleotide', 'Modified_mutations', 'Modified_effective_depth',
'Untreated_mutations', 'Untreated_effective_depth', 'R1']]
data_1 = data_1.drop_duplicates()
data_1['gene']=data_1['transcript_id'].map(lambda x:x.split('.')[0])
data=data_1
data_nuc=pd.DataFrame(columns={'gene','transcript_id','5UTR','CDS','3UTR'})
j=0
for i in gene_list:
data_ = data.loc[data['gene'] == i, :]
transcript_id = list(data_['transcript_id'])[0]
data_location_ = data_location.loc[data_location['transcript_id'] == transcript_id, :]
data_gff_ = data_gff.loc[data_gff['transcript_id'] == transcript_id, :]
if len(data_gff_) == 0:
continue
else:
if list(data_gff_['strand'])[0] == '+':
five_UTR_CDS = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'strat'])[0]
three_UTR_CDS = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'end'])[0]
else:
five_UTR_CDS = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'end'])[0]
three_UTR_CDS = list(data_gff_.loc[data_gff_['location'] == 'CDS', 'strat'])[0]
Nuc = list(data_['Nucleotide'])[0]
location = np.array(list(data_location_['location'])[0].split(',')).astype('int')
five_UTR_CDS_ = np.where(location == five_UTR_CDS)[0][0]
three_UTR_CDS_ = np.where(location == three_UTR_CDS)[0][0]
data_5UTR_ = Nuc[:five_UTR_CDS_]
data_CDS_ = Nuc[five_UTR_CDS_:three_UTR_CDS_ + 1]
data_3UTR_ = Nuc[three_UTR_CDS_ + 1:]
data_nuc.loc[j, 'gene'] = i
data_nuc.loc[j, 'transcript_id'] = transcript_id
data_nuc.loc[j, '5UTR'] = data_5UTR_
data_nuc.loc[j, '5UTR_len'] = len(data_5UTR_)
data_nuc.loc[j, 'CDS'] = data_CDS_
data_nuc.loc[j, 'CDS_len'] = len(data_CDS_)
data_nuc.loc[j, '3UTR'] = data_3UTR_
data_nuc.loc[j, '3UTR_len'] = len(data_3UTR_)
j = j + 1
print(j)
data_nuc=data_nuc[['gene','transcript_id','5UTR_len','CDS_len','3UTR_len','5UTR','CDS','3UTR']]
print('ALL',len(data_nuc))
n=15
data_nuc_=data_nuc.loc[data_nuc['5UTR_len']>n]
data_nuc_.index=range(len(data_nuc_))
print('5UTR',len(data_nuc_))
f = open('/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/merge_5UTR.fasta', 'w')
for i in range(len(data_nuc_)):
transcript = data_nuc_.loc[i, 'transcript_id']
nucleotide= data_nuc_.loc[i, '5UTR']
f.write('>' + transcript+'\n')
f.write(nucleotide + '\n')
f.close()
data_nuc_=data_nuc.loc[data_nuc['3UTR_len']>n]
data_nuc_.index=range(len(data_nuc_))
print('3UTR',len(data_nuc_))
f = open('/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/merge_3UTR.fasta', 'w')
for i in range(len(data_nuc_)):
transcript = data_nuc_.loc[i, 'transcript_id']
nucleotide= data_nuc_.loc[i, '3UTR']
f.write('>' + transcript+'\n')
f.write(nucleotide + '\n')
f.close()
output：
merge_5UTR.fasta：汇总的所有5‘UTR的序列信息
merge_3UTR.fasta：汇总的所有3’UTR的序列信息
3.3) 提取TE变化区域的3‘UTR、5’UTR序列信息
请自行根据Part II 得到的TE上/下调的gene list，在merge_5UTR.fasta和merge_3UTR.fasta中筛选，最终得到的即为TE变化基因的5'UTR，3‘UTR序列信息。
output：
merge_TE_up_5UTR.fasta/merge_TE_down_5UTR.fasta：汇总的TE上/下调的所有基因5’UTR的fasta信息。参考：/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/TE_DOWN/merge_TE_down_5UTR.fasta；/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/TE_UP/merge_TE_up_5UTR.fasta
merge_TE_up_3UTR.fasta/merge_TE_down_3UTR.fasta：汇总的TE上/下调的所有基因3’UTR的fasta信息。参考：/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/TE_DOWN/merge_TE_down_3UTR.fasta；/data/TA_QUIZ_RNA_regulation/data/gene_list/motif/TE_UP/merge_TE_up_3UTR.fasta
MEME在线motif分析
我们利用MEME在线版进行motif分析，地址为http://meme-suite.org/tools/meme。 使用MEME的Motif discovery工具，预测输入序列上的motif信息
我们上传merge_TE_up_5UTR.fasta/merge_TE_down_5UTR.fasta后提交运算。
在得到的结果中我们根据Evalue进行筛选，Sites为找到的motif在几个结构改变区域中出现过，由此可以计算出该motif的覆盖率。我们选择覆盖率较高的motif，这里取大于4%。
4) 更多思考和分析
翻译效率（TE）发生变化的基因的3‘UTR和5'UTR分别富集到的motif是什么？有什么不同？有没有之前的研究的支持？
教程中利用 MEME分析的是 sequence motif，我们能否利用一些预测structure motif 的软件和方法探索一些影响和调控翻译的structure motif？
