抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

单细胞转录组分析

单细胞转录组学习笔记

发布于：2023年6月13日

字数：2.6k 字

时长：10 分钟
次浏览

最近因为学校网络信息中心升级防火墙，导致集群无法访问公网，我搭建的反向代理服务器也暂时无法使用（真的想吐槽学校的网络管理员，三个星期了还没能解决集群的网络问题）……近期要进行中期答辩，先用向日葵远程一下实验室闲置的电脑应急，顺便把最近的学习笔记补上。

这篇笔记主要记录下单细胞组学的学习，以及单细胞转录组（Single-cell RNA-sequencing，scRNA-seq）的分析流程。

1. 单细胞测序发展

我们做植物高通量测序做的最多的是RNA-seq，比较不同组织、不同时期或者不同处理下同一个组织基因的差异表达。我们做RNA-seq是建立在对一个组织的细胞进行测序的基础上，而组织是几类细胞的集合，因此我们得到的基因表达量是所有细胞的平均值。

单细胞测序可以在单个细胞的水平上构建细胞图谱，让我们更深入了解植物组织的细胞类型，获取每个细胞的转录本信息，研究细胞的发育动态（比如细胞如何进行分化）等等。

Nature杂志每年都会总结每个领域最有价值的年度技术，2018年是单细胞转录组技术，2019年是单细胞多组学技术，2020年是空间转录组技术，可以看出带有单细胞和空间分辨率转录组技术应用的重要性。

上图是单细胞技术在植物中的研究历程，最后一个2022年的Stereo-seq也就是华大自研的时空组学技术，据华大发表在cell上的文章Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays描述，Stereo-seq的技术参数优于当前其他空间转录组技术，对空间异质性的描述更为灵敏和直观。

对于Stereo-seq我个人的理解是补上了部分空间转录组测序无法做到单个细胞分辨率的短板，真正做到对动植物的组织或者器官中的所有细胞进行转录组信息分析，鉴定不同类型细胞的差异基因表达，构建空间图谱来分辨不同细胞的发育轨迹。感兴趣可以看下面华大的这篇scStereo-seq技术用于拟南芥的文章：The single-cell stereo-seq reveals region-specific cell subtypes and transcriptome profiling in Arabidopsis leaves - ScienceDirect

2. 单细胞测序平台

  • 10x Genomics Chromium 微流控芯片技术获得单细胞反应体系，并在传统文库构建的基础上引入标签，通过追溯标签序列将众多mRNA、表面蛋白信息定位回原来的单个细胞。捕获效率最高达65%

  • BD Rhapsody™ Single-Cell Analysis System 能为单细胞中每个转录本标记特异性分子标签，实现单细胞水平上基因表达谱的绝对定量。捕获效率最高达80%

  • Illumina® Bio-Rad® Single-Cell Sequencing Solution 捕获效率低，仅为3%，但测序成本相对较低

  • ICELL8 Single-Cell System 捕获效率为30%，成本相对较低

  • C1™ 单细胞全自动制备系统 通量低、成本高、周期慢、对试验人员要求高，操作较繁琐和困难

3. 单细胞测序分析流程

质控部分和RNA-seq没有区别，去掉低质量的读序和测序接头。10X Genomics平台测的单细胞数据需要通过Cell Ranger回比到参考基因组，示例用法如下：

1
2
3
4
5
6
7
8
 cellranger count --id=sample345
    --transcriptome=/opt/refdata-cellranger-GRCh38-3.0.0
    --fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path
    --sample=mysample
# --id  文件名
# --transcriptome 参考基因组
# --fastqs 转录组数据所在的路径
# --sample 指定要使用的样本

这里主要记录下细胞过滤到亚群定义和标记基因筛选所要用到的软件Seurat4。因为我自己没有这方面的实验数据，仅仅只是尝试跑下流程，以下分析流程和代码来自：

单细胞转录组|Seurat 4.0 使用指南 - 知乎 (zhihu.com)

《Seurat 4 R包源码解析》 总目录 | 单细胞转录组分析标准流程 - 知乎 (zhihu.com)

示例数据来自10X Genomics免费提供的外周血单核细胞(PBMC)数据集：
https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

将示例数据上传至集群，解压，得到如下文件：

1
2
3
├── barcodes.tsv（10X Genomics测序的与细胞相关的barcode信息，用于标识细胞）
├── genes.tsv（基因信息，包括了基因组数据库人类登记号和基因名称）
└── matrix.mtx（基因表达量矩阵）
  • 这三个文件实际上要通过软件Cell Ranger对10X Genomics单细胞测序的下机数据进行比对基因组，统计捕获的细胞数、测序量得到。详细可以参考官方Cell Ranger软件的用法和结果分析，这里用的是官方整理后的数据，只进行下一步Seurat分析的演示。

以所在目录为工作目录，在linux中键入R进入交互命令行，运行Seurat。

3.1 导入数据和初步过滤

1
2
3
4
5
6
7
8
9
10
11
12
# 安装和加载Seurat和dplyr（用于数据清洗）R包
BiocManager::install("Seurat") 
BiocManager::install("dplyr")
library(dplyr)
library(Seurat)

# 读入数据并初步筛选
pbmc.data <- Read10X(data.dir="/path/to/yourfile")

# 创建Seurat对象，过滤检测少于200个基因的细胞，和少于3个细胞检测出的基因
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

可以看到初步过滤后，现在的一个数据集包含了2700个细胞的13714个基因。

3.2 细胞过滤

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 新增一列percent.mt用于统计映射到线粒体基因组的reads百分比
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 绘制细胞的三种特征，nFeature_RNA（每个细胞测到的特异基因数目）、nCount_RNA（每个细胞测到所有基因的表达量之和）、percent.mt的小提琴图
png("fig01_cell_feature.png")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()

# 细胞特征间的相关性绘图
png("fig02_cell_feature_correlation.png", width=1000, height=400)  
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
dev.off()

# 根据特殊基因的数目以及线粒体基因比例过滤细胞，这里选取200 < nFeature_RNA < 2500 和percent.mt < 5的数据
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) 
pbmc

fig01_cell_feature.png如下所示：

这里可以看到有的细胞检测到的特异基因数特别多，可能是因为多个细胞被同时捕获了，而特异基因数特别少的可能是低质量细胞。而线粒体基因比例特别高的细胞，有可能是一些快要死亡的细胞。这些低质量的细胞需要在这一步被过滤。

fig02_cell_feature_correlation.png如下所示：

这两个图分别表示了每个细胞测到所有基因的表达量之和与线粒体基因组百分比呈负相关，每个细胞测到所有基因的表达量之和与每个细胞测到的特异基因数目呈正相关。

这一步过滤后数据集剩下2638个细胞，共13714个基因。

3.3 细胞群体聚类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 表达水平标准化
## normalization.method  标准化方法，默认为“LogNormalize”
## scale.factor  比例因子，用以计算标准化值，默认为“10000”
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# 鉴定mark基因（也就是特征基因）
## selection.method  mark基因筛选方法，“vst”为通过方差与均值比进行筛选
## nfeatures  要筛选的mark基因的数目
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(pbmc), 10)

# 展示mark基因
png("fig03_10marker_genes.png",width=1000,height=400)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
dev.off()

fig03_10marker_genes.png如下所示：

这里筛选出2000个mark基因用于后续的下游分析，并且展示区分能力最强的前10个mark基因

1
2
3
4
5
6
7
8
9
10
11
12
# Mark基因权重标准化
## ScaleData 缩放基因的表达，给予基因同等的权重，使高表达基因不占据主导地位
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes) 

# PCA主成分分析
## RunPCA 计算特征值，在细胞间具有高度表达差异的基因有助于区分不同类型的细胞
## 采用DimPlot方法可视化
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
png("fig04_PCA.png")
DimPlot(pbmc, reduction = "pca")
dev.off()

fig04_PCA.png如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
# 评估主成分维度
## num.replicate  重采样次数
## dims  维度范围
pbmc <- JackStraw(pbmc, num.replicate = 100)  
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
## 评估主成分维度的方法1
png("fig05_PC_importance_01.png",width=700,height=400)
JackStrawPlot(pbmc, dims = 1:15)
dev.off()
## 评估主成分维度的方法2
png("fig05_PC_importance_02.png",width=700,height=400)
ElbowPlot(pbmc)
dev.off()

R包 Seurat 函数 JackStraw 、ScoreJackStraw进行重采样测试，评估用以进行细胞聚类的主成分维度。

fig05_PC_importance_01.png如下所示：

fig05_PC_importance_02.png如下所示：

上图是比较不同主成分（PC）的P值分布与均匀分布（虚线）；显著的主成分具有较低的P值（位于虚线上方）；似乎10-12个主成分后，主成分的重要性急剧下降。

下图是肘部法则以确定最佳主成分数量，9-10个主成分附近有个拐点，表明大部分真实信号是在前10个pc中捕获的。

所以这里选择10个主成分用于后续分析。

1
2
3
4
5
6
7
8
# 细胞群体聚簇
pbmc <- FindNeighbors(pbmc, dims = 1:10) 
pbmc <- FindClusters(pbmc, resolution = 0.5) 
## 绘图
pbmc <- RunUMAP(pbmc, dims = 1:10) 
png("fig06_cluster_UMAP.png",width=700,height=400) 
DimPlot(pbmc, reduction = "umap")
dev.off()

fig06_cluster_UMAP.png如下所示

可以看到，根据mark基因的表达，所有细胞最终被分为了8种类型。

做完以上分析之后，还可以根据自己的研究目标新型不同的下游分析，具体来说可以筛选亚群的标记基因（与其他类别细胞的差异表达基因）、进行细胞的发育轨迹分析、做不同品种间细胞类型的比较等等。

本文由Phantom创作，转载请注明出处

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：http://www.shelven.com/2023/06/13/a.html

更新于：2023年6月16日

Seurat

欢迎小伙伴们留言评论~