【原】轉(zhuǎn)錄組的高級(jí)分析前該如何標(biāo)準(zhǔn)化數(shù)據(jù)？

健明 2021-07-14

展開全文

我們?cè)诒局芡扑土藘善P(guān)于TCGA數(shù)據(jù)的使用，其中伸出我的小腳,，將TCGA輕輕絆倒，然后叉腰哈哈笑一文詳細(xì)描述了TCGA數(shù)據(jù)從下載到分析的全過程。在制作表達(dá)譜進(jìn)行下游WGCNA和GSEA分析時(shí),，數(shù)據(jù)標(biāo)準(zhǔn)化的工具選擇留下深坑,，今日作答。

1背景知識(shí)

一般的轉(zhuǎn)錄組數(shù)據(jù)處理流程是：

測(cè)序數(shù)據(jù)是100 bp的單端read,，用Rsubread比對(duì)到mouse reference genome(mm10), 然后使用featureCounts統(tǒng)計(jì)每個(gè)基因的count數(shù),。然后用TMM進(jìn)行標(biāo)準(zhǔn)化，轉(zhuǎn)換成log2 counts per million.最后用limma包對(duì)每個(gè)樣本每個(gè)基因的平均表達(dá)值以觀察水平權(quán)重的線性模型進(jìn)行擬合,，并用T檢驗(yàn)找到不同群體的差異表達(dá)基因,。以FDR + log2-fold-change對(duì)基因排序。參考文獻(xiàn)：A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1

以前是用基因芯片得到樣本各個(gè)基因的表達(dá)量,，服從正態(tài)分布,，但是RNA-Seq，它的抽樣過程是離散的,，結(jié)果是reads count是矩陣,，服從泊松分布，樣本間的差`異是服從負(fù)二向分布,。

這篇文章中對(duì)reads count的基因表達(dá)矩陣做的是TMM轉(zhuǎn)換,，trimmed mean of M values，被包裝到了edgeR這個(gè)R包里面,，是2010年提出的方法,，理論上是優(yōu)于RPKM: reads per kilobase per million mapped 這種normalization方法的。但是目前主流其實(shí)是DESeq2包的rlog和方差齊性轉(zhuǎn)換,，統(tǒng)計(jì)學(xué)原理不一樣,。

2 rlog和方差齊性轉(zhuǎn)換區(qū)別

許多常見的多維數(shù)據(jù)探索性分析的統(tǒng)計(jì)分析方法，例如聚類和主成分分析要求,，在那些同方差性的數(shù)據(jù)表現(xiàn)良好,。所謂的同方差性就是雖然平均值不同，但是方差相同,。

但是對(duì)于RNA-Seq count數(shù)據(jù)而言,，當(dāng)均值增加時(shí)，方差期望也會(huì)提高,。也就說直接對(duì)count matrix或標(biāo)準(zhǔn)化count（根據(jù)測(cè)序深度調(diào)整）做PCA分析,，由于高count在不同樣本間的絕對(duì)差值大，也就會(huì)對(duì)結(jié)果有很大影響,。簡(jiǎn)單粗暴的方法就是對(duì)count matrix取log后加1,。這個(gè)1也是約定俗成，看經(jīng)驗(yàn)了,。

隨便舉個(gè)栗子看下效果：

lambda <- 10^seq(from = -1, to = 2, length = 1000)
cts <- matrix(rpois(1000*100, lambda), ncol = 100)
library(vsn)
meanSdPlot(cts, ranks = FALSE)

log.cts.one <- log2(cts + 1)
meanSdPlot(log.cts.one, ranks = FALSE)

DESeq2為count數(shù)據(jù)提供了兩類變換方法,，使得不同均值的方差趨于穩(wěn)定：regularized-logarithm transformation or rlog(Love, Huber, and Anders 2014)和variance stabilizing transformation(VST)(Anders and Huber 2010)用于處理含有色散平均趨勢(shì)負(fù)二項(xiàng)數(shù)據(jù),。

2.1 到底用啥

數(shù)據(jù)集小于30 -> rlog，大數(shù)據(jù)集 -> VST,。

還有這個(gè)處理過程不是用于差異檢驗(yàn)的,，在DESeq分析中會(huì)自動(dòng)選擇最合適的所以你更不需要糾結(jié)了。只是想需要轉(zhuǎn)錄組的表達(dá)矩陣做PCA,，WGCNA,CLUSTERING等分析才用得到,。

3 測(cè)試數(shù)據(jù)

suppressPackageStartupMessages(library(airway))
suppressPackageStartupMessages(library(DESeq))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(edgeR))
suppressPackageStartupMessages(library(pasilla))
data(pasillaGenes)
data(airway)
exprSet=counts(pasillaGenes)
group_list=pData(pasillaGenes)[,2]
geneLists=row.names(exprSet)
keepGene=rowSums(edgeR::cpm(exprSet)>0) >=2
table(keepGene);dim(exprSet)

keepGene
FALSE TRUE
3545 10925
[1] 14470 7

dim(exprSet[keepGene,])

[1] 10925 7

exprSet=exprSet[keepGene,]
rownames(exprSet)=geneLists[keepGene]
(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list) )

group_list
treated1fb treated
treated2fb treated
ttreated3fb treated
tuntreated1fb untreated
tuntreated2fb untreated
tuntreated3fb untreated
tuntreated4fb untreated