|
本系列推送旨在带领生信零基础的科研人一起掌握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。 
使用手册：https:///packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf
发表文章（截至2023年11月3日被引25275次）：https://academic./nar/article/43/7/e47/2414268?ref=https://
edgeR：
使用手册： https://www./packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf 发表文章（截至2023年11月3日被引33074次）： https://academic./bioinformatics/article/26/1/139/182458 https:///packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
发表文章（截至2023年11月3日被引56577次）：https://genomebiology./articles/10.1186/s13059-014-0550-8?ref=https://四、以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./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
|
