转录组的高级分析前该如何标准化数据？

我们在本周推送了两篇关于TCGA数据的使用， 其中伸出我的小脚，将TCGA轻轻绊倒，然后叉腰哈哈笑 一文详细描述了TCGA数据从下载到分析的全过程。在制作表达谱进行下游WGCNA和GSEA分析时，数据标准化的工具选择留下深坑，今日作答。

1背景知识

一般的转录组数据处理流程是：

测序数据是100 bp的单端read，用Rsubread比对到mouse reference genome(mm10), 然后使用featureCounts统计每个基因的count数。然后用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

以前是用基因芯片得到样本各个基因的表达量，服从正态分布，但是RNA-Seq，它的抽样过程是离散的，结果是reads count是矩阵，服从泊松分布，样本间的差`异是服从负二向分布。

这篇文章中对reads count的基因表达矩阵做的是TMM转换，trimmed mean of M values，被包装到了edgeR这个R包里面，是2010年提出的方法，理论上是优于RPKM: reads per kilobase per million mapped 这种normalization方法的。但是目前主流其实是DESeq2包的rlog和方差齐性转换，统计学原理不一样。

2 rlog和方差齐性转换区别

许多常见的多维数据探索性分析的统计分析方法，例如聚类和主成分分析要求，在那些同方差性的数据表现良好。所谓的同方差性就是虽然平均值不同，但是方差相同。

但是对于RNA-Seq count数据而言，当均值增加时，方差期望也会提高。也就说直接对count matrix或标准化count（根据测序深度调整）做PCA分析，由于高count在不同样本间的绝对差值大，也就会对结果有很大影响。简单粗暴的方法就是对count matrix取log后加1。这个1也是约定俗成，看经验了。

随便举个栗子看下效果：

1. lambda <- 10^seq(from = -1, to = 2, length = 1000)

3. cts <- matrix(rpois(1000*100, lambda), ncol = 100)

5. library(vsn)

7. meanSdPlot(cts, ranks = FALSE)

mark

1. log.cts.one <- log2(cts + 1)

3. meanSdPlot(log.cts.one, ranks = FALSE)

mark

DESeq2为count数据提供了两类变换方法，使得不同均值的方差趋于稳定：regularized-logarithm transformation or rlog(Love, Huber, and Anders 2014)和variance stabilizing transformation(VST)(Anders and Huber 2010)用于处理含有色散平均趋势负二项数据。

2.1 到底用啥

数据集小于30 -> rlog，大数据集 -> VST。

还有这个处理过程不是用于差异检验的，在DESeq分析中会自动选择最合适的所以你更不需要纠结了。只是想需要转录组的表达矩阵做PCA，WGCNA,CLUSTERING等分析才用得到。

3 测试数据

1. suppressPackageStartupMessages(library(airway))

3. suppressPackageStartupMessages(library(DESeq))

5. suppressPackageStartupMessages(library(DESeq2))

7. suppressPackageStartupMessages(library(edgeR))

9. suppressPackageStartupMessages(library(pasilla))  

11. data(pasillaGenes)

13. data(airway)

15. exprSet=counts(pasillaGenes)

17. group_list=pData(pasillaGenes)[,2]

19. geneLists=row.names(exprSet)

21. keepGene=rowSums(edgeR::cpm(exprSet)>0) >=2

23. table(keepGene);dim(exprSet)

keepGene

FALSE TRUE

3545 10925

[1] 14470 7

1. dim(exprSet[keepGene,])

[1] 10925 7

1. exprSet=exprSet[keepGene,]

3. rownames(exprSet)=geneLists[keepGene]

5. (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

1. dds <- DESeqDataSetFromMatrix(countData = exprSet,

3.                              colData = colData,

5.                              design = ~ group_list)

7. dds

4 normalization对比

1. library("dplyr")

3. library("ggplot2")

5. rld <- rlog(dds, blind = FALSE)

7. head(assay(rld), 3)

1. vsd <- vst(dds, blind = FALSE)

3. head(assay(vsd), 3)

1. dds <- estimateSizeFactors(dds)

4. df <- bind_rows(

6.  as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%

8.    mutate(transformation = "log2(x + 1)"),

10.  as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"),

12.  as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"))

15. colnames(df)[1:2] <- c("x", "y")  

18. ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +

20.  coord_fixed() + facet_grid( . ~ transformation)

结果就是转换后更加集中了