scanpy踩坑实录

源代码 · 2024-9-9 08:23:44

一样平常瞎掰

迩来在分析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 scadata = sc.read_h5ad('pbmc3k_scanpy.h5ad')adataAnnData object with n_obs × n_vars = 2700 × 32738 var: 'gene_ids'adata.count = adata.Xadata.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), :]adataView 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.countTraceback (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')adataAnnData object with n_obs × n_vars = 2700 × 32738 var: 'gene_ids'adata.raw = adataadata.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), :]adataView 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.XTraceback (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数据质控

scanpy踩坑实录

所属分类: 问答交流

新帖推荐: 30日

推荐作品