时隔三个月，我带着新的单细胞笔记又来了。这是单细胞基础分析流程的最后一个板块，把这一部分搞懂之后，你就可以独立地完成一套标准的单细胞基础分析了，感兴趣的同学可以学起来呀！（当然，我们的笔记顺序并不是按照正常流程顺序来写的，所以我们后续也会专门做一篇整理性质的推送，为大家理清学习顺序的，这个大家不用担心！）
话不多说（那个那个，我，啊，这里的我是指小蛮要，咱们这篇文章是我们超级棒棒的小杨想去看极光写的，酷毙了！那我就按照惯例多说几句好不啦哈哈哈哈哈，谢谢大家！），我们直接进入正题！
-
如果小伙伴们有需求的话，可以加入我们的交流群：一定要知道 | 永久免费的环境友好型生信学习交流群又双叒叕来啦！| 伴随不定期群友好物分享！在这里，你可以稍有克制地畅所欲言！
-
超级建议大家在入群前或入群后可以看一下这个：干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！让你不虚此看！
-
如果有需要个性化定制分析服务的小伙伴，可以看看这里：你要的个性化生信分析服务今天正式开启啦！定制你的专属解决方案！全程1v1答疑！！绝对包你满意！
为什么要做多样本整合？
上一篇笔记我们讲了单细胞测序下机数据的比对质控流程，获取了相应的结果并对其中的核心内容进行了解读，详见单细胞基地 | Cell Ranger 助你轻松搞定上游分析 | 10X Genomics 平台下机数据比对质控流程。我们在文中提到了过滤矩阵文件夹filtered_feature_bc_matrix里面的信息已经满足后续分析的要求，可以直接进入seurat流程开展分析。但是呢，实际应用远没有理论情况那么理想，我们的过滤矩阵终究只是一个来自于单个样本的矩阵，无论是平行样本增加单细胞数据的生物学重复（理论上可以不用增加生物学重复，因为增加生物学重复说实在的就是在增加细胞数，但是现在的单细胞文章连细胞数都不卷的话，文章的竞争力也会大打折扣，可能发出来别人都不想相信），还是对照样本进行不同处理组间的比较分析，都涉及到多样本的情况，在标准工作流程下识别多个数据集中存在的细胞群可能会出现问题，因此多样本间的整合就很有必要。那么，“这种多个样本数据集的分析应该如何实现呢？”，这就是我们今天推送的主题–单细胞多样本数据整合。我将分成以下步骤进行阐述：准备工作、整合前样本的预处理、Seurat.v4多样本数据整合、Seurat.v5多样本数据整合。
0 准备工作
为了满足本篇笔记分析的需求，我从10X官网上下载了一组小鼠脑的单细胞测序下机数据，并按照上一篇笔记的教程单细胞基地 | Cell Ranger 助你轻松搞定上游分析 | 10X Genomics 平台下机数据比对质控流程使用cellranger v8.0.0和小鼠的GRCm39版本的参考基因集对其进行了比对定量工作。
# 下载数据并解压
# 下载数据
$ cat download.sh
# mouse_neuron_10k_1
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/4.0.0/SC3_v3_NextGem_SI_Neuron_10K/SC3_v3_NextGem_SI_Neuron_10K_fastqs.tar
# mouse_neuron_10k_2
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/4.0.0/SC3_v3_NextGem_DI_Neuron_10K/SC3_v3_NextGem_DI_Neuron_10K_fastqs.tar
# mouse_neuron_10k_3
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.0/neuron_10k_v3/neuron_10k_v3_fastqs.tar
# mouse_neuron_5k_1
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.2/5k_neuron_v3/5k_neuron_v3_fastqs.tar
# mouse_neuron_5k_2
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.2/5k_neuron_v3_nextgem/5k_neuron_v3_nextgem_fastqs.tar
# 解压
for line in $(ls ./*.tar)
do
tar -xvf ${line}
done
批量比对定量
# 先创建一个与样品编号对应的分组编号表格
$ cat sample_list.txt
SC3_v3_NextGem_SI_Neuron_10K mouse_neuron_1
SC3_v3_NextGem_DI_Neuron_10K mouse_neuron_2
neuron_10k_v3 mouse_neuron_3
5k_neuron_v3 mouse_neuron_4
5k_neuron_v3_nextgem mouse_neuron_5
# 编写批量比对定量脚本
$ cat run_qc.sh
reference=/workdir/BIO/yangqi/02.project/01.note/01.quality_control/01.genome/refdata-gex-GRCm39-2024-A
fq_path=/workdir/BIO/yangqi/02.project/01.note/01.quality_control/02.rawdata
cat sample_list.txt | while read col1 col2
do
cellranger count --id=$col2 \ # 输出文件的名字，即分析时的分组编号
--transcriptome=${reference} \
--fastqs=${fq_path}/${col1}_fastqs/ \ # 需要质控的下机fastq文件路径
--sample=$col1 \
# sample是输入文件的名字，即样品编号，通常是两种，一种是SRR开头的编号，一种是10X的文件名去掉最后的_fastqs
--create-bam=true \
--localcores=20 \
--localmem=200
done
结果简单展示
同时呢，考虑到有的同学自身电脑配置不支持比对定量所需的配置要求，我也把比对定量的结果文件（去除了.bam文件的）上传到了百度网盘，链接我设置了长期有效，大家可以根据需求下载，链接如下：
链接：https://pan.baidu.com/s/1oblJ47shAMTH5bqMrX7CAg?pwd=bten
提取码：bten
我们提供的文件信息如下图所示：
注：我们提供了除.bam文件和与矩阵等效的.h5文件之外的其他文件，感兴趣的同学可以结合我们以前的内容单细胞基地 | Cell Ranger 助你轻松搞定上游分析 | 10X Genomics 平台下机数据比对质控流程来理解输出文件内容，相信你会收获更多。
好，做完这个准备工作，就让我们进入正题。
1 整合前单样本的预处理
预处理流程主要完成两件事，分别是去除非细胞mRNA的污染(选做) 和 去除可能的双胞，这两个处理是针对单个样本进行的处理。其中呢，去除非细胞mRNA的污染(选做) 并不算是一个主流的预处理方法，并且也不会对最后结果产生很大的影响，相反，他会把矩阵中的整数counts变成小数counts，所以大家根据需要酌情使用。我这里只是见到这一个过程做一下简单介绍，感兴趣的同学可以去看看，没有需要的同学我们直接跳转 1.2 去双胞 过程。
1.1 SoupX去除非细胞mRNA的污染(选做)
参考：https://github.com/constantAmateur/SoupX
作用：用于估计和去除基于液滴的单细胞RNA-seq数据（10X的Chromium系统的油包水结构就是这种基于液滴的单细胞）中的非细胞mRNA污染。这一步做不做对后续分析影响不大，并且主流教程中也并没有过多涉及，大家可根据情况选择。
安装：我这边使用mamba进行安装，同时也需要安装一些配套使用的包 (其中seurat是V4的)，安装命令如下：
mamba install conda-forge::r-soupx
mamba install r::r-seurat
mamba install conda-forge::r-seuratobject=4.1.4
mamba install bioconda::bioconductor-dropletutils
代码解析：官网给的代码在实际运行过程中，直接运行load10X并不顺利，所以我结合了github官网上这个Issues45下面的作者回复，对代码进行了修改和整理，大家直接看下面的内容即可：
# 加载需要的R包
library(SoupX)
library(Seurat)
library(DropletUtils)
# 加载原始矩阵和过滤矩阵。作者给的命名，tod是cellranger结果中的原始矩阵；toc是过滤矩阵
tod <- Seurat::Read10X('/workdir/BIO/yangqi/02.project/01.note/01.quality_control/03.run_qc/mouse_neuron_1/outs/raw_feature_bc_matrix', gene.column = 2)
toc <- Seurat::Read10X('/workdir/BIO/yangqi/02.project/01.note/01.quality_control/03.run_qc/mouse_neuron_1/outs/filtered_feature_bc_matrix', gene.column = 2)
注：这里的gene.column为什么要等于2呢？(下面这点解释不用操作)
这是因为在cellranger的输出矩阵中features.tsv.gz里面包含了三列信息，第一列是gene ID，第二列是gene symbol，第三列是"Gene Expression"这样的字样，我们直接使用第二列的gene symbol信息，就不用后面再转换了，更方便。我随便打开一个矩阵的features.tsv，如下所示：
$head features.tsv
# ENSMUSG00000051951 Xkr4 Gene Expression
# ENSMUSG00000089699 Gm1992 Gene Expression
# ENSMUSG00000102331 Gm19938 Gene Expression
# ENSMUSG00000102343 Gm37381 Gene Expression
# ENSMUSG00000025900 Rp1 Gene Expression
# ENSMUSG00000025902 Sox17 Gene Expression
# ENSMUSG00000104238 Gm37587 Gene Expression
# ENSMUSG00000104328 Gm37323 Gene Expression
# ENSMUSG00000033845 Mrpl15 Gene Expression
# ENSMUSG00000120403 A930006A01Rik Gene Expression
好的，我们继续回到soupx上来，作者提到了需要使用setClusters这个函数添加一下聚类信息，这个咱们没有，需要简单做一下，固定参数即可，这些就不细讲啥意思了。
# 创建对象、标准化、降维、聚类、分群
object <- CreateSeuratObject(toc)
object <- SCTransform(object)
object <- RunPCA(object, features = VariableFeatures(object), npcs = 50, verbose = F)
object <- FindNeighbors(object, dims = 1:20)
object <- FindClusters(object, resolution = 0.5)
# 开始运行SoupX相关步骤
# 创建一个 SoupChannel 对象，该对象包含单个样本污染评估相关的所有内容（也就是原始矩阵和过滤矩阵）。
sc = SoupChannel(tod, toc)
# 提供聚类分群信息
sc = setClusters(sc, object@meta.data$seurat_clusters)
# autoEstCont实现污染自动估计
sc = autoEstCont(sc)
# adjustCounts调整产生新的矩阵
out = adjustCounts(sc)
# 按照10X的格式在当前路径下输出调整后的矩阵
DropletUtils:::write10xCounts("soupX_matrix", out, version="3")
我们可以看到，经过SoupX处理之后确实发生了一些变化，本来矩阵中的 整数counts 都变成了 小数形式，数值经过了略微的调整，但是变化很小，所以这一步做不做影响都不会很大。
> str(sc)
List of 5
$ toc :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. ..@ i : int [1:46803316] 0 8 11 14 18 23 25 26 35 41 ...
.. ..@ p : int [1:12442] 0 4717 10065 14643 15306 21127 21967 26714 30156 34737 ...
.. ..@ Dim : int [1:2] 33696 12441
.. ..@ Dimnames:List of 2
.. .. ..$ : chr [1:33696] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
.. .. ..$ : chr [1:12441] "AAACCAAAGCCTGTCC-1" "AAACCAAAGCTATCGA-1" "AAACCATTCACCACTT-1" "AAACCATTCGATTAAC-1" ...
.. ..@ x : num [1:46803316] 6 1 3 2 1 1 8 1 1 1 ...
.. ..@ factors : list()
> str(out)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:46803316] 0 8 11 14 18 23 25 26 35 41 ...
..@ p : int [1:12442] 0 4717 10065 14643 15306 21127 21967 26714 30156 34737 ...
..@ Dim : int [1:2] 33696 12441
..@ Dimnames:List of 2
.. ..$ : chr [1:33696] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
.. ..$ : chr [1:12441] "AAACCAAAGCCTGTCC-1" "AAACCAAAGCTATCGA-1" "AAACCATTCACCACTT-1" "AAACCATTCGATTAAC-1" ...
..@ x : num [1:46803316] 5.981 0.984 2.973 1.989 0.986 ...
..@ factors : list()
1.2 Doubletfinder去除双胞
参考：https://github.com/chris-mcginnis-ucsf/DoubletFinder
作用：用于预测并去除数据中双胞。
安装：Doubletfinder的数据结构随着seurat的版本更新进行了更新，更新的分界线是2.0.4，在2.0.4以前的版本，都是适配于seurat v4的；而在2.0.4及以后的版本都是依赖于seurat v5版本实现的。不过大家不要担心，我们这里会为大家准备两个不同版本的代码，大家可以按需选择，不过更推荐大家使用新版，及时适应seurat v5的数据结构。
DoubletFinder v2.0.3 及配套seurat V4安装
R命令行安装：DoubletFinder
remotes::install_github('https://github.com/ekernf01/DoubletFinder', force = T)
shell命令行安装下述软件
mamba install r::r-seurat
mamba install conda-forge::r-seuratobject=4.1.4
mamba install bioconda::bioconductor-dropletutils
DoubletFinder v2.0.4 及配套seurat V5安装
R命令行安装：DoubletFinder
remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')
shell命令行安装下述软件
mamba install conda-forge::r-seurat
mamba install bioconda::bioconductor-dropletutils
代码解析：
首先呢，我们先要输入数据，关于输入的数据，官网提了两点要求：
- Do not apply DoubletFinder to aggregated scRNA-seq data representing multiple distinct samples (e.g., multiple 10X lanes).也就是说，DoubletFinder是对单个样本处理的，不要对多样本数据进行操作。
- Ensure that input data is cleared of low-quality cell clusters.确保输入数据中不含低质量细胞簇。也就是说，我们需要提前对数据进行过滤。
那么，我们就先按照这两个要求对数据简单进行一个过滤处理吧。
# 加载需要的R包
library(Seurat)
library(cowplot)
library(DropletUtils)
library(DoubletFinder)
# 读入单细胞矩阵数据，并按照官网要求进行过滤预处理；我们先读取 filtered_feature_bc_matrix
object <- Read10X('your/path/to/filtered_feature_bc_matrix', gene.column = 2)
# 创建seurat对象
object <- CreateSeuratObject(counts = object.data, project = "brain_1", min.cells = 3, min.features = 200)
# 计算线粒体基因百分比，0值补成非0的极小数
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "ATP6|ATP8|COX1|COX2|COX3|CYTB|ND1|ND2|ND3|ND4|ND4L|ND5|ND6")
if (sum(object@meta.data$percent.mt) == 0){
object@meta.data$percent.mt[1] <- 0.000001
}
# 按照相对宽泛的指标过滤
# 这里呢，我简单的画了三个变量的小提琴图，评估了一下过滤指标，暂定范围如下：
# 200 < nFeature_RNA(单个细胞总基因数) < 6000;
# nCount_RNA(单个细胞总表达量数) < 20000;
# percent.mt(单个细胞内线粒体基因比例) < 5
# 注：这个是脑的数据，基因数相对来说会多一些，常规组织的话要根据组织特性设置。
object <- subset(object, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 20000 & percent.mt < 5)
接着呢，官网让我们进行一下seurat的预处理流程，提供了经典的NormalizeData标准化和SCT标准化两种方法，大家按需任选其一即可，区别不大。如果选择困难的话，我这边选的是SCT。如果你用的是三步标准化的话，记得后面每一步的sct参数的值设置成FALSE。
# 进行seurat流程的预处理(三步标准化和SCT标准化都可以，选其一即可，区别不大)
# 三步标准化流程
object <- NormalizeData(object)
object <- FindVariableFeatures(object, selection.method = "vst", nfeatures = 2000)
object <- ScaleData(object)
object <- RunPCA(object)
object <- RunUMAP(object, dims = 1:20)
# SCT标准化流程
object <- SCTransform(object)
object <- RunPCA(object)
object <- RunUMAP(object, dims = 1:20)
# 为什么比官网多了聚类分群这两步呢，因为后面需要给一个注释信息，我们用分群结果代替
object <- FindNeighbors(object, reduction = "pca", dims = 1:20)
object <- FindClusters(object, resolution = 0.5)
然后，就进入到DoubletFinder的特定步骤了，我们简单解析一下，就是他需要一个pK值，在这个pK值下能保证DoubletFinder的最佳预测准确性。这个值怎么选呢，官网提供了一个BCmvn的指标(Mean-variance normalized bimodality coefficient (BCmvn) achieves this goal, featuring a single, easily-discernible maximum at pK values that optimize AUC.)，叫什么均值方差归一化双峰系数，在这个BCmvn值最大的时候对应的pK值就是我们需要的。明白了这一点，我们就开始整理和修改代码吧。
# pK值确认(这里官网提供了有真实值参考和无真实值参考两种代码，我们这里使用无真实值参考的)
sweep.res.list_object <- paramSweep_v3(object, PCs = 1:20, sct = TRUE)
sweep.stats_object <- summarizeSweep(sweep.res.list_object, GT = FALSE)
bcmvn_object <- find.pK(sweep.stats_object)
dev.off()
这里find.pK函数会出一张图，看的是不同pK值下的BCmvn值，我们这里简单展示一下。横坐标是pK值，纵坐标是BCmvn值。
但是呢，我们是需要取这个合适的pK值，所以还是看对应的BCmvn值的dataframe来选择。
# 查看一下这个bcmvn_object
> head(bcmvn_object)
# ParamID pK MeanBC VarBC BCmetric
# 1 1 0.005 0.7574656 0.002022751 374.47304
# 2 2 0.01 0.7309493 0.002013150 363.08745
# 3 3 0.02 0.7107600 0.003004541 236.56195
# 4 4 0.03 0.6654097 0.004069111 163.52705
# 5 5 0.04 0.6253637 0.006920117 90.36895
# 6 6 0.05 0.6108826 0.015017486 40.67809
# 我们选择BCmvn对应列BCmetric的最大值对应的pK值，代码如下：
pK_value <- as.numeric(as.character(bcmvn_object$pK[bcmvn_object$BCmetric == max(bcmvn_object$BCmetric)]))
# 注：官网上说pN对DoubletFinder的影响不大，那么我们就用默认的0.25
# 同型双胞比例评估
# 这里的第一行就是将聚类分群信息当作注释信息，用于后续分析
annotations <- object@meta.data$SCT_snn_res.0.5
homotypic.prop <- modelHomotypic(annotations)
nExp_poi <- round(0.075*nrow(object@meta.data)) # 按照官网的标准，设定为7.5%
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop)) # 调整校准nExp_poi
# 使用不同的分类严格性运行DoubletFinder
# 构造reuse.pANN参数的值pANN_value
pANN_value <- paste0("pANN_0.25_", pK_value, '_', nExp_poi)
# 分别用正常值nExp_poi和矫正值nExp_poi.adj进行双胞计算
object <- doubletFinder_v3(object, PCs = 1:30, pN = 0.25, pK = pK_value, nExp = nExp_poi, reuse.pANN = FALSE, sct = TRUE)
object <- doubletFinder_v3(object, PCs = 1:30, pN = 0.25, pK = pK_value, nExp = nExp_poi.adj, reuse.pANN = pANN_value, sct = TRUE)
到这里，官网上的基础教程就做完了。接下来我们补充一个知识点：
我们可以看到在代码结束之后官网上还有这么一张图，那么这个图上的 “Doublet-Low Confidence” 和 “Doublet-High Confidence” 是什么意思呢？什么是高可信度，什么是低可信度呢？
我这样讲，咱们最后在上面运行了两次"doubletFinder_v3"函数，两组函数的差别主要是 “nExp” 和 “reuse.pANN” 两个参数的值，其中呢，用 “nExp_poi.adj” 值确定的双胞是高可信度的双胞，而使用 “nExp_poi” 值确定的双胞包含了高可信度和低可信度的双胞。这个不重要，不过大家学东西最好学透。
那我们接着就尝试把这些都区分一下吧。前面说过了，这个不重要，所以大家当作一个理解和熟悉meta.data的练习题吧。
# 我们先查看一下meta.data的信息，可以看出来meta.data分成了11列，我们关注10和11列
> head(object@meta.data)
orig.ident nCount_RNA nFeature_RNA percent.mt nCount_SCT
AAACCTGCACCAGTTA-1 D1_1 13672 3417 2.4575775 8262
AAACCTGCACCATGTA-1 D1_1 21279 3981 2.3826308 6867
AAACCTGGTGAAATCA-1 D1_1 12481 3365 3.4372246 8507
AAACCTGTCAGTCCCT-1 D1_1 25707 5562 2.1278251 7718
AAACGGGAGGACCACA-1 D1_1 8549 1423 0.9240847 8094
AAACGGGTCAAACCAC-1 D1_1 55970 6155 2.2583527 5830
nFeature_SCT SCT_snn_res.0.5 seurat_clusters
AAACCTGCACCAGTTA-1 3245 7 7
AAACCTGCACCATGTA-1 2253 10 10
AAACCTGGTGAAATCA-1 3298 16 16
AAACCTGTCAGTCCCT-1 3204 13 13
AAACGGGAGGACCACA-1 1422 3 3
AAACGGGTCAAACCAC-1 2364 2 2
pANN_0.25_0.03_183 DF.classifications_0.25_0.03_183
AAACCTGCACCAGTTA-1 0.08247423 Singlet
AAACCTGCACCATGTA-1 0.20618557 Singlet
AAACCTGGTGAAATCA-1 0.23711340 Singlet
AAACCTGTCAGTCCCT-1 0.36082474 Doublet
AAACGGGAGGACCACA-1 0.16494845 Singlet
AAACGGGTCAAACCAC-1 0.09278351 Singlet
DF.classifications_0.25_0.03_171
AAACCTGCACCAGTTA-1 Singlet
AAACCTGCACCATGTA-1 Singlet
AAACCTGGTGAAATCA-1 Singlet
AAACCTGTCAGTCCCT-1 Doublet
AAACGGGAGGACCACA-1 Singlet
AAACGGGTCAAACCAC-1 Singlet
# 我们先看一下10和11列的不同信息的数量情况，能看出来，高可信度的双胞171个，低可信度的183-171=12个
> table(object@meta.data$DF.classifications_0.25_0.03_183)
Doublet Singlet
183 2254
> table(object@meta.data$DF.classifications_0.25_0.03_171)
Doublet Singlet
171 2266
# 那么，接下来就让我们分别把他们区分出来吧
# 首先，把第10列赋值给新一列"DF_high.low"用来存储高低可信度的双胞。
object@meta.data[,"Doublet"] <- object@meta.data[,10]
# 接着，"DF_high.low"这一列中是"Doublet"，但是第11列中是"Singlet"的判定为"Doublet-Low Confidence"，我们给他命名为"Doublet_low"
object@meta.data$Doublet[which(object@meta.data$Doublet == "Doublet" & object@meta.data[,11] == "Singlet")] <- "Doublet_low"
# 然后，"DF_high.low"这一列中还是"Doublet"就是"Doublet_high"了
object@meta.data$Doublet[which(object@meta.data$Doublet == "Doublet")] <- "Doublet_high"
# 大家感兴趣的话，也可以画一下类似上面的图，我用这个2000多细胞的数据画出来是这样的。可能效果不是很好，大家将就着看！
pdf(paste0(opts$sample, "_dimplot.pdf"), width=20, height=10)
plot1 <- DimPlot(object, reduction = "umap", group.by = "Doublet", pt.size = 0.5)
plot2 <- VlnPlot(object, group.by = "Doublet", features = c("nCount_RNA", "nFeature_RNA"), pt.size = 0, ncol = 2)
plot1 + plot2
dev.off()
# 这样，我们就可以去除双胞，只取单胞了
object_Singlet <- subset(x = object, subset = Doublet == "Singlet")
# 同样的，输出为10X的标准格式。
DropletUtils:::write10xCounts("singlet_matrix", object_Singlet@assays$RNA@counts, version="3")
最后，我们编写一下批量运行的脚本：(注意：脚本里面的参数大家可以根据自己的数据自行设置。)
$ cat run_doubletfinder.sh
input_path=/workdir/BIO/yangqi/02.project/01.note/03.DoubletFinder/01.data
output_path=/workdir/BIO/yangqi/02.project/01.note/03.DoubletFinder/02.run_doubletfinder
R_path=/workdir/BIO/yangqi/01.software/mamba/mambaforge/envs/seurat.v4/bin
for line in $(ls ${input_path}/)
do
${R_path}/Rscript ./run_doubletfinder.R -i ${input_path}/${line} -s ${line} -f 6000 -c 20000 -m 5 -o ${output_path}/${line}
done
关于DoubletFinder v2.0.4（即适配于seurat v5）的使用
1.改了什么？
详见DoubletFinder官网Updates部分的(11/21/2023)这一行:
Made compatible with Seurat v5 and removed ‘_v3’ flag from relevant function names.
翻译一下就是：改成了适应seurat v5的版本，以及修改了部分函数的命名。
Seurat v5数据结构的变化，详细信息可以见这篇笔记单细胞基地 | Seurat v5 基础分析流程及其改动 | v5 和 v4 到底有啥不一样嘞！
函数名字修改就是"paramSweep_v3"变成了"paramSweep"; “doubletFinder_v3"变成了"doubletFinder”
2.使用的时候应该注意什么？
注意两件事：
① 配置个新环境（按照上面安装部分的 v2.0.4 的命令来）；
② 脚本换一下（按照下面的来，跟上面的解析类似，这里就不逐行解析了）。
# 加载需要的R包
library(Matrix)
library(Seurat)
library(cowplot)
library(DropletUtils)
library(DoubletFinder)
# 读入单细胞矩阵数据，创建seurat对象，并按照要求进行过滤预处理
object <- Read10X('your/path/to/filtered_feature_bc_matrix', gene.column = 2)
object <- CreateSeuratObject(counts = object.data, project = "brain_1", min.cells = 3, min.features = 200)
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "^mt-")
if (sum(object@meta.data$percent.mt) == 0){
object@meta.data$percent.mt[1] <- 0.000001
}
object <- subset(object, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 20000 & percent.mt < 5)
# 进行seurat流程的预处理
object <- SCTransform(object)
object <- RunPCA(object)
object <- RunUMAP(object, dims = 1:20)
object <- FindNeighbors(object, reduction = "pca", dims = 1:20)
object <- FindClusters(object, resolution = 0.5)
# pK值确认
sweep.res.list_object <- paramSweep(object, PCs = 1:20, sct = TRUE)
sweep.stats_object <- summarizeSweep(sweep.res.list_object, GT = FALSE)
bcmvn_object <- find.pK(sweep.stats_object)
dev.off()
pK_value <- as.numeric(as.character(bcmvn_object$pK[bcmvn_object$BCmetric == max(bcmvn_object$BCmetric)]))
# 同型双胞比例评估
annotations <- object@meta.data$SCT_snn_res.0.5
homotypic.prop <- modelHomotypic(annotations)
nExp_poi <- round(0.075*nrow(object@meta.data))
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
# 找双胞
pANN_value <- paste0("pANN_0.25_", pK_value, '_', nExp_poi)
object <- doubletFinder(object, PCs = 1:20, pN = 0.25, pK = pK_value, nExp = nExp_poi, reuse.pANN = FALSE, sct=TRUE)
object <- doubletFinder(object, PCs = 1:20, pN = 0.25, pK = pK_value, nExp = nExp_poi.adj, reuse.pANN = pANN_value, sct=TRUE)
object@meta.data[,"Doublet"] <- object@meta.data[,10]
object@meta.data$Doublet[which(object@meta.data$Doublet == "Doublet" & object@meta.data[,11] == "Singlet")] <- "Doublet_low"
object@meta.data$Doublet[which(object@meta.data$Doublet == "Doublet")] <- "Doublet_high"
object_Singlet <- subset(x = object, subset = Doublet == "Singlet")
# 存储结果
# v5的数据结构需要对其进行特定的转换
mtx <- as.matrix(object_Singlet@assays$RNA@layers$counts)
rownames(mtx) <- rownames(object_Singlet@assays$RNA@features)
colnames(mtx) <- rownames(object_Singlet@assays$RNA@cells)
mtx <- Matrix(mtx, sparse = TRUE)
DropletUtils:::write10xCounts("singlet_matrix", mtx, version="3")
上面我们做完了预处理的相关工作，接下来，我们就开始进行多样本的整合吧。由于v4和v5是有区别的，所以我们会分别进行解析，大家按需学习。
2 Seurat.v4多样本数据整合
这里我们先提前介绍一下我们会使用的整合方法，在标准化部分呢，我们会尝试NormalizeData标准化 和 SCTransform标准化 (想要了解两者的区别，可以看曾老师的推送SCTransform真的能完美替代Seurat早期的3个函数吗的全文和评论区；还有这个结论)；关于整合方法的话，我们会使用seurat官网的 CCA锚点整合 和 harmony整合 两种方法；至于其他的内容，我们暂时不做涉及。
2.0 准备工作
首先，我们先将去除双胞后的结果链接到新路径下
for line in {1..5}; do ln -s /workdir/BIO/yangqi/02.project/01.note/03.DoubletFinder/02.run_doubletfinder/mouse_neuron_${line}/singlet_matrix ./mouse_neuron_${line}; done
接着，我们要读取每个样本的数据，计算相关的基础指标 (注意这里就不需要再过滤了)。
# 载入需要的包
library(Seurat)
library(cowplot)
# 获取数据路径
matrix_list <- list.files("/workdir/BIO/yangqi/02.project/01.note/04.Integration/00.data")
Subset_Cells.list = vector()
# 批量处理矩阵
for (index in 1:length(matrix_list)) {
matrix <- matrix_list[index]
matrix_path <- paste0("/workdir/BIO/yangqi/02.project/01.note/04.Integration/00.data", '/', matrix)
object.data <- Read10X(matrix_path, gene.column = 2)
object <- CreateSeuratObject(counts = object.data, project = matrix, min.cells = 3, min.features = 200)
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "^mt-")
if (sum(object@meta.data$percent.mt) == 0){
object@meta.data$percent.mt[1] <- 0.000001
}
Subset_Cells.list <- append(Subset_Cells.list, object)
}
上面这个一个循环做了什么事情呢？我用一句话帮大家解释一下，就是批量读取了去双胞后产生的矩阵，然后用他们创建了seurat对象，并重新计算了过滤指标（但是不过滤），存储在一个list列表Subset_Cells.list里。现在让我们来查看一下这个 Subset_Cells.list 长什么样。
> Subset_Cells.list
[[1]]
An object of class Seurat
21544 features across 4772 samples within 1 assay
Active assay: RNA (21544 features, 0 variable features)
[[2]]
An object of class Seurat
21568 features across 4587 samples within 1 assay
Active assay: RNA (21568 features, 0 variable features)
[[3]]
An object of class Seurat
21297 features across 6467 samples within 1 assay
Active assay: RNA (21297 features, 0 variable features)
[[4]]
An object of class Seurat
20808 features across 3672 samples within 1 assay
Active assay: RNA (20808 features, 0 variable features)
[[5]]
An object of class Seurat
21002 features across 3222 samples within 1 assay
Active assay: RNA (21002 features, 0 variable features)
很明显可以看到，这个列表里面存储了5个seurat对象，分别对应我们上面的五个样本。
准备工作完成后，我们接着该怎么处理呢？无非就是两种整合思路，seurat官网的锚点整合和harmony整合，我们逐个进行解析。
2.1 锚点整合
参考： seurat v4官网的整合教程：https://satijalab.org/seurat/archive/v4.3/integration_introduction
代码详解：
seurat官网上提示我们在整合之前，需要对每个矩阵做标准化处理，标准化有两种思路，分别是seurat流程中的NormalizeData标准化 和 SCTransform标准化 ，我们分别介绍一下：
2.1.1 NormalizeData标准化 + 锚点整合CCA
# 对单个样本数据逐一进行标准化和找高变基因
Subset_Cells.list <- lapply(X = Subset_Cells.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# 寻找不同数据集中的都有的高变基因，用于数据集的整合
features <- SelectIntegrationFeatures(object.list = Subset_Cells.list)
# 用FindIntegrationAnchors函数识别锚点并使用IntegrateData函数根据这些锚点将多个数据集整合在一起
object.anchors <- FindIntegrationAnchors(object.list = Subset_Cells.list, anchor.features = features)
object.combined <- IntegrateData(anchorset = object.anchors)
# 对整合后的数据进行分析，看整合的效果如何
# 我们先切换到"integrated"这个矩阵中
DefaultAssay(immune.combined) <- "integrated"
# 然后跑seurat的标准流程进行计算，里面的部分参数值需要根据自身数据情况自主设置，我这里图省事，用了官网现成的。
object.combined <- ScaleData(object.combined, verbose = FALSE)
object.combined <- RunPCA(object.combined, npcs = 30, verbose = FALSE)
object.combined <- RunUMAP(object.combined, reduction = "pca", dims = 1:30)
object.combined <- FindNeighbors(object.combined, reduction = "pca", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# 我们按照样本分组可视化一下整合效果
pdf("Dimplot_group.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储一下整合的结果
saveRDS(object.combined, file = "integration_Normal_CCA.rds")
我们来看一下效果，至少从单样本UMAP图形状和细胞分群情况上来看挺好的。关于基因的表达情况，我们还未可知，不过我想不会太差。
2.1.2 SCTransform标准化 + 锚点整合CCA
这个部分唯一的区别就是标准化方法变成了SCTransform标准化，其他都不变，我就直接展示代码：
# 批量SCT标准化
Subset_Cells.list <- lapply(X = Subset_Cells.list, FUN = SCTransform)
# CCA锚点整合，这里的标准化方法参数记得换成SCT
features <- SelectIntegrationFeatures(object.list = Subset_Cells.list, nfeatures = 3000)
Subset_Cells.list <- PrepSCTIntegration(object.list = Subset_Cells.list, anchor.features = features)
object.anchors <- FindIntegrationAnchors(object.list = Subset_Cells.list, normalization.method = "SCT", anchor.features = features)
object.combined <- IntegrateData(anchorset = object.anchors, normalization.method = "SCT")
# seurat常规流程，降维聚类分群
object.combined <- RunPCA(object.combined, npcs = 30, verbose = FALSE)
object.combined <- RunUMAP(object.combined, reduction = "pca", dims = 1:30)
object.combined <- FindNeighbors(object.combined, reduction = "pca", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# 可视化结果
pdf("Dimplot_group_SCT_CCA.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_SCT_CCA.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_SCT_CCA.rds")
结果如下图所示，差别不是很大：
2.2 Harmony整合
参考： https://github.com/immunogenomics/harmony
安装：
注： 这里安装不要使用mamba，否则会导致harmony运行报错。详细信息见https://github.com/immunogenomics/harmony/issues/169
devtools::install_github("eddelbuettel/harmony", force = TRUE)
代码解析：
首先呢，我们依然可以使用准备工作中处理好的 Subset_Cells.list，同样地分为两种标准化方法，我们分别展示代码：
2.2.1 NormalizeData标准化 + Harmony整合
# 导入harmony包
library(harmony)
# 将不同的样本数据merge在一块
object.combined <- merge(Subset_Cells.list[[1]], y = Subset_Cells.list[-1], add.cell.ids = names(Subset_Cells.list), project = "neuron")
# seurat常规流程
object.combined <- NormalizeData(object = object.combined, normalization.method = "LogNormalize", scale.factor = 10000)
object.combined <- FindVariableFeatures(object = object.combined, selection.method = "vst", nfeatures = 2000)
object.combined <- ScaleData(object = object.combined, verbose = FALSE)
object.combined <- RunPCA(object = object.combined, npcs = 100, verbose = FALSE)
# 跑harmony整合
object.combined <- RunHarmony(object.combined, group.by.vars = "orig.ident")
# 降维，聚类
object.combined <- RunUMAP(object.combined, reduction = "harmony", dims = 1:30)
object.combined <- FindNeighbors(object.combined, reduction = "harmony", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# 画图展示整合效果
pdf("Dimplot_group_Normal_Harmony.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_Normal_Harmony.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_Normal_Harmony.rds")
结果如下图所示，整合效果也不错：
2.1.2 SCTransform标准化 + Harmony整合
# 导入harmony包
library(harmony)
# 将不同的样本数据merge在一块
object.combined <- merge(Subset_Cells.list[[1]], y = Subset_Cells.list[-1], merge.data = TRUE)
# SCT标准化
object.combined <- SCTransform(object.combined, assay = 'RNA', return.only.var.genes = F, verbose = F)
# PCA降维
object.combined <- RunPCA(object.combined, npcs = 30, verbose = FALSE)
# 跑harmony整合
object.combined <- RunHarmony(object.combined, group.by.vars = "orig.ident", assay.use = 'SCT')
# 降维，聚类
object.combined <- RunUMAP(object.combined, reduction = "harmony")
object.combined <- FindNeighbors(object.combined, reduction = "pca", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# 画图展示整合效果
pdf("Dimplot_group_SCT_Harmony.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_SCT_Harmony.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_SCT_Harmony.rds")
结果看起来也整合到一块了
3 Seurat.v5多样本数据整合
Seurat v5部分的整合呢，官方把目前主流的一些方法也集成了进去，这些在官网上也有详细的说明，👉指路，概括起来说呢，基本上是以下五种：
① 基于锚点的CCA整合 (method=CCAIntegration)
# 核心代码
obj <- IntegrateLayers(object = obj, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca", verbose = FALSE)
② 基于锚点的RPCA整合 (method=RPCAIntegration)：更快、更保守（更少校正）
# 核心代码
obj <- IntegrateLayers(object = obj, method = RPCAIntegration, orig.reduction = "pca", new.reduction = "integrated.rpca", verbose = FALSE)
③ Harmony整合 (method=HarmonyIntegration)
# 核心代码
obj <- IntegrateLayers(object = obj, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony", verbose = FALSE)
④ FastMNN整合 (method= FastMNNIntegration)
# 核心代码
obj <- IntegrateLayers(object = obj, method = FastMNNIntegration, new.reduction = "integrated.mnn", verbose = FALSE)
⑤ scVI整合 (method=scVIIntegration)
obj <- IntegrateLayers(object = obj, method = scVIIntegration, new.reduction = "integrated.scvi", conda_env = "../miniconda3/envs/scvi-env", verbose = FALSE)
不过这里呢，我们还是主要为大家展示主流的CCA和harmony，剩下的整合方法大家可以自主探索，这些方法没有太多的优劣之分，主要还是适合自己的就是最好的，大家酌情选择，核心代码一换，基本上都是适用的。
3.0 准备工作
首先呢，大家需要对seurat.v5的数据结构有一定的认识，我们在这篇笔记单细胞基地 | Seurat v5 基础分析流程及其改动 | v5 和 v4 到底有啥不一样嘞！做了简单的整理，大家可以先去补充一下相关知识。在这里我们简单地向大家说明一下seurat.v5数据集列表的构建，因为这是我们整合过程中的关键。
首先呢，我们按照 “2.0 准备工作” 先批量读进矩阵创建列表Subset_Cells.list
# 载入需要的包
library(Seurat)
library(cowplot)
# 批量读取矩阵，创建seurat对象列表Subset_Cells.list
matrix_list <- list.files("/workdir/BIO/yangqi/02.project/01.note/04.Integration/00.data")
Subset_Cells.list = vector()
for (index in 1:length(matrix_list)) {
matrix <- matrix_list[index]
matrix_path <- paste0("/workdir/BIO/yangqi/02.project/01.note/04.Integration/00.data", '/', matrix)
object.data <- Read10X(matrix_path, gene.column = 2)
object <- CreateSeuratObject(counts = object.data, project = matrix, min.cells = 3, min.features = 200)
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "^mt-")
if (sum(object@meta.data$percent.mt) == 0){
object@meta.data$percent.mt[1] <- 0.000001
}
Subset_Cells.list <- append(Subset_Cells.list, object)
}
接着，就是重点操作，我们要如何把 Subset_Cells.list 改造成适合于seurat.v5的结构呢？
我觉得官网上的示例数据的数据结构我们需要先学习一下，我们先看看seurat.v5官网上的示例数据的结构，然后按照这个结构将我们的数据改造一下吧。
# 加载需要的R包
library(Seurat)
library(SeuratData)
library(patchwork)
# 下载ifnb这个数据集
InstallData("ifnb")
# 加载ifnb数据集
ifnb <- LoadData("ifnb")
# 查看一下此时ifnb的数据结构，可以看到是一个counts矩阵，一个data矩阵
> str(ifnb)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
.. .. .. ..@ layers :List of 2
.. .. .. .. ..$ counts:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:9787436] 20 27 37 64 65 83 87 131 139 175 ...
.. .. .. .. .. .. ..@ p : int [1:14000] 0 877 1590 2440 3549 4183 4740 5720 6301 7181 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 13999
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:9787436] 1 1 1 1 1 2 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ data :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:9787436] 20 27 37 64 65 83 87 131 139 175 ...
.. .. .. .. .. .. ..@ p : int [1:14000] 0 877 1590 2440 3549 4183 4740 5720 6301 7181 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 13999
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:9787436] 1 1 1 1 1 2 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. ..@ cells :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:13999, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:13999] "AAACATACATTTCC.1" "AAACATACCAGAAA.1" "AAACATACCTCGCT.1" "AAACATACCTGGTA.1" ...
.. .. .. .. .. .. .. ..$ : chr [1:2] "counts" "data"
.. .. .. .. .. ..$ dim : int [1:2] 13999 2
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:13999] "AAACATACATTTCC.1" "AAACATACCAGAAA.1" "AAACATACCTCGCT.1" "AAACATACCTGGTA.1" ...
.. .. .. .. .. .. ..$ : chr [1:2] "counts" "data"
.. .. .. ..@ features :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:14053, 1:2] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:14053] "AL627309.1" "RP11-206L10.2" "LINC00115" "NOC2L" ...
.. .. .. .. .. .. .. ..$ : chr [1:2] "counts" "data"
.. .. .. .. .. ..$ dim : int [1:2] 14053 2
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:14053] "AL627309.1" "RP11-206L10.2" "LINC00115" "NOC2L" ...
.. .. .. .. .. .. ..$ : chr [1:2] "counts" "data"
.. .. .. ..@ default : int 1
.. .. .. ..@ assay.orig: chr(0)
.. .. .. ..@ meta.data :'data.frame': 14053 obs. of 0 variables
.. .. .. ..@ misc : list()
.. .. .. ..@ key : chr "rna_"
我们看到官网上依据ifnb$stim信息，把ifnb拆分了。那么我们在拆分前先看一下ifnb$stim里面是什么信息。
> head(ifnb)
orig.ident nCount_RNA nFeature_RNA stim seurat_annotations
AAACATACATTTCC.1 IMMUNE_CTRL 3017 877 CTRL CD14 Mono
AAACATACCAGAAA.1 IMMUNE_CTRL 2481 713 CTRL CD14 Mono
AAACATACCTCGCT.1 IMMUNE_CTRL 3420 850 CTRL CD14 Mono
AAACATACCTGGTA.1 IMMUNE_CTRL 3156 1109 CTRL pDC
AAACATACGATGAA.1 IMMUNE_CTRL 1868 634 CTRL CD4 Memory T
AAACATACGGCATT.1 IMMUNE_CTRL 1581 557 CTRL CD14 Mono
AAACATACTGCGTA.1 IMMUNE_CTRL 2747 980 CTRL T activated
AAACATACTGCTGA.1 IMMUNE_CTRL 1341 581 CTRL CD4 Naive T
AAACATTGAGTGTC.1 IMMUNE_CTRL 2155 880 CTRL CD8 T
AAACATTGCTTCGC.1 IMMUNE_CTRL 2536 669 CTRL CD14 Mono
> table(ifnb$stim)
CTRL STIM
6548 7451
很明显这个数据集由两个样本构成的，分别是CTRL和STIM，这就是拆分的依据，我们按照官网的命令拆分一下吧。
ifnb[["RNA"]] <- split(ifnb[["RNA"]], f = ifnb$stim)
ifnb
# 查看一下此时ifnb的数据结构，可以看到是两个counts矩阵，两个data矩阵
# 说明ifnb被根据样本信息拆分成了两份。
str(ifnb)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
.. .. .. ..@ layers :List of 4
.. .. .. .. ..$ counts.CTRL:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:4626015] 20 27 37 64 65 83 87 131 139 175 ...
.. .. .. .. .. .. ..@ p : int [1:6549] 0 877 1590 2440 3549 4183 4740 5720 6301 7181 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 6548
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:4626015] 1 1 1 1 1 2 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.STIM:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:5161421] 7 10 20 25 27 79 87 113 176 194 ...
.. .. .. .. .. .. ..@ p : int [1:7452] 0 588 1380 1964 2696 3242 4558 5140 5702 6210 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 7451
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:5161421] 15 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ data.CTRL :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:4626015] 20 27 37 64 65 83 87 131 139 175 ...
.. .. .. .. .. .. ..@ p : int [1:6549] 0 877 1590 2440 3549 4183 4740 5720 6301 7181 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 6548
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:4626015] 1 1 1 1 1 2 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ data.STIM :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:5161421] 7 10 20 25 27 79 87 113 176 194 ...
.. .. .. .. .. .. ..@ p : int [1:7452] 0 588 1380 1964 2696 3242 4558 5140 5702 6210 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 14053 7451
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:5161421] 15 1 1 1 1 1 1 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. ..@ cells :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:13999, 1:4] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:13999] "AAACATACATTTCC.1" "AAACATACCAGAAA.1" "AAACATACCTCGCT.1" "AAACATACCTGGTA.1" ...
.. .. .. .. .. .. .. ..$ : chr [1:4] "counts.CTRL" "counts.STIM" "data.CTRL" "data.STIM"
.. .. .. .. .. ..$ dim : int [1:2] 13999 4
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:13999] "AAACATACATTTCC.1" "AAACATACCAGAAA.1" "AAACATACCTCGCT.1" "AAACATACCTGGTA.1" ...
.. .. .. .. .. .. ..$ : chr [1:4] "counts.CTRL" "counts.STIM" "data.CTRL" "data.STIM"
.. .. .. ..@ features :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:14053, 1:4] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:14053] "AL627309.1" "RP11-206L10.2" "LINC00115" "NOC2L" ...
.. .. .. .. .. .. .. ..$ : chr [1:4] "counts.CTRL" "counts.STIM" "data.CTRL" "data.STIM"
.. .. .. .. .. ..$ dim : int [1:2] 14053 4
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:14053] "AL627309.1" "RP11-206L10.2" "LINC00115" "NOC2L" ...
.. .. .. .. .. .. ..$ : chr [1:4] "counts.CTRL" "counts.STIM" "data.CTRL" "data.STIM"
.. .. .. ..@ default : int 2
.. .. .. ..@ assay.orig: chr(0)
.. .. .. ..@ meta.data :'data.frame': 14053 obs. of 0 variables
.. .. .. ..@ misc : list()
.. .. .. ..@ key : chr "rna_"
看完了官网的数据结构，我想我们也需要一个类似的数据结构，我们尝试构建一个吧。我找了一些笔记后，发现这两篇笔记seurat v5直播，一键完成五种数据整合：harmony，CCA，RPCA,FastMNN,scVI、单细胞专题 | 14 Seurat v5版本融合多样本的单细胞数据以及一些报错问题中提到的merge这个函数可以实现这一目的，我们操作一下看看。
# 将不同的对象merge进一个对象中
> object.combined <- merge(Subset_Cells.list[[1]], y = Subset_Cells.list[-1], add.cell.ids = names(Subset_Cells.list), project = "neuron")
# 查看一下数据结构
> str(object.combined)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
.. .. .. ..@ layers :List of 5
.. .. .. .. ..$ counts.mouse_neuron_1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:17520797] 0 7 10 11 13 16 22 28 31 36 ...
.. .. .. .. .. .. ..@ p : int [1:4829] 0 5132 9976 10443 14442 19667 24059 27643 32582 37329 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 23476 4828
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:17520797] 15 2 3 1 2 8 17 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.mouse_neuron_2:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:17579967] 0 6 9 10 12 15 21 27 30 35 ...
.. .. .. .. .. .. ..@ p : int [1:4635] 0 5132 9974 13920 19103 23468 27006 31962 36700 41116 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 23449 4634
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:17579967] 15 3 2 1 1 9 15 2 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.mouse_neuron_3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:24534828] 0 5 18 35 41 52 58 84 90 92 ...
.. .. .. .. .. .. ..@ p : int [1:6590] 0 2889 7701 10902 14465 19720 24087 27646 30798 35197 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 22697 6589
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:24534828] 6 1 1 1 4 6 2 4 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.mouse_neuron_4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:14452170] 0 1 5 8 11 14 18 20 21 28 ...
.. .. .. .. .. .. ..@ p : int [1:3742] 0 4155 8053 13249 17362 21735 25762 28913 34211 34878 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 22104 3741
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:14452170] 17 1 1 1 1 1 3 4 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. .. ..$ counts.mouse_neuron_5:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:12662480] 0 2 5 8 11 14 29 33 35 40 ...
.. .. .. .. .. .. ..@ p : int [1:3295] 0 4431 9270 11211 15271 19271 24435 28371 31396 33928 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 22275 3294
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:12662480] 6 2 2 2 2 1 2 1 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. ..@ cells :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:23086, 1:5] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:23086] "AAACCCACAAGAATAC-1_1" "AAACCCACAATAAGGT-1_1" "AAACCCACACCATAAC-1_1" "AAACCCACAGCCCACA-1_1" ...
.. .. .. .. .. .. .. ..$ : chr [1:5] "counts.mouse_neuron_1" "counts.mouse_neuron_2" "counts.mouse_neuron_3" "counts.mouse_neuron_4" ...
.. .. .. .. .. ..$ dim : int [1:2] 23086 5
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:23086] "AAACCCACAAGAATAC-1_1" "AAACCCACAATAAGGT-1_1" "AAACCCACACCATAAC-1_1" "AAACCCACAGCCCACA-1_1" ...
.. .. .. .. .. .. ..$ : chr [1:5] "counts.mouse_neuron_1" "counts.mouse_neuron_2" "counts.mouse_neuron_3" "counts.mouse_neuron_4" ...
.. .. .. ..@ features :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:24771, 1:5] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:24771] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
.. .. .. .. .. .. .. ..$ : chr [1:5] "counts.mouse_neuron_1" "counts.mouse_neuron_2" "counts.mouse_neuron_3" "counts.mouse_neuron_4" ...
.. .. .. .. .. ..$ dim : int [1:2] 24771 5
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:24771] "Xkr4" "Gm1992" "Gm19938" "Gm37381" ...
.. .. .. .. .. .. ..$ : chr [1:5] "counts.mouse_neuron_1" "counts.mouse_neuron_2" "counts.mouse_neuron_3" "counts.mouse_neuron_4" ...
.. .. .. ..@ default : int 5
.. .. .. ..@ assay.orig: chr(0)
.. .. .. ..@ meta.data :'data.frame': 24771 obs. of 0 variables
.. .. .. ..@ misc : list()
.. .. .. ..@ key : chr "rna_"
刚好是我们需要的数据格式，非常好，我们接着就开始按照官网的方式开始整合吧。
3.1 锚点整合
Seurat.v5也是两种标准化方法，我们逐个展示。
3.1.1 NormalizeData标准化 + 锚点整合CCA
# 首先是seurat标准流程
object.combined <- NormalizeData(object.combined)
object.combined <- FindVariableFeatures(object.combined)
object.combined <- ScaleData(object.combined)
object.combined <- RunPCA(object.combined)
# 接着，运行CCA整合，这一步是v5特有的。
object.combined <- IntegrateLayers(object = object.combined, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca", verbose = FALSE)
# 使用JoinLayers()函数将矩阵融合，这一步也是v5特有的。
# 这一步拆开看就是把object.combined@assays$RNA@layers$counts.mouse_neuron_1(2/3/4/5)
# 融合成一个object.combined@assays$RNA@layers$counts
object.combined[["RNA"]] <- JoinLayers(object.combined[["RNA"]])
# 聚类、分群
object.combined <- FindNeighbors(object.combined, reduction = "integrated.cca", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# UMAP降维
object.combined <- RunUMAP(object.combined, dims = 1:30, reduction = "integrated.cca")
# 画图展示整合效果
pdf("Dimplot_group_Normal_CCA.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_Normal_CCA.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_Normal_CCA.rds")
整合结果如下：
3.1.1 SCTransform标准化 + 锚点整合CCA
# SCT标准化
object.combined <- SCTransform(object.combined)
object.combined <- RunPCA(object.combined)
# 运行CCA整合
object.combined <- IntegrateLayers(object.combined, method = CCAIntegration, orig.reduction = "SCT", verbose = FALSE)
# 聚类、分群
object.combined <- FindNeighbors(object.combined, reduction = "integrated.dr", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
# UMAP降维
object.combined <- RunUMAP(object.combined, dims = 1:30, reduction = "integrated.dr")
# 画图展示整合效果
pdf("Dimplot_group_SCT_CCA.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "umap", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "umap", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_SCT_CCA.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "umap", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_SCT_CCA.rds")
整合结果如下：
3.2 Harmony整合
安装：
这个是可以用mamba直接安装的，并且安装出来时1.2.0版本的，我看官网都没有这个版本，不太清楚原因，但是能用，命令如下：
mamba install bioconda::r-harmony
3.2.1 NormalizeData标准化 + Harmony整合
# 导入harmony包
library(harmony)
# seurat标准流程
object.combined <- NormalizeData(object.combined)
object.combined <- FindVariableFeatures(object.combined)
object.combined <- ScaleData(object.combined)
object.combined <- RunPCA(object.combined)
# 运行harmony整合
object.combined <- IntegrateLayers(object.combined, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony", verbose = FALSE)
# 使用JoinLayers()函数将矩阵融合
object.combined[["RNA"]] <- JoinLayers(object.combined[["RNA"]])
# 聚类、分群、降维
object.combined <- FindNeighbors(object.combined, reduction = "harmony", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
object.combined <- RunUMAP(object.combined, reduction = "harmony", dims = 1:30, reduction.name = "harmony")
# 画图展示整合效果
pdf("Dimplot_group_Normal_Harmony.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "harmony", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "harmony", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_Normal_Harmony.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "harmony", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_Normal_Harmony.rds")
结果如下：
3.2.1 SCTransform标准化 + Harmony整合
# 导入harmony包
library(harmony)
# SCT标准化
object.combined <- SCTransform(object.combined)
object.combined <- RunPCA(object.combined)
# 运行harmony整合
object.combined <- IntegrateLayers(object.combined, method = HarmonyIntegration, normalization.method = "SCT", orig.reduction = "pca", new.reduction = "harmony",verbose = FALSE)
# 聚类、分群、降维
object.combined <- FindNeighbors(object.combined, reduction = "harmony", dims = 1:30)
object.combined <- FindClusters(object.combined, resolution = 0.5)
object.combined <- RunUMAP(object.combined, dims = 1:30, reduction = "harmony", reduction.name = "harmony")
# 画图展示整合效果
pdf("Dimplot_group_SCT_Harmony.pdf", width = 20, height = 10)
p1 <- DimPlot(object.combined, reduction = "harmony", group.by = "orig.ident", pt.size = 1.0)
p2 <- DimPlot(object.combined, reduction = "harmony", label = TRUE, pt.size = 1.0)
p1 + p2
dev.off()
pdf("Dimplot_split_SCT_Harmony.pdf", width = 50, height = 10)
DimPlot(object.combined, reduction = "harmony", split.by = "orig.ident", ncol = 5, pt.size = 1.0)
dev.off()
# 存储结果
saveRDS(object.combined, file = "integration_SCT_Harmony.rds")
结果如下：
最后呢，我拿这个seurat.v5的 SCT + Harmony 的整合结果画了一些脑中经典marker基因的表达情况，基本上分布都挺特异的，进一步说明整合效果应该是挺好的，不过一些细节肯定是还需要再把控一下的。
OK! 那么到这里本次的笔记基本上也就结束了，本次笔记原理的部分介绍得可能相对较少，主要是为了理清楚流程，在拿到比对定量的矩阵数据后我们应该做哪些工作，以得到满足seurat标准流程分析的结果。在本教程中，我们为了展示整合后较好地消除了批次效果，用了固定参数对整合后的数据集进行了简单分析，大家在分析自己的数据时需要结合肘形图和clustertree的细胞分群树来确定相关参数，我们这里就不做扩展了。
另外呢，在这里提前做一个预告，为了方便大家的使用呢，我们将会在下一期更新一个合集系列，将我们的单细胞推文排排顺序（虽然只有六七篇，哈哈哈）。同时呢，我们也将会理清楚我们已有的笔记，将他们封装整理成系统的脚本，大家按照要求配置好环境后，只需要向脚本传参即可实现相关的操作，大家可以尽情期待一下。
最后的最后，如果大家觉得我们的推文对你的单细胞学习有帮助的话，欢迎大家点赞、关注、转发，让更多的人能够看到！
文末碎碎念
那今天的分享就到这里啦！我们下期再见哟！
最后顺便给自己推荐一下嘿嘿嘿！
如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！
蟹蟹你们的喜欢和支持！！！
扫一扫
1117
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
