单细胞表达矩阵格式转换

作者：ahworld
链接：单细胞表达矩阵格式转换
来源：微信公众号-seqyuan
著作权归作者所有，任何形式的转载都请联系作者。

10x Genomics官方软件CellRanger输出的单细胞表达矩阵有三个文件:

• barcodes.tsv
• genes.tsv
• matrix.mtx

barcodes.tsv存储的是标识细胞的barcodes列表，格式如下：

AAACCTGAGCATCATC-1
AAACCTGAGCTAACTC-1
AAACCTGAGCTAGTGG-1
AAACCTGCACATTAGC-1
AAACCTGCACTGTTAG-1
AAACCTGCATAGTAAG-1

genes.tsv存储的是基因列表，共两列tab分割，第一列为gene.ids，第二列为gene.symbol：

ENSG00000243485 RP11-34P13.3
ENSG00000237613 FAM138A
ENSG00000186092 OR4F5
ENSG00000238009 RP11-34P13.7
ENSG00000239945 RP11-34P13.8
ENSG00000239906 RP11-34P13.14

matrix.mtx存储的是基因在各细胞（barcode）中的表达count：

• 第一列：gene在genes.tsv中的行号
• 第二列：barcode在barcodes.tsv中的行号
• 第三列：基因表达count

前两行固定，第三行统计了对应列数值之和。

%%MatrixMarket matrix coordinate integer general
%
33694 8381 11788294
33665 1 5
33663 1 5
33662 1 13
33661 1 1
33660 1 3

今天遇到一个问题：我下载了一个基因表达矩阵GSM3270887_countTable_colonCreMin.txt.gz，作为测试使用，这个表达矩阵为标准的Matrix，行名为gene.symbol，列名为barcode，如下：

AAACCTGAGCGGATCA	AAACCTGAGCTCAACT	AAACCTGCACTTAACG	AAACCTGCAGCGTCCA
Xkr4	0	0	0
Gm1992	0	0	0
Gm37381	0	0	0
Rp1	0	0	0
Rp1.1	0	0	0
Sox17	0	0	0
Gm37323	0	0	0
Mrpl15	0	1	0
Lypla1	0	0	1

因为是测试使用，这种N * N格式的 Matrix在读取速度上远逊于CellRanger的矩阵格式，所以我想把这个矩阵格式转换为CellRanger三个文件样式的矩阵。我的解决方案参考了biostars.org下面问题的答案，并做了修改。

Question: Storing a gene expression matrix in a matrix.mtx

用R读入数据

library(Matrix)
colon.data <- read.csv(file='GSM3270887_countTable_colonCreMin.txt.gz', sep="\t", header = T, row.names = 1)
colon.data <- Matrix(as.matrix(colon.data), sparse=T)

ngenes <- nrow(colon.data)
psedu_gene.ids <- paste0("ENSG0000", seq_len(ngenes))

耗时记录

> system.time(colon.data <- read.csv(file='GSM3270887_countTable_colonCreMin.txt.gz', sep="\t", header = T, row.names = 1))
   user  system elapsed 
 62.657   6.444  86.127

解决方案 1

writeMM(obj = colon.data, file="./matrix.mtx")
write.table(data.frame(psedu_gene.ids,rownames(colon.data)), file="./genes.tsv", 
            col.names=F,row.names = F, sep = "\t", quote=FALSE)
write(x = colnames(colon.data), file = "./barcodes.tsv")

解决方案 2

BiocManager::install("DropletUtils")
library(DropletUtils)

write10xCounts(path=getwd(), colon.data, gene.id=psedu_gene.ids,
               gene.symbol=rownames(colon.data), barcodes=colnames(colon.data))

推荐第一种解决方案，第二种解决方案需要安装额外的包，而且在输出路径参数上有些问题。

CellRanger矩阵读入时间

library(Seurat)
> system.time(Read10X(data.dir = "./"))
   user  system elapsed 
  7.294   0.488   8.027