不管在RNA-seq還是ChIP-seq分析中骆莹，都涉及到一個(gè)概念萨蚕，就是標(biāo)準(zhǔn)化。
作為一個(gè)剛?cè)腴T的小白蹄胰，總是聽王院長(zhǎng)和其他熟悉生信分析的同學(xué)們強(qiáng)調(diào)標(biāo)準(zhǔn)化的重要性岳遥。

為什么說(shuō)它重要呢？
因?yàn)闇y(cè)序讀出來(lái)的reads落在一個(gè)基因內(nèi)的counts數(shù)目取決于基因長(zhǎng)度和測(cè)序深度裕寨。如果一個(gè)基因越長(zhǎng)浩蓉，測(cè)序深度越高，那么落在其內(nèi)部的read counts數(shù)目就會(huì)相對(duì)越多宾袜。所以在多個(gè)樣本（不同細(xì)胞系條件捻艳、不同組織、不同器官庆猫、不同個(gè)體等等等）中比較某一個(gè)基因的表達(dá)量认轨，如果不進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化，比較結(jié)果是沒有意義的月培。所以你拿到比對(duì)出來(lái)的read count是不能直接使用的｀易帧（https://www.cnblogs.com/beckygogogo/p/9270698.html）

既然這么重要，所以這個(gè)標(biāo)準(zhǔn)化到底是個(gè)什么鬼杉畜？
簡(jiǎn)單地說(shuō)：標(biāo)準(zhǔn)化的對(duì)象就是基因長(zhǎng)度與測(cè)序深度纪蜒。（https://vip.biotrainee.com/d/63-rpkm-fpkm-rpm-tpm）我在網(wǎng)上搜索了很多關(guān)于標(biāo)準(zhǔn)化的教程，這時(shí)我才發(fā)現(xiàn)學(xué)好數(shù)學(xué)的重要性此叠。纯续。。我是一看到公式就頭大懊鹪ｂ怼！茸歧！不管多么簡(jiǎn)單的公式兔魂，都會(huì)讓我瞬間上頭。举娩。析校。

那么目前有幾種標(biāo)準(zhǔn)化的方法？
目前有這么幾種標(biāo)準(zhǔn)化方法：RPKM,FPKM,RPM,TPM

First of all：
在使用這幾種方法之前铜涉，你都需要拿到一個(gè)東西：就是read count矩陣智玻。
什么是read count? 答：高通量測(cè)序中比對(duì)到轉(zhuǎn)錄本/基因組上的reads數(shù)。（這個(gè)在我之前寫的幾個(gè)實(shí)戰(zhàn)的文章里有提到）

下面就來(lái)具體的了解一下幾種標(biāo)準(zhǔn)化方法：

(1)RPKM:
Reads Per Kilobase of exon model per Million mapped reads
(每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的reads)芙代，主要用來(lái)對(duì)單端測(cè)序（single-end RNA-seq）進(jìn)行定量的方法吊奢。
RPKM(軟件:Range) 的計(jì)算公式：
RPKM= total exon reads/ (mapped reads (Millions) × exon length(KB))
total exon reads：某個(gè)樣本mapping到特定基因的外顯子上的所有的reads。
mapped reads ( Millions ) :某個(gè)樣本的所有reads總和。
exon length( KB )：某個(gè)基因的長(zhǎng)度（外顯子的長(zhǎng)度的總和页滚，以KB為單位）召边。
優(yōu)點(diǎn)：tophat-cufflinks流程固定，應(yīng)用范圍廣裹驰。理論上隧熙，可彌補(bǔ)reads Count的缺點(diǎn)，消除樣本間和基因間差異幻林。（http://www.reibang.com/p/c25e84383ae3）

比如現(xiàn)在有3個(gè)樣品贞盯，4個(gè)不同的基因，1沪饺、2躏敢、3樣本的總reads數(shù)分別是35，45整葡，106件余。每個(gè)基因的read count如下圖所示：

乍一看上去好像read count數(shù)差異還挺大的，但是不要高興的太早遭居。蛾扇。。我們需要先進(jìn)行標(biāo)準(zhǔn)化分析魏滚。根據(jù)公式：RPKM= total exon reads/ (mapped reads (Millions) ×exon length(KB))
我們可以計(jì)算每一個(gè)樣品的RPKM镀首，比如樣品一的RPKM=10/（35×2）=1.43，其他的同理 ：

現(xiàn)在可以看出鼠次，雖然我們拿到的read count值不一樣更哄，但是經(jīng)過(guò)標(biāo)準(zhǔn)化以后，這幾個(gè)基因在3個(gè)樣品的表達(dá)量實(shí)際上是差不多的腥寇。

（2）FPKM
全稱：Fragments Per Kilobase of exon model per Million mapped fragments(每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的fragments)成翩，主要是針對(duì)pair-end測(cè)序表達(dá)量進(jìn)行計(jì)算。
FPKM (軟件：cufflinks)赦役。公式與RPKM一致的麻敌。兩種計(jì)算方法唯一的區(qū)別就是：
一個(gè)是fragment，一個(gè)是read掂摔。

如果是PE(Pair-end)測(cè)序术羔，每個(gè)fragments會(huì)有兩個(gè)reads，F(xiàn)PKM只計(jì)算這一對(duì)reads能比對(duì)到同一個(gè)轉(zhuǎn)錄本的fragments數(shù)量乙漓；而RPKM計(jì)算的是可以比對(duì)到轉(zhuǎn)錄本的reads數(shù)量级历，它不管PE的兩個(gè)reads是否能比對(duì)到同一個(gè)轉(zhuǎn)錄本上。所以叭披，如果是SE(single-end)測(cè)序寥殖，那么FPKM和RPKM計(jì)算的結(jié)果將是一致的。（http://www.reibang.com/p/35e861b76486）

（3）RPM
Reads/Counts of exon model per Million mapped reads (每百萬(wàn)映射讀取的reads).
RPM的計(jì)算公式：
RPM=(total exon reads×10^6) / mapped reads (Millions)
total exon reads：某個(gè)樣本mapping到特定基因的外顯子上的所有的reads；
mapped reads (Millions) :某個(gè)樣本的所有reads總和嚼贡；
優(yōu)點(diǎn)：利于進(jìn)行樣本間比較熏纯。根據(jù)比對(duì)到基因組上的總reads count，進(jìn)行標(biāo)準(zhǔn)化粤策。即：不論比對(duì)到基因組上的總reads count是多少樟澜，都將總reads count標(biāo)準(zhǔn)化為10^6。
缺點(diǎn)：未消除exon長(zhǎng)度造成的表達(dá)差異掐场，難以進(jìn)行樣本內(nèi)exon差異表達(dá)的比較往扔。
參考文章：關(guān)于readsCount贩猎、RPKM/FPKM熊户、RPM、TPM的理解

（4）TPM
Transcripts Per Kilobase of exonmodel per Million mapped reads (每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的Transcripts)吭服，是一種優(yōu)化后的RPKM嚷堡，用于同一物種不同組織的比較。（https://zhuanlan.zhihu.com/p/37196518）
TPM (軟件：RSEM) 的計(jì)算公式：
TPMi={( Ni/Li )*1000000 } / sum( Ni/Li+……..+ Nm/Lm )

Ni：mapping到某個(gè)基因上的read數(shù)艇棕；
Li：某個(gè)基因的外顯子長(zhǎng)度的總和说铃。

在一個(gè)樣本中一個(gè)基因的TPM：先對(duì)每個(gè)基因的read數(shù)用基因的長(zhǎng)度進(jìn)行校正纤泵，之后再用校正后的這個(gè)基因read數(shù)(Ni/Li)與校正后的這個(gè)樣本的所有read數(shù)（sum(Ni/Li+……..+ Nm/Lm)）求商。由此可知，TPM概括了基因的長(zhǎng)度墙杯、表達(dá)量和基因數(shù)目。TPM可以用于同一物種不同組織間的比較钻洒，因?yàn)閟um值總是唯一的酗昼。（https://zhuanlan.zhihu.com/p/37196518）
優(yōu)點(diǎn)：首先消除exon長(zhǎng)度造成的差異，隨后消除樣本間測(cè)序總reads count不同造成的差異闺骚。
缺點(diǎn)：因?yàn)椴皇遣捎帽葘?duì)到基因組上的總reads count彩扔，所以特殊情況下不夠準(zhǔn)確。例如：某突變體對(duì)exon造成整體影響時(shí)僻爽，難以找出差異虫碉。
這里有一個(gè)例子可以比較容易理解：TPM值就是RPKM的百分比嘛！(轉(zhuǎn)自生信菜鳥團(tuán)）

這些方法的區(qū)別是什么胸梆？
FPKM/RPKM：先按列（也就是這個(gè)樣本的總read數(shù)）進(jìn)行標(biāo)化敦捧，之后再對(duì)對(duì)個(gè)基因的長(zhǎng)度進(jìn)行標(biāo)準(zhǔn)化。
TPM：先對(duì)基因長(zhǎng)度進(jìn)行標(biāo)準(zhǔn)化碰镜，之后再對(duì)列（這個(gè)時(shí)候就不再是這個(gè)樣本的總read數(shù)了）進(jìn)行標(biāo)化绞惦。這樣使得最終的TPM矩陣的每列都相同（列和都等于1），也就是說(shuō)每個(gè)樣本中的TPM的和都是一樣的洋措。這樣就會(huì)使得我們更容易去比較同一個(gè)基因在不同樣本中所占的read數(shù)的比例济蝉。而RPKM/FPKM由于最終的表達(dá)值矩陣的列和不同，故而不能直接比較同一個(gè)基因在不同樣本中所占的read數(shù)的比例。
參考文章：
RNA-Seq數(shù)據(jù)標(biāo)準(zhǔn)化方法

哪一種方法適合我的數(shù)據(jù)分析王滤？
說(shuō)了那么多贺嫂，好像每一個(gè)看起來(lái)都還可以，但是不知道在什么條件下用哪一個(gè)方法雁乡。這篇文章：RNA-Seq數(shù)據(jù)標(biāo)準(zhǔn)化方法 是這樣說(shuō)的：
1.學(xué)術(shù)界已經(jīng)不再推薦RPKM第喳、FPKM；（但實(shí)際上還是有很多人在用的踱稍。曲饱。。）那么為什么說(shuō)RPKM和FPKM是錯(cuò)的呢珠月？而又為什么目前很多人都還在用扩淀？這里的一篇文章寫的非常清楚，非常深?yuàn)W啤挎，寫了一大堆公式Ｗぷ弧！庆聘！為什么說(shuō)FPKM和RPKM都錯(cuò)了胜臊？ 這篇文章的詳細(xì)分析過(guò)程我沒看，真的是心有余而力不足啊伙判。象对。。但是結(jié)論還是看得懂的宴抚。
2.比較基因的表達(dá)豐度勒魔，例如哪個(gè)基因在哪個(gè)組織里高表達(dá)，用TPM做均一化處理酱塔；
3.不同組間比較沥邻，找差異基因，先得到read counts羊娃，然后用DESeq2或edgeR唐全，做均一化和差異基因篩選；如果對(duì)比某個(gè)基因的KO組和對(duì)照蕊玷，推薦DESeq2邮利。DESeq2和edgeR不用RPKM/FPKM或TPM做均一化，而是直接用原始的read counts做均一化處理垃帅。（這里要多說(shuō)一句：是王院長(zhǎng)說(shuō)的：如果你的RNA-seq沒有生物學(xué)重復(fù)延届，也就是說(shuō)每一種處理就一個(gè)樣品，那就不能用DESeq2來(lái)做標(biāo)準(zhǔn)化贸诚，要用DESeq做方庭。）
4.如果組間RNA成分差異較大厕吉，例如liver跟spleen比，組織特異性表達(dá)基因在里面搗亂械念；再例如常見的某個(gè)基因KO組與對(duì)照相比头朱，該基因成分在KO組里缺失。就用edgeR/DESeq2作均一化龄减。

下面是我最近在重復(fù)一篇文章的分析時(shí)遇到的問(wèn)題：沒有生物學(xué)重復(fù)的RNA-seq怎么分析项钮？？希停？
我在之前的RNA-seq分析練習(xí)的文章里用到了DESeq2進(jìn)行read count的標(biāo)準(zhǔn)化烁巫，也試過(guò)single-end和pair-end的reads分析。但是前兩天我在試著重復(fù)“Targeting super-enhancer-associated oncogenes in oesophageal squamous cell carcinoma”文章里RNA-seq的部分的時(shí)候（這篇文章的ChIP-seq分析部分見http://www.reibang.com/p/05aae63d7c52）宠能，發(fā)現(xiàn)這篇文章里的RNA-seq數(shù)據(jù)沒有生物學(xué)重復(fù)Ｑ窍丁！棍潘！
這個(gè)就比較郁悶了恃鞋，因?yàn)槲抑坝玫腄ESeq2在分析時(shí)崖媚，在構(gòu)建矩陣的時(shí)候需要輸入replicates的read count文件亦歉。據(jù)王院長(zhǎng)說(shuō)，沒有生物學(xué)重復(fù)的RNA-seq在做標(biāo)準(zhǔn)化的時(shí)候肴楷，就不能用DESeq2，而是要用DESeq赛蔫。另外我又在網(wǎng)上查閱了一些資料：
在3個(gè)生物學(xué)重復(fù)不同處理下泥张，如何篩選差異基因
簡(jiǎn)單使用DESeq2/EdgeR做差異分析
簡(jiǎn)單使用DESeq做差異分析
readbio/DESeq.R
這些文章都表示：對(duì)于沒有生物學(xué)重復(fù)的RNA-seq分析呵恢，需要用DESeq標(biāo)準(zhǔn)化。

下面是我的代碼：

#導(dǎo)入表達(dá)矩陣
>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進(jìn)行標(biāo)準(zhǔn)化
>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小時(shí)樣品與對(duì)照樣品的標(biāo)準(zhǔn)化后的數(shù)據(jù)渗钉，這里conditionA設(shè)置為0小時(shí)處理钞钙，conditionB設(shè)置為2小時(shí)處理樣品
>res_2h = nbinomTest(cds,"KYSE510_0hr", "KYSE510_2hr")
>write.table(res_2h, "KYSE510_2h_vs_0h.res", append = FALSE, quote = FALSE, sep = "\t")

#提取處理4小時(shí)樣品與對(duì)照樣品的標(biāo)準(zhǔn)化后的數(shù)據(jù)
>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小時(shí)樣品與對(duì)照樣品的標(biāo)準(zhǔn)化后的數(shù)據(jù)
>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小時(shí)樣品與對(duì)照樣品的標(biāo)準(zhǔn)化后的數(shù)據(jù)
>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

我查了網(wǎng)上對(duì)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.
以上就是無(wú)生物學(xué)重復(fù)的標(biāo)準(zhǔn)化過(guò)程。