单细胞基地 | 单细胞多样本数据整合 | 两万字大干货

文摘   2024-07-18 18:11   上海  

时隔三个月,我带着新的单细胞笔记又来了。这是单细胞基础分析流程的最后一个板块,把这一部分搞懂之后,你就可以独立地完成一套标准的单细胞基础分析了,感兴趣的同学可以学起来呀!(当然,我们的笔记顺序并不是按照正常流程顺序来写的,所以我们后续也会专门做一篇整理性质的推送,为大家理清学习顺序的,这个大家不用担心!)

话不多说(那个那个,我,啊,这里的我是指小蛮要,咱们这篇文章是我们超级棒棒的小杨写的,酷毙了!那我就按照惯例多说几句好不啦哈哈哈哈哈,谢谢大家!),我们直接进入正题!

为什么要做多样本整合?

上一篇笔记我们讲了单细胞测序下机数据的比对质控流程,获取了相应的结果并对其中的核心内容进行了解读,详见单细胞基地 | 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

代码解析:官网给的代码(https://github.com/constantAmateur/SoupX#quickstart)在实际运行过程中,直接运行load10X并不顺利,所以我结合了github官网上这个Issues45(https://github.com/constantAmateur/SoupX/issues/45)下面的作者回复,对代码进行了修改和整理,大家直接看下面的内容即可:

# 加载需要的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,如下所示:

好的,我们继续回到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

代码解析:首先呢,我们先要输入数据,关于输入的数据,官网提了两点要求:

  1. Do not apply DoubletFinder to aggregated scRNA-seq data representing multiple distinct samples (e.g., multiple 10X lanes).也就是说,DoubletFinder是对单个样本处理的,不要对多样本数据进行操作。
  2. 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()

DoubletFinder2
# 这样,我们就可以去除双胞,只取单胞了
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部分(https://github.com/chris-mcginnis-ucsf/DoubletFinder?tab=readme-ov-file#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 长什么样。

很明显可以看到,这个列表里面存储了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部分的整合呢,官方把目前主流的一些方法也集成了进去,这些在官网上也有详细的说明,👉指路(https://satijalab.org/seurat/articles/seurat5_integration#perform-streamlined-one-line-integrative-analysis),概括起来说呢,基本上是以下五种: 

① 基于锚点的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里面是什么信息。

很明显这个数据集由两个样本构成的,分别是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的细胞分群树来确定相关参数,我们这里就不做扩展了。

另外呢,在这里提前做一个预告,为了方便大家的使用呢,我们将会在下一期更新一个合集系列,将我们的单细胞推文排排顺序(虽然只有六七篇,哈哈哈)。同时呢,我们也将会理清楚我们已有的笔记,将他们封装整理成系统的脚本,大家按照要求配置好环境后,只需要向脚本传参即可实现相关的操作,大家可以尽情期待一下

最后的最后,如果大家觉得我们的推文对你的单细胞学习有帮助的话,欢迎大家点赞、关注、转发,让更多的人能够看到!

文末碎碎念

那今天的分享就到这里啦!我们下期再见哟!

最后顺便给自己推荐一下嘿嘿嘿!

如果我的分享对你有用的话,欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴!这可是我的第一原动力!

蟹蟹你们的喜欢和支持!!!

生信小白要知道
主打小白保姆级教程,因为自己淋过雨,所以想给大家撑把伞! 记录从小白到现在小灰的过程,希望以后可以成为小黑!
 最新文章