• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)单细胞数据scanpy 学习笔记

      (未测试)单细胞数据scanpy 学习笔记

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2020年2月23日
      • 评论 0评论

      Scanpy 是一个基于 Python 分析单细胞数据的软件包,内容包括预处理,可视化,聚类,拟时序分析和差异表达分析等。本文翻译自 scanpy 的官方教程 Preprocessing and clustering 3k PBMCs[1],用 scanpy 重现 Seurat 聚类教程[2] 中的绝大部分内容。

      0. scanpy 安装

      直接用 pip ,安装成功。

      pip install scanpy[louvain]

      Docker

      docker pull fastgenomics/scanpy:1.4-p368-v1-stretch-slim

      1. 载入数据

      # 下载PBMC 数据集
      ## 其实就是 Seurat 那个示例数据,之前下过就不用重复下了
      !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
      import numpy as np
      import pandas as pd
      import scanpy as sc
      /Users/baozhiwei/anaconda3/lib/python3.7/site-packages/anndata/core/anndata.py:17: FutureWarning: pandas.core.index is deprecated and will be removed in a future version.  The public classes are available in the top-level namespace.
        from pandas.core.index import RangeIndex
      # verbosity 的取值表示测试结果显示的详细程度,数字越大越详细
      ## errors (0), warnings (1), info (2), hints (3)
      sc.settings.verbosity = 3
      # 输出版本号
      sc.logging.print_versions()
      # set_figure_params 设置图片的分辨率/大小以及其他样式
      sc.settings.set_figure_params(dpi=80)
      scanpy==1.4.5.post3 anndata==0.6.22.post1 umap==0.3.10 numpy==1.18.1 scipy==1.4.1 pandas==1.0.1 scikit-learn==0.22.1 statsmodels==0.11.0 python-igraph==0.7.1
      # 设置结果文件保存路径
      results_file = './pbmc3k.h5ad'
      # 导入 10X 数据
      adata = sc.read_10x_mtx(
          './data/filtered_gene_bc_matrices/hg19/',  # 包含有 `.mtx` 文件的目录
          var_names='gene_symbols',                  # 用 gene symbols 作为变量名 (variables-axis index)
          cache=True)                                # 使用缓存文件加快读取
      ... writing an h5ad cache file to speedup reading next time

      scanpy 的结果文件是 AnnData 对象,意为 annotated data。AnnData 用了一些广义的单词来描述「细胞」和「基因」:将细胞称为观察值 observations ,将基因称为变量 variables 。AnnData 包括四个可以存储信息的区域:

      •adata.X 存储 count matrix,数据类型为稀疏矩阵 scipy.sparse.csr.csr_matrix

      •adata.obs 存储关于 obervations(cells) 的 metadata,数据类型为 dataframe

      •adata.var 存储关于 variables(genes) 的 metadata,数据类型为 dataframe

      •AnnData.uns 存储后续附加的其他非结构化信息

      •adata.obs_names 和 adata.var_names index

      细胞名和基因名可分别通过 adata.obs_names 和 adata.var_names 查看。 AnnData 对象可以像 dataframe 一样进行切片操作,例如,adata_subset = adata[:, list_of_gene_names]。

      # 索引去重,若上一步中使用 `var_names='gene_ids'` 则这一步非必须进行
      # 其实也可以直接用 pandas 判断索引是否重复
      # adata.var_names.is_unique
      adata.var_names_make_unique()

      adata

      AnnData object with n_obs × n_vars = 2700 × 32738 
          var: 'gene_ids'

      2. 数据预处理

      可视化所有细胞中计数最多的基因。

      sc.pl.highest_expr_genes(adata, n_top=20, )
      normalizing counts per cell
          finished (0:00:00)

      数据初筛

      # 保留至少在三个细胞中表达的基因,保留至少包含 200 个基因的细胞
      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

      标注线粒体基因并计算每个细胞中的线粒体基因比例。

      mito_genes = adata.var_names.str.startswith('MT-')
      # 计算每个细胞中的线粒体基因比例
      # 使用`.A1` 将 numpy.matrix 转为一维数组 ndarray
      adata.obs['percent_mito'] = np.sum(
          adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
      # 把每个细胞的基因计数添加到 adata 中
      adata.obs['n_counts'] = adata.X.sum(axis=1).A1

      方法二

      还可以使用 pp.calculate_qc_metrics 来计算每个细胞中线粒体基因的数量,以及其他一系列 QC 数据。关于更多 preprocessing 的内容可参见文档[3]。

      # 标记线粒体基因
      adata.var['mito'] = mito_genes
      calculate_qc_metrics 函数会返回由两个 dataframe 组成的元组,一个是细胞的 QC 矩阵,另一个是基因的 QC 矩阵。
      qc = sc.pp.calculate_qc_metrics(adata,qc_vars=['mito'])
      cell_qc_dataframe = qc[0]
      gene_qc_dataframe = qc[1]
      cell_qc_dataframe.iloc[:4, :4]
      gene_qc_dataframe.iloc[:4, :4]

      用小提琴图可视化

      sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'],
                   jitter=0.4, multi_panel=True)

      可视化特征之间的关系。

      sc.pl.scatter(adata, x='n_counts', y='percent_mito')
      sc.pl.scatter(adata, x='n_counts', y='n_genes')
      adata
      AnnData object with n_obs × n_vars = 2700 × 13714 
          obs: 'n_genes', 'percent_mito', 'n_counts'
          var: 'gene_ids', 'n_cells'

      过滤线粒体基因比例 > 5% 和基因总数 >2500 的细胞。

      adata = adata[adata.obs.n_genes < 2500, :]
      adata = adata[adata.obs.percent_mito < 0.05, :]
      adata
      AnnData object with n_obs × n_vars = 2638 × 13714 
          obs: 'n_genes', 'percent_mito', 'n_counts'
          var: 'gene_ids', 'n_cells'

      标准化数据

      sc.pp.normalize_total(adata, target_sum=1e4)
      normalizing counts per cell
          finished (0:00:01)

      进行自然对数转换

      sc.pp.log1p(adata)

      可以将 AnnData 对象的 .raw 属性设置为经归一化和对数化的原始基因表达值,供之后的可视化分析使用。

      adata.raw = adata

      选择差异基因

      sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
      extracting highly variable genes
          finished (0:00:00)
      --> added
          'highly_variable', boolean vector (adata.var)
          'means', float vector (adata.var)
          'dispersions', float vector (adata.var)
          'dispersions_norm', float vector (adata.var)

      这一步在 adata.var 中添加了四列内容(highly_variable,means,dispersions,dispersions_norm)。

      sc.pl.highly_variable_genes(adata)

      取出高度差异的基因。

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

      校正细胞基因计数和线粒体基因比例的影响。

      sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])
      regressing out ['n_counts', 'percent_mito']
          sparse input is densified and may lead to high memory use
          finished (0:00:04)

      数据缩放

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

      3. PCA

      sc.tl.pca(adata, svd_solver='arpack')
      # svd_solver 指定奇异值分解 SVD 的方法
      computing PCA with n_comps = 50
          on highly variable genes
          finished (0:00:00)

      绘制碎石图,确定数据维度。

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

      保存结果。

      adata.write(results_file)
      adata
      AnnData object with n_obs × n_vars = 2638 × 1838
          obs: 'n_genes', 'percent_mito', 'n_counts'
          var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
          uns: 'log1p', 'pca'
          obsm: 'X_pca'
          varm: 'PCs'

      4. 细胞聚类

      为了重现 Seurat 的结果,我们使用下面的参数。

      sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
      computing neighbors
          using 'X_pca' with n_pcs = 40
          finished: added to `.uns['neighbors']`
          'distances', distances for each pair of neighbors
          'connectivities', weighted adjacency matrix (0:00:01)
      sc.tl.leiden(adata)
      running Leiden clustering
          finished: found 8 clusters and added
          'leiden', the cluster labels (adata.obs, categorical) (0:00:00)

      非线性降维 UMAP

      sc.tl.umap(adata)
      computing UMAP
          finished: added
          'X_umap', UMAP coordinates (adata.obsm) (0:00:04)

      绘制聚类图

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

      保存结果。

      adata.write(results_file)

      5. 找 marker gene

      接下来计算每个 cluster 中高度差异基因的排名。默认,若之前已初始化 AnnData 的 .raw 属性,则使用会该属性。

      t 检验

      sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
      sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
      ranking genes
          finished: added to `.uns['rank_genes_groups']`
          'names', sorted np.recarray to be indexed by group ids
          'scores', sorted np.recarray to be indexed by group ids
          'logfoldchanges', sorted np.recarray to be indexed by group ids
          'pvals', sorted np.recarray to be indexed by group ids
          'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:00)
      sc.settings.verbosity = 2  # reduce the verbosity

      Wilcoxon 秩和检验

      与 t 检验的结果非常相似。我们建议在文章中使用后者,参见 Sonison&Robinson(2018)[4]。

      sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
      sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
      ranking genes
          finished (0:00:02)

      保存结果

      adata.write(results_file)

      逻辑回归

      或者,也可以根据 Natranos et al. (2018)[5] 的建议使用逻辑回归对基因进行排名。本质区别在于,这里使用多元变量,而传统的差异检验是单变量。更多细节可参考 Clark et al. (2014)[6] 。

      sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
      sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
      ranking genes
          finished (0:00:03)

      除 IL7R 仅在 t 检验的结果中发现,以及仅在其他两种检验方法中发现的 FCER1A以外,其他标记基因均可通过所有检验方法得到。

      Louvain GroupMarkersCell Type
      0IL7RCD4 T cells
      1CD14, LYZCD14+ Monocytes
      2MS4A1B cells
      3CD8ACD8 T cells
      4GNLY, NKG7NK cells
      5FCGR3A, MS4A7FCGR3A+ Monocytes
      6FCER1A, CST3Dendritic Cells
      7PPBPMegakaryocytes

      定义一个标记基因的列表,方便之后的分析。

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

      重新加载保存有 Wilcoxon 秩和检验结果的对象。

      显示每个 cluster 中排名前五的基因。

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

      显示每个基因在每一 cluster 中的对应 p 值。

      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)

      cluster 与 cluster 之间进行比较。

      sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
      sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)
      ranking genes
          finished (0:00:01)

      还可用 sc.pl.rank_genes_groups_violin 进一步获得信息更丰富的图。

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

      重新载入对象,将单一 cluster 与其他剩余的 cluster 进行比较。

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

      如果要比较某个基因在不同组之间的差异,可使用以下代码。

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

      6. 标记细胞类型

      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='', frameon=False, save='.pdf')
      WARNING: saving figure to file figures/umap.pdf

      可视化 marker genes。

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

      绘制不同细胞亚群中多个基因的小提琴图。

      ax = sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90)
      adata
      AnnData object with n_obs × n_vars = 2638 × 1838
          obs: 'n_genes', 'percent_mito', 'n_counts', 'leiden'
          var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
          uns: 'leiden', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'
          obsm: 'X_pca', 'X_umap'
          varm: 'PCs'

      保存数据。

      adata.write(results_file, compression='gzip')

      引用链接

      [1]  https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html
      [2]  http://satijalab.org/seurat/pbmc3k_tutorial.html
      [3]  https://icb-scanpy.readthedocs-hosted.com/en/stable/api/index.html#module-scanpy.pp
      [4]  https://doi.org/10.1038/nmeth.4612
      [5]  https://doi.org/10.1101/258566
      [6]  https://doi.org/10.1186/1471-2105-15-79

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      fasta、fastq、vcf、bam、bed、gtf--多种生信格式的R语言读取
      2020年2月23日

      下一篇文章

      (未测试)跟着豆豆一起探索可变剪切分析
      2020年2月23日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?