BioinformaticsRNA-seq生物信息学与算法

Bioconductor分析RNA-seq数据

2017-10-27  本文已影响623人  王诗翔

参考学习《R语言与Bioconductor生物信息学应用》第六章

前言

Y叔的公众号biobabble发过一篇【听说你想学R？】,七月份发的但昨天推送了一次，所以我看到了，看到了对《R语言与Bioconductor生物信息学应用》这本书的强力吐槽，而我发的这篇笔记，连同上次发的Bioconductor分析基因芯片数据都是来自于对该书内容的提炼和学习。所以呢我觉得有必要在这里说几句。

坦白说，这本书确实有很多问题，我自己讲讲几点吧：首先，它有点过时了，公众号评论就有人说基因芯片分析过时了，我个人觉得不客观，但确实这本书有些过时，由此产生了一系列问题，特别要提到的就是代码的可重复性，我个人在运行书中一些代码时很多时候会不work，然后我会自己思考怎么让它work，实在不行就放弃；第二点，它不适合R入门的人，特别是你去看它的第一章，云里雾里，好歹我也摸了半年R了啊，好吧，我就直接跳过了；第三点，它的流程不完整，就是不是很连续的都能让你从头到尾的go through下去，自然心里会感觉不爽快。如此种种。

接着你心里已经有准备听我下面要将的“but”的跳转了，我为什么会学它并做笔记，乃至分享出来？首先，在我需要学习芯片和RNAseq分析的时候身边刚好有这本书，我也不知道是谁的，好吧，那就拿起来看看，发现正是我需要的，所以我看它，第一章看不太懂，没意思，我也不想看，直接跳过从第二章看起，到现在整本书基本看了大半，看过的代码都尝试着去运行过，确实有所收获，所以我会写前言，算是对这本书的客观评价吧；其次，我想谈谈有哪些收获，我本人可以算是有编程基础的，算不得菜鸟，但是对于基因芯片的基础也好，RNAseq乃至基因组分析流程、背景等等可以说是菜鸟的不能再菜鸟，这本书给了我对芯片数据来源、处理流程的一些基本认知，其实这在一些国内资源上是找不到的。有一点我心里非常的不服气，为什么我听说中国做生信非常厉害的人很多，找得到的中文资料却很少？为什么百度其他搜索引擎很强大，一涉及科研领域就非常之垃圾？这也是我挺佩服生信技能树或类似的这样的团队以及相关个人，当我们在喷一些书籍很垃圾，而实际它确实有很多问题的时候，我们能不能贡献自己的力量呢？几个月的学习里我深知自己才能有限，所知甚少，所以不断模仿和记录。我把这些笔记陈列出来并不是它写的有多好，多值得模仿，而是它能够给予我们新的知识，又能够在我们忘记时方便查询。学习必然是一个探寻和思索的过程，技能的掌握它不是一本书可以带给你的，特别是一本技术类书籍，它给了你一个看似可行的方案，你要实际去操作它，然后心里给予评价，在你不确定时，需要多方面整理实践不同的解决方案，然后找到自己的出路，建立自己对该某个问题处理的完整体系。

这篇笔记并不会带你真正学好RNA-seq的分析，至少我看完之后没有，但它确实可以补一些知识的模块。它不适应入门R，也不适应完全模仿做具体的分析，而是适合你在掌握R之后，你在做测序分析之前想了解的一些知识。当你知道它非你所需，你可以完全不看它。

整体而言，这本书非常短，整体评价偏差，但国内在这个方面学习恐怕没有比这好的中文书籍吧？（所以我建议多看网络资料，这也是我在交互学习的，比如生信媛公众号文章目录）我希望那些厉害的人物（教授级人物们）能够多拓展一些中文科研的视界（提升国内人员对生信的整体认知与分析能力，加速学习周期），我也会持续记录这样一些知识，与生信技能树里面的小伙伴一起从不同的研究方向，角度去拓展形形色色的基础知识与理论。我再次强调，我专注于笔记的目的除了自身学习以外，是当你在面对一些概念或者问题的疑惑时，你能键入百度搜索后快速地链接到本文，并从中找到可执行的方案或者帮助你理解，而不是完完全全整篇的通读。而当你确实是需要对所有的知识点有学习的需求时，你再选择读它，不仅仅是这一篇笔记或博文。不要浪费自己的时间，也要耐心地投入自己的时间。

本笔记涉及一些公式，简书暂时还不支持Latex显示公式，更好地阅读请移步博文http://www.flypeom.site/bioinformatics/r/2017/10/27/RNAseq-data-analysis/

如果想了解测序基本原理和知识，查看我整理的测序与平台的几个链接。

使用的学习数据：NCBI SRA (Sequence Read Archive)数据库，数据集编号SRA091277

使用的是菊花转录组样品，分析过程包括原始数据获取、数据清理、质量控制、转录组拼接、转录本定量、标准化和表达差异分析等过程。

样品名称 样品描述 RUN编号 测序长度
T1 处理组（脱水处理3小时） SRR921340 100bp
T2 处理组（脱水处理3小时） SRR921341 100bp
T3 处理组（脱水处理3小时） SRR921342 100bp
T1-1 处理组（脱水处理3小时） SRR921346 51bp
T2-1 处理组（脱水处理3小时） SRR921344 51bp
T3-1 处理组（脱水处理3小时） SRR921345 51bp
CK1 对照组（不做任何处理） SRR921321 100bp
CK2 对照组（不做任何处理） SRR921322 100bp
CK3 对照组（不做任何处理） SRR921324 100bp
CK1-1 对照组（不做任何处理） SRR921336 51bp
CK2-1 对照组（不做任何处理） SRR921337 51bp
CK3-1 对照组（不做任何处理） SRR921338 51bp

高通量测序基础知识

（这里只记录书中重要的知识点并加以理解）

基于第二代测序建立起来的基因组测序、RNA-seq和Small RNA-seq等应用，都由样本收集、文库制备和测序三个过程组成，不同之处在于样品收集和文库制备。

rna-seq.jpg

第二代测序仪（Illumina）测的序列，无论来自DNA-seq文库还是RNA-seq文库，从左到右依次分为3个区域：5'接头（Adapter）区、目标序列区和3'接头区。当多个样品在一个泳道（Lane）中同时测序时，我们可以使用多样品（Multiplex）技术，具体而言是给每个样品分配一个不重复的条形码（Barcode），其实质是一个6-8位的DNA序列，测序后可以通过这个序列将不同的样品分开。单端测序（Single end）指仅从正向测序；双端测序（Paired end）指先从正向测序，然后从反向测序。Barcode则根据另一引物“Sb”独立测序得到。

理论上，由于制备的RNA-seq文库插入长度的峰值常常为200或300bp，所以测序应该只得到文库中目标序列5'端开始的部分（这就是常说的read了，以前总搞不懂~）。但是呢，文库中会有少量目标序列不到测序长度（比如测100bp实际目标序列只有几十），那么测序就可能会测到3'端接头序列，这就是所谓的接头污染。数据预处理时，如果发现接头序列过多，一般是RNA-seq文库插入长度没有控制好；如果出现大量全长的3'接头，一般是接头过量，导致了大量接头自连（Self ligation）。

实际应用中，估计测序深度使用更多的是达到质量标准的有效数据量，而不是原始数据量。当RNA-seq有效数据比例过低时，无法检测一些低丰度的转录本，要考虑重新测序。

测序深度，也叫乘数，指每个碱基被测序的平均次数，是用来衡量测序量的首要参数。测序覆盖度，也叫覆盖率，指被测序到的碱基占全基因组大小的比率。假如用Illumina 2000测序仪完成一次人类基因组（3G大小）单端测序，即可得到300G数据（假设全部是有效数据），估计的测序深度即为100倍（300G/3G），常见表示为100X。将所有读段比对到人类基因组，如果发现只有2.7G的碱基至少有1个读段覆盖到，其实际测序深度为111X（300G/2.7G），测序覆盖度为90%（2.7G/3G）。

不同的测序目的要使用不同的测序策略。如DNA组装使用较多的是2X100bp或更长的双端测序；RNA-seq使用较多的是100bp或更长的单端链特异性测序；small RNA-seq多用50bp单端测序。

从测序得到的读段组装成目标基因组或者转录组的基因策略是比对和拼接，比对是把读段定位到参考基因组或者转录组上，然后再拼接成连续序列；拼接也胶从头组装（Denovo assembly），是在没有参考基因组或者转录组前提下，根据读段之间的重叠区，把所有读段拼接起来，直接获得基因组或者转录组（参考文章）。

转录组比对常用的软件有BWA、Bowtie和Tophat；拼接常用的软件是Trinity。

基因组和转录组组装的不同点：基因组组装希望尽量获得唯一或较少的组装结果，即一致性序列（Consensus sequence）。一致性序列上并不是每个位点都只有一种碱基，它实际上只代表该位点出现频率最高的碱基，存在两种以上碱基的位点叫做杂合位点。注意，过分追求一致性会导致过拼接，即来自不同基因的相似序列会被误拼接到一起。基因组和转录组组装可以用一个非常重要的指标N50来评价，即将所有组装后的序列按照长度从大到小排列，累加值接近所有序列长度总和一半时的那个位置对应的序列长度。N50越大，组装的结果越好，类似的有N90。

测序的质量分数

Phred分数

测序中常用错误概率$P_e$（Error probability）来表示每个核苷酸测量的准确性，还可以赋予一个数值来更简便地表示这个意思，叫做测序质量分数（Quality score）。因为这个分数最开始通过Phred软件从测序仪生成的色谱图中得到的，所以也叫做Phred分数（$Q_{Phred}$）。Phred分数的取值范围是0到93，可以表示很宽的误差范围，即从1（完全错误）到非常低的错误率$10^{-93}$。Phred分数是最基本的质量分数，其他的质量计分标准都来自Phred分数。
$$
Q_{Phred} = -10\times log_{10}P_e
$$

$Q_{Phred}$ $P_{e}$ Base call accuracy
0 1 0
10 0.1 0.9
20 0.01 0.99
30 0.001 0.999
40 0.0001 0.9999
50 0.00001 0.99999

Sanger分数（Phred+33）

Phred分数包括2位数字，还需要用空格分隔，不方便阅读，又要占用大量存储空间，实际上文件中不采用它。为了在文件中方便地表示质量，常常将Phred分数加上33（从33到126变化，ASCII码正好覆盖了可打印区），并用其ASCII码值对应的字符表示，这就是Sanger分数。Sanger分数常用于FASTQ格式的文件。

Illumina/Solexa分数（Phred+64）

分数之间的转换公式：
$$
Q_{Solexa} = -10 \times log_{10}(\frac{P_e}{1-P_e}) \

Q_{Solexa} = 10 \times log_{10}(10^{\frac{Q_{Phred}}{10}}- 1) \

Q_{Phred} = 10 \times log_{10}(10^{\frac{Q_{Solexa}}{10}}- 1) \
$$
Solexa分数的取值范围是-5到62，它在FASTQ文件中需要加上64并转换为相应的ASCII码值（59到126）对应的字符来表示质量。2006年，Illumina公司收购Solexa公司后继续沿用其标准。显著Illmina采用新的标准，采用了Phred分数（范围0-62）加64的质量分数。

Sanger分数（Phred+33）和Illumina分数（Phred+64）是当前应用最为普遍的质量分数系统。

以Phred=20（即常见的Q20标准）为例，其Sanger分数为53，对应数字5；其Illumina分数为84，对应字母T。

Bioconductor中的ShortRead包提供了SolexaQuality和PhredQuality函数分别生成Illumina分数和Sanger分数。

source("http://Bioconductor.org/biocLite.R")
biocLite("ShortRead")
library(ShortRead)
Q=20
PhredQuality(as.integer(Q))
SolexaQuality(as.integer(Q))

> Q=20
> PhredQuality(as.integer(Q))
  A PhredQuality instance of length 1
    width seq
[1]     1 5
> SolexaQuality(as.integer(Q))
  A SolexaQuality instance of length 1
    width seq
[1]     1 T

高通量测序文件格式

FASTQ格式

FASTQ格式是序列文件中常见的一种，它一般包括四部分：第一部分是由“@”开始，后面跟着序列的描述信息（对于高通量数据，这里是读段的名称），这点跟FASTA格式是一样的（起头的符号不一样）；第二部分是DNA序列；第三部分是由“+”号开始，后面或者是读段的名称，或者为空；第四部分是DNA序列上每个碱基的质量分数，每个质量分数对应一个DNA碱基。

Bioconductor中的ShortRead包提供了quality函数可以自动识别FASTQ文件中的质量分数的种类。

我随便写了一个fastq的demo文件，内容

@HWUSI-EAS100R:123:COEPYACXX:6:73:941:1973#0/1
GATTTGGGGTTCAAAGCAGTATCGATCAAAATAGTAAAATCCATTTGTTCAACTCACAGTTT
+ HWUSI-EAS100R:123:***********************
!"************************************************************

进行一些操作：

library(ShortRead)
# 读入FASTQ文件
reads <- readFastq("./demo.fastq")
# 得到质量分数的类型
score_sys <- data.class(quality(reads))
# 得到质量分数
qual <- quality(quality(reads)) # 这里好像一个quality函数就够了，还是尊重原文吧
# 质量分数转为16进制表示
myqual_mat <- charToRaw(as.character(unlist(qual)))

# 如果是Phred+64分数表示系统
if(score_sys=="SFastqQuality"){
  # 显示分数系统类型
  cat("The quality score system is Phred+64", "\n")
  # 输出原始分数值
  strtoi(myqual_mat, 16L)-64
}
# 如果是Phred+33分数表示系统
if(score_sys=="FastqQuality"){
  # 显示分数系统类型
  cat("The quality score system is Phred+33", "\n")
  # 输出原始分数值
  strtoi(myqual_mat, 16L)-33
}


The quality score system is Phred+33 
 [1] 0 1 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9

NCBI中的FASTQ与SRA格式

NCBI的Sequence Read Archive （SRA）数据库，接受FASTQ格式的高通量数据上传，并将分数标准从开始的Illumina分数转换成了Sanger分数。而且，NCBI在读段名称前面增加了数据在SRA库的编号和版本，后面增加了读段的长度。SRA数据库为了节省存储空间，将FASTQ文件压缩为二进制的SRA格式进行保存。用户如果下载SRA格式数据，可以使用工具软件fastq-dump将数据从SRA格式转回FASTQ格式。

QUAL格式文件

Solid测序仪产生分离的序列文件（CSFASTA格式）和质量文件（QUAL格式），两者必须成对出现。QUAL文件采用Phred分数，而且行必须与FASTA文件中的行一一对应。CSFASTA文件与FASTA格式看似相同，但实际不同，Solid不是用核苷酸残基表示序列数据，而是采用了颜色空间的表示方法。Solid的序列文件如果和质量文件合并，可以产生CSFASTAQ格式的文件，也可以根据颜色编码转为真正的FASTAQ格式的文件。

Solid的结果文件转为标准FASTQ格式的文件需要注意两个问题：第一是第1个碱基由于来自测序引物，必须删除；第二就是由于Solid颜色空间的编码是前后依赖的，一旦错一个，会导致后面连续错误，一般都将参考基因组反转为颜色空间编码再进行比对等分析，而不主张将Solid的结果文件直接转换为FASTQ文件。

关于常见的生信数据文件格式，参见生信常见数据格式。

RNA-seq技术的特点

RNA-seq对芯片的优势

RNA-seq检测基因表达主要在7个方面比基因芯片有优势。

基因芯片 RNA-seq
参考序列 需要 不需要
动态范围 小 大
背景噪声 大 小
受降解影响 大 小
序列变异 无法检测 可以检测
转录组方向 不能确定 能确定
可重复性 一般 高

RNA-seq存在的问题

下面看一下如何计算由测序误差引起的Barcode的错误分配，假设Barcode（Barcode唯一确定样本）长度为6个碱基，每个Barcode两两之间两个碱基不同，所有的Barcode都用满，同时假设错误的发生符合二项分布，那么只要2个碱基错误，就会发生一次错误分配，在Illumina测序仪每个碱基的平均错误率0.5%的前提下，下面例子可以计算出一个泳道的错误错误分配概率。

> p = 0.0005
> sum(sapply(2:6, FUN=function(k) choose(6,k)*p^k*(1-p)^(6-k)))
[1] 3.745003e-06

也就是说，在一个泳道，每百万读段就会有370个读段分配错误。如果还考虑DNA簇混合和跳跃PCR引起的Barcode错误分配，这个数值还要高很多。这种分配错误对一般的转录组分析没有影响，但是对一些高灵敏度的突变检测项目影响很大。

Illumine 2000测序仪中，一次运行（Run）可以使用2个流动槽（Flow cell），每个流动槽包括8个泳道（Lane），一个泳道包含2个面（Surface），每个面还有3个条（Swath）也叫列（Column），每一列由16个小区（Tile）组成，后者又由大量DNA簇（Cluster）组成。

Illumine 2000测序仪每次运行（单端测序）理论上可以产生大约30亿个DNA簇，每个DNA簇理论上可以产生一条读段（Read），如果测序长度为100bp，一次运行可以得到3G个读段，其原始数据量为300G个碱基。

RNA-seq数据预处理

RNA-seq数据预处理与基因芯片预处理的目的一致，都是要得到基因表达数据，这里确切说是转录本。但它们的实现细节和方法则有很大不同，下面逐一介绍。

质量控制

数据质量分析报告可以调用ShortRead包中的qa函数得到，另一个常用的工具是FastQC。

下面是一个用qa函数生成质量分析报告的实例，示例数据来自SRA数据库的RNA-seq数据SRR921344。

library(BiocInstaller)
# 安装包
biocLite("SRAdb")
biocLite("R.utils")
# 导入包
library(ShortRead)
library(SRAdb)
library(R.utils)
# download data
db_dir <- "~/Projects/coGECP-pro/db/"
sra_dbname <- paste0(db_dir, 'SRAmetadb.sqlite')    
sra_con <- dbConnect( dbDriver("SQLite"), sra_dbname )
getFASTQfile("SRR921344", sra_con, destDir = "~/Projects/coGECP-pro/fast-format/", srcType = 'ftp', ascpCMD = NULL )
dbDisconnect( sra_con )

# 解压后改名
gunzip("../fast-format/SRR921344.fastq.gz", destname="../fast-format/T2-1.fastq")

# 需要分析的数据文件名称
fastqfile="T2-1.fastq"
# 得到质量分析的记过
qa <- qa(dirPath = "../fast-format/", pattern=fastqfile, type="fastq")

# 输出html格式的分析报告
report(qa, dest="qcReport", type="html")

执行完毕后会得到名为qcReport的结果目录，主要包括质量分析的报告文件“index.html”和文件夹“image”，image文件夹下包括所有的统计信息图示。

比较重要的图示信息有前20个高频出现的读段统计信息，这对于确定接头或者其他污染很重要。FastQC软件带有一个文件包括了常用的Illumina接头序列，会把高频出现的序列比对到这些接头序列，并给予提示；Biocondutor则需要编程比对到NCBI UniVec数据库（http://www.ncbi.nlm.nih.gov/tools/vecscreen/univec）来确定接头序列。

接下来是四种碱基的逐点质量图，该图横坐标是测序的循环数，对应测序时5'端开始的每个位点，纵坐标是所有该位点测量的碱基的质量分数平均数。

质量控制中最重要的一个图是逐点质量图，它与四种碱基逐点质量图的区别在于不区分碱基种类，给出每个位点的测序质量分数的平均数、中位数和上下四分位线。

通过质量控制，可以确定当前样品的数据是否应该保留进入下一步分析或者丢弃。对于通过质量控制的数据，还要进行读段处理，清理后的数据才是实际分析中使用的数据。

读段清理

读段清理主要去掉读段中多个“N”碱基，读段两端的低质量区域（质量分数少于Q20），读段3'端可能混入的接头序列，还有可能污染进来的rRNA和病毒序列。

转录组组装

转录组组装包括从头组装和基于参考序列的组装。从头组装最常用的软件是Trinity，其优点是：不依赖参考序列；能较好地重建变异的、可变剪切的或者来自染色体重组的转录本。从头组装的缺点是会消耗大量的内存资源，测序深度的需求也很高，对测序错误很敏感，高相似度的转录本可能会被误拼到一起。基于参考序列的组装常使用TopHat和Cufflinks（主要用于转录本的识别、定量、标准化与差异分析）两种软件的组合来完成。其优点是：内存需求小；污染影响小，因为污染读段不能比对到参考序列；灵敏度高，所需测序深度低，能检测低丰度的转录本；组装的转录本序列更完整；可以增加参考基因组中的转录本注释。基于参考序列的组装的缺点：严重依赖参考序列及其注释信息等。

转录组定量和标准化

有参的转录组定量，需要利用TopHat比对软件将所有读段比对到参考基因组上，然后由Cuffinks软件完成定量；无参的转录组定量需要利用Bowtie或BWA比对软件将所有读段比对到组装得到的转录组上直接计数。

单端测序读段的比对比较简单，双端测序的质量不一致，往往是反向一端测序质量低，如果按照同样的标准要求两端测序的读段都比对上，会丢失很多比对结果。一般采取的方式是两端读段分别依据不同的标准（例如正向允许错配一个，反向允许错配两个）做单端比对，然后根据两端对齐后中间距离抽取成对的比对结果。

转录组定量得到的基因表达矩阵是简单计数得来的，因此称作原始计数（Raw counts）。有参的转录组定量，可以用Bioconductor的GenomicRanges/Rsamtools软件包中的summarizeOverlaps函数实现。输入的数据为比对后得到的Bam或者Sam文件，经过基因水平或者外显子水平的计数，可以直接输出为某种预定义对象，便于下游软件（如edgeR包和DESeq包）继续处理。

下面例子使用summarizeOverlaps函数从Bam文件中获取原始计数，并输出为edgeR包和DESeq包中的数据对象。

require(BiocInstaller)
# biocLite("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
# biocLite("DESeq")
require(Rsamtools)
require(DESeq)
require(edgeR)
require(pasillaBamSubset)
library(GenomicAlignments)

#此处数据已经不能用了，所以从其他包里弄了个bam格式文件试试
# bamfile <- system.file("extdata", "ex1.bam", package="Rsamtools")
# bf1 <- BamFileList(bamfile, index=character())
# features <- GRanges(seqnames = c())

# from GenomicRanges::GenomicRangesHOWTOs   
reads <- c(untrt1=untreated1_chr4(),
             untrt3=untreated3_chr4())
# single-end reads
# paired-end reads

library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
exbygene <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, "gene")
se <- summarizeOverlaps(exbygene, reads, mode="IntersectionNotEmpty")

# back 
deseq <- newCountDataSet(assays(se)$counts, rownames(colData(se)))
edger <- DGEList(assays(se)$counts, group=rownames(colData(se)))

类似基因芯片，RNA-seq定量后的数据需要标准化，使得所有的样品具有可比性。最常见的一个指标是RPKM（Reads Per Kilo base per Million reads），即每百万读段中来自某一个基因每千碱基长度的读段数目。具体计算时使用比对到某个基因的读段个数除以比对到基因组或者转录组的所有读段个数（以百万为单位），再除以基因或转录本的长度（以KB为单位）。基因表达差异分析的常用软件edgeR软件和DESeq都自带了数据标准化功能，可以直接处理用原始计数表示的基因表达矩阵，RPKM表示的基因表达矩阵的使用主要是为了方便用户直接观察数据。

RNA-seq数据分析

RNA-seq数据与表达谱芯片数据在基因表达差异的显著性分析流程基本相同，不同的地方只在确定差异表达基因方面。

基因表达差异的显著性分析

表达差异分析只对比不同样本之间的同一个转录本，所以不需要考虑转录本长度，只考虑总读段数。一个最简单思想就是，样本测序得到的总读段数（实际上是可以比对到转录组的总读段数）越多，则每个基因分配到的读段越多。因此最简单的标准化因子就是总读段数，用总读段数作标准化的前提是大部分基因的表达是非显著变化的，这与基因芯片中的基本假设相同。但是实际工作中发现很多情况下总读段数主要是一小部分大量表达的基因贡献的。Bullard等（2010）在比较了几种标准化方法的基础上发现在每个泳道内使用非零计数分布的上四分位数（Q75%）作为标准化因子是一种更稳健的选择，总体表现是所研究方法中最优的。

Bioconductor的edgeR包和DESeq包分别提供了上四分位数和中位数来作为标准化因子，就是出于这个思想。

edgeR提供了三种标准化算法，分别是M值加权截断均值法（Weighted trimmed mean of M-values, TMM），相对对数表达值法（Relative log expression, RLE）和上四分位法（Upperquartile），其中TMM是默认设定。这些标准化方法大同小异，其基本思想就是去除表达值较大的少数基因的影响，而保留大部分没有显著变化的基因。

由于基因芯片检测杂交的荧光强度信号是连续值，往往假设它符合正态分布；而RNA-seq测量的是离散值，最简单的假设就是二项分布。由于RNA-seq读段数量非常多，而且一条读段映射到一个给定基因的概率足够小，在实际计算中，二项分布常用它的极限形式泊松分布来替代。泊松分布的一个性质是其方差等于均值，但是当有生物学重复时，RNA-seq数据会表现出比泊松分布期望的更高的变异性，对相当多的基因来说方差可能超过均值，这种现象叫过离散。对过离散数据，基于泊松分布假设的分析容易低估不同生物学重复带来的取样误差而得到过多的假阳性的差异表达基因。为了允许额外的变异，一个自然的想法就是给均值加上一个散度参数，以使方差可以大于均值。于是作为泊松分布的推广，又引入了负二项分布（NB）来作为基本假设，负二项分布是当前基因表达的显著性分析中最常用假设。

现在从百度百科摘取二项分布与负二项分布的公式，以加强理解。

二项分布，即重复n次的贝努利试验，用$\xi$表示随机试验的结果。如果事件发生的概率是p,不发生的概率是1-p，那么N次独立重复试验中发生k次的概率是：
$$
P(\xi=k)=C_n^k * p^k *(1-p)^{n-k}
$$
期望：$E_\xi = np$

方差： $D_\xi = npq$

而负二项分布的公式（概率密度）为：
$$
f(k; r, p)=P_r(x=k) = C_{k+r+1}{r-1}*pk*(1-p)^{n-k}
$$
它表示，已知一个事件在贝努利试验中每次的出现概率是p，在一连串贝努利试验中，一件事刚好在第r+k次试验出现第r次的概率。

写作：X ~ NB(r; p)

基于负二项分布的edgeR、DESeq包是当前最主要的分析程序。在edgeR中，对于任一样品i中的任一个基因g，假设它的分布符合负二项分布。
$$
Y_{gi} = NB(M_ip_{gj}, \phi_g)
$$
其中$M_i$是样品i中的读段总数（实际中是可以比对上的读段总数）；$\phi_g$就是基因g的散度；$p_{gj}$是基因g在某个实验条件j下或者分组j中的相对丰度。第g个基因在某个实验条件j下或分组j中，NB分布的均值为$\mu_{gj}=p_gjM_i$，方差为$\mu_{gj}(1+\mu_{gj}\phi_g)$，对于表达差异分析，需要估计的是散度$\phi_g$，当它趋于0时，负二项分布退化为泊松分布，这时方差退化为第一项$\mu_{gj}$，一般认为来自技术重复，方差第2项$\mu^2_{gj}\phi_g$来自生物学重复。

RNA-seq数据分析往往只有很小的样本量，为每个基因估计一个散度非常困难。比较好的策略是允许不同的基因有不同的个体散度，而这些个体散度的估计可以用一些合适的统计方法借助基因间的信息来改进。相对于edgeR，DESeq默认设置采取了最保守的估计策略，即选取每个基因的经验散度和拟合得到的散度趋势线取值中最大的作为最终的散度估计值，因此DESeq往往选出更少的差异表达基因（为什么呢？）。DESeq由于可以利用同一个样本基因间的数据估计散度，而不一定需要重复样本来计算，因此可以直接用于无重复实验的表达差异分析。

后面一些实例代码不好重复，不再描述。

上一篇下一篇

猜你喜欢

热点阅读