最近一直在用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
的数据。