用VSCode Jupyter 学习Scanpy——预处理与分群

作者: 米妮爱分享 | 来源:发表于2020-12-19 18:29 被阅读0次

预处理与分群的步骤大体和Ｒ包　Ｓｅｕｒａｔ（https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html）类似哦，
数据的下载和目录的创建可以参考一下代码
[1]

# !mkdir data
# !wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz
# !cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz
# !mkdir write

引入包
[2]

import pandas as pd
import scanpy as sc
import scanpy as sc

设置和查看当前环境
[3]

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')

image.png

指定结果路径
[4]

results_file = 'write/pbmc3k.h5ad'  # the file that will store the analysis results

读入单细胞测序文件为　AnnData 对象，它包括许多注释和代表鼠的slots,它有自己的hdf5格式

读入数据
[5]

adata = sc.read_10x_mtx(
    'data/filtered_gene_bc_matrices/hg19/',  # the directory with the `.mtx` file
    var_names='gene_symbols',                # use gene symbols for the variable names (variables-axis index)
    cache=True)                              # write a cache file for faster subsequent reading

[6]

adata.var_names_make_unique()  # this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`

[7]

adata

image.png

预处理

显示所有细胞中每个单个细胞中计数最高的基因。
[8]

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

normalizing counts per cell finished (0:00:00)

image.png

基本的数据过滤
[9]

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

filtered out 19024 genes that are detected in less than 3 cells

我们还需要收集一些有关线粒体的信息，因为这些信息对于质量控制很重要
高比例线粒体基因表明细胞质量较差（https://master.bioconductor.org/packages/release/workflows/html/simpleSingleCell.html#examining-gene-level-metrics），这可能是因为细胞穿孔后细胞质RNA丢失。线粒体大于单个转录物分子，不太可能通过细胞膜逸出。
使用 pp.calculate_qc_metrics，我们可以非常有效地计算许多指标
[10]

adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

一些计算质量测度的小提琴图：

count 数中表达的基因数量
每个细胞的总count数
线粒体基因count数的百分比
[11]

sc.pl.violin(adata, ['n_genes_by_counts'])
sc.pl.violin(adata, ['total_counts'])
sc.pl.violin(adata, ['pct_counts_mt'])

image.png

(按官网给的代码画出来是横着的，为什么呢？有一样的小伙伴留言，不知道是不是环境包版本不同的原因)

需要删除表达太多线粒体基因或总数过多的细胞。
[12]

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')

image.png

实际上是通过 AnnData对象的slots来进行过滤的。
[13]

adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]

总count数归一化（文库大小正确）数据矩阵 X 为每个细胞最多10,000reads，因此细胞的count数就可比。
[14]

sc.pp.normalize_total(adata, target_sum=1e4)

normalizing counts per cell
finished (0:00:00)

对数据进行对数转换
[15]

sc.pp.log1p(adata)

鉴定高度表达变化的基因
[16]

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

image.png

[17]

sc.pl.highly_variable_genes(adata)

image.png

将.AnnData对象的.raw 属性设置为归一化和对数化的原始基因表达，以供之后在差异测试和基因表达的可视化中使用。相当于是冻结了AnnData对象的状态。

可以通过调用.raw.to_adata() 来获取AnnData对象.raw中的对象
[18]

adata.raw = adata

如果在后面没有继续使用进行校正数据做 sc.pp.regress_out和sc.pp.scale 缩放它，也可以完全不用使用.raw。

先前的高可变基因检测结果将作为注释存储在.var.highly_variabl中，并会被PCA ，sc.pp.neighbors 以及随后的流形/图形工具自动检测。在这种情况下，不需要执行下面的过滤步骤。

[19]

adata = adata[:, adata.var.highly_variable]

淘汰每个细胞的总计数和表达的线粒体基因百分比的影响。将数据缩放为单位方差。
[20]

sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])

image.png

将每个基因缩放至单位方差，剪辑超过标准偏差10的值
[21]

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

主成分分析

通过运行主成分分析（PCA）来降低数据的维数，该分析可以揭示变化的主轴并对数据进行去噪。
[22]

sc.tl.pca(adata, svd_solver='arpack')

image.png

我们可以在PCA坐标中绘制散点图
[23]

sc.pl.pca(adata, color='CST3')

image.png

查看单个PC对数据总方差的贡献。这为我们提供了有关为计算细胞的群关系时（例如，在聚类函数sc.tl.louvain()或tSNE sc.tl.tsne() 中使用）应考虑的PC数量。根据经验，通常可以粗略估计PC的数量。
[24]

sc.pl.pca_variance_ratio(adata, log=True)

image.png

保存结果
[25]

adata.write(results_file)

[26]

adata

image.png

计算邻域图

使用数据矩阵的PCA表示来数据矩阵的聚类图。可以在此处简单地使用默认值。为了再现此数据Seurat的结果，我们采用以下值。
[27]

sc.pp.neighbors(adata, n_neighbors=175, n_pcs=40)

image.png

(版本不一样，使用同样参数结果也会不一样呢，n_neighbors指的是每个点的邻近点的数量，neighbors的个数越多，聚类数会越少。需要自己调整参数对比seurat 的分类结果参考是否正确，我也是调了很多次参数才细胞分成类似的哦)
pc的数量依赖于上面所做图，一般是选在拐点处的的主成分，只需要一个粗略值，不同的pc数量所产生的聚类形状也不同

嵌入邻域图

建议使用UMAP将图形嵌入二维中（[McInnes et al，2018] https://arxiv.org/abs/1802.03426。它比tSNE更忠实于流形的全局连通性，即它可以更好地保存轨迹。在某些情况下，可能会发现群集断开连接和类似的连接冲突。通常可以通过运行以下命令来补救它们：

tl.paga(adata)
pl.paga(adata, plot=False)  # remove `plot=False` if you want to see the coarse-grained graph
tl.umap(adata, init_pos='paga')

[28]

sc.tl.umap(adata)

image.png

[29]

sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

image.png

当设置adata的.raw属性时，先前的图显示了“raw”（标准化，对数化但未校正）的基因表达。您还可以通过明确声明不想使用来绘制缩放和校正后的基因表达。
[30]

sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)

image.png

聚类邻域图

与Seurat和其他一样，scanpy 推荐Traag 等人（2018）的Leiden图聚类方法（基于优化模块化的社区检测）。请注意，Leiden聚类直接将聚类邻域图细胞，这已在上一节中进行了计算。
[31]

sc.tl.leiden(adata)

image.png

绘制群集，这与Seurat的结果非常吻合
[32]

sc.pl.umap(adata, color=['leiden', 'CST3', 'NKG7'])

image.png

保存结果
[33]

adata.write(results_file)

寻找基因Marker

让我们计算每个簇中高度差异化基因的排名。为此，默认情况下AnnData属性.raw 被使用如果之前已被初始化。最简单，最快的方法是t检验。
[34]

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

image.png

[35]

sc.settings.verbosity = 2  # reduce the verbosity

Wilcoxon rank-sum (Mann-Whitney-U)
的结果的测试是非常相似的,建议在出版物中使用后者，例如，参见Sonison & Robinson (2018).
。您可能还会考虑功能更强大的差异测试程序包，例如MAST，limma，DESeq2，对于python，则是最新的diffxpy。
[36]

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

image.png

保存结果
[37]

adata.write(results_file)

或者，还可以使用逻辑回归对基因进行排名。例如，Natranos et al. (2018).
已经提出了这一点。本质上的区别在于，在这里，使用多变量方法，而传统的差异检验是单变量的， Clark et al. (2014)有更多细节。
[38]

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

image.png

除 IL7R，其仅由t检验发现，FCER1A，仅由其他两个方法发现，其它的所有的标记基因在所有的方法中发现。

Louvain Group	Markers	Cell Type
0	IL7R	CD4 T cells
1	CD14, LYZ	CD14+ Monocytes
2	MS4A1	B cells
3	CD8A	CD8 T cells
4	FCGR3A, MS4A7	FCGR3A+ Monocytes
5	GNLY, NKG7	NK cells
6	FCER1A, CST3	Dendritic Cells
7	PPBP	Megakaryocytes

定义标记基因的列表，以供参考
[39]

marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
                'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']

重新加载已保存的Wilcoxon Rank-Sum 测试结果的对象。
在数据框中显示每个聚类0、1，...，7的前10个排名最高的基因
[40]

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

[41]

image.png

获取带有评分和组的表。
[42]

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)

image.png

与单个群集进行比较
[43]

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

image.png

如果想要某个组的更详细的视图，使用 sc.pl.rank_genes_groups_violin。
[44]

sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)

image.png

与其余组进行比较，重新加载对象计算差异表达式
[45]

adata = sc.read(results_file)

[46]

sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)

image.png

如果要跨组比较某个基因，可以使用
[47]

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

image.png

实际标记细胞类型

[48]

new_cluster_names = [
    'CD4 T', 'CD14 Monocytes',
    'B', 'CD8 T',
    'FCGR3A Monocytes','NK',
    'Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)

[49]

sc.pl.umap(adata, color='leiden',size=50, legend_loc='on data',legend_fontsize=6, title='', frameon=False, save='.pdf')

image.png

对比官网，基本一样了

注释了细胞类型可视化标记基因。
[50]

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

image.png

还有一个非常紧凑的小提琴图
[51]

sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90);

image.png

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

adata

image.png

[53]

adata.write(results_file, compression='gzip')  # `compression='gzip'` saves disk space, but slows down writing and subsequent reading

使用下来可以大致了解该文件h5ls，该文件有很多选项-有关更多详细信息，请参见此处。文件格式将来可能仍会受到进一步优化。但是，所有读取功能都将保持向后兼容。

如果要与只想使用它进行可视化的人共享该文件，则减小文件大小的一种简单方法是删除密集缩放和校正的数据矩阵。该文件仍包含的可视化文件使用的原始数据为adata.raw。
[54]

adata.raw.to_adata().write('./write/pbmc3k_withoutX.h5ad')

[55]
如果要导出到“ csv”，则可以使用以下选项：

# Export single fields of the annotation of observations
# adata.obs[['n_counts', 'louvain_groups']].to_csv(
#     './write/pbmc3k_corrected_louvain_groups.csv')


# Export single columns of the multidimensional annotation
# adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv(
#     './write/pbmc3k_corrected_X_pca.csv')

# Or export everything except the data using `.write_csvs`.
# Set `skip_data=False` if you also want to export the data.
# adata.write_csvs(results_file[:-5], )

网友评论

本文标题：用VSCode Jupyter 学习Scanpy——预处理与分群

本文链接：https://www.haomeiwen.com/subject/rafcnktx.html

延伸阅读

深度阅读

您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

用VSCode Jupyter 学习Scanpy——预处理与分群

预处理

主成分分析

计算邻域图

嵌入邻域图

聚类邻域图

寻找基因Marker

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读