单细胞专题 | 1.单细胞测序(10×genomics技术)的原理 单细胞专题 | 2.如何开始单细胞RNASeq数据分析 单细胞专题 | 3.单细胞转录组的上游分析-从BCL到FASTQ 单细胞专题 | 4.单细胞转录组的上游分析-从SRA到FASTQ 单细胞专题 | 5.单细胞转录组的上游分析-从FASTQ到count矩阵
1.数据读入 Cell Ranger生成的主要表格文件主要包括3个文件,分别是barcodes, features, matrix三个表格文件,其中barcodes文件存储的是细胞的信息,features存储的是基因的名称信息,matrix存储的是表达量的信息。还有一种数据是作者在GEO数据库直接提供表达矩阵(csv或txt) (1).读入csv文件的表达矩阵构建Seurat对象 Seurat需要的输入信息为表达量矩阵,矩阵行为基因,列为细胞。使用Seurat提供的Read10X函数可以很方便的将10x结果读入到R矩阵中。使用CreateSeuratObject生成Seurat对象,后续分析都是在该对象上进行操作。
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
# Load data
dir='./'
Sys.time()
raw.data <- read.csv(paste(dir,"Data_input/csv_files/S01_datafinal.csv", sep=""), header=T, row.names = 1)
Sys.time()
dim(raw.data)
raw.data[1:4,1:4]
head(colnames(raw.data))
# Load metadata
metadata <- read.csv(paste(dir,"Data_input/csv_files/S01_metacells.csv", sep=""), row.names=1, header=T)
head(metadata)
# Find ERCC's, compute the percent ERCC, and drop them from the raw data.
erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = TRUE)
percent.ercc <- Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
fivenum(percent.ercc)
ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE)
raw.data <- raw.data[-ercc.index,]
dim(raw.data)
有了表达矩阵,直接使用 CreateSeuratObject 函数即可,然后慢慢添加这个表达矩阵的一些其它外部属性,全部代码如下:
# Create the Seurat object with all the data (unfiltered)
main_tiss <- CreateSeuratObject(counts = raw.data)
# add rownames to metadta
row.names(metadata) <- metadata$cell_id
# add metadata to Seurat object
main_tiss <- AddMetaData(object = main_tiss, metadata = metadata)
main_tiss <- AddMetaData(object = main_tiss, percent.ercc, col.name = "percent.ercc")
# Head to check
head(main_tiss@meta.data)
# Calculate percent ribosomal genes and add to metadata
ribo.genes <- grep(pattern = "^RP[SL][[:digit:]]", x = rownames(x = main_tiss@assays$RNA@data), value = TRUE)
percent.ribo <- Matrix::colSums(main_tiss@assays$RNA@counts[ribo.genes, ])/Matrix::colSums(main_tiss@assays$RNA@data)
fivenum(percent.ribo)
main_tiss <- AddMetaData(object = main_tiss, metadata = percent.ribo, col.name = "percent.ribo")
main_tiss
# Filter cells so that remaining cells have nGenes >= 500 and nReads >= 50000
main_tiss_filtered <- subset(x=main_tiss, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
main_tiss_filtered
我们得到了main_tiss_filtered这个变量,是一个Seurat对象。
(2).读入10X标准输出的3个文件和融合多个样本数据
这三个文件指的是:barcodes.tsv, features.tsv, matrix.mtx。这个情况就比较好处理了,barcodes.tsv就是 cell id,features.tsv就是 gene id,matrix.mtx就是计数 counts 矩阵。
例如:
############### 10X标准输出的3个文件 ##############
# 列出当前目录下所有开头是GSM的文件
fs <- list.files('./data/GSE106273_RAW/','^GSM')
# 然后获取四个样本信息
library(stringr)
samples <- str_split(fs,'_', simplify = T)[,1]
# 设置一个循环,对每个样本信息做同样的事:
#(1)找到包含这个样本的文件(用grepl)
# (2)设置对应的目录名(str_split+paste)然后创建目录(用dir.create)
# (3)将文件放到对应目录(采用的是file.rename)并重命名文件
setwd("data/GSE106273_RAW/")
library(R.utils)
lapply(unique(samples),function(x){
y <- fs[grepl(x,fs)]
folder <- paste(str_split(y[1],'_',simplify = T)[,2:3],
collapse = '')
# 创建文件夹
dir.create(folder, recursive = T)
# 文件重命名:
file.rename(y[1],file.path(folder,"barcodes.tsv.gz"))
file.rename(y[2],file.path(folder,"genes.tsv.gz"))
file.rename(y[3],file.path(folder,"matrix.mtx.gz"))
# 解压,其实不解压也行
sapply(paste0(folder,c("/barcodes.tsv.gz","/genes.tsv.gz","/matrix.mtx.gz")),
gunzip)
})
# 批量读取成10X对象:
# 因为Read10X函数需要对目录进行操作,所以先把目录名提取出来
folders=list.files('./',pattern='[12]$')
folders
# [1] "G1" "G2" "L1" "L2" "NP1" "NP2" "PI1" "PI2"
# 然后使用lapply进行循环(之前出过一期apply系列函数教程,可以查阅一下,
# lapply是对列表或向量进行循环,而apply是对数据框或矩阵操作)
library(Seurat)
sceList <- lapply(folders,function(folder){
CreateSeuratObject(counts = Read10X(folder),
project = folder )
})
# 此时的sceList是一个包含8个10X对象的集合,下一步需要将其合并
sceList
# 合并:
sce_big <- merge(sceList[[1]],
y = c(sceList[[2]],sceList[[3]],sceList[[4]],
sceList[[5]],sceList[[6]],
sceList[[7]],sceList[[8]]),
add.cell.ids = folders,
project = "mouse8")
table(sce_big$orig.ident)
# G_1 G_2 L_1 L_2 NP_1 NP_2 PI_1 PI_2
# 2915 3106 5906 3697 2249 2127 1500 4306
# 保存:
save(sce_big,file = 'sce_big.Rdata') # 保存的数据
再如:
###### step1:导入数据 ######
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(data.table)
samples=list.files('GSE129139_RAW//')
samples
sceList = lapply(samples,function(pro){
#pro=samples[1]
folder=file.path('GSE129139_RAW',pro )
print(pro)
print(folder)
print(list.files(folder))
print( Sys.time() )
sce=CreateSeuratObject(counts = Read10X(folder),
project = pro )
print( Sys.time() )
return(sce)
})
sceList
samples
sce.all = merge(sceList[[1]], y = sceList[-1],
add.cell.ids = samples,
merge.data = TRUE)
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident)
library(stringr)
table(sce.all$orig.ident)
(3).读入h5格式文件
############### h5格式文件 ##############
# 读取单个
sce <- Read10X_h5(filename = "GSM4107899_LH16.3814_raw_gene_bc_matrices_h5.h5")
sce <- CreateSeuratObject(counts = sce)
# 批量读取并合并:
library(hdf5r)
setwd("../GSE138433_RAW/")
files <- list.files('./')
files
names <- str_split(files,'_', simplify = T)[,2]
# 批量读取:
sceList <- lapply(files, function(files){
CreateSeuratObject(counts = Read10X_h5(files),
project = names)
})
# 合并:
sce_big <- merge(sceList[[1]],
y = c(sceList[[2]],sceList[[3]],sceList[[4]],
sceList[[5]],sceList[[6]]),
add.cell.ids = names)
table(sce_big$orig.ident)
# LH16.3814 LH17.3222 LH17.3554 LH17.364 LH17.530 LH18.277
# 737280 737280 737280 737280 737280 737280
# 保存:
save(sce_big,file = 'sce_big.Rdata') # 保存的数据
(4). h5ad格式
需要安装SeuratDisk包,先将后h5ad格式转换为h5seurat格式,再使用LoadH5Seurat()函数读取Seurat对象。
############### h5ad格式文件 ##############
setwd("../GSE153643_RAW/")
# remotes::install_github("mojaveazure/seurat-disk")
library(SeuratDisk)
# 先将后`h5ad`格式转换为`h5seurat`格式:
Convert("GSM4648565_liver_raw_counts.h5ad", "h5seurat",
overwrite = TRUE, assay = "RNA")
# 再使用`LoadH5Seurat()`函数读取Seurat对象。
sce <- LoadH5Seurat("GSM4648565_liver_raw_counts.h5seurat")
2.Seurat对象
Seurat对象可以把它比作一个大容器,几乎存储了一切项目相关信息,包括每个细胞的barcodes,所有定量的基因和每个细胞的UMI矩阵。我们后面的大量分析包括对数据的降维、聚类分群、注释、等都是可以写入到Seurat对象来保存。
本文分享自 MedBioInfoCloud 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!