scanpy怎么分开读取GEO数据库的10X单细胞3个文件
scanpy怎么分开读取GEO数据库的10X单细胞3个文件
之前在在单细胞天地教程：表达矩阵逆转为10X的标准输出3个文件，详细介绍过 10X文件的3个标准文件。比如SRR7722939数据集里面，文件
barcodes.tsv和genes.tsv，就是表达矩阵的行名和列名：
jmzengdeMacBook-Pro:SRR7722939 jmzeng$ head barcodes.tsv
AAACCTGAGCGAAGGG-1
AAACCTGAGGTCATCT-1
AAACCTGAGTCCTCCT-1
AAACCTGCACCAGCAC-1
AAACCTGGTAACGTTC-1
AAACCTGGTAAGGATT-1
AAACCTGGTTGTCGCG-1
AAACCTGTCCTGCCAT-1
AAACGGGAGTCATCCA-1
AAACGGGCATGGATGG-1
jmzengdeMacBook-Pro:SRR7722939 jmzeng$ head genes.tsv
hg38_ENSG00000243485 hg38_RP11-34P13.3
hg38_ENSG00000237613 hg38_FAM138A
hg38_ENSG00000186092 hg38_OR4F5
hg38_ENSG00000238009 hg38_RP11-34P13.7
hg38_ENSG00000239945 hg38_RP11-34P13.8
hg38_ENSG00000239906 hg38_RP11-34P13.14
hg38_ENSG00000241599 hg38_RP11-34P13.9
hg38_ENSG00000279928 hg38_FO538757.3
hg38_ENSG00000279457 hg38_FO538757.2
hg38_ENSG00000228463 hg38_AP006222.2
但是matrix.mtx，就稍微复杂一点，仔细看：
jmzengdeMacBook-Pro:SRR7722939 jmzeng$ head matrix.mtx
%%MatrixMarket matrix coordinate integer general
%
33694 2049 1878957
28 1 1
55 1 2
59 1 1
60 1 1
62 1 1
78 1 2
111 1 1
如果你关注这3个文件的行数：
2049 barcodes.tsv
33694 genes.tsv
1878960 matrix.mtx
就会发现，matrix.mtx文件里面的33694、2049、1878957数值，分别是细胞数量，基因数量，以及有表达量的值的数量（全部的值应该是33694X2049接近7000万，但是有值的仅仅是不到200万，所以单细胞矩阵里面只有3%左右的值大于0 ）。每个10X样本都是走流程拿到10x单细胞转录组数据的3个文件的表达矩阵。
但是，值得注意的是上面的3个文件，其实是10x这个公司的单细胞的输出文件，其实是cellranger流程的，这个时候如果是跑的老版本cellranger，就是genes.tsv，但是如果是新版本，就是（gz文件），详见Cell Ranger 3.0 VS 2.0做了哪些改动。
可以看到早期版本的文件是不压缩，后期新版本都是压缩的，所以如果碰到了压缩了的早期版本文件，就会冲突，在Python的scanpy里面。比如我们 进入网站GSM3148575，如下图所示：
可知这是一个10X数据文件，很标准的单细胞3个文件，但是它被压缩了，理论上应该是新版，而它文件名里面是genes.tsv，又意味着是旧版，就冲突了！我们直接看网页描述吧：
由上图可知，这是一个使用Cell Ranger V2处理后的10X数据文件，下载该数据。如果你直接使用scanpy的read_10x_mtx()函数进行读取的话，会读取失败，因为Cell Ranger V3将输出文件压缩以减少磁盘空间（gz文件），详见Cell Ranger 3.0 VS 2.0做了哪些改动 ,scanpy会把它认为这是一个使用Cell Ranger V3处理后的10X数据文件，会使用读取Cell Ranger V3处理后的10X数据文件的方式进行读取，然而这其实是一个“伪Cell Ranger V3”的10X文件，是一个真正的Cell Ranger V2处理后的10X数据文件.
那么要怎么成功读取这个文件呢?
我这里演示一下，使用Python的scanpy分别读取3个文件后，再组合起来：
代码如下：
import os
import pandas as pd
from scipy.io import mmread
import anndata
#10x文件存放的路径
path='C:/Users/data/GSM3148575/'
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)
#创建scanpy对象
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()
查看scanpy对象,代码如下：
print(adata)
结果如下：
AnnData object with n_obs × n_vars = 7096 × 33694
var: 'gene_ids'