scanpy踩坑实录

日常瞎掰

  最近在分析scRNA-seq的数据，三个样本数据不大，用Seurat分析真的很香！奈何又要下载公共数据的样本，这样一捣腾，Seuart倒显得有些力不从心了，何也？没成想公共数据提供的矩阵里面有48万多个细胞，加上自有的两万多个细胞，共有50多万个细胞，这数据量用Seuart处理多少有点费劲了，怎么办？所以，很自然地就想到了scanpy，python处理大数据的能力比起R还是没话说的。可是，习惯了Seuart，突然转到较为陌生的scanpy，不知不觉就花了好多时间在踩坑上面，特此记录一下。

踩坑1

  scanpy对象的.X属性里面存储的是表达矩阵，可是这个矩阵会随着处理步骤一直变化，如normalize、log1p步骤后，里面就变成标准化的log矩阵了，就不再是原始count值了，并不会像Seuart那样会分别存储原始矩阵和处理后的矩阵。多样本的对象合并前最好不要做任何数据处理步骤，否则就要注意了。像本人在合并自有数据和公共数据时，先处理了自有数据随后合并公共数据，然后就遇到了下图中的问题。自己还在网上各种查找解决问题的办法，还别说真有遇到同样问题的情况，不过给出的答案也没法解决。

  当时并没有反应过来，现在想来这个坑还是自己挖的，原因就在于，在合并自有和公共数据前，先处理了自有数据，所以里面的矩阵值是标准化的log值，合并公共数据后再对数据做normalize、log1p时产生了很多nan值，从而引发这个错误。想想真是无语至极，此处省略一万字的内心活动。
  不过，话说回来，踩坑也并不完全是坏事，这样以后合并不同的样本对象时，就知道怎么做了。如果再遇到同样的问题，首先就会知道应该先检查对象矩阵里面是不是原始count值。

踩坑2

  scanpy对象里面保存原始count矩阵，其实对于这个需求scanpy已经提前想到了，所以scanpy对象预留了.raw属性来存储。当然，这个属性只能保存scanpy对象不能直接只保存表达矩阵，如果想只保存矩阵可以自定义一个属性来存储。那么，问题来了，自定义属性其实并不安全，数据处理过程中有可能丢失。这到底怎么回事呢？看看下面一段代码：

import scanpy as sc

adata = sc.read_h5ad('pbmc3k_scanpy.h5ad')
adata
AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'

adata.count = adata.X
adata.count
<2700x32738 sparse matrix of type '<class 'numpy.float32'>'
        with 2286884 stored elements in Compressed Sparse Row format>

adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)

sc.pp.filter_cells(adata, min_genes = 200)
sc.pp.filter_genes(adata, min_cells = 3)
adata = adata[(adata.obs.total_counts > 500) & (adata.obs.n_genes > 300) & (adata.obs['pct_counts_mt'] <= 5), :]
adata
View of AnnData object with n_obs × n_vars = 2633 × 13714
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'

adata.count
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'AnnData' object has no attribute 'count'

  从上面的代码可以知道，如果对象采用类似R数据框形式的方式筛选子集会丢失掉自定义属性。所以，用自定义属性方式存储原始矩阵时，得要保证后面没有这样方式的取子集操作，否则属性丢了都不知道。
  当然，我们可以用.raw属性来保存原始的对象，这样原始矩阵也就会一同被保留下来。而且这个方式不会因为子集操作而丢失，是不是就可以随意操作了呢？看看下面的代码：

adata = sc.read_h5ad('pbmc3k_scanpy.h5ad')
adata
AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'

adata.raw = adata
adata.raw.count
<2700x32738 sparse matrix of type '<class 'numpy.float32'>'
        with 2286884 stored elements in Compressed Sparse Row format>

adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)

sc.pp.filter_cells(adata, min_genes = 200)
sc.pp.filter_genes(adata, min_cells = 3)
adata = adata[(adata.obs.total_counts > 500) & (adata.obs.n_genes > 300) & (adata.obs['pct_counts_mt'] <= 5), :]
adata
View of AnnData object with n_obs × n_vars = 2633 × 13714
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'

adata.raw.X
<2633x32738 sparse matrix of type '<class 'numpy.float32'>'
        with 2254830 stored elements in Compressed Sparse Row format>

adata.X = adata.raw.X
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "~/.local/lib/python3.6/site-packages/anndata/_core/anndata.py", line 696, in X
    f"Data matrix has wrong shape {value.shape}, "
ValueError: Data matrix has wrong shape (2633, 32738), need to be (2633, 13714).

  从上面的代码可以看到，.raw.X确实保存了原始矩阵，而且也没有被取子集操作影响，可以说没有一点变化还是保留了所有细胞。当然，使用这种方式还是有好处的，方便后续直接使用，因为scanpy有些步骤可以指定使用原始矩阵，默认就是从这个属性里面读取的。
  所以，从上面两种方式可以得知，想要在对象里面保留原始矩阵，需要注意时间节点，不然原始矩阵与正在使用的矩阵维度不同，也没法直接使用原始矩阵。

结束语

  scanpy的处理数据的过程基本与Seuart保持一致，完全可以类比，如果熟悉Seurat的使用，那么使用scanpy的学习成本就会降低不少。一般按照常规分析流程，也不会遇到什么问题。如果进行一些个性操作，就要注意scanpy的特性了，否则不出意外的话就会出现意外了。在不熟悉的情况下，来一些非常规操作就出现了意外，说不定就会像咱一样晕乎乎的走一些弯路。因此，既然要进行数据分析了，势必会进行这样那样的操作，所以分析前还是有必要先看看文档认识一下scanpy，正所谓磨刀不误砍柴工嘛，这样做也许后续可以避免不少的麻烦。

往期回顾

差异基因密度分布
R绘图配色总结
saddleplot | A/B compartments
双曲线火山图一键拿捏
ChIP-seq数据质控