本系列推送旨在带领生信零基础的科研人一起掌握Bulk RNA-seq数据分析,同时为其他Bulk组学和单细胞(核)转录组测序的数据分析奠定基础。
往期回顾:
第3期. Bulk RNA-seq | 第3期. 基因ID转换,一键搞定
今天我们来接触一下Bulk RNA-seq数据分析过程中的关键一环——差异分析!后文将分享以下4个方面:
一、什么是差异分析
为了回答这个问题,我们来看下面一个例子。下面的数据框就是样本*基因的count矩阵(列为样本名,行为基因名,中间的数字为count[可以简单理解为某基因在测序时被测到的次数]),这个矩阵可以展示每个基因在每个样本中的count,但是不能直接体现每个基因在组间的差异(包括fold change,p value, adjusted p value等),而差异分析就是实现上述转化的过程。
差异分析前:
差异分析后:
二、为什么要进行差异分析?
简单来说:差异分析是其他多种分析的基础。比如:做富集分析的时候(比如GO或KEGG),我们是想看所有差异基因富集的通路,所以前提是得到组间差异基因的list;再比如做火山图的时候,需要知道每个基因在两组间的fold change和adjusted p value。
既然差异分析如此重要,那有哪些方法可以实现差异分析呢?
三、差异分析的三巨头是哪些以及有什么区别?
到目前为止,Bulk RNA-seq的差异分析主要涉及三种R包(又称为差异分析的三巨头):limma, edgeR, DESeq2。
edgeR:
使用手册:
https://www.bioconductor.org/packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
发表文章(截至2023年11月3日被引33074次):
https://academic.oup.com/bioinformatics/article/26/1/139/182458
四、以DESeq2为例演示全过程
篇幅有限,本文仅演示基于DESeq2的差异分析全过程(基于counts进行分析,不能用tpm、fpkm等归一化后的数据,想获得练习数据,可在公众号输入:Bulk RNA-seq练习数据2)。
install.packages('R.utils')
#BiocManager::install('DESeq2')
library(R.utils)
library(tidyverse)
library(DESeq2)
# 2.1 Raw count读取 ----
data <- readxl::read_xlsx('./1.数据/GSE46224_raw_counts_GRCh38.p13_NCBI.xlsx') # #路径需要根据文件的保存位置对应更改
# 2.2 注释信息下载 (用于后续GeneID转化为Symbol)----
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
apath <- paste(urld, "type=rnaseq_counts", "file=Human.GRCh38.p13.annot.tsv.gz", sep="&") #如果是小鼠的,需要相应更改
annot <- data.table::fread(apath, header=T, quote="", stringsAsFactors=F, data.table=F) #快速读取大型数据,但此步耗时较长
rownames(annot) <- annot$GeneID #加上行名,行名按照GeneID这一列
annot_2 <- annot[,c('GeneID','Symbol')]
3. 基因注释-GeneID转Gene Symbol
data_anno <- data %>%
as.data.frame() %>%
dplyr::inner_join(annot_2, by = "GeneID") %>% #data和annot_2针对GeneID这一列取交集
na.omit() %>% #去除NA
stats::aggregate(. ~ Symbol, data =., FUN = mean) %>% #对于多个ID对应同一个symbol,对这些ID取均值
tibble::column_to_rownames("Symbol") %>% #symbol这一列变成行名
dplyr::select(-"GeneID") %>% #删除"GeneID"列
round() # 对于小数,需要取整(虽然counts都是整数,但是前面有一步是取了平均值,所以可能出现小数)
4. 选择需要分析的基因,一般只分析protein_coding RNA,lncRNA, miRNA。本案例选择protein-coding。
gene_selected <- rownames(data_anno[which(annot$GeneType %in% c('protein-coding')),])
data_anno <- data_anno[gene_selected,]
5. 创建分组信息(meta data)
meta <- data.frame(sample = colnames(data_anno),
group = c(rep('NF',8),
rep('ICM',8)))
6. DESeq2
#6.1 过滤基因,仅保留count值足够大的基因,默认在70%的样本中有表达的基因,为后续差异分析降级假阴性率 ----
meta$group <- as.factor(meta$group)
##6.1.1 从计数数据创建DESeq2数据集
dds <- DESeqDataSetFromMatrix(countData = data_anno,
colData = meta,
design = ~ group)
##6.1.2 得到keep(默认在70%的样本中有表达的基因,为后续差异分析降级假阴性率)
###(1)构建a DGEList object
dge <- DGEList(counts = data_anno)
#DGEList函数用于汇总和比较数据集之间的差异指标,具体用法可以问ChatGPT,但是这个里面的内容,怎么打开呢?, group = as.factor(meta_selected$group)
###(2)构建design matrix
group <- as.factor(meta$group)
design <- model.matrix(~0+group) #创建模型矩阵
colnames(design) <- levels(group) # 设置列变量的名称
###(3)根据a DGEList object和design matrix,筛选基因
keep <- filterByExpr(dge, design)
#6.2 过滤后的count矩阵,后续两种差异分析均需要基于此文件进行分析 ----
data_anno <- data_anno[keep,]
dds <- dds[keep,]
#6.3 差异分析 ----
dds <- DESeq(dds,quiet = F)
res <- results(dds,contrast=c("group", "ICM", "NF")) #指定提取为exp/ctr结果(病理-对照)
summary(res)
resOrdered <- res[order(res$padj),] #order根据padj从小到大排序结果
tempDEG <- as.data.frame(resOrdered) #把结果放进数据框
DEseq2_DEG <- na.omit(tempDEG) #删除NA行
DEseq2_DEG$gene <- rownames(DEseq2_DEG) #添加一列,gene
联系客服