|
本文首发于生信菜鸟团公众号，直达链接是https://mp.weixin.qq.com/s/NaZ5kz3ew2O01cFEnK8sXg 今天接到浙江大学的学徒求助，他在学习 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 )
因为这是一个小众需求，所以就不需要耗费时间去纠结代码运行效率了，反正运行起来是没有问题的。 单细胞转录组数据分析系列教程目录如下：
|
