最近一直在用scanpy分析单细胞数据，真的是越用越顺手，但常在河边走哪有不湿鞋的呢？这不一不小心就给自己埋下了一个难以察觉的潜在问题。为什么这么说呢？因为这个问题通常情况下不会遇到，而且就算一不小心入坑，大部分情况下也不会受任何影响，仅在某些特殊需求时出现问题。
当确定了细胞类型后，下一步肯定要看看各细胞类型间的差异基因，如果做了两种条件的处理，也会看看各细胞类型在两个条件下的差异基因。差异基因做完了，接着会选择一些关注的基因展示一下具体的表达情况，不管是出于验证处理过程正确与否还是作为示例展示，看差异基因的具体表达情况都是最直接有效的方式，类似如下图。
通过这样的小提琴图，可以很直观地看出基因在不同条件下的表达分布情况。为了做出这样的图，结果意外出现了。用了好长时间才反应过来，发现问题所在。
import scanpy as sc
import matplotlib.pyplot as plt
adata = sc.read_h5ad('human_heart.h5ad')
adata = adata[adata.obs['sample'].isin(['ctrl1','ctrl2']),:]
adata.layers['counts'] = adata.X
adata.X = adata.layers['counts'].toarray()
原始read count存储在对象的counts属性层里面，从中提取想要的子集出来分析，然后用read count替换X矩阵的数据以便从头分析。也许，眼尖的朋友在这里已经发现了一些端倪。后续接着就是一套标准处理：
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
adata.raw = adata
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pp.scale(adata, max_value=10)
sc.pp.pca(adata, n_comps=100)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=50)
sc.tl.leiden(adata, resolution=0.5)
sc.tl.umap(adata)
细胞类型鉴定过程忽略，这里以Endothelial为例，来做ctrl1和ctrl2两者的差异基因分析。到这里，一切都很顺利没有任何问题。
ad_endo = adata[adata.obs.celltype=='Endothelial',:]
sc.tl.rank_genes_groups(ad_endo , groupby="sample", groups=['ctrl1'], reference='ctrl2', pts=0.1, method="wilcoxon")
做完差异分析，想看看基因的表达值在两个条件下的具体情况时，前面埋的雷在这里爆了。抛出了下面的错误：
sc.pl.rank_genes_groups_violin(ad_endo, groups='ctrl1', n_genes=8)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/software/lib/python3.6/site-packages/scanpy/plotting/_tools/__init__.py", line 913, in rank_genes_groups_violin
df[g] = X_col
File "/home/software/lib/python3.6/site-packages/pandas/core/frame.py", line 3044, in __setitem__
self._set_item(key, value)
File "/home/software/lib/python3.6/site-packages/pandas/core/frame.py", line 3119, in _set_item
self._ensure_valid_index(value)
File "/home/software/lib/python3.6/site-packages/pandas/core/frame.py", line 3170, in _ensure_valid_index
value = Series(value)
File "/home/software/lib/python3.6/site-packages/pandas/core/series.py", line 327, in __init__
data = sanitize_array(data, index, dtype, copy, raise_cast_failure=True)
File "/home/software/lib/python3.6/site-packages/pandas/core/construction.py", line 496, in sanitize_array
raise Exception("Data must be 1-dimensional")
Exception: Data must be 1-dimensional
可以看出，这个潜在的问题不影响前面的整体分析流程，仅是在绘制差异基因小提琴图时有问题。从这点来说，其实也不算什么问题，只不过影响了绘图的便捷性。既然，问题出现了，还是要确定一下原因，以便确认问题影响的范围。
一开始挺蒙的，毕竟，一连串的标准过程都正确执行了，如果有问题早该出现了。经过一番思索和排除还是不知问题所在。但事出反常必有妖，还是要搞清楚问题，于是就继续倒腾。最终，在浪费了一天的时间后，搞清楚啥时候埋的雷了。
有些不起眼的反常，可能就是问题的罪魁祸首。所以，没事不要瞎搞，还是按照官方路线走，不然，表面虽然看似风平浪静，但实际已埋下祸根，不知后面啥时候就会爆出来。
就像这个问题，一个格式引发的异常，愣是没法联想到，为何？因为前面一系列的操作都毫无影响，有时候这种情况更危险，内部已经不正常，外面却没有任何表现，一副正常有序的景象，让人不易察觉。
正常情况下，adata对象里面的矩阵为基于numpy array的sparse matrix格式：
<3467x33538 sparse matrix of type '<class 'numpy.float32'>'
with 100688614 stored elements in Compressed Sparse Row format>
回过头看看，自己生成的adata对象里矩阵的格式，是普通的numpy array：
ad_endo.X
array([[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ],
[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ],
[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ],
...,
[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ],
[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ],
[-0.00851632, 0. , 0. , ..., -0.06525466,
0. , 0. ]], dtype=float32)
现在一切都明朗了，原来是前面导出原始read count矩阵的时候，不知当时咋回事在替换的同时加了一个格式转换，为后面的错误埋下了伏笔：
# 错误方式
adata.X = adata.layers['counts'].toarray()
# 正确方式
adata.X = adata.layers['counts']
原因找到了，纠正起来也就方便了。可以选择从头再分析一遍彻底解决隐患，或者仅是验证一下可以简单点，既然不影响分析结果，修改一下矩阵格式即可。
from scipy.sparse import csr_matrix
ad_endo.X = csr_matrix(ad_endo.X)
ad_endo.X
<3467x33538 sparse matrix of type '<class 'numpy.float32'>'
with 100688614 stored elements in Compressed Sparse Row format>
此时，修改后再绘图就可以顺利出图了：
sc.pl.rank_genes_groups_violin(ad_endo, groups='ctrl1', n_genes=8)
不过，这样虽然可以画出图，但是看起来有点怪怪的，因为这里绘图使用scale的数据。
网友评论