|
本系列推送旨在帶領(lǐng)生信零基礎(chǔ)的科研人一起掌握Bulk RNA-seq數(shù)據(jù)分析，同時為其他Bulk組學(xué)和單細(xì)胞（核）轉(zhuǎn)錄組測序的數(shù)據(jù)分析奠定基礎(chǔ)。 往期回顧： 第3期. Bulk RNA-seq | 第3期. 基因ID轉(zhuǎn)換，一鍵搞定 今天我們來接觸一下Bulk RNA-seq數(shù)據(jù)分析過程中的關(guān)鍵一環(huán)——差異分析！后文將分享以下4個方面： 三、差異分析的三巨頭是哪些以及有什么區(qū)別？一、什么是差異分析 為了回答這個問題，我們來看下面一個例子。下面的數(shù)據(jù)框就是樣本*基因的count矩陣（列為樣本名，行為基因名，中間的數(shù)字為count[可以簡單理解為某基因在測序時被測到的次數(shù)]），這個矩陣可以展示每個基因在每個樣本中的count，但是不能直接體現(xiàn)每個基因在組間的差異（包括fold change，p value， adjusted p value等），而差異分析就是實現(xiàn)上述轉(zhuǎn)化的過程。
差異分析前：
差異分析后： 
二、為什么要進行差異分析？ 簡單來說：差異分析是其他多種分析的基礎(chǔ)。比如：做富集分析的時候（比如GO或KEGG），我們是想看所有差異基因富集的通路，所以前提是得到組間差異基因的list；再比如做火山圖的時候，需要知道每個基因在兩組間的fold change和adjusted p value。 既然差異分析如此重要，那有哪些方法可以實現(xiàn)差異分析呢？ 三、差異分析的三巨頭是哪些以及有什么區(qū)別？ 到目前為止，Bulk RNA-seq的差異分析主要涉及三種R包（又稱為差異分析的三巨頭）：limma, edgeR, DESeq2。 
下面先提供一下3種R包的官網(wǎng)使用說明：使用手冊：https:///packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf
發(fā)表文章（截至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 發(fā)表文章（截至2023年11月3日被引33074次）： https://academic./bioinformatics/article/26/1/139/182458 https:///packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
發(fā)表文章（截至2023年11月3日被引56577次）：https://genomebiology./articles/10.1186/s13059-014-0550-8?ref=https://四、以DESeq2為例演示全過程 篇幅有限，本文僅演示基于DESeq2的差異分析全過程（基于counts進行分析，不能用tpm、fpkm等歸一化后的數(shù)據(jù)，想獲得練習(xí)數(shù)據(jù)，可在公眾號輸入：Bulk RNA-seq練習(xí)數(shù)據(jù)2）。 install.packages('R.utils')
#BiocManager::install('DESeq2')
library(R.utils)
library(tidyverse)
library(DESeq2)
# 2.1 Raw count讀取 ----
data <- readxl::read_xlsx('./1.數(shù)據(jù)/GSE46224_raw_counts_GRCh38.p13_NCBI.xlsx') # #路徑需要根據(jù)文件的保存位置對應(yīng)更改
# 2.2 注釋信息下載 (用于后續(xù)GeneID轉(zhuǎn)化為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="&") #如果是小鼠的，需要相應(yīng)更改
annot <- data.table::fread(apath, header=T, quote="", stringsAsFactors=F, data.table=F) #快速讀取大型數(shù)據(jù)，但此步耗時較長
rownames(annot) <- annot$GeneID #加上行名，行名按照GeneID這一列
annot_2 <- annot[,c('GeneID','Symbol')]
3. 基因注釋-GeneID轉(zhuǎn)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對應(yīng)同一個symbol，對這些ID取均值
tibble::column_to_rownames("Symbol") %>% #symbol這一列變成行名
dplyr::select(-"GeneID") %>% #刪除"GeneID"列
round() # 對于小數(shù)，需要取整(雖然counts都是整數(shù)，但是前面有一步是取了平均值，所以可能出現(xiàn)小數(shù))
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. 創(chuàng)建分組信息(meta data) meta <- data.frame(sample = colnames(data_anno),
group = c(rep('NF',8),
rep('ICM',8)))
6. DESeq2 #6.1 過濾基因，僅保留count值足夠大的基因，默認(rèn)在70%的樣本中有表達(dá)的基因,為后續(xù)差異分析降級假陰性率 ----
meta$group <- as.factor(meta$group)
##6.1.1 從計數(shù)數(shù)據(jù)創(chuàng)建DESeq2數(shù)據(jù)集
dds <- DESeqDataSetFromMatrix(countData = data_anno,
colData = meta,
design = ~ group)
##6.1.2 得到keep(默認(rèn)在70%的樣本中有表達(dá)的基因,為后續(xù)差異分析降級假陰性率)
###（1）構(gòu)建a DGEList object
dge <- DGEList(counts = data_anno)
#DGEList函數(shù)用于匯總和比較數(shù)據(jù)集之間的差異指標(biāo),具體用法可以問ChatGPT,但是這個里面的內(nèi)容，怎么打開呢？, group = as.factor(meta_selected$group)
###（2）構(gòu)建design matrix
group <- as.factor(meta$group)
design <- model.matrix(~0+group) #創(chuàng)建模型矩陣
colnames(design) <- levels(group) # 設(shè)置列變量的名稱
###（3）根據(jù)a DGEList object和design matrix，篩選基因
keep <- filterByExpr(dge, design)
#6.2 過濾后的count矩陣，后續(xù)兩種差異分析均需要基于此文件進行分析 ----
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結(jié)果(病理-對照)
summary(res)
resOrdered <- res[order(res$padj),] #order根據(jù)padj從小到大排序結(jié)果
tempDEG <- as.data.frame(resOrdered) #把結(jié)果放進數(shù)據(jù)框
DEseq2_DEG <- na.omit(tempDEG) #刪除NA行
DEseq2_DEG$gene <- rownames(DEseq2_DEG) #添加一列,gene
|
