如何进行批量差异分析并绘制其火山图及拼图
2023-02-19 12:23:54 时间
缘起
上周，曾老师给了我一个8个样本8个组别的转录组数据，即每组只有一个样本的转录组数据。我的处理方式是从中抽取两组进行差异分析，与文中描述的显著差异基因数目以及文中指出的差异基因进行比较，看看分析结果是否大致相同。但是，其实我有些没有理解到老师的意思。老师的初衷是想同原文一样批量绘制「同个部位」两两组别间的差异分析结果，看看其差异基因数量的分布，然后进行比较。在看到我理解偏差后，老师还飞快地给我提供了单样本批量差异分析的脚本。「因此，本周我们主要重点展示如何批量进行单样本差异分析以及批量绘制火山图并拼图」。
探究
同样地，今天，我们使用的转录组数据集还是2019年发表在Diabetes杂志上，文献名称为Sarm1 Gene Defificiency Attenuates Diabetic Peripheral Neuropathy in Mice。该数据集由8个样本组成，每个样本代表一个分组。
转录组数据集介绍
该数据集提交在ENA官网，其PRJ项目号是PRJNA540413。若有小伙伴想从上游开始处理，具体数据下载链接见 https://www.ebi.ac.uk/ena/browser/view/PRJNA540413。该数据集有8个组别，但是8个组别中分为两个部位，SC与SN，处理中含基因敲除鼠与WT鼠，以及STZ处理与药物处理。
重点强调，SC与SN代表「两个部位」哈，文中展示的图也是按照两个部位分为上半部分（SC部位4组4副图）与下半部分（SN组4组4副图）
正式分析
1.清空环境，加载R包
rm(list = ls())
library(edgeR)
library(DESeq2)
library(FactoMineR)
library(factoextra)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(stringr)
library(stringi)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(VennDiagram)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(EnhancedVolcano)
2.读取数据，获得基因表达矩阵
rawcount <- read.table("./all_trim.id.txt",header = T,sep="\t",row.names = 1)
colnames(rawcount)
## [1] "Chr" "Start" "End"
## [4] "Strand" "Length" "SRR8988890.sort.bam"
## [7] "SRR8988891.sort.bam" "SRR8988892.sort.bam" "SRR8988893.sort.bam"
## [10] "SRR8988894.sort.bam" "SRR8988895.sort.bam" "SRR8988896.sort.bam"
## [13] "SRR8988897.sort.bam"
rawcount[1:4,1:4]
## Chr
## ENSMUSG00000102693.1 chr1
## ENSMUSG00000064842.1 chr1
## ENSMUSG00000051951.5 chr1;chr1;chr1;chr1;chr1;chr1;chr1
## ENSMUSG00000102851.1 chr1
## Start
## ENSMUSG00000102693.1 3073253
## ENSMUSG00000064842.1 3102016
## ENSMUSG00000051951.5 3205901;3206523;3213439;3213609;3214482;3421702;3670552
## ENSMUSG00000102851.1 3252757
## End
## ENSMUSG00000102693.1 3074322
## ENSMUSG00000064842.1 3102125
## ENSMUSG00000051951.5 3207317;3207317;3215632;3216344;3216968;3421901;3671498
## ENSMUSG00000102851.1 3253236
## Strand
## ENSMUSG00000102693.1 +
## ENSMUSG00000064842.1 +
## ENSMUSG00000051951.5 -;-;-;-;-;-;-
## ENSMUSG00000102851.1 +
rawcount=rawcount[,6:13]
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
filter_count <- rawcount[keep,] #获得filter_count矩阵
express_cpm <- log2(cpm(filter_count)+ 1)
express_cpm[1:4,1:4] #获得cpm矩阵
## SRR8988890.sort.bam SRR8988891.sort.bam
## ENSMUSG00000051951.5 1.5845466 1.6597803
## ENSMUSG00000102331.1 0.5182788 0.6228564
## ENSMUSG00000025900.13 0.2167244 0.2165807
## ENSMUSG00000025902.13 3.7177481 3.7461175
## SRR8988892.sort.bam SRR8988893.sort.bam
## ENSMUSG00000051951.5 2.0250141 1.7654425
## ENSMUSG00000102331.1 0.4685779 0.8318200
## ENSMUSG00000025900.13 0.2196636 0.5655548
## ENSMUSG00000025902.13 3.9206494 4.4710626
a=str_split(rownames(filter_count),"\\.",simplify = T)[,1]
rownames(filter_count)=a
rawcount=filter_count
##转换ID
id2symbol <- bitr(rownames(rawcount),
fromType = "ENSEMBL",
toType = "SYMBOL",
OrgDb = org.Mm.eg.db)
rawcount <- cbind(GeneID=rownames(rawcount),rawcount) #多一列换成基因ID
rawcount <- merge(id2symbol,rawcount,
by.x="ENSEMBL",by.y="GeneID",all.y=T)
rawcount=rawcount[!duplicated(rawcount$SYMBOL),]
rawcount <- drop_na(rawcount) #去除含有NA值的行
rownames(rawcount)=rawcount$SYMBOL
rawcount=rawcount[,3:10]
head(rawcount)
## SRR8988890.sort.bam SRR8988891.sort.bam SRR8988892.sort.bam
## Gnai3 732 601 592
## Cdc45 75 79 47
## H19 28 59 15
## Scml2 6 12 10
## Apoh 16 19 15
## Narf 996 801 857
## SRR8988893.sort.bam SRR8988894.sort.bam SRR8988895.sort.bam
## Gnai3 647 2715 1778
## Cdc45 56 57 50
## H19 41 1700 2338
## Scml2 10 35 26
## Apoh 10 30 41
## Narf 916 598 933
## SRR8988896.sort.bam SRR8988897.sort.bam
## Gnai3 2214 2644
## Cdc45 104 63
## H19 1749 1273
## Scml2 35 31
## Apoh 7 17
## Narf 927 761
3.自定义针对单样本的差异分析函数
（统计单样本之间差异分析上下调显著基因的数目）
deg_single <- function(exprSet){
#加载包
# exprSet= rawcount[,c(1,3)]
exprSet=exprSet[rowSums(exprSet)>1,]
group_list = c('A','B')
library(edgeR)
#设置分组信息
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.1#设置bcv为0.1
et <- exactTest(exprSet, dispersion=bcv^2)
DEG_edgeR=as.data.frame(topTags(et, n = nrow(exprSet$counts)))
head(DEG_edgeR)
#筛选上下调，设定阈值
fc_cutoff <- 2
pvalue <- 0.01
DEG_edgeR$regulated <- "normal"
loc_up <- intersect(which(DEG_edgeR$logFC>log2(fc_cutoff)),
which(DEG_edgeR$PValue<pvalue))
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
which(DEG_edgeR$PValue<pvalue))
DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"
return(table(DEG_edgeR$regulated))
}
4. 1-4号（SC）样本差异分析后，显著上下调基因的统计
注意哈，前4个样本是针对SC的，也就是文章图的上半部分。区别于文章图有两点，一是可视化的方式，二是文章展示的是4个组别间表达量取lg值绘制散点图，然而我们用火山图展示了组别间两两差异分析的12个差异分析结果（除了未进行自身之间的差异分析外，我们都进行比较了哈）
B1=lapply(1:4, function(i){
lapply(1:4, function(j){
if(i != j){
deg_single(rawcount[,c(i,j)])
}
})
})
B1
## [[1]]
## [[1]][[1]]
## NULL
##
## [[1]][[2]]
##
## down normal up
## 110 17454 217
##
## [[1]][[3]]
##
## down normal up
## 173 17457 92
##
## [[1]][[4]]
##
## down normal up
## 62 17382 306
##
##
## [[2]]
## [[2]][[1]]
##
## down normal up
## 217 17454 110
##
## [[2]][[2]]
## NULL
##
## [[2]][[3]]
##
## down normal up
## 260 17431 39
##
## [[2]][[4]]
##
## down normal up
## 191 17237 361
##
##
## [[3]]
## [[3]][[1]]
##
## down normal up
## 92 17457 173
##
## [[3]][[2]]
##
## down normal up
## 39 17431 260
##
## [[3]][[3]]
## NULL
##
## [[3]][[4]]
##
## down normal up
## 39 17236 453
##
##
## [[4]]
## [[4]][[1]]
##
## down normal up
## 306 17382 62
##
## [[4]][[2]]
##
## down normal up
## 361 17237 191
##
## [[4]][[3]]
##
## down normal up
## 453 17236 39
##
## [[4]][[4]]
## NULL
5.自定义差异分析结果绘制火山图的函数
deg_EnhanceV <- function(exprSet){
#加载包
# exprSet= rawcount[,c(1,3)]
exprSet=exprSet[rowSums(exprSet)>1,]
group_list = c('A','B')
library(edgeR)
#设置分组信息
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.1#设置bcv为0.1
et <- exactTest(exprSet, dispersion=bcv^2)
DEG_edgeR=as.data.frame(topTags(et, n = nrow(exprSet$counts)))
head(DEG_edgeR)
write.csv(DEG_edgeR, 'single-sample-deg-result.csv', quote = FALSE)
#筛选上下调，设定阈值
fc_cutoff <- 2
pvalue <- 0.01
DEG_edgeR$regulated <- "normal"
loc_up <- intersect(which(DEG_edgeR$logFC>log2(fc_cutoff)),
which(DEG_edgeR$PValue<pvalue))
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
which(DEG_edgeR$PValue<pvalue))
DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"
EnhanceV=list()
EnhanceV=print(as.ggplot(EnhancedVolcano(DEG_edgeR,
lab =rownames(DEG_edgeR),
x = 'logFC',
y = 'FDR')))
# b=list(EnhanceV,table(DEG_edgeR$regulated))
# return(b)
return(EnhanceV)
}
6. 1-4号（SC）样本批量绘制火山图并且拼图
C1=lapply(1:4, function(i){
lapply(1:4, function(j){
if(i != j){
C= deg_EnhanceV(rawcount[,c(i,j)])
}
})
})
# 注意批量运行完，会展示12副图的结果，下面显展示第一幅图的差异分析结果，12副图的结果以拼图的形式展现
C=C1
library(cowplot)
b=plot_grid(plotlist = c(C[[1]],C[[2]],C[[3]],C[[4]]), align = "h",
nrow = 4, labels=LETTERS[1:16])
b
7. 5-8号（SN）样本差异分析后，显著上下调基因统计
B5=lapply(5:8, function(i){
lapply(5:8, function(j){
if(i != j){
deg_single(rawcount[,c(i,j)])
}
})
})
B5
## [[1]]
## [[1]][[1]]
## NULL
##
## [[1]][[2]]
##
## down normal up
## 945 16161 608
##
## [[1]][[3]]
##
## down normal up
## 211 17041 492
##
## [[1]][[4]]
##
## down normal up
## 440 16834 451
##
##
## [[2]]
## [[2]][[1]]
##
## down normal up
## 608 16161 945
##
## [[2]][[2]]
## NULL
##
## [[2]][[3]]
##
## down normal up
## 242 16966 499
##
## [[2]][[4]]
##
## down normal up
## 584 16107 1025
##
##
## [[3]]
## [[3]][[1]]
##
## down normal up
## 492 17041 211
##
## [[3]][[2]]
##
## down normal up
## 499 16966 242
##
## [[3]][[3]]
## NULL
##
## [[3]][[4]]
##
## down normal up
## 500 16947 293
##
##
## [[4]]
## [[4]][[1]]
##
## down normal up
## 451 16834 440
##
## [[4]][[2]]
##
## down normal up
## 1025 16107 584
##
## [[4]][[3]]
##
## down normal up
## 293 16947 500
##
## [[4]][[4]]
## NULL
8. 5-8号（SN）样本的批量绘制火山图并且拼图
C5=lapply(5:8, function(i){
lapply(5:8, function(j){
if(i != j){
C= deg_EnhanceV(rawcount[,c(i,j)])
}
})
})
# 注意批量运行完，会展示12副图的结果，下面显展示第一幅图的差异分析结果，12副图的结果同样在下面的代码中以拼图的形式展现
C=C5
library(cowplot)
b=plot_grid(plotlist = c(C[[1]],C[[2]],C[[3]],C[[4]]), align = "h",
nrow = 4, labels=LETTERS[5:21]) #此处标记咱们就从第五个字母E开始吧b以上演示了批量差异分析以及批量绘制火山图及拼图的方式，感兴趣的小伙伴们可以尝试一下。12副火山图的拼图其实有点大，如果不清楚的话，自身探究时可以将图片保存大些即可。
值得注意的是：原文是对两组之间的lg值，绘制散点图；而不是像我们一样两两组合进行差异分析；我们在获得了差异分析的结果之后，如果有余力的话，其实也可以向作者一样进行两两组合绘制散点图探索下，感兴趣的小伙伴们可以自身尝试下哈。
相关文章
- Jgit的使用笔记
- 利用Github Action实现Tornadofx/JavaFx打包
- 叹息！GitHub Trending 即将成为历史！
- 微软软了？开源社区讨论炸锅，GitHub CEO 亲自来答
- GitHub Trending 列表频现重复项，前后端都没去重？
- Photoshop Elements 2021版本软件安装教程（mac+windows全版本都有）
- (ps全版本)Photoshop 2020的安装与破解教程（mac+windows全版本都有）
- (ps全版本)Photoshop cc2018的安装与破解教程（mac+windows全版本，包括2023
- 环境搭建：Oracle GoldenGate 大数据迁移到 Redshift/Flat file/Flume/Kafka测试流程
- 每个开发人员都要掌握的：最小 Linux 基础课
- 来撸羊毛了！Windows 环境下 Hexo 博客搭建，并部署到 GitHub Pages
- 超实用！手把手入门 MongoDB：这些坑点请一定远离
- 【GitHub日报】22-10-09 zustand、neovim、webtorrent、express 等4款App今日上新
- 【GitHub日报】22-10-10 brew、minio、vite、seaweedfs、dbeaver 等8款App今日上新
- 【GitHub日报】22-10-11 cobra、grafana、vue、ToolJet、redwood 等13款App今日上新
- Photoshop 2018 下载及安装教程（mac+windows全版本都有，包括最新的2023）
- Photoshop 2017 下载及安装教程（mac+windows全版本都有，包括最新的2023）
- Photoshop 2020 下载及安装教程（mac+windows全版本都有，包括最新的2023）
- Photoshop 2023 资源免费下载（mac+windows全版本都有，包括最新的2023）
- 最新版本Photoshop CC2018软件安装教程（mac+windows全版本都有，包括2023