ixxmu
commented
1 year ago Rscript脚本解放双手：以CCA整合分析为例 by 生信随笔
随着我深入学习生信分析，我发现生信真的是学无止境：从R语言编程的学习，到常规统计和画图、芯片数据，NGS多组学上下游包括转录组、基因组、表观遗传和蛋白组学等，再到单细胞多组学和空间转录组，还有我从未涉猎的宏基因组和影像组学等等。
从编程学习的角度来说，除了需要熟练掌握R语言以外，我发现还需要学习Linux操作和shell命令，如果需要分析单细胞/空转大数据的话，还需要学习Python。 从单细胞多组学和空转学习的角度来说，除了常规的标准分析性，我们还需要学习各种个性化分析（本质上是学习各种R包和Python包，包括拟时序、转录因子分析和inferCNV等等），这些分析的难易基本取决于前期R语言、Linux和Python的掌握情况，编程能力越强，这些个性化分析代码层面的理解越容易。特别是空转的学习，运行资源消耗大，需要掌握的Python包特别多。 从知识架构习的角度来说，除了上述代码编程层面的知识体系建立，还需要建立生物信息学领域的常规知识体系，例如测序，多组学，Bulk和单细胞测序等这些基本概念是什么及其原理。最重要的是第三层次的知识体系的学习，围绕本领域的背景知识和本领域关注的生物学问题。我深有体会的是，如果只掌握了前述2个层次的知识体系的话，对于我们做科学研究，仍然是不够的，可以说没有第三层次的内核作为支撑，大概率会被吐槽为“程序员”。此外，如果是从事干湿结合的话，有必要了解和学习各种湿实验及其大致原理和用途（可以回答哪些层面的科学问题？）。
坦率的说，我从2020年8月到现在2023年8月学习了3年的生信分析，我感觉我才刚刚入门生信这个领域。。。越学越发现自己的渺小。但是我也鼓励自己，同时我也想鼓励各位：生信的学习本身是一个学习曲线比较陡的学科，基础入门难，但是咬咬牙坚持下去，后续的问题基本都是举一反三的。就比如我现在在学习Python单细胞分析，R语言单细胞分析其实是我的一个舒适区，我是极不想去学习Python的单细胞分析体系，但我又深知Python的实用性和价值。所以这一年多，我几次学习Python，几次悻悻而归。最近因为分析需要，咬一咬牙坚持，突然发现Python也没有像20年那会学习R语言、生信和Linux的难（从0学起的难度，现在回想起来还是很痛苦）。总体来说，因为有前期的生信基础，Python的学习曲线的陡坡会比想象中的低很多。
最后，我再分享最近的一个感悟（鸡汤）：生物信息学的学习是永无止境的，对于未来从事生信这行也包括从事湿实验这行的人来说（现在的湿实验基本离不开干实验分析），很可能是一辈子的学习，所以不要觉得目前自己能作出点啥东西就觉得满足了；同时，面对自己的无知，也不要恐慌和自责，毕竟闻道有先后，术业有专攻，哪里不会从哪里学起，一定要相信自己的学习能力。
言归正传，以下是正文。
本期帖子分享一个Linux服务器下Rscript传参脚本，用于单细胞批次的整合分析。其实熟悉我的朋友都知道，关于Linux服务器下的脚本，我已做了多次的分享，例如：
举一反三，Shell脚本可以挂后台运行，R语言和Python的当然也是可以：
在熟练掌握这些批处理脚本之后，有一个疑问随之而来：如何像Linux下运行那些软件一样，进行参数的传递？这里我以单细胞整合分析的CCA为例写了一个传参R语言脚本：
cat >Rscript_seurat_CCA.R
首先在Linux命令行下运行上述代码，即创建一个空白的Rscript_seurat_CCA.R脚本，然后复制粘贴下述代码：
##### 参数接收
library(argparse)
parser = ArgumentParser()
parser$add_argument("--seurat_dir", help="The Dir for input seurat file.",required =T)
parser$add_argument("--sc_format", help="The format of input file",required =T,
choices = c('rds','qs',"RDS"),default = "rds")
parser$add_argument("--normalize_method",
help="Using normalizeData by default. SCT is optional",
choices = c('SCT','Standard'),default = "Standard")
parser$add_argument("--cluster", choices = c("T","F"),default = "F",
help="Whether clustering needs to be performed? TRUE or FALSE")
parser$add_argument("--save_anchors", choices = c("T","F"), default = "F",
help="Whether save anchors? If save please provide path.")
parser$add_argument("--outdir", help="Outdir",required =T)
parser$add_argument("--batchID", help="Batch ID",required =T)
args <- parser$parse_args()
input.data = args$seurat_dir
sc_format = args$sc_format
save_anchors = args$save_anchors
cluster_perform = args$cluster
outdir = args$outdir
batchID = args$batchID
normalize_method = args$normalize_method
##### 默认参数的一些判断
if(cluster_perform=='F'){
cluster_perform = FALSE
}else if(save_anchors=='T'){
cluster_perform = TRUE
}
if(normalize_method=='Standard'){
normalize_method = "Standard"
warning("Normalize method: Standard analysis by NormalizeData")
}else if(save_anchors=='SCT'){
normalize_method = "SCT"
warning("Normalize method: SCT")
}
if(save_anchors=='F'){
save_anchors = FALSE
}else if(save_anchors=='T'){
save_anchors = TRUE
}
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = TRUE)
print(paste("Folder", outdir, "created."))
}
######## Run CCA整合去批次
library(Seurat)
library(dplyr)
library(patchwork)
library(readr)
library(ggplot2)
library(RColorBrewer)
library(future)
library(cowplot)
library(stringr)
library(qs)
# change the current plan to access parallelization
plan("multisession", workers = 2)
plan()
#设置可用的内存
options(future.globals.maxSize = 10 * 1024^3)
#### Step0. 加载数据集
if(sc_format=="rds"|sc_format=="RDS"){
seurat.data = read_rds(file = input.data)
}else if (sc_format=="qs"){
seurat.data = qread(file = input.data)
} else {
stop("For now, input data only support seurat object for RDS or qs files!")
}
if(unique(table(seurat.data@meta.data[,batchID])<200)[1]){
stop("CCA requires a minimum of 200 cells in each batch for integration!")
}
#### Step1. split, anchors and Integrate
seurat_list <- SplitObject(seurat.data, split.by = batchID)
if(normalize_method=="SCT"){
##SCT
seurat_list <- lapply(X = seurat_list, FUN = SCTransform, method = "glmGamPoi")
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = seurat_list, nfeatures = 3000)
seurat_list <- PrepSCTIntegration(object.list = seurat_list, anchor.features = features)
seurat_anchors <- FindIntegrationAnchors(object.list = seurat_list, normalization.method = "SCT",
anchor.features = features)
if(save_anchors){
saveRDS(seurat_anchors,file = paste0(outdir,"/Step1.Seurat_anchors_CCA_SCT.rds"))
}
# 利用anchors进行整合
seurat_int <- IntegrateData(anchorset = seurat_anchors, normalization.method = "SCT")
}else if(normalize_method=="Standard")
{
# normalize and identify variable features for each dataset independently
seurat_list <- lapply(X = seurat_list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
return(x)
})
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = seurat_list)
#找到整合的锚点anchors
seurat_anchors <- FindIntegrationAnchors(object.list = seurat_list,
anchor.features = features, dims = 1:30)
if(save_anchors){
saveRDS(seurat_anchors,file = paste0(outdir,"/Step1.Seurat_anchors_CCA_Standard.rds"))
}
# 利用anchors进行整合
seurat_int <- IntegrateData(anchorset = seurat_anchors, dims = 1:30)
names(seurat_int@assays)
seurat_int@active.assay
}
#### Step2. 是否选择继续降维聚类？
if(cluster_perform){
#移除不用的对象
rm(seurat_list)
rm(seurat_anchors)
### 2.1. 降维聚类
n.pcs = 30
if(normalize_method=="Standard"){
seurat_int <- ScaleData(seurat_int, verbose = FALSE)
}
seurat_int <- seurat_int %>%
RunPCA(npcs = n.pcs, verbose = FALSE) %>%
FindNeighbors(dims = 1:n.pcs) %>%
RunUMAP(dims = 1:n.pcs)
for (res in c(0.05,0.1,0.3,0.5,0.8,1,1.2,1.4)){
print(paste0("FindClusters resolution: ",res))
seurat_int <- FindClusters(seurat_int,resolution = res, algorithm = 1)
}
#umap可视化
cluster_umap <- wrap_plots(ncol = 4,
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.0.05", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.0.1", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.0.3", label = T)& NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.0.5", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.0.8", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.1", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.1.2", label = T) & NoAxes(),
DimPlot(seurat_int, reduction = "umap", group.by = "integrated_snn_res.1.4", label = T)& NoAxes()
)
cluster_umap
ggsave(cluster_umap,filename = paste0(outdir,"/Step2.After_inter.cluster_UMAP_CCA.pdf"),
width = 20, height = 14)
}
#### Step3. save IntegrateData
if(normalize_method=="SCT"){
saveRDS(seurat_int,file = paste0(outdir,"/Seurat_IntegrateData_SCT_Standard.rds"))
}else if(normalize_method=="Standard"){
saveRDS(seurat_int,file = paste0(outdir,"/Seurat_IntegrateData_CCA_Standard.rds"))
}
基于写好的Rscript_seurat_CCA.R脚本，然后写一个Run_CCA.sh脚本，批量运行：
测试数据在：https://pan.baidu.com/s/1IPZ3G8Jle_UpV3S6ZRkS3Q?pwd=1234 提取码：1234
参数的话，看看上面的help内容，我都标注了：
cat >Run_CCA.sh
Rscript Rscript_seurat_CCA.R --seurat_dir ./CCA_test.rds --sc_format RDS --normalize_method Standard --cluster T --save_anchors T --outdir ./ --batchID stim
## Run
nohup bash Run_CCA.sh 1>Run_seurat.CCA.log 2>&1 &
运行完毕之后，目标目录下会输出结果，可以用于后续的继续分析：
上述脚本旨在一个抛砖引玉的作用，具体一些细节用户可以自行修改。另外如果有发现任何错误，欢迎大家后台留言进行批评指正。
传参脚本基本上是一劳永逸的，对于运行大样本的单细胞数据分析大有裨益。其实举一反三，CCA的脚本可以这样写，RPCA、Harmony和其他整合方法都可以写成一个传参R脚本，包括：
也有用python做过BBKNN和harmony：
另外注意，脚本的形式比较适合Linux服务器下运行。如果你所在的课题组没有服务器的话，我比较推荐两家租赁的云服务器：
生信技能树多人共享没有root权限，但是价格便宜：144线程640Gb内存服务器共享一年仍然是仅需800 西柚云服务器支持独享，有root权限可以做很多不一样的探索：独享服务器，满足你对生信分析的所有幻想 我向西柚云拿到了200的无门槛现金优惠，新用户注册即可获取，注册链接：https://www.xiyoucloud.net/aff/GKRYOBCA
最后，由衷希望大家都能找到属于自己的前进方向！
- END -
https://mp.weixin.qq.com/s/l9mwAQB_4lDvMhdTlBodTQ