scanpy读取10x单细胞数据

做单细胞或空间组课题时经常会需要导入文献中的单细胞数据作为参考，市面上最常见的格式又以10x genomics为主要代表，通常包括barcodes.tsv.gz、features.tsv.gz（或者genes.tsv.gz）、matrix.mtx.gz三种格式文件。

在面对数据读取问题时，R语言Seurat包有Read10X函数，Python中scanpy包则对应scanpy.read_10x_mtx()函数。官网的导入教程scanpy.read_10x_mtx — scanpy但在具体使用时，这两个函数对于单细胞格式和命名要求都很高。在scanpy中尝试将所有文件去除前缀后仍然找不到文件后，索性决定单独读取各个要素，生成AnnData格式的数据文件即可。

首先是导入包和数据路径

import os
import pandas as pd
from scipy.io import mmread
import anndata
#10x文件存放的路径
path='/data/work/scRNA_load/2020.05.Dev_cell/'
genes=None
cell=None
mtx=None
for name in os.listdir(path):
    file= os.path.join(path, name)
    if 'barcodes.tsv' in name:
        cell = pd.read_csv(file, header=None)
    elif 'genes.tsv' in name:
        genes = pd.read_csv(file, sep='\t', header=None)
    elif 'matrix.mtx' in name:
        mtx = mmread(file)

查看mtx矩阵的数据结构

mtx

<33694x42834 sparse matrix of type '<class 'numpy.int64'>'
	with 76381528 stored elements in COOrdinate format>

对应我们的数据，42834个细胞和33694个基因，这里表达矩阵的顺序反了，AnnData格式对应的应该是“细胞obs*基因vars”。建立正确的AnnData矩阵

adata = anndata.AnnData(mtx)
adata=adata.T #转置
adata.obs.index = pd.Index(cell[0])
adata.var.index = pd.Index(genes[1])
adata.var['gene_ids'] = genes[0].to_list()
adata.var_names = genes[1].to_list()
adata.obs_names = cell[0].to_list()

这里是genes变量的前五行

0     ENSG00000243485     RP11-34P13.3
1     ENSG00000237613     FAM138A
2     ENSG00000186092     OR4F5
3     ENSG00000238009     RP11-34P13.7
4     ENSG00000239945     RP11-34P13.8

其中0，1，2，3，4是行索引，ENSG00000243485 等表示基因的唯一标识符（Ensembl ID），RP11-34P13.3 等，表示基因的符号名称（Gene Symbol）。所以我们在这里将基因的唯一标识符列转化成AnnData中的‘gene_ids’,将基因名称作为索引。

最后需要将稀疏矩阵转成密集矩阵，否则无法保存结果文件。

adata.X = adata.X.A

顺利保存.h5ad的原始AnnData文件

output_file = '/data/work/scRNA_load/2020.05.Dev_cell/GSE116106_retina_raw.h5ad'
adata.write_h5ad(output_file)

后续读取时正常使用read_h5ad即可。