獲取本章節(jié)數據和代碼：關注微信公眾號：小杜的生信筆記（ID：Du_Bioinformatics），回復關鍵詞：limma差異分析
------------------------
基于R語言進行差異分析的包有很多個毯辅，比如我自己常用的有DESeq2、limma嚼贡、edge等等弥姻。我們本期的專題是進行各種包差異分析的專題。
本專題是使用limma包差異分析慢叨。
1.數據準備
差異分析是兩兩數據集間的比較盖矫，一個是對照組樣本(CK)丽惭，一個是處理組樣本(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>
注：數據樣本根據自己實驗數據決定
2.數據分析
2.1 導入所需的包
## 導入R包
library(limma)
library(dplyr)
2.2 加載數據
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設置成一個model matrix
colnames(list) <- c("CK", "Treat")
df.fit <- lmFit(df, list) ## 數據與list進行匹配
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
到這部分炼彪，差異分析就結束了吐根。對Ｕ摺辐马！ 沒錯就是這么簡簡單單...................
- 后面的部分就是提取差異結果，和繪制火山圖。
2.5 導出差異結果
## 導出所有的差異結果
nrDEG = na.omit(tempOutput) ## 去掉數據中有NA的行或列
diffsig <- nrDEG
write.csv(diffsig, "all.limmaOut.csv")
2.6 篩選差異基因
差異基因基因的篩選喜爷，我們一般是使用P值和LogFC篩選冗疮，常用的篩選標準P<0.05,|LogFC| > 1，這是最常規(guī)的篩選標準檩帐，如果你的數據差異較大术幔，也可以更改P值和LogFC的大小。
## 我們使用|logFC| > 0.5湃密，padj < 0.05（矯正后P值）
foldChange = 0.5
padj = 0.05
## 篩選出所有差異基因的結果
All_diffSig <- diffsig[(diffsig$adj.P.Val < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),]
#---------------------
> dim(All_diffSig)
[1] 0 6
## 我們發(fā)現竟然沒有差異基因诅挑，這是應該我這邊的數據是隨機的結果，如果你的數據有這樣的問題泛源，你需要在仔細檢查一下哦拔妥。
## 我們?yōu)榱讼旅娴牟僮髡＿M行，我們選用的P值（未矯正）進行篩選达箍。
All_diffSig <- diffsig[(diffsig$P.Value < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),]
write.csv(All_diffSig, "all.diffsig.csv") ##輸出差異基因數據集
#-----------------
> dim(All_diffSig)
[1] 335 6
## 共有335個差異基因
7）篩選上調和下調的基因
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")
到這部分没龙，我們差異分析就全部結束了。已經拿到差異文件缎玫，及上調和下調的基因文件硬纤。
3. 繪制火山圖
在常規(guī)的差異分析之后，我們需要進行差異數據的可視化赃磨，火山圖和差異基因熱圖是最常用的可視化圖形筝家。
## 導入R包
library(ggplot2)
library(ggrepel)
## 繪制火山圖
## 進行分類別
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)
# 開始繪圖
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. 繪制差異基因表達熱圖
注:本章節(jié)，我們是只使用count值進行熱圖繪制邻辉，后續(xù)的分子中肛鹏，我們會對其進行優(yōu)化
## 導入R包
library(pheatmap)
## 提取差異基因的表達量
DEG_id <- read.csv("all.diffsig.csv", header = T) # 讀取差異基因的文件
head(DEG_id)
## 匹配差異基因的表達量
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é)數據和代碼：關注微信公眾號：小杜的生信筆記（ID：Du_Bioinformatics），回復關鍵詞：limma差異分析
“小杜的生信筆記” 公眾號恩沛、知乎在扰、簡書、B站平臺雷客，主要發(fā)表或收錄生物信息學的教程芒珠，以及基于R的分析和可視化（包括數據分析，圖形繪制等）搅裙；分享感興趣的文獻和學習資料皱卓！
