表达矩阵逆转为10X的标准输出3个文件

2020 年 3 月 27 日
筆記

今天接到浙江大学的学徒求助，他在学习 TooManyCellsR 包和 too-many-cells 软件的过程中遇到了一个很有趣的问题，就是这个软件的输入必须是 cellranger 的三个结果文件，matrix.mtx，barcodes.tsv 和 genes.tsv。而有些公共数据并不会提供3个数据，比如: SE117988_raw.expMatrix_PBMC.csv.gz , 就是 10x的表达矩阵。

我们会使用下面的代码来读取这个表达矩阵，进行Seurat分析。

library(data.table)  # 这个表达矩阵其实是10x的，不过可以演示  a=fread('GSE117988_raw.expMatrix_PBMC.csv.gz',header = TRUE)  length(a$V1)  length(unique(a$V1))  hg=a$V1  dat=a[,2:ncol(a)]  rownames(dat)=hg  hg[grepl('^MT-',hg)]  colnames(dat)  rownames(dat)  meta=as.data.frame(colnames(dat))  colnames(meta)=c('cell name')  rownames(meta)=colnames(dat)  head(meta)  ## 前面大量的代码，都是数据预处理  library(Seurat)  dat[1:4,1:4]  class(dat)  # 重点是构建 Seurat对象  pbmc  <- CreateSeuratObject(counts = dat,                              meta.data = meta,                              min.cells = 3, min.features = 200,project = '10x_PBMC')    pbmc

这样的例子，就是作者不提供cellranger 的三个结果文件，matrix.mtx，barcodes.tsv 和 genes.tsv，不过大部分其它数据集，比如 GSE128033 和 GSE135893，你随便下载其中一个，就能看到每个样本都是走流程拿到10x单细胞转录组数据的3个文件的表达矩阵。

2.2M Mar  8  2019 GSM3660655_SC94IPFUP_barcodes.tsv.gz  259K Mar  8  2019 GSM3660655_SC94IPFUP_genes.tsv.gz   26M Mar  8  2019 GSM3660655_SC94IPFUP_matrix.mtx.gz    2.2M Mar  8  2019 GSM3660656_SC95IPFLOW_barcodes.tsv.gz  259K Mar  8  2019 GSM3660656_SC95IPFLOW_genes.tsv.gz   31M Mar  8  2019 GSM3660656_SC95IPFLOW_matrix.mtx.gz

下游处理的时候，一定要保证这3个文件同时存在，而且在同一个文件夹下面。

示例代码是：

rm(list=ls())  options(stringsAsFactors = F)  library(Seurat)  sce1 <- CreateSeuratObject(Read10X('../10x-results/WT/'),                            "wt")

重点就是 Read10X 函数读取文件夹路径，比如：../10x-results/WT/ ，保证文件夹下面有3个文件。

因为10x单细胞转录组表达矩阵里面的0值非常多，所以换成3个文件存储更节省空间。

本质上仍然是一个表达矩阵而已，如果你都有了表达矩阵，就没必要去想那3个文件了。

自己制作那3个文件也不是不可以

极特殊情况下，比如使用 too-many-cells 软件，这个软件的输入必须是 cellranger 的三个结果文件，matrix.mtx，barcodes.tsv 和 genes.tsv，不得不进行格式转换了。

首先需要解析3个文件的规律

前两个文件比较好理解，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  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

但是matrix.mtx，就稍微复杂一点，仔细看

    2049 barcodes.tsv     33694 genes.tsv   1878960 matrix.mtx

读取这3个文件，进入R里面：

rm(list=ls())  options(stringsAsFactors = F)  library(Seurat)  sce <- CreateSeuratObject(Read10X('SRR7722939/'),                             "SRR7722939")    ct=GetAssayData(object = sce, assay = "RNA", slot = "counts")  ct=as.data.frame(ct)  fivenum(apply(ct,2,function(x) sum(x>0)))  table(ct>0)

一些简单的摸索，如下：

> fivenum(apply(ct,2,function(x) sum(x>0)))  AGTGAGGCAAGCTGTT TGGTTAGAGGTGCACA CACAGTATCCACGAAT CCTATTAGTCCCGACA TCAGCTCCATGTCTCC                57              744              874             1020             3777  > table(ct>0)       FALSE     TRUE  67160049  1878957  > ct[1:5,1:4]                    AAACCTGAGCGAAGGG AAACCTGAGGTCATCT AAACCTGAGTCCTCCT AAACCTGCACCAGCAC  hg38-RP11-34P13.3                0                0                0                0  hg38-FAM138A                     0                0                0                0  hg38-OR4F5                       0                0                0                0  hg38-RP11-34P13.7                0                0                0                0  hg38-RP11-34P13.8                0                0                0                0  >

经过对比，可以看到matrix.mtx 文件记录的就是 1878957 个非0值，表达矩阵的行名和列名是有顺序的。制作barcodes.tsv 和 genes.tsv，的代码非常简单：

write.table(data.frame(rownames(ct),rownames(ct)),file = 'genes.tsv',              quote = F,sep = 't',              col.names = F,row.names = F)  write.table(colnames(ct),file = 'barcodes.tsv',quote = F,              col.names = F,row.names = F)

而matrix.mtx 文件是3列，第一列是行号，第二列是列号，第三列是基因表达量，而且仅仅是列出有表达量的基因即可，也就是说67160049 个0值都删掉，仅仅是显示 1878957 个非0值即可。这就是为什么10X公司采取这个方式来存储表达矩阵了，的确是非常的压缩空间啊！

原理很简单，但是代码运行速度很考验编程底层能力

首先写一个头信息

file="matrix.mtx"  sink(file)  cat("%%MatrixMarket matrix coordinate integer generaln")  cat("%n")  cat(paste(nrow(ct),ncol(ct),sum(ct>0),"n"))  sink()

再写入表达量信息

tmp=ct[1:5,1:4]  tmp  tmp=do.call(rbind,lapply(1:ncol(ct),function(i){    return(data.frame(row=1:nrow(ct),                      col=i,                      exp=ct[,i]))  }) )  tmp=tmp[tmp$exp>0,]  head(tmp)  write.table(tmp,file = 'matrix.mtx',quote = F,              col.names = F,row.names = F,append = T )

因为这是一个小众需求，所以就不需要耗费时间去纠结代码运行效率了，反正运行起来是没有问题的。如果你也正好有这个需求，学习一下，然后打赏吧！