如何直接用Seurat读取GEO中的单细胞测序表达矩阵

1

 常见的单细胞count matrix

  • Cell Ranger生成的raw count

    Cell Ranger (v3.0)中生成的文件除了bam文件外主要就是如下的三个表格文件(Seurat 的示例文件，2700个pbmc细胞单细胞测序）：

    我们可以利用head命令检查数据三个表格的内容。

    Barcodes通俗来讲就是每个细胞的代码，组成就是ATCG四个碱基排列组合成的不同的14个碱基组合；

    Gene.tsv或者features.tsv一般是基因的ensembl ID 和symbol

    matrix.mtx说白了就是每个细胞不同基因的表达矩阵，我们利用分别检查文件的开头和结尾：

这里我们可以发现其实就是2700个细胞不同基因的表达（第一列是基因的ID，用于与genes.tsv对应转换；第二列则是细胞的编号，匹配barcodes.tsv；第三列则是基因的表达量TPM）（没有表达的基因不做记录）这三组表格组合成。理解这三个表格组成后我们也不难发现，缺一不可的是matrx.mtx文件，而genes.tsv则一般是用于注释的基因组通用文件；而如果缺失barcodes.tsv的话，则可以根据matrix判断细胞数量自己“人为构建出”相应数量不同的barcode表格或者利用samtools从bam文件获取。当我们把这三个文件后存在一个独立文件夹后可以直接利用Seurat （v3.0）的Read10X()命令读取并构建成行名称为基因名，列名称为barcode序列(基因名x细胞）的表达矩阵（也就是SeuratObject)进行后续分析。如果我们只想从这三个表格直接整合成一个（基因名x细胞）的表达矩阵，可以利用以下代码完成：

library(Matrix)
matrix_dir = "~/filtered_feature_bc_matrix/hg19/"   ##根据实际文件夹进行修改
barcode.path <- paste0(matrix_dir, "barcodes.tsv")
features.path <- paste0(matrix_dir, "genes.tsv")
matrix.path <- paste0(matrix_dir, "matrix.mtx")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                          header = FALSE,
                          stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                          header = FALSE,
                          stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

  • 从公共数据库中获取的count matrix

    拿我们常见的GEO数据库为例，如果是上传到GEO数据的数据必须要上传处理后的数据（https://www.ncbi.nlm.nih.gov/geo/info/seq.html），这一方面方便其他研究人员直接更快速的获取或者验证最初的高通量测序，减少了下载SRA粗数据并进行重新比对的时间。

    一般来讲这些数据往往是整合好的一个count matrix，比如最新上传的一组造血干细胞单细胞测序数据（A 3D Atlas of Hematopoietic Stem and Progenitor Cell Expansion by Multi-dimensional RNA-Seq Analysis）（GSE120503），我们看到的处理后数据是单个文件，如下图所示：

    解压后我们得到只有一个叫做“GSM3402061_zebrafish_HSC_counts_change.merge.txt”的文件，而不是Cell Ranger输出的三个文件。

    我们检查一下文件的内容：

    其实这就是我们在上一步整合出的（基因 x 细胞）的表达矩阵，那么如果我们想直接利用Seurat导入这个表达矩阵进行后续分析该如何做呢？

2

 Count matrix导入Seur

对于上述的表达矩阵，我们不能直接使用Seurat的Read10X()函数进行读取，但是要进行后续分析我们可以直接把这个表达矩阵变成SeuratObject。这是一个R读取表格的基本操作：

setwd("/test/")  ##注意工作目录
library(Seurat)  ##version 3.0
library(dplyr)
new_counts <- read.table(file="/test/GSM3402061_zebrafish_HSC_counts_change.merge.txt")
head(new_counts)
mydata <- CreateSeuratObject(counts = new_counts, min.cells = 3, project = "mydata_scRNAseq")

通过以上两种操作我们就可以完成Cell Ranger产出数据与SeuratObject之间的互相转换。而利用这种简单的几行命令，我们可以较快的从他人上传好的数据中获取我们所需的信息（当然这需要我们充分相信合作者或者数据上传人对于数据处理的数据质量），节省了大量下载和处理数据的时间。

赞 (0)

相关推荐