最近，我们一直在围绕着RNA-Seq的相关技术在进行讨论，我们聊过了RNA-seq的建库polyA+与rRNA-的2种建库策略（第一个是利用特异性和ployA结合的磁珠把含有ployA尾巴的RNA“抓下来”，第二种是把rRNA降解掉，对剩下的RNA进行建库。可见BBQ(生物信息基础问题31)：RNA-Seq建库用哪种策略？），我们也聊过了测序深度与测序仪通量的问题(可见BBQ(生物信息基础问题32)：1个RNA-Seq样本到底要测多少序列？)，那么今天我们来聊一下链特异性RNA-Seq的建库问题。
1.首先我们来思考一下RNA-seq最开始的设计思路是什么？
（1.0版本）我们都知道，在RNA-Seq测的也是DNA，所以在RNA-Seq建库的时候，需要把RNA反转成cDNA。
RNA ---> cDNA
cDNA + adapter ---> RCR 扩增
done
上机测序
但是，这样子做可能会存在一个问题，就是我们不能区分测出来的到底是正链还是负链！那么后期分析的时候就会有坑了！在计算FPKM的时候就不准了！
有没有什么方法可以标记一下链的信息呢？链特异性的RNA-Seq建库就可以！
如图所示，红色的reads表示mapping到了genome的foward strand，蓝色的reads表示mapping到了genome的reverse strand。对于这个gene来说，它本身存在于foward strand，链特意建库的F图就完全能够说明这个现象；但是G图是普通的RNA-Seq建库，也就是不区分链的信息，因此在这个foward strand gene附近的reads有mapping到foward strand的，也有mapping到reverse strand的情况。
图1：F为链特异性建库，G为非特异性建库（Dmitri Parkhomchuk et al., NAR, 2009）
2.链特异性建库的优势与劣势
是在于它能够处理一些gene overlap比较复杂的情况。我们都知道，几乎所有高等生物的gene在genome中的分布都是非均匀的，而且一般都是没有链的偏好性。
用链特异性的建库方法，是能够更加准确进行gene定量的。
至于链特异性建库的劣势，大概有2点吧：1个是贵，1个是操作复杂对于珍贵样品（比如人体组织样品）一旦建库不成功就game over了。
3.链特异性建库的方法有几种，最常用的是哪种?
至于链特异性的建库方法，现在有若干种，但是最常用的其实是2种，1种是RNA ligation methods；1种是dUTP method。
图2.1：RNA ligation （Joshua Z Levin et al., Nature Methods, 2010）
这种放大就是对mRNA先加上一个3'adapter然后然后进行片段筛选（or not）再加上5'adapter.再进行扩增。
除此之外，还有其他的建库方法，如SMART建库等等。
图2.2：其他的RNA ligation method（Joshua Z Levin et al., Nature Methods, 2010）
第二种方法是 dUTP method
图2.3：dUTP RNA library method（Joshua Z Levin et al., Nature Methods, 2010）
这种方法是先合成一条cDNA，在后续的过程中，把dTTP换成dUTP，这样DNA中原来T的位置就全都变成了U，随后两边加上adapter序列。到这里，是存在两种情况的。
第1种情况：5' adapter -> gene -> 3' adapter
第2种情况：5' adapter -> gene 反向互补 -> 3' adapter
但是，无论是这2中情况中的哪一种，gene的反相互链中都原来的T都被替换成了U，这时候我们使用USER酶就可以特异性降解带有U的那条链。因此这时候就只能保存第1种情况了。因此也就可以进行链特异性的测序。
在原文献中，还有一BS转换建库，这里不是很推荐，因为在BS转换的过程中损失量是很大的。有兴趣的同学可以去看看相关文献。
图2.4 BS 建库（Joshua Z Levin et al., Nature Methods, 2010）
4.提问环节
4.1 对于普通的RNA-Seq建库来说，衡量gene表达量，对于foward gene来说，如果有互补链的reads，在计算FPKM的时候需不需要把互补链的reads都去掉？如果是链特异性的RNA-Seq建库呢？为什么？
个人认为，如果是普通的RNA-seq的建库，它没有区分正链负链是无法把互补链的reads去掉的。因为分不清楚，而且在进行FPKM计算的时候，结果可能会偏高！如果是链特异的RNA-seq建库的话，肯定是需要去掉互补链的。两个gene重叠，一个在foward strand，一个在reverse strand是很常见的显现，在计算的时候，算一条链就行了。
4.2 第2个问题，请结合我们的Bioconductor教程来进行。请尝试分析human genome gene 的overlap的情况。
可以利用GenomicRanges这个包进行操作，没下载的首先先安装这个包，再下载一下人类的基因组数据，利用这个包里的reduce函数看一下结果
library(GenomicRanges) ## 包的导入
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("EnsDb.Hsapiens.v75", version = "3.8") ###到这一步是进行人类参考基因组的下载
library(EnsDb.Hsapiens.v75 ) ##参考基因组的导入
ensembl.hg38 =EnsDb.Hsapiens.v7 ##赋值一下
ensembl.hg38.gene=genes(ensembl.hg38)##进行把人类基因取出来
length(ensembl.hg38.gene)##看一下长度
[1] 64102
reduce(ensembl.hg38.gene) ##利用reduce函数看一下结果
length(reduce(ensembl.hg38.gene))
[1] 53110
Ref:
1:生物信息学100个基础问题 —— 第33题 什么是链特异性的RNA-Seq？
2:生物信息学100个基础问题 —— 番外3: R与Bioconductor的入门教学
网友评论