表達(dá)矩陣逆轉(zhuǎn)為10X的標(biāo)準(zhǔn)輸出3個文件 | 生信菜鳥團(tuán)

頭頭了不起 2020-07-13

展開全文

本文首發(fā)于生信菜鳥團(tuán)公眾號，直達(dá)鏈接是https://mp.weixin.qq.com/s/NaZ5kz3ew2O01cFEnK8sXg

今天接到浙江大學(xué)的學(xué)徒求助，他在學(xué)習(xí) TooManyCellsR 包和 too-many-cells 軟件的過程中遇到了一個很有趣的問題，就是這個軟件的輸入必須是 cellranger 的三個結(jié)果文件，matrix.mtx，barcodes.tsv 和 genes.tsv。而有些公共數(shù)據(jù)并不會提供3個數(shù)據(jù)，比如: SE117988_raw.expMatrix_PBMC.csv.gz , 就是 10x的表達(dá)矩陣。

我們會使用下面的代碼來讀取這個表達(dá)矩陣，進(jìn)行Seurat分析。

library(data.table)
# 這個表達(dá)矩陣其實是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)
## 前面大量的代碼，都是數(shù)據(jù)預(yù)處理
library(Seurat)
dat[1:4,1:4]
class(dat)
# 重點是構(gòu)建 Seurat對象
pbmc <- CreateSeuratObject(counts = dat,
 meta.data = meta,
 min.cells = 3, min.features = 200,project = '10x_PBMC')

pbmc

這樣的例子，就是作者不提供cellranger 的三個結(jié)果文件，matrix.mtx，barcodes.tsv 和 genes.tsv，不過大部分其它數(shù)據(jù)集，比如 GSE128033 和 GSE135893，你隨便下載其中一個，就能看到每個樣本都是走流程拿到10x單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的3個文件的表達(dá)矩陣。

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 函數(shù)讀取文件夾路徑，比如：../10x-results/WT/ ，保證文件夾下面有3個文件。

因為10x單細(xì)胞轉(zhuǎn)錄組表達(dá)矩陣?yán)锩娴?值非常多，所以換成3個文件存儲更節(jié)省空間。

本質(zhì)上仍然是一個表達(dá)矩陣而已，如果你都有了表達(dá)矩陣，就沒必要去想那3個文件了。

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

極特殊情況下，比如使用 too-many-cells 軟件，這個軟件的輸入必須是 cellranger 的三個結(jié)果文件，matrix.mtx，barcodes.tsv 和 genes.tsv，不得不進(jìn)行格式轉(zhuǎn)換了。

首先需要解析3個文件的規(guī)律

前兩個文件比較好理解，barcodes.tsv 和 genes.tsv，就是表達(dá)矩陣的行名和列名：

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，就稍微復(fù)雜一點，仔細(xì)看：

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

讀取這3個文件，進(jìn)入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
>

經(jīng)過對比，可以看到matrix.mtx 文件記錄的就是 1878957 個非0值，表達(dá)矩陣的行名和列名是有順序的。制作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列，第一列是行號，第二列是列號，第三列是基因表達(dá)量，而且僅僅是列出有表達(dá)量的基因即可，也就是說67160049 個0值都刪掉，僅僅是顯示 1878957 個非0值即可。這就是為什么10X公司采取這個方式來存儲表達(dá)矩陣了，的確是非常的壓縮空間??！

原理很簡單，但是代碼運行速度很考驗編程底層能力

首先寫一個頭信息

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()

再寫入表達(dá)量信息

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 )