我们阅读量破万的综述:RNA-seq这十年(3万字长文综述) 给粉丝朋友们带来了很多理解上的挑战,所以我们开辟专栏慢慢介绍其中的一些概念性的问题,上一期: 箱线图的生物学含义
这一讲我们来说一下limma/voom,edgeR,DESeq2,转录组差异分析的三大R包! 差异分析的第一步是要构建符合不同模型的R对象,主要包括两部分的信息:表达矩阵和分组信息。 这次主要讨论一下limma/voom,edgeR,DESeq2是转录组差异分析的三大R包的表达矩阵和分组矩阵构建,主要针对二分组转录组数据的差异分析。 一、limma和edgeR包的表达矩阵和分组信息1.limma和edgeR包DEGList对象的构建 limma和edgeR包都是由一个研究团队开发,方法之间互相继承。edgeR是专门针对转录组数据开发的,limma包最早是用来进行芯片数据的差异分析,对转录组数据差异分析的功能是后来添加的,表达矩阵的构建方法直接使用edgeR包中的DGEList函数。 DEGList函数的参数示例: DGEList(counts = matrix(0, 0, 0), lib.size = colSums(counts), norm.factors = rep(1,ncol(counts)), samples = NULL, group = NULL, genes = NULL, remove.zeros = FALSE)
使用airway中的转录组表达矩阵来演示 # BiocManager::install(c("airway", "edgeR")) > library(airway) > data(airway) # 获取基因counts矩阵 > exprSet <- assay(airway) > exprSet[1:6,1:6] SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 ENSG00000000003 679 448 873 408 1138 1047 ENSG00000000005 0 0 0 0 0 0 ENSG00000000419 467 515 621 365 587 799 ENSG00000000457 260 211 263 164 245 331 ENSG00000000460 60 55 40 35 78 63 ENSG00000000938 0 0 2 0 1 0 exprSet <- assay(airway) # 获取分组信息 > group_list <- colData(airway)$dex > group_list [1] untrt trt untrt trt untrt trt untrt trt Levels: trt untrt
使用 DEGList函数构建limma和edgeR包需要的输入矩阵: > dge <- edgeR::DGEList(counts=exprSet) > dge An object of class "DGEList" $counts SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 ENSG00000000003 679 448 873 408 1138 1047 770 ENSG00000000005 0 0 0 0 0 0 0 ENSG00000000419 467 515 621 365 587 799 417 ENSG00000000457 260 211 263 164 245 331 233 ENSG00000000460 60 55 40 35 78 63 76 SRR1039521 ENSG00000000003 572 ENSG00000000005 0 ENSG00000000419 508 ENSG00000000457 229 ENSG00000000460 60 64097 more rows ...
$samples group lib.size norm.factors SRR1039508 1 20637971 1 SRR1039509 1 18809481 1 SRR1039512 1 25348649 1 SRR1039513 1 15163415 1 SRR1039516 1 24448408 1 SRR1039517 1 30818215 1 SRR1039520 1 19126151 1 SRR1039521 1 21164133 1
可以看到包含了counts矩阵和一些其他用于差异分析要使用的信息,还可以把分组信息添加进来。 > DEG <- edgeR::DGEList(counts=exprSet,group=factor(group_list)) > DEG An object of class "DGEList" $counts SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 ENSG00000000003 679 448 873 408 1138 1047 770 ENSG00000000005 0 0 0 0 0 0 0 ENSG00000000419 467 515 621 365 587 799 417 ENSG00000000457 260 211 263 164 245 331 233 ENSG00000000460 60 55 40 35 78 63 76 SRR1039521 ENSG00000000003 572 ENSG00000000005 0 ENSG00000000419 508 ENSG00000000457 229 ENSG00000000460 60 64097 more rows ...
$samples group lib.size norm.factors SRR1039508 untrt 20637971 1 SRR1039509 trt 18809481 1 SRR1039512 untrt 25348649 1 SRR1039513 trt 15163415 1 SRR1039516 untrt 24448408 1 SRR1039517 trt 30818215 1 SRR1039520 untrt 19126151 1 SRR1039521 trt 21164133 1
这些使用方法适用于绝大多数limma和edgeR包差异分析的表达矩阵构建。 2.limma和edgeR包分组矩阵的设置 limma和edgeR的假设都是数据符合正态分布,构建线性模型。 使用model.matrix函数构建分组信息的矩阵,就是将分组信息二值化,用0和1构成的矩阵来代表不同的分组信息。 > design <- model.matrix(~0+factor(group_list)) > colnames(design) <- levels(factor(group_list)) > rownames(design) <- colnames(exprSet) > design trt untrt SRR1039508 0 1 SRR1039509 1 0 SRR1039512 0 1 SRR1039513 1 0 SRR1039516 0 1 SRR1039517 1 0 SRR1039520 0 1 SRR1039521 1 0 attr(,"assign") [1] 1 1 attr(,"contrasts") attr(,"contrasts")$`factor(group_list)` [1] "contr.treatment"
需要注意的一点是下面两种构建模型矩阵的区别 > design <- model.matrix(~factor(group_list)) > design (Intercept) factor(group_list)untrt 1 1 1 2 1 0 3 1 1 4 1 0 5 1 1 6 1 0 7 1 1 8 1 0 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$`factor(group_list)` [1] "contr.treatment"
> design <- model.matrix(~0+factor(group_list)) > design factor(group_list)trt factor(group_list)untrt 1 0 1 2 1 0 3 0 1 4 1 0 5 0 1 6 1 0 7 0 1 8 1 0 attr(,"assign") [1] 1 1 attr(,"contrasts") attr(,"contrasts")$`factor(group_list)` [1] "contr.treatment"
第一种方法是将第一列的分组信息作为线性模型的截距,第二列开始依次与第一列比较,通过coef参数可以把差异分析结果依次提取出来。 第二种方法,仅仅是分组信息而已,需要通过makeContrasts函数来制作差异比较矩阵控制。 > # 通过makeContrasts设置需要进行对比的分组 > comp='trt-untrt' > cont.matrix <- makeContrasts(contrasts=c(comp),levels = design) > cont.matrix Contrasts Levels trt-untrt trt 1 untrt -1
二、DESeq2包的表达矩阵和分组信息的构建1.DESeq2包输入文件:DESeqDataSet对象的制作 构建DESeqDataSet函数的参数示例: DESeqDataSet(se, design, ignoreRank = FALSE)
DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE, ignoreRank = FALSE, ...)
DESeqDataSetFromHTSeqCount(sampleTable, directory = ".", design, ignoreRank = FALSE, ...)
DESeqDataSetFromTximport(txi, colData, design, ...)
DESeqDataSet使用示例:从SummarizedExperiment流程中产生的数据导入 > library(airway) > data(airway) > ddsSE <- DESeq2::DESeqDataSet(airway, design = ~ cell + dex) > ddsSE class: DESeqDataSet dim: 64102 8 metadata(2): '' version assays(1): counts rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99 rowData names(0): colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521 colData names(9): SampleName cell ... Sample BioSample >
DESeqDataSetFromMatrix使用示例:从count矩阵中构建DESeqDataSet: > colData <- data.frame(row.names=colnames(exprSet),group_list=group_list) > dds <- DESeq2::DESeqDataSetFromMatrix(countData = exprSet,colData = colData, design = ~ group_list) > dds class: DESeqDataSet dim: 64102 8 metadata(1): version assays(1): counts rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99 rowData names(0): colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521 colData names(1): group_list
DESeqDataSetFromHTSeqCount使用示例:从HTSeqCount中导入数据 # 指定表达矩阵文件所在的文件夹 > directory <- system.file("extdata", package="pasilla", + mustWork=TRUE) > > directory [1] "/Library/Frameworks/R.framework/Versions/3.6/Resources/library/pasilla/extdata" > sampleFiles <- grep("treated",list.files(directory),value=TRUE) > sampleFiles [1] "treated1fb.txt" "treated2fb.txt" "treated3fb.txt" "untreated1fb.txt" [5] "untreated2fb.txt" "untreated3fb.txt" "untreated4fb.txt" > sampleCondition <- sub("(.*treated).*","\\1",sampleFiles) > sampleCondition [1] "treated" "treated" "treated" "untreated" "untreated" "untreated" "untreated" # 构建包含表达矩阵文件信息的数据框 > sampleTable <- data.frame(sampleName = sampleFiles, + fileName = sampleFiles, + condition = sampleCondition) # 导入为DESeqDataSet对象 > ddsHTSeq <- DESeq2::DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, + directory = directory, + design= ~ condition) > ddsHTSeq class: DESeqDataSet dim: 70463 7 metadata(1): version assays(1): counts rownames(70463): FBgn0000003:001 FBgn0000008:001 ... FBgn0261575:001 FBgn0261575:002 rowData names(0): colnames(7): treated1fb.txt treated2fb.txt ... untreated3fb.txt untreated4fb.txt colData names(1): condition
DESeqDataSetFromTximport使用示例:通过Tximport导入不基于比对的基因定量矩阵,主要是以下四个常用软件 Salmon (Patro et al. 2017) Sailfish (Patro, Mount, and Kingsford 2014) kallisto (Bray et al. 2016) RSEM (Li and Dewey 2011)
# 加载R包及示例数据 > library("tximport") > library("readr") > library("tximportData") > dir <- system.file("extdata", package="tximportData") # 设置分组信息 > group_list <- factor(rep(c("A","B"),each=3)) # 获取表达矩阵所在的文件夹,salmon的结果为例 files <- list.files(file.path(dir,"salmon"),pattern = "*quant.sf.gz", recursive = TRUE) full_files_path <- file.path(dir,"salmon",files) # 读入转录本和基因对应列表 > tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv")) Parsed with column specification: cols( TXNAME = col_character(), GENEID = col_character() ) > head(tx2gene) # A tibble: 6 x 2 TXNAME GENEID <chr> <chr> 1 ENST00000456328.2 ENSG00000223972.5 2 ENST00000450305.2 ENSG00000223972.5 3 ENST00000473358.1 ENSG00000243485.5 4 ENST00000469289.1 ENSG00000243485.5 5 ENST00000607096.1 ENSG00000284332.1 6 ENST00000606857.1 ENSG00000268020.3 # txi倒入需要两个参数:表达矩阵所在路径和基因转录本对应的列表 > txi <- tximport(full_files_path, type="salmon", tx2gene=tx2gene) reading in files with read_tsv 1 2 3 4 5 6 summarizing abundance summarizing counts summarizing length > sampleName<- c("ERR188021", "ERR188088", "ERR188288","ERR188297", "ERR188329", "ERR188356") > colData <- cbind(sampleName, group_list) > ddsTxi <- DESeq2::DESeqDataSetFromTximport(txi, + colData = colData, + design = ~ group_list) using counts and average transcript lengths from tximport > ddsTxi class: DESeqDataSet dim: 58288 6 metadata(1): version assays(2): counts avgTxLength rownames(58288): ENSG00000000003.14 ENSG00000000005.5 ... ENSG00000284747.1 ENSG00000284748.1 rowData names(0): colnames: NULL colData names(2): sampleName group_list
2.DESeq2分组信息的设置 DESeq2的差异分析的分组信息设置比较简单,主要通过resuls函数实现 results(object, contrast, name, lfcThreshold = 0, altHypothesis = c("greaterAbs", "lessAbs", "greater", "less"), listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE, alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun, format = c("DataFrame", "GRanges", "GRangesList"), test, addMLE = FALSE, tidy = FALSE, parallel = FALSE, BPPARAM = bpparam(), minmu = 0.5)
resultsNames(object)
removeResults(object)
使用示例: results(dds, contrast=c("condition","C","B"))
dds代表DESeq2得到了差异分析的结果,contrast的输入一个长度为3的向量,与上面构建DESeqDataSet时输入的分组信息对应。 # 如输入的分组信息是如下的因子向量 > group_list [1] A A A B B B Levels: A B # 提取A和B差异分析结果的示例如下,A代表对照组,B代表处理组,注意先后顺序,与edgeR正好相反 results(dds, contrast=c("group_list","B","A"))
三、总结limma,edgeR,DESeq2三大包基本是做转录组差异分析的金标准,大多数转录组的文章都是用这三个R包进行差异分析。 edgeR差异分析速度快,得到的基因数目比较多,假阳性高(实际不差异结果差异)。 DESeq2差异分析速度慢,得到的基因数目比较少,假阴性高(实际差异结果不差异)。 需要注意的是制作分组信息的因子向量是,因子水平的前后顺序,在R的很多模型中,默认将因子向量的第一个水平看作对照组。 四、假如是多个分组呢比如,大家都知道,TCGA的乳腺癌可以分成PAM50的5类,那么差异分析就复杂了,大家可以拿我3年前的WGCNA的教程做例子,下面是分组信息啦 这个时候有两个策略来做差异分析,当然,分组比较多的时候,差异分析并不是最好的策略啦,WGCNA等其它算法更好! 策略1:在分组信息里面挑选 参考我GitHub代码, https://github.com/jmzeng1314/my-R/tree/master/10-RNA-seq-3-groups group_list cont.matrix=makeContrasts(contrasts=c('treat_12-control','treat_2-control'),levels = design) fit2=contrasts.fit(fit,cont.matrix) fit2=eBayes(fit2)
tempOutput = topTable(fit2, coef='treat_12-control', n=Inf) DEG_treat_12_limma_voom = na.omit(tempOutput) write.csv(DEG_treat_12_limma_voom,"DEG_treat_12_limma_voom.csv",quote = F)
tempOutput = topTable(fit2, coef='treat_2-control', n=Inf) DEG_treat_2_limma_voom = na.omit(tempOutput) write.csv(DEG_treat_2_limma_voom,"DEG_treat_2_limma_voom.csv",quote = F)
在提取差异分析结果的时候,需要指定是哪个组和哪个组在进行比较。 值得一提的是, 我的GitHub简直就是宝藏,我上面提到的3年前的WGCNA的教程做例子,最近看到有两个文章就拿同样的数据代码和图片发了一个4分,一个5分的文章!!! 你懂得! 策略2:提取子矩阵和子分组信息 这个很容易理解了,把表达矩阵根据自己想要进行的两两比对来筛选即可,这样就可以多次做差异分析啦,而且保证每次都只有两个分组。 参考资料http://www./1514.html http://www./255.html http://www./1533.html https://ucdavis-bioinformatics-training./2018-June-RNA-Seq-Workshop/thursday/DE.html http://www./packages/release/bioc/html/limma.html http://www./packages/release/bioc/html/edgeR.html http://www./packages/release/bioc/html/DESeq2.html https:///bioc/DESeq2/man/DESeqDataSet.html https:///bioc/DESeq2/man/results.html https:///bioc/limma/man/11RNAseq.html
|