获取本章节数据和代码：关注微信公众号：小杜的生信筆記（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，这是最常规的筛选标准，如果你的数据差异较大，也可以更改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
## 我们发现竟然没有差异基因，这是应该我这边的数据是随机的结果，如果你的数据有这样的问题，你需要在仔细检查一下哦。
## 我们为了下面的操作正常进行，我们选用的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. 绘制火山图
在常规的差异分析之后，我们需要进行差异数据的可视化，火山图和差异基因热图是最常用的可视化图形。
## 导入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. 绘制差异基因表达热图
注:本章节，我们是只使用count值进行热图绘制，后续的分子中，我们会对其进行优化
## 导入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
获取本章节数据和代码：关注微信公众号：小杜的生信筆記（ID：Du_Bioinformatics），回复关键词：limma差异分析
“小杜的生信筆記” 公众号、知乎、简书、B站平台，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！
网友评论