limma差異表達分析
本篇筆記的內容是在R語言中利用limma包進行差異表達分析，主要針對轉錄組測序得到的基因表達數(shù)據進行下游分析斋扰，并將分析結果可視化渡八，繪制火山圖和熱圖
[TOC]
基因表達差異分析是我們做轉錄組最關鍵根本的一步，不管哪種差異分析传货，其本質都是廣義線性模型屎鳍，limma也是廣義線性模型的一種，其對每個gene的表達量擬合一個線性方程问裕。
limma包是2015年發(fā)表在Nucleic Acids Resarch一個做差異分析的工具逮壁，目前引用次數(shù)高達七千多次，最流行的差異分析軟件之一就是limma僻澎。
環(huán)境部署與安裝
- 安裝limma包
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("limma")
- 安裝ggVolcano包
# install.packages("devtools")
devtools::install_github("BioSenior/ggVolcano")
- 安裝ggplot2包
#直接安裝tidyverse貌踏，一勞永逸（推薦,數(shù)據分析大禮包）
install.packages("tidyverse")
#直接安裝ggplot2
install.packages("ggplot2")
#從Github上安裝最新的版本十饥，先安裝devtools(如果沒安裝的話)
devtools::install_github("tidyverse/ggplot2")
輸入數(shù)據準備
-
樣本信息表：第一列為樣本名稱ID窟勃，第二列為樣本的分組（CK、HT）逗堵，
sampleinfo.csv格式如下這一步需要注意：每個分組下至少有兩個生物學重復秉氧，即至少4行數(shù)據，如果沒有重復的話在后續(xù)貝葉斯檢驗會出錯蜒秤，原因是該模型利用統(tǒng)計假設檢驗汁咏，多個重復能夠評估變異性，而如果僅有一組數(shù)據作媚，則無法檢驗攘滩。
sample group A01 CK A02 HT … …
-
表達矩陣：每行是一個基因，每列是一個樣本纸泡，表達數(shù)據為TPM漂问，
data.csv格式如下A01 A02 … gene1 xx xx … gene2 xx xx … … … … …
差異表達分析過程
準備環(huán)節(jié)
載入R包，設置參數(shù)，其中job變量用于項目輸出文件前綴標識蚤假，可以自定義修改栏饮。
setwd("D:/LABdata")
options(stringsAsFactors = F)
rm(list=ls()) #清空變量
job <- "test" #設定項目名稱
library(limma)
library(ggplot2) #用于繪制火山圖
library(ggVolcano)
數(shù)據導入
導入樣本信息和表達量數(shù)據，然后進行刪除表達量之和為0基因磷仰、log2化袍嬉、替換異常值等步驟，得到原始數(shù)據矩陣灶平。
# 輸入表達矩陣和分組文件 -------------------------------------------------------------
expr_data<-read.csv("data.csv",header = T,row.names = 1,sep = ",") #輸入文件TPM原始值伺通，行名是基因，列名是樣本
expr_data <- expr_data[which(rowSums(expr_data)!=0),] #刪除表達量為0的基因
expr_data = log2(expr_data) #log化處理
expr_data[expr_data == -Inf] = 0 #將log化后的負無窮值替換為0
group<-read.csv("sampleinfo.csv",header = T,row.names = 1,sep = ",") #輸入文件逢享，樣本信息表泵殴，包含分組信息
構建分組矩陣
根據樣本的分組信息，構建分組矩陣拼苍，最終得到的design矩陣由0和1構成笑诅，為斜對角矩陣。
# 構建分組矩陣--design ---------------------------------------------------------
design <- model.matrix(~0+factor(group$group))
colnames(design) <- levels(factor(group$group))
rownames(design) <- colnames(expr_data)
構建比較矩陣
設置樣本的比較方式疮鲫，這里為CK對照比HT處理吆你，該步驟生成的文件為1和-1構成的矩陣。
# #構建比較矩陣——contrast -------------------------------------------------------
contrast.matrix <- makeContrasts(CK-HT,levels = design) #根據實際的樣本分組修改俊犯，這里對照組CK妇多，處理組HT
線性混合模擬
該步驟是limma包的核心步驟，首先使用lmFit函數(shù)進行非線性最小二乘法分析燕侠，然后用經驗貝葉斯調整t-test中方差部分者祖，得到差異表達結果。
# #線性擬合模型構建 ---------------------------------------------------------------
fit <- lmFit(expr_data,design) #非線性最小二乘法
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)#用經驗貝葉斯調整t-test中方差的部分
DEG <- topTable(fit2, coef = 1,n = Inf,sort.by="logFC")
DEG <- na.omit(DEG)
最終生成的DEG文件包含以下幾列信息：
> colnames(DEG)
[1] "logFC" "AveExpr" "t" "P.Value" "adj.P.Val" "B"
[7] "regulate" "Genes"
差異基因標注
本步驟中對P.Value和logFC進行篩選绢彤，并對每個基因進行標注七问，顯示其表達變化情況。
DEG$regulate <- ifelse(DEG$P.Value > 0.05, "unchanged",
ifelse(DEG$logFC > 1, "up-regulated",
ifelse(DEG$logFC < -1, "down-regulated", "unchanged")))
結果保存
生成兩個文件茫舶，保存差異表達結果械巡。
write.table(table(DEG$regulate),file = paste0(job,"_","DEG_result_1_005.txt"),
sep = "\t",quote = F,row.names = T,col.names = T)
write.table(data.frame(gene_symbol=rownames(DEG),DEG),file = paste0(job,"_","DEG_result.txt"),
sep = "\t",quote = F,row.names = F,col.names = T)
區(qū)分上下調基因
該步驟中DEG$P.Value<0.05&abs(DEG$logFC)>1為參數(shù)，篩選了差異倍數(shù)和顯著性饶氏，并將基因按照上調和下調進行區(qū)分讥耗。
DE_1_0.05 <- DEG[DEG$P.Value<0.05&abs(DEG$logFC)>1,]
upGene_1_0.05 <- DE_1_0.05[DE_1_0.05$regulate == "up-regulated",]
downGene_1_0.05 <- DE_1_0.05[DE_1_0.05$regulate == "down-regulated",]
write.csv(upGene_1_0.05,paste0(job,"_","upGene_1_005.csv"))
write.csv(downGene_1_0.05,paste0(job,"_","downGene_1_005.csv"))
差異基因名稱提取
提取差異倍數(shù)最大的1000個基因，按照順序保存到txt文本中疹启，只保留基因名稱ID古程，用于后續(xù)驗證。
tem1 <- head(rownames(upGene_1_0.05),1000) #以logFC差異倍數(shù)從大到小為序喊崖，提取前1000個基因名稱
tem2 <- as.data.frame(tem1) # 轉化數(shù)據類型為數(shù)據框
write.table(tem2$tem,file = paste0(job,"_","upgene_head100_name.txt"),
row.names = FALSE,col.names = FALSE)
tem1 <- head(rownames(downGene_1_0.05),1000) #以logFC差異倍數(shù)從大到小為序挣磨，提取前1000個基因名稱
tem2 <- as.data.frame(tem1) # 轉化數(shù)據類型為數(shù)據框
write.table(tem2$tem,file = paste0(job,"_","downgene_head100_name.txt"),
row.names = FALSE,col.names = FALSE)
自定義篩選閾值
之前的結果均為默認設置菲宴，如果你需要修改，僅需更改下面開頭兩行參數(shù)即可趋急，運行后可以得到3個文件喝峦，分別是差異基因集、上下調過濾所得基因信息呜达。
foldChange = 2 # 自定義修改篩選參數(shù)
padj = 0.05 # 自定義修改篩選參數(shù)
All_diffSig <- DEG[(DEG$adj.P.Val < padj & (DEG$logFC > foldChange | DEG$logFC < (-foldChange))),]
#dim(All_diffSig)
write.csv(All_diffSig, paste0(job,"_","all_diffsig_filtered.csv")) ##輸出差異基因數(shù)據集
### 自定義篩選上調和下調的基因 ===================================================================
diffup <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC > foldChange)),]
write.csv(diffup, paste0(job,"_","diffup_filtered.csv"))
diffdown <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC < -foldChange)),]
write.csv(diffdown, paste0(job,"_","diffdown_filtered.csv"))
利用limma分析后谣蠢，正常情況下應該會生成下面這些數(shù)據文件，可以用于驗證過程中是否有問題查近。
image-20230222202113721
結果可視化
火山圖繪制方法一
有了上面的DEG差異分析結果眉踱，根據gene、logFC霜威、adj.P.Val三列信息即可繪制火山圖谈喳，第一種方法使用ggvolcano，繪圖代碼如下：
# ggvolcano繪制火山圖 ----------------------------------------------------------
pdf(paste0(job,"_","volcano1.pdf"),width = 10,height = 10)
ggvolcano(data = DEG,x = "logFC",y = "P.Value",output = FALSE,label = "Genes",
fills = c("#00AFBB", "#999999", "#FC4E07"),
colors = c("#00AFBB", "#999999", "#FC4E07"),
x_lab = "log2FC",
y_lab = "-Log10P.Value",
legend_position = "UR") #標簽位置為up right
dev.off()
image-20230222202352090
火山圖繪制方法二
第二種方法使用ggplot2戈泼，得到另外一種形式的火山圖婿禽，繪圖代碼如下：
# 火山圖的繪制 ------------------------------------------------------------------
DEG$Genes <- rownames(DEG)
pdf(paste0(job,"_","volcano2.pdf"),width = 7,height = 7)
ggplot(DEG,aes(x=logFC,y=-log10(P.Value)))+ #x軸logFC,y軸adj.p.value
geom_point(alpha=0.5,size=2,aes(color=regulate))+ #點的透明度，大小
ylab("-log10(P.Value)")+ #y軸的說明
scale_color_manual(values = c("blue", "grey", "red"))+ #點的顏色
geom_vline(xintercept = c(-1,1),lty=4,col ="black",lwd=0.8)+ #logFC分界線
geom_hline(yintercept=-log10(0.05),lty=4,col = "black",lwd=0.8)+ #adj.p.val分界線
theme_bw() #火山圖繪制
dev.off()
image-20230222202415207
熱圖繪制
根據基因的表達變化信息大猛，繪制熱圖并展示聚類樹扭倾，詳細代碼如下：
# 熱圖的繪制 -------------------------------------------------------------------
DEG_genes <- DEG[DEG$P.Value<0.05&abs(DEG$logFC)>1,]
DEG_gene_expr <- expr_data[rownames(DEG_genes),]
#DEG_gene_expr[is.infinite(DEG_gene_expr)] = 0
#DEG_gene_expr[DEG_gene_expr == -Inf] = 0
pdf(paste0(job,"_","pheatmap.pdf"))
pheatmap(DEG_gene_expr,
color = colorRampPalette(c("blue","white","red"))(100), #顏色
scale = "row", #歸一化的方式
border_color = NA, #線的顏色
fontsize = 10, #文字大小
show_rownames = F)
dev.off()
image-20230222202505089
參考資料：https://zhuanlan.zhihu.com/p/437712423
本文由mdnice多平臺發(fā)布
