scanpy读取数据失败记录：张泽民团队的单细胞数据GSE212890
scanpy读取数据失败记录：张泽民团队的单细胞数据GSE212890
尝试用 scanpy 读取单细胞数据 GSE212890 的时候，因为
scanpy.read_10x_mtx()函数对单细胞数据的格式和命名要求很高，记录一下在数据读取方面遇到的问题和解决办法。新手上路，有什么不足还请多多批评指正。
先来看一下数据的基本信息：
数据来源于 2023 年发表在 Cell 上的《A pan-cancer single-cell panorama of human natural killer cells》，来自张泽民大佬团队的单细胞文章。其中，数据集 GSE212890 包含 47 个样本，包含子宫内膜癌、肾癌、食管癌、甲状腺癌、胰腺癌、卵巢癌、乳腺癌、输卵管癌等 8 种癌症的肿瘤组织和癌旁正常组织。
原始测序数据可在：https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000321 申请获取。
本文利用计算门控（computational gating），通过 CD3⁻CD56⁺/KLRF1⁺ 筛选 NK 细胞（排除 T 细胞和其他免疫细胞），下图是筛选后的 NK 细胞数据：
处理过程：
这个数据看着和 10x Genomics 格式非常相似，所以我一开始的想法是改一下数据命名和格式以适配 scanpy.read_10x_mtx() 函数。10x数据通常包括 barcodes.tsv.gz、features.tsv.gz（或 genes.tsv.gz）、matrix.mtx.gz 三种文件，因此先将数据重命名，并将 csv 格式改为 tsv 格式：
gunzip GSE212890_NK.barcodes.csv.gz GSE212890_NK.genes.csv.gz
less GSE212890_NK.barcodes.csv |cut -d ',' -f 2 |sed '1d' >barcodes.tsv
less GSE212890_NK.genes.csv |cut -d ',' -f 2 |sed '1d' >features.tsv
gzip barcodes.tsv features.tsv
cp GSE212890_NK_counts.mtx.gz matrix.mtx.gz
数据分别是这样的：
barcodes.tsv.gz：
features.tsv.gz：
matrix.mtx.gz：
尝试读取数据：
import scanpy as sc
adata = sc.read_10x_mtx(
sample_dir, #存放 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz 三个文件的路径
var_names="gene_symbols" #可选 gene_ids 和 gene_symbols
)
报错 KeyError：1，经检查发现，标准的 features.tsv.gz 通常有两列，第一列是 gene_ids，第二列是 gene_symbols。当设置 var_names="gene_symbols" 时，scanpy 会尝试从 features.tsv.gz 文件中读取第二列作为基因名称，但这里的 features.tsv.gz 文件只有一列 gene_symbols，所以会报下标溢出错误。因此试着在 features.tsv.gz 第一列前面加一个制表符 \t 再次尝试运行：
less GSE212890_NK.genes.csv |cut -d ',' -f 2 |sed '1d' |sed 's/^/\t/' >features.tsv
gzip features.tsv
再次报错：
这次报错的信息是 ValueError: Length of passed value for var_names is 28991, but this AnnData has shape: (28991, 11963)
这个错误的意思是：我们传递给 AnnData 的 var_names（即基因名）有28991个，但是 AnnData 矩阵的形状是（28991, 11963），这个错误乍一看很奇怪啊，貌似感觉是维度不匹配的问题。再检查一下：
features.tsv.gz 有 28991 行，表示 28991 个基因
barcodes.tsv.gz 有 11963 行，表示 11963 个细胞
这里看着和 AnnData 的维度是一致的，检查一下 mtx 文件的格式：
先来了解一下标准的 mtx 文件的格式：
- 以
%符号开头的行是注释行，提供关于文件的元数据。 - 第一个非注释行是最关键的一行，这一行包含三个整数：第一个是矩阵的行数，代表基因的数量；第二个是矩阵的列数，代表细胞的数量；第三个是矩阵中非零元素的种数。
- 之后的所有行都包含实际的矩阵数据，每一行代表矩阵中的一个非零元素，每一行的三个值分别代表该元素所在的行索引；该元素所在的列索引；该元素的值（代表该基因在该细胞中的表达量）。
可以看到这里的 matrix.mtx.gz 数据应该是第一列变成了细胞数，第二列变成了基因数，和正常的 mtx 数据相比反了过来，所以导致 AnnData 和 var_names 不匹配：
因此把这两列调换一下：
zcat matrix.mtx.gz |awk '/^%/ {print; next} {print $2,$1,$3}' |gzip >matrix_transposed.mtx.gz
mv matrix_transposed.mtx.gz matrix.mtx.gz
重新跑一遍
还在报错 KeyError：2，实在不知道怎么办了，还是把三个文件分开读取尝试一下吧：
import gzip
import pandas as pd
import scipy.io
from scipy.sparse import csr_matrix
import anndata as ad
# 读取细胞信息
with gzip.open(os.path.join(sample_dir,'barcodes.tsv.gz'), 'rt') as f:
barcodes = [line.strip() for line in f]
# 读取基因信息
with gzip.open(os.path.join(sample_dir,'features.tsv.gz'), 'rt') as f:
features = [line.strip() for line in f]
# 读取 mtx 数据
with gzip.open(os.path.join(sample_dir,'matrix.mtx.gz'), 'rb') as f:
mtx = scipy.io.mmread(f).tocsr()
# 把三个文件合并在一起
adata = ad.AnnData(
X=mtx.T,
obs=pd.DataFrame(index=barcodes),
var=pd.DataFrame(index=features)
)
# 处理重复基因名
adata.var_names_make_unique()
# 验证结果
print(adata)
AnnData object with n_obs × n_vars = ×
顺便读取一下 metadata 信息
metadata = pd.read_csv(os.path.join(sample_dir,'GSE212890_NK_metadata.csv.gz'), compression = 'gzip', index_col = )
adata.obs = metadata
print(adata)
AnnData object with n_obs × n_vars = ×
obs: 'sampleID', 'cellID', 'meta_patientID', 'datasets', 'meta_histology', 'meta_platform', 'meta_tissue', 'Majortype', 'celltype'
这样就没问题啦~
腾讯云开发者
