单细胞测序简介

单细胞测序是指在单个细胞水平上进行测序，单细胞转录组测序（single cell RNA Seq，scENA-seq）是指对于单个细胞水平上将mRNA反转录扩增后进行高通量测序的技术。

单细胞测序通过在单个细胞水平上进行测序，解决了用组织样本无法获得不同细胞间的异质性信息或样本量太少无法进行常规测序的难题，为研究人员研究单个细胞的行为，机制等提供了新的方向。

下面展示了单细胞测序技术的特点：

单细胞测序可以获得单个细胞基因的表达，与普通转录组测序相比，普通转录组测序是把组织的总RNA提取出来，进行建库检测，获得的是一个平均的基因表达；单细胞测序会先将单个细胞从组织中解离出来，测量每个细胞的基因表达情况；

获得单细胞测序数据后，我们的任务是对其做数据挖掘工作，通常，单细胞测序数据挖掘流程通常为：

单细胞测序数据下载；
质控分析和数据过滤；
对数据降维：线性降维保留主成分（PCA），非线性降维t-SNE保留更具代表性的信息；
对降维后的数据进行聚类，将单个细胞归属到各自的簇内；
为每一个簇找到具有显著性质的基因（Marker基因）；
根据经验对细胞进行注释；
细胞轨迹分析：分析哪些细胞先出现，哪些细胞后出现，得到细胞的分化过程；

数据预处理

我将分析 10X Genomics 平台免费提供的外周血单核细胞 (PBMC3k) 数据集。数据为 2700 个单细胞在 Illumina NextSeq 500 上测序；

数据集下载地址：pbmc3k_filtered_gene_bc_matrices

下载后需解压得到以下目录：

|- data|- filtered_gene_bc_matrices|- hg19|- barcodes.tsv|- genes.tsv|- matrix.mtx

导入相关包：

import numpy as np
import pandas as pd
import scanpy as sc

进行一些设置：

sc.settings.verbosity=3 # verbosity: errors (0), warnings (1), info (2), hints(提示) (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80,facecolor='white')
result_file='./write/pbmc3k.h5ad' # 用于缓存分析结果

读入PBMC3k数据集：

adata=sc.read_10x_mtx('./data/filtered_gene_bc_matrices/hg19/', # .mtx文件的所在目录var_names='gene_symbols', # 使用gene_symbols作为AnnData的特征名称cache=True) # 写入缓存便于快速读写adata.var_names_make_unique() # 如果在 sc.read_10x_mtx 中使用 var_names='gene_ids'，这是不必要的
print(adata)
"""
AnnData object with n_obs × n_vars = 2700 × 32738var: 'gene_ids'
"""

针对AnnData的特征（基因），统计每个基因在所有细胞中的分布，并绘制这些基因分布的箱型图：

sc.pl.highest_expr_genes(adata,n_top=20)

使用两个工具：sc.pp.filter_cells进行细胞的过滤，sc.pp.filter_genes进行基因的过滤。这两个工具的说明可回顾Scanpy(一)AnnData数据结构与一些API用法介绍；过滤过程如下：

sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=3)

上面过程为基因过滤（Basic filtering），下面我需要收集关于线粒体基因（mitochondrial genes）的信息，这些信息对质量控制（Quality Control，QC）很重要。

我们需要过滤包含线粒体基因和基因过多的细胞，因为：

检测出高比例的线粒体基因，表明细胞的测量质量较差；
表达基因过多可能是由于一个测量油滴包裹了多个细胞，从而检测出比正常情况要多的基因数，因此要过滤这些细胞。

首先，查看数据的特征量（AnnData.var）：

print(adata.var)
"""gene_ids  n_cells
AL627309.1     ENSG00000237683        9
...                        ...      ...
SRSF10-1       ENSG00000215699       69
"""

对于var这个Dataframe的第0列，也就是index（基因名称），也被称为var_names，现在我们要为var中的基因做标记，判断哪些基因是线粒体相关的：

# startswith()方法用于检查字符串是否是以指定子字符串开头, 如果是则返回 True, 否则返回 False
adata.var['mt']=adata.var_names.str.startswith('MT-') # var新增一列'mt', 这是用于标记基因是否是线粒体相关的一个mask

使用 pp.calculate_qc_metrics，我们可以计算许多QC指标：

sc.pp.calculate_qc_metrics(adata,qc_vars=['mt'],percent_top=None,log1p=False,inplace=True)

扩展部分：QC与sc.pp.calculate_qc_metrics

在单细胞分析中，我们通常要对某些基因计算数据指标，根据指标过滤数据，以便保证后续计算的质量，这就是质量控制；

sc.pp.calculate_qc_metrics是用于计算QC指标的函数：

sc.pp.calculate_qc_metrics(data, qc_vars=(), percent_top=(50, 100, 200, 500), inplace=False, log1p=True)

参数说明：

data：Anndata；
qc_vars：用于标识要控制的特征（基因），布尔型元素，用于作为mask使用；
percent_top：计算与常出现基因的比，percent_top=[50] 计算与第 50 个最常出现基因的比例，None则不计算；
inplace：决定是否将计算指标添加到var和obs中；
log1p：设置为False可以跳过转换到log1p空间的过程；log1p即log(1+number)，用于压缩数据并确保结果是一个正数；

实例如下：

import seaborn as sns
import matplotlib.pyplot as pltpbmc=sc.datasets.pbmc3k() # 自动下载到./data/pbmc3k_raw.h5ad# 标记线粒体基因, 并保存为 mask
pbmc.var['mito']=pbmc.var_names.str.startswith('MT-')
print(pbmc.var)
"""gene_ids   mito
index
MIR1302-10     ENSG00000243485  False
...                        ...    ...
AC002321.1     ENSG00000215611  False
"""
print(pbmc.var.columns)
"""
Index(['gene_ids', 'mito'], dtype='object')
"""
print(pbmc.obs.columns)
"""
Index([], dtype='object')此时obs只有index, obs的index为细胞的名称
"""# 计算线粒体基因的 QC 指标
sc.pp.calculate_qc_metrics(pbmc,qc_vars=['mito'],inplace=True)
print(pbmc.var.columns) # 关于基因轴上的QC指标
"""
Index(['gene_ids', 'mito', 'n_cells_by_counts', 'mean_counts','log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts','log1p_total_counts'],dtype='object')
"""
print(pbmc.obs.columns) # 关于细胞轴上的QC指标
"""
Index(['n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts','log1p_total_counts', 'pct_counts_in_top_50_genes','pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes','pct_counts_in_top_500_genes', 'total_counts_mito','log1p_total_counts_mito', 'pct_counts_mito'],dtype='object')
"""# 对obs中的计算指标可视化
sns.jointplot(data=pbmc.obs,x="log1p_total_counts",y="log1p_total_counts_mito",kind="hex")
plt.show()
sns.histplot(pbmc.obs["pct_counts_mito"])
plt.show()

total_counts记录了对每个细胞，所包含基因的总数；

total_counts_mito记录了对每个细胞，所包含控制基因（qc_vars中列出的基因）的总数；

pct_counts_mito记录了对每个细胞，含线粒体的基因总数；（其中 ‘mito’ 来自var中的index）

现在，我们用小提琴图可视化两个QC指标：

# jitter 表示抖动
sc.pl.violin(adata, ['pct_counts_mt'], jitter=0.4)
sc.pl.violin(adata, ['n_genes_by_counts'], jitter=0.4)

从上面两个指标可以看出，大部分细胞的线粒体基因数在5以下，大部分细胞的基因总数在2500以下，所以我们要过滤掉那些异常的细胞：

# 去除线粒体基因过多或基因总数过多的细胞
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]

然后对数据标准化，并转换到 log1p 空间：

# 标准化
sc.pp.normalize_total(adata, target_sum=1e4)# 压缩到logp1空间
sc.pp.log1p(adata)

确定高变基因并过滤数据：

# 确定高变基因
sc.pp.highly_variable_genes(adata)
"""
extracting highly variable genesfinished (0:00:00)
--> added'highly_variable', boolean vector (adata.var) # 布尔型mask'means', float vector (adata.var)'dispersions', float vector (adata.var)'dispersions_norm', float vector (adata.var)
"""# 根据高变基因过滤数据
adata=adata[:,adata.var.highly_variable]

对数据进行回归并校正数据：

"""
sc.pp.regress_out(data, keys)
keys:要回归的data.obs中的数据列
"""
# 回归 adata.obs 中的列 (columns)
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])

将数据adata.X缩放到单位方差和零均值，对于缩放后的数据，在值为10处截断：

sc.pp.scale(adata, max_value=10)

降维与聚类

我们可以用PCA进行降维：

# PCA降维并可视化
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color='CST3')# 写入缓存
adata.write(result_file)

sc.tl.pca
使用PCA对数据进行变换；

sc.pl.pca

sc.pl.pca(adata,color=None)

注意参数color，我们可以传递adata中的'ann1' or ['ann1', 'ann2']，这些ann可以是基因，可以是聚类结果，所以我们可以根据基因标记颜色，也可以根据聚类算法结果标记颜色；

tl用于做变换，pl用于可视化（只取前2个维度），调用pl前需要先找到tl计算得到的数据空间；

比如上图，PCA将高维数据adata.X聚类到2维空间后得到的只是一些平面下的散点，但我们可以根据每个散点（细胞）中包含基因CST3的数量为散点标记颜色。

我们也可以用UMAP进行降维并可视化，我们现在使用的数据是经过PCA变换后的数据（保留了所有成分），使用这个数据并用UMAP方法做变换：

# UMAP 降维得到 embedding 并可视化
sc.pp.neighbors(adata,n_neighbors=10,n_pcs=40)
sc.tl.umap(adata)
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

现在，利用经过UMAP变换后的数据，使用 leiden 方法进行聚类：

# 用 leiden 方法聚类
sc.tl.leiden(adata)
"""
running Leiden clusteringfinished: found 8 clusters and added'leiden', the cluster labels (adata.obs, categorical) (0:00:00)
"""
print(adata.obs)
"""n_genes  n_genes_by_counts  ...  pct_counts_mt  leiden
AAACATACAACCAC-1      781                779  ...       3.017776       0
AAACATTGAGCTAC-1     1352               1352  ...       3.793596       2
...                   ...                ...  ...            ...     ...
TTTGCATGCCTCAC-1      724                723  ...       0.806452       0
"""# 用 UMAP 可视化聚类结果
sc.pl.umap(adata, color='leiden')# 写入结果
adata.write(result_file)

找到Marker基因

让我们计算每个簇中高度差异基因的排名。最简单和最快的方法是 t 检验。

可以使用：

scanpy.tl.rank_genes_groups(adata, groupby, groups='all', reference='rest', n_genes=None, method=None)groupby：分组的依据groups：组的子集，例如 ['g1', 'g2', 'g3']，对于所有组，为'all'（默认）。reference：如果是 'rest'，则将每个组与该组其余部分的并集进行比较。如果是组标识符，则相对于reference进行比较。n_genes：出现在返回表中的基因数量。默认为所有基因。method：计算方法

为之前leiden聚类的每个簇计算基因排名，用当前簇与其他簇进行比较获得排名：

sc.tl.rank_genes_groups(adata,'leiden',method='t-test')
# 每个簇只可视化前25个基因
sc.pl.rank_genes_groups(adata,n_genes=25,sharey=False)

也可以用别的方法做检验：

sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

不同的方法得到的基因排名结果是不同的；

在dataframe中显示每个集群 0, 1, …, 7 中排名最高的 10 个基因：

pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(10)

得到一张包含基因排名和分数的表：

result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame({group + '_' + key[:1]: result[key][group]for group in groups for key in ['names', 'pvals']}).head(5)

也可以只比较单个簇：

sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)

我们可以针对某些基因，分析它们在各个簇中的分布：

sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')

我们可以给细胞的簇命名，注意列表内的细胞名顺序要与簇顺序对应：

new_cluster_names = ['CD4 T', 'CD14 Monocytes','B', 'CD8 T','NK', 'FCGR3A Monocytes','Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='mark cell type',frameon=False)

基于注释后的细胞，再可视化marker基因，假设经过前面基因排名的选择，我们已经找到了marker基因：

marker_genes = ['CD79A', 'MS4A1','LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1','FCGR3A', 'FCER1A', 'CST3', 'PPBP']

可视化如下，这是这些基因在各类细胞中的分布：

sc.pl.dotplot(adata, marker_genes, groupby='leiden')

也可以用紧凑的小提琴图集合表示：

在此分析过程中，AnnData 积累了以下注释：

adata
"""
AnnData object with n_obs × n_vars = 2638 × 1838obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'uns: 'log1p', 'hvg', 'neighbors', 'umap', 'leiden', 'leiden_colors', 'rank_genes_groups'obsm: 'X_pca', 'X_umap'obsp: 'distances', 'connectivities'
"""

Scanpy(二)对PBMC3k聚类相关推荐

百度新闻的索引机制(二)：智能聚类
百度新闻的索引机制(二):智能聚类 http://net.chinabyte.com/377/2520877.shtml 转载于:https://www.cnblogs.com/cy163/archi ...
文本聚类（二）—— KMeans 聚类
目录二.KMeans 聚类 2.1 加载数据集 2.2 数据清洗 2.3 文本向量化 2.4 文本聚类 2.5 关键词展示 2.6 判定最佳聚类数参考文档在第一篇内容中,我们介绍了 LDA 主题 ...
K-means 算法实现二维数据聚类
所谓聚类分析,就是给定一个元素集合D,其中每个元素具有n个观测属性,对这些属性使用某种算法将D划分成K个子集,要求每个子集内部的元素之间相似度尽可能高,而不同子集的元素相似度尽可能低.聚类分析是一种无 ...
机器学习算法（十二）：聚类
目录 1 K的选择 1.1 肘部法则(Elbow method) 1.2 目标法则 1.3 间隔统计量 Gap Statistic 1.4 关于K值选择的改进算法--ISODATA算法 2 聚类算法的 ...
Scanpy(二)可视化函数
目录 embedding的散点图可视化基因表达量与其他变量用已知marker基因识别clusters dotplot violin plot stacked-violin plot matrixp ...
文本聚类分析算法_Kmeans 算法实现二维数据聚类
所谓聚类分析,就是给定一个元素集合D,其中每个元素具有n个观测属性,对这些属性使用某种算法将D划分成K个子集,要求每个子集内部的元素之间相似度尽可能高,而不同子集的元素相似度尽可能低.聚类分析是一种无 ...
机器学习算法（十二）：聚类（2）层次聚类 Hierarchical Clustering
目录 1 层次聚类 1.1 层次聚类的原理 1.2 两个组合数据点间的距离: 2 自底向上的合并算法 2.1 AGNES算法 (AGglomerative NESting) 2.1.1 原理 2.1. ...
新手探索NLP（十二）——文本聚类
简介聚类又称群分析,是数据挖掘的一种重要的思想,聚类(Cluster)分析是由若干模式(Pattern)组成的,通常,模式是一个度量(Measurement)的向量,或者是多维空间中的一个点.聚类分 ...
机器学习笔记（十二）：聚类
目录 1)Unsupervised learning introduction 2)K-means algorithm 3)Optimization objective 4)Random initia ...

Scanpy(二)对PBMC3k聚类

目录

单细胞测序简介

数据预处理

降维与聚类

找到Marker基因

Scanpy(二)对PBMC3k聚类相关推荐

最新文章

热门文章