獲取本章節(jié)數(shù)據(jù)和代碼：關(guān)注微信公眾號(hào)：小杜的生信筆記（ID：Du_Bioinformatics），回復(fù)關(guān)鍵詞：limma差異分析
------------------------
基于R語(yǔ)言進(jìn)行差異分析的包有很多個(gè)，比如我自己常用的有DESeq2、limma、edge等等。我們本期的專(zhuān)題是進(jìn)行各種包差異分析的專(zhuān)題。
本專(zhuān)題是使用limma包差異分析。
1.數(shù)據(jù)準(zhǔn)備
差異分析是兩兩數(shù)據(jù)集間的比較，一個(gè)是對(duì)照組樣本(CK)，一個(gè)是處理組樣本(Treat)。
CK01 CK02 CK03 T1 T2 T3
gene0001 56 46 91 36 19 28
gene0002 0 0 5 0 31 8
gene0003 4 0 6 2 2 0
gene0004 257 254 235 241 162 183
gene0005 30 17 22 27 9 75
gene0006 503 565 490 369 426 480
gene0007 201 82 62 296 206 189
gene0008 56 6 12 108 62 0
gene0009 6 10 9 13 8 14
gene0010 19 0 0 27 0 0
gene0011 0 0 0 0 0 0
<center>樣本格式
</center>
注：數(shù)據(jù)樣本根據(jù)自己實(shí)驗(yàn)數(shù)據(jù)決定
2.數(shù)據(jù)分析
2.1 導(dǎo)入所需的包
## 導(dǎo)入R包
library(limma)
library(dplyr)
2.2 加載數(shù)據(jù)
df <- read.table("Counts_data.txt", header = T, sep = "\t", row.names = 1, check.names = F)
head(df)
image.png
2.3 樣本信息注釋
list <- c(rep("Treat", 66), rep("CK",148)) %>% factor(., levels = c("CK", "Treat"), ordered = F)
##--------------
> head(list)
CK Treat
1 0 1
2 0 1
3 0 1
4 0 1
5 0 1
6 0 1
.............
list <- model.matrix(~factor(list)+0) #把group設(shè)置成一個(gè)model matrix
colnames(list) <- c("CK", "Treat")
df.fit <- lmFit(df, list) ## 數(shù)據(jù)與list進(jìn)行匹配
2.4 差異分析
df.matrix <- makeContrasts(tumor - normal, levels = list)
fit <- contrasts.fit(df.fit, df.matrix)
fit <- eBayes(fit)
tempOutput <- topTable(fit,n = Inf, adjust = "fdr")
> head(tempOutput)
logFC AveExpr t P.Value adj.P.Val B
gene11796 -1.2620803 5.551402 -8.392278 5.886136e-15 7.083964e-11 23.25729
gene7364 -0.9825962 7.471963 -7.598215 8.629186e-13 5.192613e-09 18.51600
gene11454 -1.3204341 6.238318 -7.368543 3.472935e-12 1.301549e-08 17.19354
gene11689 -1.0769861 4.967290 -7.331950 4.325880e-12 1.301549e-08 16.98503
gene11227 -1.2035217 5.925234 -7.263016 6.531456e-12 1.572122e-08 16.59390
gene3259 -2.3695741 9.467290 -6.962632 3.829470e-11 5.760959e-08 14.91576
到這部分，差異分析就結(jié)束了。對(duì)??！ 沒(méi)錯(cuò)就是這么簡(jiǎn)簡(jiǎn)單單...................
- 后面的部分就是提取差異結(jié)果，和繪制火山圖。
2.5 導(dǎo)出差異結(jié)果
## 導(dǎo)出所有的差異結(jié)果
nrDEG = na.omit(tempOutput) ## 去掉數(shù)據(jù)中有NA的行或列
diffsig <- nrDEG
write.csv(diffsig, "all.limmaOut.csv")
2.6 篩選差異基因
差異基因基因的篩選，我們一般是使用P值和LogFC篩選，常用的篩選標(biāo)準(zhǔn)P<0.05,|LogFC| > 1，這是最常規(guī)的篩選標(biāo)準(zhǔn)，如果你的數(shù)據(jù)差異較大，也可以更改P值和LogFC的大小。
## 我們使用|logFC| > 0.5，padj < 0.05（矯正后P值）
foldChange = 0.5
padj = 0.05
## 篩選出所有差異基因的結(jié)果
All_diffSig <- diffsig[(diffsig$adj.P.Val < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),]
#---------------------
> dim(All_diffSig)
[1] 0 6
## 我們發(fā)現(xiàn)竟然沒(méi)有差異基因，這是應(yīng)該我這邊的數(shù)據(jù)是隨機(jī)的結(jié)果，如果你的數(shù)據(jù)有這樣的問(wèn)題，你需要在仔細(xì)檢查一下哦。
## 我們?yōu)榱讼旅娴牟僮髡＿M(jìn)行，我們選用的P值（未矯正）進(jìn)行篩選。
All_diffSig <- diffsig[(diffsig$P.Value < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),]
write.csv(All_diffSig, "all.diffsig.csv") ##輸出差異基因數(shù)據(jù)集
#-----------------
> dim(All_diffSig)
[1] 335 6
## 共有335個(gè)差異基因
7）篩選上調(diào)和下調(diào)的基因
diffup <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig$logFC > foldChange)),]
write.csv(diffup, "diffup.csv")
#
diffdown <- All_diffSig[(All_diffSig$P.Value < padj & (All_diffSig < -foldChange)),]
write.csv(diffdown, "diffdown.csv")
到這部分，我們差異分析就全部結(jié)束了。已經(jīng)拿到差異文件，及上調(diào)和下調(diào)的基因文件。
3. 繪制火山圖
在常規(guī)的差異分析之后，我們需要進(jìn)行差異數(shù)據(jù)的可視化，火山圖和差異基因熱圖是最常用的可視化圖形。
## 導(dǎo)入R包
library(ggplot2)
library(ggrepel)
## 繪制火山圖
## 進(jìn)行分類(lèi)別
logFC <- diffsig$logFC
deg.padj <- diffsig$P.Value
data <- data.frame(logFC = logFC, padj = deg.padj)
data$group[(data$padj > 0.05 | data$padj == "NA") | (data$logFC < foldChange) & data$logFC > -foldChange] <- "Not"
data$group[(data$padj <= 0.05 & data$logFC > 1)] <- "Up"
data$group[(data$padj <= 0.05 & data$logFC < -1)] <- "Down"
x_lim <- max(logFC,-logFC)
# 開(kāi)始繪圖
pdf('volcano.pdf',width = 7,height = 6.5) ## 輸出文件
label = subset(diffsig,P.Value <0.05 & abs(logFC) > 0.5)
label1 = rownames(label)
colnames(diffsig)[1] = 'log2FC'
Significant=ifelse((diffsig$P.Value < 0.05 & abs(diffsig$log2FC)> 0.5), ifelse(diffsig$log2FC > 0.5,"Up","Down"), "Not")
ggplot(diffsig, aes(log2FC, -log10(P.Value)))+
geom_point(aes(col=Significant))+
scale_color_manual(values=c("#0072B5","grey","#BC3C28"))+
labs(title = " ")+
geom_vline(xintercept=c(-0.5,0.5), colour="black", linetype="dashed")+
geom_hline(yintercept = -log10(0.05),colour="black", linetype="dashed")+
theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))+
labs(x="log2(FoldChange)",y="-log10(Pvalue)")+
theme(axis.text=element_text(size=13),axis.title=element_text(size=13))+
str(diffsig, max.level = c(-1, 1))+theme_bw()
dev.off()
4. 繪制差異基因表達(dá)熱圖
注:本章節(jié)，我們是只使用count值進(jìn)行熱圖繪制，后續(xù)的分子中，我們會(huì)對(duì)其進(jìn)行優(yōu)化
## 導(dǎo)入R包
library(pheatmap)
## 提取差異基因的表達(dá)量
DEG_id <- read.csv("all.diffsig.csv", header = T) # 讀取差異基因的文件
head(DEG_id)
## 匹配差異基因的表達(dá)量
DEG_id <- unique(DEG_id$X)
DEG_exp <- df[DEG_id,]
hmexp <- na.omit(DEG_exp)
## 樣本注釋信息
annotation_col <- data.frame(Group = factor(c(rep("Treat", 66), rep("CK",148))))
rownames(annotation_col) <- colnames(hmexp)
## 繪制熱圖
pdf(file = "heatmap02.pdf", height = 8, width = 12)
pheatmap(hmexp,
annotation_col = annotation_col,
color = colorRampPalette(c("green","black","red"))(50),
cluster_cols = F,
show_rownames = F,
show_colnames = F,
scale = "row", ## none, row, column
fontsize = 12,
fontsize_row = 12,
fontsize_col = 6,
border = FALSE)
print(p)
dev.off()
image.png
獲取本章節(jié)數(shù)據(jù)和代碼：關(guān)注微信公眾號(hào)：小杜的生信筆記（ID：Du_Bioinformatics），回復(fù)關(guān)鍵詞：limma差異分析
“小杜的生信筆記” 公眾號(hào)、知乎、簡(jiǎn)書(shū)、B站平臺(tái)，主要發(fā)表或收錄生物信息學(xué)的教程，以及基于R的分析和可視化（包括數(shù)據(jù)分析，圖形繪制等）；分享感興趣的文獻(xiàn)和學(xué)習(xí)資料！
