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

生信菜鸟团

发布于 2020-03-27 13:57:37

4.5K0

发布于 2020-03-27 13:57:37

文章被收录于专栏：生信菜鸟团

今天接到浙江大学的学徒求助，他在学习 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 general\n")
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 )

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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-03-16，如有侵权请联系 cloudcommunity@tencent.com 删除

typescript

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

typescript

登录后参与评论

0 条评论

热度