我們在本周推送了兩篇關于TCGA數(shù)據(jù)的使用， 其中伸出我的小腳，將TCGA輕輕絆倒，然后叉腰哈哈笑 一文詳細描述了TCGA數(shù)據(jù)從下載到分析的全過程。在制作表達譜進行下游WGCNA和GSEA分析時，數(shù)據(jù)標準化的工具選擇留下深坑，今日作答。
1背景知識
一般的轉錄組數(shù)據(jù)處理流程是：
測序數(shù)據(jù)是100 bp的單端read，用Rsubread比對到mouse reference genome(mm10), 然后使用featureCounts統(tǒng)計每個基因的count數(shù)。然后用TMM進行標準化，轉換成log2 counts per million.最后用limma包對每個樣本每個基因的平均表達值以觀察水平權重的線性模型進行擬合，并用T檢驗找到不同群體的差異表達基因。以FDR + log2-fold-change對基因排序。 參考文獻：A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1
以前是用基因芯片得到樣本各個基因的表達量，服從正態(tài)分布，但是RNA-Seq，它的抽樣過程是離散的，結果是reads count是矩陣，服從泊松分布，樣本間的差`異是服從負二向分布。
這篇文章中對reads count的基因表達矩陣做的是TMM轉換，trimmed mean of M values，被包裝到了edgeR這個R包里面，是2010年提出的方法，理論上是優(yōu)于RPKM: reads per kilobase per million mapped 這種normalization方法的。但是目前主流其實是DESeq2包的rlog和方差齊性轉換，統(tǒng)計學原理不一樣。
2 rlog和方差齊性轉換區(qū)別
許多常見的多維數(shù)據(jù)探索性分析的統(tǒng)計分析方法，例如聚類和主成分分析要求，在那些同方差性的數(shù)據(jù)表現(xiàn)良好。所謂的同方差性就是雖然平均值不同，但是方差相同。
但是對于RNA-Seq count數(shù)據(jù)而言，當均值增加時，方差期望也會提高。也就說直接對count matrix或標準化count（根據(jù)測序深度調(diào)整）做PCA分析，由于高count在不同樣本間的絕對差值大，也就會對結果有很大影響。簡單粗暴的方法就是對count matrix取log后加1。這個1也是約定俗成，看經(jīng)驗了。
隨便舉個栗子看下效果：
lambda <- 10 ^ seq ( from = - 1 , to = 2 , length = 1000 )
cts <- matrix ( rpois ( 1000 * 100 , lambda ), ncol = 100 )
library ( vsn )
meanSdPlot ( cts , ranks = FALSE )
mark log . cts . one <- log2 ( cts + 1 )
meanSdPlot ( log . cts . one , ranks = FALSE )
mark 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ù)據(jù)。
2.1 到底用啥
數(shù)據(jù)集小于30 -> rlog，大數(shù)據(jù)集 -> VST。
還有這個處理過程不是用于差異檢驗的，在DESeq分析中會自動選擇最合適的所以你更不需要糾結了。只是想需要轉錄組的表達矩陣做PCA，WGCNA,CLUSTERING等分析才用得到。
3 測試數(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
dds <- DESeqDataSetFromMatrix ( countData = exprSet ,
colData = colData ,
design = ~ group_list )
dds
4 normalization對比 library ( "dplyr" )
library ( "ggplot2" )
rld <- rlog ( dds , blind = FALSE )
head ( assay ( rld ), 3 )
vsd <- vst ( dds , blind = FALSE )
head ( assay ( vsd ), 3 )
dds <- estimateSizeFactors ( dds )
df <- bind_rows (
as_data_frame ( log2 ( counts ( dds , normalized = TRUE )[, 1 : 2 ]+ 1 )) %>%
mutate ( transformation = "log2(x + 1)" ),
as_data_frame ( assay ( rld )[, 1 : 2 ]) %>% mutate ( transformation = "rlog" ),
as_data_frame ( assay ( vsd )[, 1 : 2 ]) %>% mutate ( transformation = "vst" ))
colnames ( df )[ 1 : 2 ] <- c ( "x" , "y" )
ggplot ( df , aes ( x = x , y = y )) + geom_hex ( bins = 80 ) +
coord_fixed () + facet_grid ( . ~ transformation )
結果就是轉換后更加集中了
