转录组那些事儿 Part I

刘小泽写于18.8.19
之前写过三次关于转录组的实战小文，但是说实话，的确当时还不理解其中的含义，只是想跑个流程，更别说脚本优化了。就像刚买来一部心仪的电子产品，只想着拆开包装，嗅一下新机的味道，看看做工怎样，手感如何，但是它的内涵却不曾理解。后来学了其他的组学，接触久了，发现不深入了解，自己是无论如何都发掘不了它的潜力。

转录组火起来的原因主要是它能结合高通量测序，快速准确地识别转录本，进行表达定量，当然这也是它的核心功能。一般常见的转录组分析是找差异基因、协同变化基因、标记基因、融合基因、新转录本、可变剪切。结合R语言进行可视化、功能注释、网络分析。它既可以单枪匹马，也可以为别的组学打辅助。

好的结果离不开设计

差异来源？

• 样本差异：选取基因表达水平差异明显的不同组织细胞
• 处理差异：处理组和对照组相比，可能由于物理方法（物理损伤、照射等）、化学方法（药品刺激、抑制剂使用等）、生物方法（细菌、病毒侵染等）
• 剂量差异：做药物实验时，需要设计处理组各个药物剂量梯度，来验证药物的作用范围和效力
• 时间差异：时间不同，结果不同，找到特定时间点的影响结果或者分析某个发育现象与时间的关系

能干什么？

1. 比较mRNA水平
   • 同物种同组织不同处理：研究不同条件下基因的表达差异
   • 同物种异组织：不同组织中基因的表达差异
   • 同组织异物种：基因进化上的关系
   • 以上三种可以加上时间变化：研究不同发育/用药时期基因表达差异
2. 基因互作：大量样本建立基因的网络关系，找出通路，发现功能
3. 表达模式：大量样本进行分类，发现与性状相关的基因，对样本进行预测

要测什么？

• mRNA：最常见的转录组测序，建库一般选200-300bp的片段，150或125PE测序

• microRNA：将microRNA分离出来直接单独测序

• IncRNA：长链非编码RNA，有正向、反向转录，要进行链特异性建库

  关于链特异性建库：作用就是测序过程保留转录本的方向信息，让我们知道转录本是来自正义链还是反义链。方便后来区分不同的IncRNA类型以及它的定位，可以更准确获得基因结构和表达信息

怎么提取？

原核生物大部分是核糖体RNA（rRNA），它的mRNA只占据了1-5%。要测它的mRNA，需要先提取纯化。

• 提取：大多数动植物组织样品，使用Trizol试剂即可；多糖含量丰富的植物，可以用多糖多酚试剂盒；脂肪组织可以用QIAGEN的RNeasy lipidmini kit ；
• 纯化：真核生物纯化mRNA，是利用它的3‘端polyA，采用oligoT磁珠将其富集纯化。但是原核没有polyA，因此需要先去除total RNA中的rRNA，需要使用去rRNA试剂盒（Ribo-Zero或KAPA试剂盒），另外对于要测物种IncRNA的实验，如果有适用的试剂盒就用，否则不适用就会影响下游数据质量

检测是否合格的指标：RNA总量、RIN值、OD260/280以及真核28S/18S、原核23S/16S。RIN值越高，28S/18S越接近2表示提取的RNA完整性越好。RIN值高于6.5可以做建库准备，太低影响准确度。有一些昆虫或者水生动物没有28S条带，因此RIN值不能作为参考，但是18S的前基线平稳即可

重复、深度、长度？

重复

• 生物学重复：不同的生物样本做同样的实验
• 技术重复：一个生物样本测定多次

一般生物学重复要保持3以上，另外重复之间的Spearman相关系数要大于0.9（遗传背景不一致的相关系数要大于0.8）

另外，日常公司所说的“样本数量”=生物处理数*重复数，比如你有对照和处理组，各有三个重复，那么就是6个样本，当然测序分析的费用也是按样本收取

有一篇文章用酵母做过实验，doi: 10.1261/rna.053959.115，结果发现，随着重复的增加，找到差异基因越多；要筛到90%以上的差异表达基因，需要30个重复；其实实际分析，也不需要这么多的差异基因，使用合适的软件如edgeR或Deseq，可以控制假阳性率，即使样本重复数比较低，筛出的差异基因可信度还是比较高的。两个结论：生物学重复至少6个；对于每个实验处理要找到大部分（大概是80%以上）差异基因，至少12个生物学重复

根据文章https://doi.org/10.1186/s13059-016-0881-8当重复数为3时，差异倍数（Fold Change，FC）为1.5的基因只能找到43%。另外差异倍数较大的（FC>4）一般都能被覆盖到

1.png

深度

指的是测序得到的总碱基数与待测基因组大小的比值。深度越大，得到的reads条数越多，碱基越多。鉴定表达量中等深度即可，PE 150的reads数20M，测6G数据

结论就是：要为了找差异，花同样钱，多测样本好过加大深度

长度

多测一些长片段可以提高比对效率和转录本识别率，意思就是目标越大越易寻找，尤其对于基因组注释不好或者没有参考基因组的物种，双端测序加上长reads，会增加结果准确性

批次效应？

也许你见过它的分身"batch effect"。它是怎么回事呢？

不同测序平台的数据，同一个平台不同时间或者不同lane上产生的数据，同一样本不同时期，不同试剂做的同一样本等等，这些条件下产生的数据都是批次效应。简单说，就是你的数据量比较大时就容易出现。

是什么？

最常见的就是公司给你测的时候，不是放在一个测序仪上，对照组和实验组分开放置，比如对照是敏感组，实验组是抗性组，先测了抗性组，后来才测了敏感组，结果确实分析出了许多的差异基因。但差异基因是准的吗？会不会有可能是由于敏感组后来测的时候又发生了一些变化呢？说不准，但的确这里上机测序的时间成了干扰因素。要减少批次效应，一大方案就是选择支持更多样本的测序仪，例如NovaSeq一次建库就能容纳96*4个样本

尤其在分析公共数据库时，整合多个不同测序平台数据一起分析差异，这时很容易引入批次效应

如何检测？

怀疑哪个因素产生了干扰，就把它标记出来，比如怀疑时间产生了影响；然后对差异基因聚类分析，看看与时间前后是否相关，若相关就存在批次效应

如何校正？

详细内容可以看这本书https://genomicsclass.github.io/book/第10章

2.png

分析模式

3.png

• 基因组比对：有参考基因组，想分析新转录本

  注释信息不是很完善，或者想找一些非编码RNA

  一般步骤：测序reads比对到基因组=》基于比对结果组装转录本=〉基因/转录本表达定量=》差异或富集分析

• 转录组比对：有参考基因组，分析已知转录本

  参考基因组注释较完善，如人、小鼠等模式生物，带着明确的目的去分析已知基因在样本中的表达。这种模式最简单、快捷

  一般步骤：测序reads与转录本比对=》转录本定量=〉差异或富集分析

• 转录组组装后比对：没有参考基因组，或者有组装质量不好的，需要自己组装转录本（应用场景少，不适合入门）

  一般步骤：测序reads进行De novo组装=》reads比对到组装的转录本=〉转录本表达定量=》差异或富集分析

好的结果离不开分析

一、数据预处理

• 原始数据：Illumina测序仪下机的数据通常为Bcl格式，然后公司使用Bcl2Fastq软件，根据Index序列分割转换成每个样品的Fastq文件，用户拿到的就是fastq格式的原始数据。

• 质控：使用fastqc，查看碱基质量、接头情况、GC含量、序列长度、重复序列等

• 过滤：一般需要去掉低质量碱基或者未识别碱基（N）太多的reads；另外如果测序文库的插入片段太短，比如insert size=50，但采用PE 150测序，read1和read2就会测到接头，所谓的“测通“就是这意思。此时需要去掉接头序列 4.png

  采用Trim Galore或者trimmomatic进行过滤；又或者采用国内的fastp软件，直接将质控和过滤整合，输出文件就是clean data

二、比对

虽然双端测序一次测了头和尾，但是并不能将整个mRNA转录本覆盖。我们结果得到了几百万条reads，要是知道哪条reads来自哪个转录本就能有的放矢，下面计算reads表达量也就知道了某个转录本或者基因的表达量。
将测序reads与参考转录本/基因的比较匹配过程就是比对，可以说mapping 或者alignment

由于DNA转录得到mRNA时将内含子切除，因此mRNA反转录得到的cDNA不一定十分完美的还原回原基因组，相当大一部分会被分开。因此原来为基因组比对设计的软件如bwa可能效果会下降，可以采用专门为转录组开发的比对软件Hisat2、STAR，可以找到剪切位点。当然，如果只为了寻找差异基因，可以用bowtie2以及更快的非比对软件salmon、sailfish、jellyfish。

文章A survey of best practices for RNA-seq data analysis中做了几款主流的比对软件Tophat、Hisat2、STAR评测，结果发现：

• Hisat2比对回去的拼接点比较少，但是找的拼接点成功率高（也就是说，它比对的量少但质优）。
• STAR比对到基因组上唯一的reads数最多，对于双端reads，比对不上的STAR就移除，不会选择妥协只比对单端；它的稳定性最好，体现在处理较长reads和较短reads的结果不会波动太大；STAR容忍性比较好，容易接受错配碱基和soft-clipping（没比对上的不去除，只标记出来），只为帮更多reads“找到家”
• 速度方面，Hisat2比STAR快2.5倍，比上一任tophat快约100倍

看似简单的比对过程，就是帮150bp的reads找到家，其中可能还要让reads付出点“被分割”的代价。但是， 基因组有多大？人类的是3G，也就是30亿碱基，一个150bp对于整个基因组来说，简直不值一提，要从头一个一个比对吗？姑且这样可以，那么我们有多少reads？一般6G数据，150PE，会有20Mreads（=60亿/150/2），也就是2000万条reads。这该怎么办？怎样保证高效和低错误率？

reads是测序仪决定的，是固定的，就是这么多，就是这么长。那么我们只能从参考基因组入手，怎么让他找的快？这里就用到了一个算法——BWT（Burrows–Wheeler_transform），他其实是一个压缩技术，将原来文本转换成一个相似的文本，转换后让相同的字符位置连续【https://blog.csdn.net/luanzheng_365/article/details/78575429】 。使用这个算法，将基因组变成了一个索引index，而我们要查找的序列就是索引中的一个子集，这样比对的任务就不再是将reads从头到尾和基因组去比较，而是转换成了把子集reads和索引index去比较，做到了有的放矢。

比对完我们需要的是bam文件，然后使用bam还可以做其他一些比对统计，或者导入IGV查看

5.png

三、定量（三个水平）

• 基因水平

  常用htseq-count、featureCounts、bedtools(multicov)、Qualimap、GenomicRanges，它们根据read和基因位置的重叠区域判断read的归属。 6.png

  以htseq-count为例，它默认采用union方式进行统计哪些reads分配到哪个基因上。从图上看，软件对前几种都容易判断，但是后三种出现了多比对现象(multi-mapping reads)，这时各个定量软件就出现了差别，htseq-count选择无视这种情况，但是Qualimap选择将geneA、B都算上。这个软件性能不错，但就是速度慢。

  如果有许多样本等待处理，那么featureCounts或许是不错的选择。featureCounts被整合到了subread中，它对于多重比对的reads并不像htseq一样全部丢弃，而是根据比对的不同区域大小比较，最终选择排除、全部或部分计算

  每个样品进行计数后，都是一个个分散的文件，需要将他们合并成一个表达矩阵，行为基因名，列为样品名，中间是计数结果。对于这个矩阵matrix，后期分析需要再标准化【一般产生偏差的因素主要是：基因长度、测序深度、GC含量、测序仪系统误差】。标准化的方式有：RPKM（单端测序用的多）、FPKM（目前主流）、TMM、TPM。也有的软件会自动进行标准化。

  另外，有的软件需要标准化后的矩阵，有的不需要（如DESeq2）

• 外显子水平

  在基因水平之上，又分析的差异的外显子，使用DEXSeq的dexseq_prepare_annotation.py脚本。另外需要提供无重叠的外显子区域gtf文件

• 转录本水平

  基于比对的一般采用Stringtie、eXpress；不比对直接得计数结果的kallisto、sailfish、salmon

明天第二部分是用R语言进行下游的特别分析，包括可视化、差异基因筛选、富集分析等，另外还有实战脚本奉上

---

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!