从单细胞表达矩阵构建Anndata对象

之前的文章整理了读取单细胞10X标准数据的流程，后面发现更多的单细胞数据只是一个单纯的表达矩阵，文件后缀可能是csv.gz/matrix.csv.gz/txt.gz等等，原理基本是一致的。

这里记录一下如何快速有效地读取一个单细胞表达矩阵，以及简单的过程检查。

这里以 Yamagata M, Yan W, Sanes JR. A cell atlas of the chick retina based on single-cell transcriptomics. Elife 2021 Jan 4;10. PMID: 33393903 文章的数据为例

导入需要的包

import pandas as pd
import numpy as np
import scanpy as sc
import anndata as ad
from anndata import AnnData
import gzip
import os
import sys

为了快速了解这些压缩矩阵的数据结构，我们可以在读取前先查看前五行

file_path = '/data/work/retina_sc/chick/GSE159107_E12chick_count.matrix.csv.gz'

with gzip.open(file_path, 'rt') as file:
    for i in range(5):
        line = file.readline()
        print(line)

下拉查看发现第一行是单细胞标识

从第二行开始是基因名和计数，即首列为基因名，首行为细胞名，且每个元素以逗号分隔

于是我们正式读取压缩矩阵，file_path是下载的压缩包存放的位置

file_path = '/data/work/retina_sc/chick/GSE159107_E12chick_count.matrix.csv.gz'
adata_df = pd.read_csv(file_path, sep=',', index_col=0)
adata_df = adata_df[~adata_df.index.duplicated(keep='first')]
adata = sc.AnnData(X=adata_df.values.T, obs=pd.DataFrame(index=adata_df.columns), var=pd.DataFrame(index=adata_df.index))
adata

最后检查一下anndata对象的index是否正确

adata.var.index

adata.obs.index

var索引是基因名，obs观测值索引是细胞标识，没啥问题就可以保存一下转化之后的原始数据

adata.write('/data/work/retina_sc/chick/GSE159107_E12chick.h5ad')

另外两个批次地数据可以重复上述操作，最后如果想要合并多个批次数据

dataset = AnnData.concatenate(*[sc.read_h5ad(os.path.join("/data/work/retina_sc/chick", f)) 
                                for f in sorted(os.listdir("/data/work/retina_sc/chick")) if f.endswith("h5ad")])
dataset

观测值obs中会多一个列叫'batch'，三个批次数据默认标记为'0''1''2'