单细胞数据的导入与质控 - Seurat

######################################

  • 题目：单细胞数据的导入与质控 - Seurat
  • 语言：R
  • Author: YQ
    ######################################

step0.设置工作环境

  • 方法一：新建repository，在该repo内写R.scripts

  • 方法二：setwd()


setwd('/Users/apple/Desktop/Rsave/')
getwd()

创建以下目录

single_cell_rnaseq/
├── data
├── results
└── figures

step1.导入数据，创建计数矩阵

目前我碰到的有这几种情况；

- 情况一：三个文件

三个文件指的是“barcodes.tsv","features.tsv","matrix.mtx";
这个情况就比较好处理了，barcodes.tsv就是cell id，features.tsv就是gene id，matrix.mtx就是计数counts矩阵

##示例 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136916

# 设置默认路径
matrix_dir = '/Users/apple/Desktop/Rwork/Glasslung/' 

# 按路径读取三个文件
barcode.path<-paste0(matrix_dir,"barcodes.tsv")
genes.path<-paste0(matrix_dir,"features.tsv")
matrix.path<-paste0(matrix_dir,"matrix.mtx")

# readMM读取数据的values,columns,index
zebrafish.data <- readMM(file = matrix.path) ##mac上不能读压缩文件
gene.names = read.delim(genes.path,header = FALSE, stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,header = FALSE,stringsAsFactors = FALSE)
colnames(zebrafish.data) = barcode.names$V1
rownames(zebrafish.data) = gene.names$V2 ##把示例中的V1改成V2

# check矩阵
zebrafish.data[1:6, 1:6] ##check矩阵
dim(zebrafish.data) ##check矩阵

- 情况二：直接给了计数矩阵的csv

单个csv矩阵读入

# 读取csv格式的矩阵
data1 <- read.csv("/Users/quyue/Desktop/gene_bc_matrices.csv", header = T, row.names = 1)

# csv矩阵转换成数据框
datan = data.frame(data1)

# 数据框转换成稀疏矩阵matrix
dataan <- as(as.matrix(datan), "dgCMatrix")

多个csv矩阵读入


# 设置默认路径
matrix_dir = '/Users/quyue/Desktop/bioinfo/bonemarrowniche/0526data/'  # 为了避免重复输入默认路径，赋值给matrix_dir

# 读取多个csv
data1 <- read.csv(paste0(matrix_dir, "GSM3714747_chow1_filtered_gene_bc_matrices.csv"), header = T,row.names= 1)
data2 <- read.csv(paste0(matrix_dir, "GSM3714748_chow2_filtered_gene_bc_matrices.csv"), header = T,row.names= 1)
data3 <- read.csv(paste0(matrix_dir, "GSM3714749_chow3_filtered_gene_bc_matrices.csv"), header = T,row.names= 1)

# csv合并成数据框
datan = data.frame(data1,data2,data3)

# 数据框转换成稀疏矩阵matrix
dataan <- as(as.matrix(datan), "dgCMatrix")

- 情况三：直接给了计数矩阵的txt

单个txt矩阵读入

# 读取txt格式的矩阵
data1 <- read.table("/Users/quyue/Desktop/GSM2915579_5FU-cntrl-col23.txt", header = T, row.names = 1)

# txt矩阵转换成数据框
datan = data.frame(data1)

# 数据框装转换成稀疏矩阵matrix
dataan <- as(as.matrix(datan),"dgCMatrix")

多个txt矩阵读入


# 设置默认路径
matrix_dir = "/Users/quyue/Desktop/"

# 读取多个txt格式的矩阵
data1 = read.table(paste0(matrix_dir,"GSM2915579_5FU-cntrl-col23-1.txt"), header = T, row.names = 1)

data2 = read.table(paste0(matrix_dir,"GSM2915579_5FU-cntrl-col23-2.txt"), header = T, row.names = 1)

data3 = read.table(paste0(matrix_dir,"GSM2915579_5FU-cntrl-col23-3.txt"), header = T, row.names = 1)

dataa = cbind(data1, data2, data3) #行数不同，会生成NaN值

# txt矩阵转换成数据框
datan = data.frame(dataa)

# 数据框转换成稀疏矩阵matrix
dataan <- as(as.matrix(datan),"dgCMatrix") # 此步骤是否必要有待验证

- 情况四：xslx文件

其实就是用R读取excel文件，把xslx转化成csv再读取

# 加载openxslx这个R包
library('openxslx')

# 读入xslx格式矩阵
a <- read.xslx("test.xslx",sheet = 1) # 文件名+sheet序号

# 把xslx格式矩阵转换成csv格式矩阵
# 获取excel中工作簿的名称
sheetnames<-getSheetNames('test.xlsx')

#把每个工作薄的数据按照'工作薄名称.csv'的名称存储
for(iin(1:length(sheetnames))){

write.table(read.xlsx('/Users/quyue/Desktop/test.xlsx',sheet=i),paste(sheetnames[i],'.csv',sep=''),row.names=F,sep=',')

} # 有待验证，不是我自己写的

# 读取csv格式矩阵
data1 <- read.csv("/Users/quyue/Desktop/test.csv", header = T, row.names = 1)

# 把csv格式矩阵转换成数据框
datan = data.frame(data1)

# 把数据框转换成稀疏矩阵matrix
dataan <- as(as.matrix(datan), "dgCMatrix")

- 情况五：直接给了已经构建好seurat对象的R.data

# 直接load R.data
load(file = "/Users/quyue/Desktop/bioinfo/bonemarrowniche/RNAMagnetDataBundle/NicheData10x.rda")
    > Loading required package: Seurat

# load完以后，看到底是什么数据类型，并查看数据
NicheData10x
    > An object of class Seurat # 已经是seurat对象了
    16701 features across 7497 samples within 1 assay 
    Active assay: RNA (16701 features, 2872 variable features)
     2 dimensional reductions calculated: pca, tsne

View(NicheData10x) # 显示seurat的具体参数

step2.构建矩阵为seurat对象

找到一个解释Seurat中所有函数的网站，比R自带的要清楚很多 rdrr

工作流程

  • generation of the count matrix 生成计数矩阵
    格式化读取、分离样本、映射(mapping)、定量(quantification);

  • quality control of the raw counts 原始矩阵质量控制
    过滤掉质量差的细胞

  • clustering of filtering counts 过滤后计数的聚类
    将转录活动相似的细胞归为一类（细胞类型=不同的聚类）

  • Marker identification 标记识别
    识别每个细胞群的基因标记（marker）

  • Optional downstream steps 其他可选的下游步骤

2.1 稀疏矩阵到seurat对象

2.1.1 使用ReadMM()或者Read10x()读取单个样本数据
  • 函数用法说明：

1.readMM():

input:"xx.mtx"

备注：readMM is the function of Matrix packages, it changes the standard matrix into sparse matrix.Of note,features.tsv and barcodes.tsv should be library first, and then combine sparse matrix、features.tsv and barcodes.tsv form a counts matrix with cell id and gene id.（详见step1情况一）

2.read10X():

input:"raw_feature_bc_matrix"(files from CellRanger output)

备注：read10X id the function of Seurat,its input from CellRanger output(10X genomics data的专用软件Cell Ranger的output会有一个out目录，"raw_feature_bc_matrix"在该目录中)

library(Seurat)

# 按1中方法完成数据读入到稀疏矩阵输出的过程
# sparse matrix稀疏矩阵
ctrl_counts <- Read10x(data.dir = "data/ctrl_raw_feature_bc_matrix")


# 将稀疏矩阵转换成一个Seurat
ctrl <- CreateSeuratObject(counts = ctrl_counts,
                          min.features = 100)

# 查看Seurat对象的元数据
head(ctrl@meta.data)

## 备注：元数据 meta.data代表什么？
orig.ident # 通常包含所知的样品名字，默认为“SeuratProject”
nCount_RNA # 每个细胞的UMI数目
nFeature_RNA # 每个细胞所检测到的基因数目
2.1.2 使用for循环读入多个样本数据
  • R中for 循环的语法结构
for (variable in input){
    command1
    command2
    command3
}
  • R中for循环实现读入多个样本数据
## 目的：for loop iterate 两个样本的文件
## 思路：执行三个命令 (1) 使用Read10x()读取计数数据；（2）使用CreateSeuratObject()从读取的数据中创建Seurat对象；（3）根据样本名将Seurat对象赋值给一个新的变量（好处是不会覆盖上一次迭代中创建的Seurat对象）

mtx_dir = "/Users/quyue/Desktop/"

for (files in c("ctrl_raw_feature_bc_matrix","stim_raw_feature_bc_matrix")){
    seurat_data <- Read10X(data.dir = paste0(mtx_dir, file)) # (1)
    seurat_obj <- CreateSeuratObject(counts = seurat_data,
                                    min.features = 100,
                                    project = file) #(2)
    assign(file, seurat_obj) #(3)
}

2.2 查看seurat对象元数据&merge

再完成2.1已经创建好了Seurat对象以后，快速查看元数据以了解其大概；

## 查看新Seurat对象的元数据
head(ctrl_raw_feature_bc_matrix@meta.data)
head(stim_raw_feature_bc_matrix@meta.data)

接着，需要把多个seurat对象合并成一个，这样比较方便我们同时为每个样品进行质控，且方便我们更容易地比较所有样本的数据质量，用merge()函数执行该操作；

## 创建一个合并的Seurat对象
## 用add.cell.id参数将样品特异的前缀添加到每个细胞ID，用以区分不同样本中的同个细胞ID
merged_seurat <- merge(x = ctrl_raw_feature_bc_matrix,
                       y = stim_raw_feature_bc_matrix,
                       add.cell.id = c("ctrl","stim"))

## 查看合并对象是否有对应的样本前缀
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)

step3. 进行质控流程（QC and pre-processing）

3.1 熟悉质控指标

通常根据以下四个指标筛选：

# 每个细胞的read counts数目（UMI指的是molecular ID）
nCount_RNA  # (> 1000)

# 每个细胞所检测到的基因数目
nFeature_RNA # (> 1000)

# 线粒体基因比例 mitoRatio
mitoRatio  # (< 5%)

# 单位UMI(read count=molecular id)检测到的gene数目的比例(这个数据让我们知道数据的复杂性)
GenesPerUMI

3.2 标准的质控分析流程 - Standard pre-processing workflow

3.2.1 计算质控指标
  • 需要计算的主要是mitoRatio和GenesPerUMI（因为nCount_RNA和nFeature_RNA在meta.data中）

1.计算mitoRatio

主要是利用的Seurat的PercentageFeatureSet()功能，这个函数将使用一个模式(pattern)搜索基因标识符，对于每一列（细胞），它将选取特征基因的计数之和，除以所有基因的计数之和

## count mitoRatio with PercentageFeatureSet() 并添加 mito Ratio这一项到 meta.data当中

# 方法一
## The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 方法二
pbmc$mitoRatio <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc$mitoRatio <- pbmc@meta.data$mitoRatio / 100

2.计算log10GenesPerUMI（option）

GenesPerUMI = nFeature_RNA / nCount_RNA
这个指标代表的意思是 “单位read counts的gene数目的比例” 反映的是数据的复杂度；常用于评估ctrl和stim基因的整体复杂性

## count log10GenesPerUMI and add it to seurat object metadata
pbmc$log10GenesPerUMI <- log10(pbmc$nFeature_RNA)/log10(pbmc$nCount_RNA)
3.2.2 可视化质控指标
## 一般可视化以下三个参数即可："nFeature_RNA", "nCount_RNA", "percent.mt"
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size=0.05)

## 可视化log10GenesPerUMI（option）
# 通过可视化每一个UMI检测到的基因数来可视化基因表达的整体复杂性
metadata %>%
    ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
    geom_density(alpha = 0.2) +
    theme_classic() +
    geom_vline(xintercept = 0.8)
3.2.3 筛选 Filtering（option - 因为一般给的是经过筛选过后的gene mtx）
3.2.3.1 细胞水平的筛选（Cell-level filtering）
  • 参考阈值
nUMI > 500
nGene > 250
log10GenesPerUMI > 0.8
mitoRatio < 0.2

细胞水平的过滤, 主要是用subset()函数过滤

# 使用选择的阈值筛掉低质量读写 - 这些将会随实验而改变
filtered_pbmc <- subset(x = pbmc, subset= (nUMI >= 500) & 
                                          (nGene >= 250) & 
                                          (log10GenesPerUMI > 0.80) & 
                                          (mitoRatio < 0.20))
3.2.3.2 基因水平的筛选（Gene-level filtering）

在我们的数据中，会有许多基因的计数为零。这些基因会显著降低一个细胞的平均表达量，因此我们将从数据中删除它们。

主要删除两类基因，
（1）删除所有细胞中零表达的基因；

（2）如果一个基因只在少数几个细胞中表达，那么它的意义并不是特别大，因为它仍然会把没有表达它的其他所有细胞的平均表达量降下来。对于我们的数据，我们选择只保留在10个或更多细胞中表达的基因。

# 提取计数
counts <- GetAssayData(object = filtered_pbmc, slot = "counts")

# 根据在每个细胞的计数是否大于0为每个基因输出一个逻辑向量
nonzero <- counts > 0

# 将所有TRUE值相加，如果每个基因的TRUE值超过10个，则返回TRUE。
keep_genes <- Matrix::rowSums(nonzero) >= 10

# 仅保留那些在10个以上细胞中表达的基因
filtered_counts <- counts[keep_genes, ]

# 重新赋值给经过过滤的Seurat对象
filtered_pbmc <- CreateSeuratObject(filtered_counts, meta.data = filtered_pbmc@meta.data)
3.2.4 重新评估质量指标 (Re-assess QC metrics)

在进行过滤后，建议回过头来再看看指标，确保你的数据符合你的预期，对下游分析有好处。

3.2.5 保存筛选过的细胞 (Saving filtered cells)

我们将保存筛选后的细胞对象用于聚类 (clustering) 和标记识别 (marker identification)，将过滤以后的Seurat数据创建为.RData对象

# 创建.RData对象以方便在任何时候导入
save(filtered_pbmc, file="data/pbmc_seurat_filtered.RData")
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,491评论 5赞 459
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 81,856评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 141,745评论 0赞 319
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,196评论 1赞 263
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,073评论 4赞 355
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,112评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,531评论 3赞 381
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,215评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,485评论 1赞 290
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,578评论 2赞 309
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,356评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,215评论 3赞 312
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,583评论 3赞 299
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,898评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,497评论 2赞 341
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,697评论 2赞 335