欢迎来到尧图网

客户服务 关于我们

您的位置:首页 > 文旅 > 八卦 > 应审稿人要求| pseudo bulk差异分析

应审稿人要求| pseudo bulk差异分析

2025/4/19 20:10:01 来源:https://blog.csdn.net/weixin_47195452/article/details/143561351  浏览:    关键词:应审稿人要求| pseudo bulk差异分析

一、写在前面

最近有粉丝提问,收到了如下的审稿人意见:

图片

审稿人认为在单细胞测序过程中,利用findMarker通过Wilcox获得的差异基因虽然考虑到了不同组别细胞数量的不同,但是未能考虑到每组样本数量的不同。因此作者希望纳入样本水平的pseudo-bulk分析能够有助于确认两种条件下的差异基因。

首先我个人觉得审稿人自己的话有些矛盾,使用pseudo-bulk计算差异基因,岂不是无法考虑到不同组别中样本数量的差异?另外这有些吹毛求疵,在Seurat V5之前,作者甚至没有在包里集成pseudo-bulk的函数与算法(当然也可以自己提取矩阵计算)。难道能说作者发表的这几篇Cell和Nature给大家推荐的流程不好吗:


Hao, Hao, et al., Cell 2021 [Seurat v4]
Stuart, Butler, et al., Cell 2019 [Seurat v3]
Butler, et al., Nat Biotechnol 2018 [Seurat v2] Satija, Farrell, et al., Nat Biotechnol 2015 [Seurat v1]

再者说,scRNA-seqBulk RNA-Seq更加的稀疏,将前者模拟为后者参与差异计算,其实也没那么科学。当然,审稿人的观点也不是全无道理,若能够通过不同的算法得到相同的差异基因结果,的确有较高的说服力。

二、pseudo-bulk差异分析走起

测试文件可以自行下载:

图片

链接:https://pan.baidu.com/s/12dEGTJy4DnQ7gH2mbxCf-A?pwd=7qfm 

提取码:7qfm 

图片

2.1 数据载入

# 加载R包
library(Seurat)
## 载入需要的程序包:SeuratObject
## 载入需要的程序包:sp
## 
## 载入程序包:'SeuratObject'
## The following objects are masked from 'package:base':
## 
##     intersect, t
# 读取数据:
scRNA <- readRDS('test_data/T1D_scRNA.rds')# 这个数据包含24个样本:
unique(scRNA$sample)
##  [1] "D_503"  "H_120"  "H_630"  "H_3060" "D_609"  "H_727"  "H_4579" "D_504" 
##  [9] "H_3128" "H_7108" "D_502"  "D_497"  "D_506"  "H_409"  "H_6625" "D_610" 
## [17] "D_501"  "D_500"  "H_4119" "H_1334" "D_498"  "H_2928" "D_644"  "D_505"

# 包含两个组别的数据:
DimPlot(scRNA,split.by = 'Group')

图片

2.2 差异计算

(1) pseudo-bulk差异计算

### 生成拟bulk 数据 ###
bulk <- AggregateExpression(scRNA, return.seurat = T, slot = "counts", assays = "RNA",group.by = c("cell_type", "sample", "Group")# 分别填写细胞类型、样本变量、分组变量的slot名称)
## Names of identity class contain underscores ('_'), replacing with dashes ('-')
## Centering and scaling data matrix
## 
## This message is displayed once every 8 hours.

# 生成的是一个新的Seurat对象
bulk
## An object of class Seurat 
## 41056 features across 345 samples within 1 assay 
## Active assay: RNA (41056 features, 0 variable features)
##  3 layers present: counts, data, scale.data

我们可以像普通scRNA-seqSeurat对象一样,利用FindMarkers()进行差异分析,我们这里用celltype10做演示。

# 取出celltype10对应的对象:
ct10.bulk <- subset(bulk, cell_type == "celltype10")
# 改变默认分类变量:
Idents(ct10.bulk) <- "Group"# 下面的额计算依赖DESeq2,做过Bulk RNA-Seq的同学都知道:
if(!require(DESeq2))BiocManager::install('DESeq2')
## 载入需要的程序包:DESeq2
## 载入需要的程序包:S4Vectors
## 载入需要的程序包:stats4
## 载入需要的程序包:BiocGenerics
## 
## 载入程序包:'BiocGenerics'
## The following object is masked from 'package:SeuratObject':
## 
##     intersect
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
##     tapply, union, unique, unsplit, which.max, which.min
## 
## 载入程序包:'S4Vectors'
## The following object is masked from 'package:utils':
## 
##     findMatches
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 载入需要的程序包:IRanges
## 
## 载入程序包:'IRanges'
## The following object is masked from 'package:sp':
## 
##     %over%
## The following object is masked from 'package:grDevices':
## 
##     windows
## 载入需要的程序包:GenomicRanges
## 载入需要的程序包:GenomeInfoDb
## 载入需要的程序包:SummarizedExperiment
## 载入需要的程序包:MatrixGenerics
## 载入需要的程序包:matrixStats
## 
## 载入程序包:'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 载入需要的程序包:Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 载入程序包:'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## 载入程序包:'SummarizedExperiment'
## The following object is masked from 'package:Seurat':
## 
##     Assays
## The following object is masked from 'package:SeuratObject':
## 
##     Assays

# 差异计算:
bulk_deg <- FindMarkers(ct10.bulk, ident.1 = "D", ident.2 = "H", # 这样算出来的Fold Change就是D/Hslot = "counts", test.use = "DESeq2",# 这里可以选择其它算法verbose = F# 关闭进度提示)
## converting counts to integer mode
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates

head(bulk_deg)# 看一下差异列表
##                        p_val avg_log2FC pct.1 pct.2 p_val_adj
## ENSG00000047346 1.900955e-05 -1.1201537 0.917 1.000  0.780456
## ENSG00000168685 7.606182e-05 -0.5817112 1.000 1.000  1.000000
## ENSG00000131759 9.696591e-05 -2.1668759 0.667 1.000  1.000000
## ENSG00000166750 2.136829e-04 -1.5473545 0.750 0.917  1.000000
## ENSG00000163947 3.127113e-04 -0.9136378 0.917 1.000  1.000000
## ENSG00000239713 4.090397e-04 -2.2337008 0.250 0.917  1.000000

如何写循环计算所有细胞类型的差异基因,就留在这里当习题啦。

(2)细胞水平的差异计算

# 整理分组变量:
scRNA$CT_Group <- paste(scRNA$cell_type,scRNA$Group,sep = '_')# 查看新的分组变量:
unique(scRNA$CT_Group)
##  [1] "celltype12_D" "celltype3_H"  "celltype13_H" "celltype6_H"  "celltype0_D" 
##  [6] "celltype12_H" "celltype4_H"  "celltype11_D" "celltype14_H" "celltype9_H" 
## [11] "celltype11_H" "celltype2_H"  "celltype0_H"  "celltype7_H"  "celltype14_D"
## [16] "celltype1_D"  "celltype4_D"  "celltype1_H"  "celltype8_H"  "celltype3_D" 
## [21] "celltype13_D" "celltype8_D"  "celltype7_D"  "celltype5_H"  "celltype6_D" 
## [26] "celltype15_H" "celltype2_D"  "celltype5_D"  "celltype10_H" "celltype9_D" 
## [31] "celltype10_D" "celltype15_D"

# 差异计算:
cell_deg <- FindMarkers(scRNA,ident.1 = 'celltype10_D',ident.2 = 'celltype10_H' ,group.by = 'CT_Group')# 同样得到的是celltype10在D组 vs H组的结果
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session

(3)两种算法的对比

library(dplyr)
## 
## 载入程序包:'dplyr'
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following object is masked from 'package:matrixStats':
## 
##     count
## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect
## The following objects are masked from 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# 先看一下两种算法的显著差异基因数量  
bulk_sig <- filter(bulk_deg,p_val < 0.05)
nrow(bulk_sig)
## [1] 308

cell_sig <- filter(cell_deg,p_val < 0.05)
nrow(cell_sig)
## [1] 1494

可以看出,pseudo-bulk得到的差异基因数量要少很多,画一个韦恩图看看二者交集

if(!require(VennDiagram))install.packages("VennDiagram")
## 载入需要的程序包:VennDiagram
## 载入需要的程序包:grid
## 载入需要的程序包:futile.logger

venn.plot <- venn.diagram(x = list(Bulk = rownames(bulk_sig), Cell = rownames(cell_sig)
),category.names = c("Bulk DEG", "Single-Cell DEG"),filename = NULL,output = TRUE,main = "Venn Diagram of Significant Genes"
)
grid.draw(venn.plot)

图片

可以看出包含关系还是挺明显的,那我们再用交集基因的avg_log2FC做一个线性回归看看两次差异分析的相关性如何:

# 获得两次差异分析共同出现的基因:
inter_gene <- intersect(rownames(bulk_sig),rownames(cell_sig))
# 取出avg_log2FC整理为数据框
data4plot <- data.frame(Bulk = bulk_sig[inter_gene,'avg_log2FC'],Cell = cell_sig[inter_gene,'avg_log2FC']  )# 线性回归分析:
lm.model <- lm(Bulk ~ Cell,data = data4plot)
summary(lm.model)#看一下统计学参数
## 
## Call:
## lm(formula = Bulk ~ Cell, data = data4plot)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01613 -0.24964 -0.04723  0.17148  2.19351 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.17249    0.02757  -6.257 1.58e-09 ***
## Cell         0.70081    0.01959  35.767  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4107 on 263 degrees of freedom
## Multiple R-squared:  0.8295, Adjusted R-squared:  0.8288 
## F-statistic:  1279 on 1 and 263 DF,  p-value: < 2.2e-16

mypara <- coefficients(lm.model)#得到截距和斜率
a <- mypara[2]#斜率
b <- mypara[1]#截距
a <- round(a,2)#取两位有效数字
b <- round(b,2)library(ggplot2)
library(ggpubr)
## 
## 载入程序包:'ggpubr'
## The following object is masked from 'package:VennDiagram':
## 
##     rotate
# 来个散点图吧~
lmplot <- ggplot( data4plot,aes(x=Bulk, y=Cell))+geom_point(color="black")+stat_smooth(method="lm",se=TRUE)+stat_cor(data=data4plot, method = "pearson")+#加上置信区间、R值、P值ggtitle(label = paste(": y = ", a, " * x + ", b, sep = ""))+geom_rug()+#加上线性回归方程labs(x='Bulk DEG',y= 'single-cell DEG')lmplot
## `geom_smooth()` using formula = 'y ~ x'

图片

R=0.91,那么R^2就是0.83,可以看出二者的相关性还是不错的,就看能不能过审稿人这关啦。

环境信息

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.utf8 
## [2] LC_CTYPE=Chinese (Simplified)_China.utf8   
## [3] LC_MONETARY=Chinese (Simplified)_China.utf8
## [4] LC_NUMERIC=C                               
## [5] LC_TIME=Chinese (Simplified)_China.utf8    
## 
## time zone: Asia/Shanghai
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggpubr_0.6.0                ggplot2_3.5.1              
##  [3] VennDiagram_1.7.3           futile.logger_1.4.3        
##  [5] dplyr_1.1.4                 DESeq2_1.44.0              
##  [7] SummarizedExperiment_1.34.0 Biobase_2.64.0             
##  [9] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## [11] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [13] IRanges_2.38.1              S4Vectors_0.42.1           
## [15] BiocGenerics_0.50.0         Seurat_5.1.0               
## [17] SeuratObject_5.0.2          sp_2.1-4                   
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22        splines_4.4.1           later_1.3.2            
##   [4] tibble_3.2.1            polyclip_1.10-7         fastDummies_1.7.4      
##   [7] lifecycle_1.0.4         rstatix_0.7.2           globals_0.16.3         
##  [10] lattice_0.22-6          MASS_7.3-60.2           backports_1.5.0        
##  [13] magrittr_2.0.3          plotly_4.10.4           sass_0.4.9             
##  [16] rmarkdown_2.28          jquerylib_0.1.4         yaml_2.3.10            
##  [19] httpuv_1.6.15           sctransform_0.4.1       spam_2.10-0            
##  [22] spatstat.sparse_3.1-0   reticulate_1.39.0       cowplot_1.1.3          
##  [25] pbapply_1.7-2           RColorBrewer_1.1-3      abind_1.4-5            
##  [28] zlibbioc_1.50.0         Rtsne_0.17              purrr_1.0.2            
##  [31] GenomeInfoDbData_1.2.12 ggrepel_0.9.6           irlba_2.3.5.1          
##  [34] listenv_0.9.1           spatstat.utils_3.1-0    openintro_2.5.0        
##  [37] airports_0.1.0          goftest_1.2-3           RSpectra_0.16-2        
##  [40] spatstat.random_3.3-1   fitdistrplus_1.2-1      parallelly_1.38.0      
##  [43] leiden_0.4.3.1          codetools_0.2-20        DelayedArray_0.30.1    
##  [46] tidyselect_1.2.1        UCSC.utils_1.0.0        farver_2.1.2           
##  [49] spatstat.explore_3.3-2  jsonlite_1.8.8          progressr_0.14.0       
##  [52] ggridges_0.5.6          survival_3.6-4          tools_4.4.1            
##  [55] ica_1.0-3               Rcpp_1.0.13             glue_1.7.0             
##  [58] gridExtra_2.3           SparseArray_1.4.8       mgcv_1.9-1             
##  [61] xfun_0.47               withr_3.0.1             formatR_1.14           
##  [64] fastmap_1.2.0           fansi_1.0.6             digest_0.6.37          
##  [67] R6_2.5.1                mime_0.12               colorspace_2.1-1       
##  [70] scattermore_1.2         tensor_1.5              spatstat.data_3.1-2    
##  [73] utf8_1.2.4              tidyr_1.3.1             generics_0.1.3         
##  [76] data.table_1.16.0       usdata_0.3.1            httr_1.4.7             
##  [79] htmlwidgets_1.6.4       S4Arrays_1.4.1          uwot_0.2.2             
##  [82] pkgconfig_2.0.3         gtable_0.3.5            lmtest_0.9-40          
##  [85] XVector_0.44.0          htmltools_0.5.8.1       carData_3.0-5          
##  [88] dotCall64_1.1-1         scales_1.3.0            png_0.1-8              
##  [91] spatstat.univar_3.0-1   knitr_1.48              lambda.r_1.2.4         
##  [94] rstudioapi_0.16.0       tzdb_0.4.0              reshape2_1.4.4         
##  [97] nlme_3.1-164            cachem_1.1.0            zoo_1.8-12             
## [100] stringr_1.5.1           KernSmooth_2.23-24      parallel_4.4.1         
## [103] miniUI_0.1.1.1          pillar_1.9.0            vctrs_0.6.5            
## [106] RANN_2.6.2              promises_1.3.0          car_3.1-2              
## [109] xtable_1.8-4            cluster_2.1.6           evaluate_0.24.0        
## [112] readr_2.1.5             cli_3.6.3               locfit_1.5-9.10        
## [115] compiler_4.4.1          futile.options_1.0.1    rlang_1.1.4            
## [118] crayon_1.5.3            future.apply_1.11.2     ggsignif_0.6.4         
## [121] labeling_0.4.3          plyr_1.8.9              stringi_1.8.4          
## [124] viridisLite_0.4.2       deldir_2.0-4            BiocParallel_1.38.0    
## [127] munsell_0.5.1           lazyeval_0.2.2          spatstat.geom_3.3-2    
## [130] Matrix_1.7-0            RcppHNSW_0.6.0          hms_1.1.3              
## [133] patchwork_1.2.0         future_1.34.0           shiny_1.9.1            
## [136] highr_0.11              ROCR_1.0-11             broom_1.0.6            
## [139] igraph_2.0.3            bslib_0.8.0             cherryblossom_0.1.0

欢迎致谢

如果以上内容对你有帮助,欢迎在文章的Acknowledgement中加上这一段,联系客服微信可以发放奖励:

Since Biomamba and his wechat public account team produce bioinformatics tutorials and share code with annotation, we thank Biomamba for their guidance in bioinformatics and data analysis for the current study. 

欢迎在发文/毕业时向我们分享你的喜悦~

版权声明:

本网仅为发布的内容提供存储空间,不对发表、转载的内容提供任何形式的保证。凡本网注明“来源:XXX网络”的作品,均转载自其它媒体,著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处。

我们尊重并感谢每一位作者,均已注明文章来源和作者。如因作品内容、版权或其它问题,请及时与我们联系,联系邮箱:809451989@qq.com,投稿邮箱:809451989@qq.com

热搜词