PRO-seq数据分析

背景知识

大多数RNA-seq都是研究不同条件下细胞内mRNA变化。除了基因的编码区（CDS）可以转录成mRNA，基因组上的其他区域也能不同程度地转录（例如poly A,下游区域以及Enhancer），Enhancer可以产生短的且不稳定的RNA来调控转录，而这种调控的错误会引发多种疾病，因此，理解这种调控机理十分重要，然而传统RNA-seq技术在检测这种不稳定的RNA方面效率很低。
而PRO-seq技术就是对传统RNA-seq技术在这方面的改进，它可以富集并且测出刚刚被RNA聚合酶转录出来的新生RNA，并且精度达到一个碱基对。

相关文献：Nature protocol Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq)
https://www.ncbi.nlm.nih.gov/pubmed/?term=27442863
http://www.docin.com/p-1691424856.html

数据来源

文章标题：Nascent RNA sequencing reveals a dynamic global transcriptional response at genes and enhancers to the natural medicinal compound celastrol
数据来源：2017年5月23日冷泉港实验室更新的PRO-seq表达谱
实验设计：
K562细胞系在加入雷公藤红素（中药的一种）后，于0min,10min,20min,40min,60min,160min共六个时间点进行测量，每次2个重复，共计12个数据。
数据下载网址：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96869

创新亮点

传统的RNA-seq研究都是在测量mRNA的量，而通过PRO-seq，可以测出新合成的RNA，并且检测到几分钟后转录水平的变化，这可以更精确地分辨出调控的信号传导通路。

数据预处理

由于此数据原始数据sra太大，没有表达矩阵，只提供了测序序列reads在染色体上分布的位置文件(bw文件)，所以需要进行数据预处理，统计每个基因上reads的数量，作为表达量，此处调用了R语言的rtracklayer包读入bw文件，接下来使用GenomicRanges包统计每个基因上的reads数。

差异表达基因筛选

由于样本量较少，故考虑专为小样本设计的T检验方法。又因为是时间序列，不能轻易划分成两个大组，考虑到时间是一个连续性因素，设计了如下的筛选方法：
1、对0min和10min的两组四个样本进行t检验；接下来依次进行10min和20min；20min和40min；40min和60min；60min和160min的t检验，相邻两个时间点进行t检验共计5次。
2、在上述5次t检验中，如果有4次发生显著性差异，且p<0.1，则说明细胞在加入雷公藤红素后，该基因表达有着显著性变化。
经过这样的筛选后，共有19个差异表达的基因，详见diffgene.txt，第一列为EntrezID。

表达量变化图

对于上述的19个差异表达基因，随机选取几个绘制表达量变化图，纵坐标为表达量，横坐标为不同时间。
可以发现两个基因表达量都在降低，并且在40-60之间有个转录反应的峰，这与文献摘要的This transcriptional response occurred in two major waves, one within 10 minutes, and a second 40-60 minutes after treatment.相对应。

图表 1BTBD2基因

图表 2PEAR1基因

表达谱绘制

首先，对于差异表达基因绘制表达谱，先是只对基因聚类，可以看出来从左到右，颜色由红到绿；这表明随着时间增长，大部分基因的表达量都是由高变低，这与文献摘要中提到的“雷公藤红素会抑制大部分的基因转录”相吻合。

图表 3表达谱单向聚类
接下来，对表达谱双向聚类，可以发现同一时间测得的两个重复试验都能聚类到一起去。

图表 4表达谱双向聚类

分析与讨论

1、对于时间序列数据的处理，这个相邻两组t检验的模型显然还是太过于简单，
2、GEO下载下来的PRO-seq数据是有作为内参对照的spikein数据，可以利用这些内参对照数据对数据进行归一化，将预处理做得更精细。