单细胞表达矩阵格式转换
2020-04-25 本文已影响0人
seqyuan
作者：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
