几种标准化方法的学习笔记

不管在RNA-seq还是ChIP-seq分析中,都涉及到一个概念,就是标准化。
作为一个刚入门的小白,总是听王院长和其他熟悉生信分析的同学们强调标准化的重要性。

为什么说它重要呢?
因为测序读出来的reads落在一个基因内的counts数目取决于基因长度和测序深度。如果一个基因越长,测序深度越高,那么落在其内部的read counts数目就会相对越多。所以在多个样本(不同细胞系条件、不同组织、不同器官、不同个体等等等)中比较某一个基因的表达量,如果不进行数据标准化,比较结果是没有意义的。所以你拿到比对出来的read count是不能直接使用的!(https://www.cnblogs.com/beckygogogo/p/9270698.html

既然这么重要,所以这个标准化到底是个什么鬼?
简单地说:标准化的对象就是基因长度与测序深度。(https://vip.biotrainee.com/d/63-rpkm-fpkm-rpm-tpm)我在网上搜索了很多关于标准化的教程,这时我才发现学好数学的重要性。。。我是一看到公式就头大啊!!!不管多么简单的公式,都会让我瞬间上头。。。

那么目前有几种标准化的方法?
目前有这么几种标准化方法:RPKM,FPKM,RPM,TPM

First of all:
在使用这几种方法之前,你都需要拿到一个东西:就是read count矩阵。
什么是read count? 答:高通量测序中比对到转录本/基因组上的reads数。(这个在我之前写的几个实战的文章里有提到)

下面就来具体的了解一下几种标准化方法:

(1)RPKM:
Reads Per Kilobase of exon model per Million mapped reads
(每千个碱基的转录每百万映射读取的reads),主要用来对单端测序(single-end RNA-seq)进行定量的方法。
RPKM(软件:Range) 的计算公式:
RPKM= total exon reads/ (mapped reads (Millions) × exon length(KB))
total exon reads:某个样本mapping到特定基因的外显子上的所有的reads。
mapped reads ( Millions ) :某个样本的所有reads总和。
exon length( KB ):某个基因的长度(外显子的长度的总和,以KB为单位)。
优点:tophat-cufflinks流程固定,应用范围广。理论上,可弥补reads Count的缺点,消除样本间和基因间差异。(//www.greatytc.com/p/c25e84383ae3

比如现在有3个样品,4个不同的基因,1、2、3样本的总reads数分别是35,45,106。每个基因的read count如下图所示:

乍一看上去好像read count数差异还挺大的,但是不要高兴的太早。。。我们需要先进行标准化分析。根据公式:RPKM= total exon reads/ (mapped reads (Millions) ×exon length(KB))
我们可以计算每一个样品的RPKM,比如样品一的RPKM=10/(35×2)=1.43,其他的同理 :

现在可以看出,虽然我们拿到的read count值不一样,但是经过标准化以后,这几个基因在3个样品的表达量实际上是差不多的。

(2)FPKM
全称:Fragments Per Kilobase of exon model per Million mapped fragments(每千个碱基的转录每百万映射读取的fragments),主要是针对pair-end测序表达量进行计算。
FPKM (软件:cufflinks)。公式与RPKM一致的。两种计算方法唯一的区别就是:
一个是fragment,一个是read。

如果是PE(Pair-end)测序,每个fragments会有两个reads,FPKM只计算这一对reads能比对到同一个转录本的fragments数量;而RPKM计算的是可以比对到转录本的reads数量,它不管PE的两个reads是否能比对到同一个转录本上。所以,如果是SE(single-end)测序,那么FPKM和RPKM计算的结果将是一致的。(//www.greatytc.com/p/35e861b76486

(3)RPM
Reads/Counts of exon model per Million mapped reads (每百万映射读取的reads).
RPM的计算公式:
RPM=(total exon reads×10^6) / mapped reads (Millions)
total exon reads:某个样本mapping到特定基因的外显子上的所有的reads;
mapped reads (Millions) :某个样本的所有reads总和;
优点:利于进行样本间比较。根据比对到基因组上的总reads count,进行标准化。即:不论比对到基因组上的总reads count是多少,都将总reads count标准化为10^6。
缺点:未消除exon长度造成的表达差异,难以进行样本内exon差异表达的比较。
参考文章:关于readsCount、RPKM/FPKM、RPM、TPM的理解

(4)TPM
Transcripts Per Kilobase of exonmodel per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts),是一种优化后的RPKM,用于同一物种不同组织的比较。(https://zhuanlan.zhihu.com/p/37196518
TPM (软件:RSEM) 的计算公式:
TPMi={( Ni/Li )*1000000 } / sum( Ni/Li+……..+ Nm/Lm )

Ni:mapping到某个基因上的read数;
Li:某个基因的外显子长度的总和。

在一个样本中一个基因的TPM:先对每个基因的read数用基因的长度进行校正,之后再用校正后的这个基因read数(Ni/Li)与校正后的这个样本的所有read数(sum(Ni/Li+……..+ Nm/Lm))求商。由此可知,TPM概括了基因的长度、表达量和基因数目。TPM可以用于同一物种不同组织间的比较,因为sum值总是唯一的。(https://zhuanlan.zhihu.com/p/37196518
优点:首先消除exon长度造成的差异,随后消除样本间测序总reads count不同造成的差异。
缺点:因为不是采用比对到基因组上的总reads count,所以特殊情况下不够准确。例如:某突变体对exon造成整体影响时,难以找出差异。
这里有一个例子可以比较容易理解:TPM值就是RPKM的百分比嘛!(转自生信菜鸟团)

这些方法的区别是什么?
FPKM/RPKM:先按列(也就是这个样本的总read数)进行标化,之后再对对个基因的长度进行标准化。
TPM:先对基因长度进行标准化,之后再对列(这个时候就不再是这个样本的总read数了)进行标化。这样使得最终的TPM矩阵的每列都相同(列和都等于1),也就是说每个样本中的TPM的和都是一样的。这样就会使得我们更容易去比较同一个基因在不同样本中所占的read数的比例。而RPKM/FPKM由于最终的表达值矩阵的列和不同,故而不能直接比较同一个基因在不同样本中所占的read数的比例。
参考文章:
RNA-Seq数据标准化方法

哪一种方法适合我的数据分析?
说了那么多,好像每一个看起来都还可以,但是不知道在什么条件下用哪一个方法。这篇文章:RNA-Seq数据标准化方法 是这样说的:
1.学术界已经不再推荐RPKM、FPKM;(但实际上还是有很多人在用的。。。)那么为什么说RPKM和FPKM是错的呢?而又为什么目前很多人都还在用?这里的一篇文章写的非常清楚,非常深奥,写了一大堆公式!!!为什么说FPKM和RPKM都错了? 这篇文章的详细分析过程我没看,真的是心有余而力不足啊。。。但是结论还是看得懂的。
2.比较基因的表达丰度,例如哪个基因在哪个组织里高表达,用TPM做均一化处理;
3.不同组间比较,找差异基因,先得到read counts,然后用DESeq2或edgeR,做均一化和差异基因筛选;如果对比某个基因的KO组和对照,推荐DESeq2。DESeq2和edgeR不用RPKM/FPKM或TPM做均一化,而是直接用原始的read counts做均一化处理。(这里要多说一句:是王院长说的:如果你的RNA-seq没有生物学重复,也就是说每一种处理就一个样品,那就不能用DESeq2来做标准化,要用DESeq做。)
4.如果组间RNA成分差异较大,例如liver跟spleen比,组织特异性表达基因在里面捣乱;再例如常见的某个基因KO组与对照相比,该基因成分在KO组里缺失。就用edgeR/DESeq2作均一化。

下面是我最近在重复一篇文章的分析时遇到的问题:没有生物学重复的RNA-seq怎么分析???
我在之前的RNA-seq分析练习的文章里用到了DESeq2进行read count的标准化,也试过single-end和pair-end的reads分析。但是前两天我在试着重复“Targeting super-enhancer-associated oncogenes in oesophageal squamous cell carcinoma”文章里RNA-seq的部分的时候(这篇文章的ChIP-seq分析部分见//www.greatytc.com/p/05aae63d7c52),发现这篇文章里的RNA-seq数据没有生物学重复!!!
这个就比较郁闷了,因为我之前用的DESeq2在分析时,在构建矩阵的时候需要输入replicates的read count文件。据王院长说,没有生物学重复的RNA-seq在做标准化的时候,就不能用DESeq2,而是要用DESeq。另外我又在网上查阅了一些资料:
在3个生物学重复不同处理下,如何筛选差异基因
简单使用DESeq2/EdgeR做差异分析
简单使用DESeq做差异分析
readbio/DESeq.R
这些文章都表示:对于没有生物学重复的RNA-seq分析,需要用DESeq标准化。

下面是我的代码:

#导入表达矩阵
>mycounts <- read.csv("/media/yanfang/FYWD/RNA_seq/matrix/KYSE510_treatment_readcount_unnorm.csv")
>head(mycounts)
>rownames(mycounts)<-mycounts[,1]
>mycounts_withoutx<-mycounts[,-1]
>head(mycounts_withoutx)

#利用DESeq进行标准化
>library(DESeq)
>count<-data.frame(row.names=colnames(mycounts_withoutx),condition = c("KYSE510_0hr","KYSE510_2hr","KYSE510_4hr","KYSE510_6hr","KYSE510_12hr"),libType  =  c("pair-end","pair-end","pair-end","pair-end","pair-end"));
>condition = count$condition;condition;
>cds  =  newCountDataSet(mycounts_withoutx,condition);
>cds  =  estimateSizeFactors(cds);
>sizeFactors(cds);

### This step is telling the DESeq you don't have replicates. 
#Note the option sharingMode="fit-only". Remember that the default, sharingMode="maximum", takes care of
#outliers, i.e., genes with dispersion much larger than the fitted values. Without replicates, we cannot catch such
#outliers and so have to disable this functionality.
>cds  =  estimateDispersions(cds ,method="blind",sharingMode="fit-only")

#提取处理2小时样品与对照样品的标准化后的数据,这里conditionA设置为0小时处理,conditionB设置为2小时处理样品
>res_2h = nbinomTest(cds,"KYSE510_0hr", "KYSE510_2hr")
>write.table(res_2h, "KYSE510_2h_vs_0h.res", append = FALSE, quote = FALSE, sep = "\t")

#提取处理4小时样品与对照样品的标准化后的数据
>THZ1_4h_vs_0h = nbinomTest(cds,"KYSE510_0hr", "KYSE510_4hr")
>write.table(THZ1_4h_vs_0h, "KYSE510_4h_vs_0h.res", append = FALSE, quote = FALSE, sep = "\t")

#提取处理6小时样品与对照样品的标准化后的数据
>THZ1_6h_vs_0h = nbinomTest(cds,"KYSE510_0hr", "KYSE510_6hr")
>write.table(THZ1_6h_vs_0h, "KYSE510_6h_vs_0h.res", append = FALSE, quote = FALSE, sep = "\t")

#提取处理12小时样品与对照样品的标准化后的数据
>THZ1_12h_vs_0h = nbinomTest(cds,"KYSE510_0hr", "KYSE510_12hr")
>write.table(THZ1_12h_vs_0h, "KYSE510_12h_vs_0h.res", append = FALSE, quote = FALSE, sep = "\t")

输出的是csv文件,打开之后是这样的:

> a<- read.table("KYSE510_2h_vs_0h.res",sep="\t",header=T)
> head(a)
               id     baseMean    baseMeanA baseMeanB foldChange log2FoldChange      pval padj
1 ENSG00000000003 1403.3666530 1292.5892915 1514.1440  1.1714038      0.2282385 0.8127726    1
2 ENSG00000000005    0.0000000    0.0000000    0.0000         NA             NA        NA   NA
3 ENSG00000000419 1559.1391516 1717.4263314 1400.8520  0.8156693     -0.2939437 0.7601109    1
4 ENSG00000000457  507.8739093  624.6008395  391.1470  0.6262351     -0.6752237 0.4944952    1
5 ENSG00000000460  946.8086288 1018.7049871  874.9123  0.8588475     -0.2195260 0.8212624    1
6 ENSG00000000938    0.4519543    0.9039086    0.0000  0.0000000           -Inf 0.9919796    1

我查了网上对baseMean,baseMeanA和baseMeanB的解释:
(http://seqanswers.com/forums/showthread.php?t=23204)
baseMean mean normalised counts, averaged over all samples from both conditions.
baseMeanA mean normalised counts from condition A.
baseMeanB mean normalised counts from condition B.
以上就是无生物学重复的标准化过程。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,376评论 6 491
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,126评论 2 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 156,966评论 0 347
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,432评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,519评论 6 385
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,792评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,933评论 3 406
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,701评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,143评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,488评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,626评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,292评论 4 329
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,896评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,742评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,977评论 1 265
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,324评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,494评论 2 348

推荐阅读更多精彩内容