scRNA基础分析-1:安装包、导入数据、过滤质控
2020-08-25 本文已影响0人
小贝学生信
已经学习了转录组基本入门知识,现在了解、学习下单细胞转录组(scRNA)相关知识。
依旧是站在巨人的肩膀上看世界,入门学习目前主要参考 微信公众号--生信会客厅的单细胞转录组分析教程 。笔记代码,示例文件均来自推文~
scRNA基础分析-1:安装包、导入数据、过滤质控 - 简书
scRNA基础分析-2:降维与聚类 - 简书
scRNA基础分析-3:鉴定细胞类型 - 简书
scRNA基础分析-4:细胞亚类再聚类、注释 - 简书
scRNA基础分析-5:伪时间分析 - 简书
scRNA基础分析-6:富集分析 - 简书
一、安装相关R包
install.packages('Seurat', dependencies=T)
#scRNA分析的基础包
install.packages("tidyverse", dependencies=T)
install.packages("BiocManager", dependencies=T)
BiocManager::install(c("SingleR","monocle", "DESeq2", "MAST"),ask = F,update = F)
BiocManager::install(c("clusterProfiler","DOSE","pheatmap"),ask = F,update = F)
BiocManager::install(c("org.Hs.eg.db","org.Mm.eg.db","org.Rn.eg.db"),ask = F,update = F)
二、数据导入,创建seurat对象
library(Seurat)
scRNA.counts <- Read10X(data.dir = "filtered_feature_bc_matrix")
#filtered_feature_bc_matrix为一个目录,里面有同一样本的三种文件。
2-1
scRNA = CreateSeuratObject(scRNA.counts)
#后续的分析大都以此为基础
table(scRNA@meta.data$orig.ident) #查看样本的细胞数量
#1222个
三、质控,过滤不感兴趣的细胞
- 过滤标准有细胞基因数量、mRNA数量、核糖体基因等。因此,前期我们需要对所有细胞的相关概况进行初步探索
探索
- 1、细胞的基因(nFeature_RNA)、mRNA数量(nCount_RNA)
head(scRNA@meta.data)
summary(scRNA@meta.data)
3-1
- 2、细胞中核糖体基因比例
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
head(scRNA@meta.data)
#计算红细胞比例,这一步没有看懂;红细胞基因比例?
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB.genes, rownames(scRNA@assays$RNA))
HB.genes <- rownames(scRNA@assays$RNA)[HB_m]
HB.genes <- HB.genes[!is.na(HB.genes)]
scRNA[["percent.HB"]]<-PercentageFeatureSet(scRNA, features=HB.genes)
- 上述可视化
#小提琴图
col.num <- length(levels(scRNA@active.ident))
# 1L
violin <- VlnPlot(scRNA,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
cols =rainbow(col.num),
pt.size = 0.01, #不需要显示点,可以设置pt.size = 0
ncol = 3)
3-2
#细胞基因数量与mRNA、核糖体基因数量的相关性
plot1 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot1|plot2
3-3
质控、过滤
- 一般根据基因数量和核糖体比例来过滤细胞
- 作者建议minGene=500,maxGene=4000,pctMT=15(线粒体基因比例低于15%)。
minGene=500
maxGene=4000
pctMT=15
scRNA <- subset(scRNA, subset = nFeature_RNA > minGene & nFeature_RNA < maxGene & percent.mt < pctMT)
col.num <- length(levels(scRNA@active.ident))
violin <-VlnPlot(scRNA,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
cols =rainbow(col.num),
pt.size = 0.1,
ncol = 3)
3-4
数据标准化与中心化(按列)
- 标准化是原始分数减去平均数然后除以标准差,差异基因分析会用到;
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
- 中心化是原始分数减去平均数,PCA主成分分析会用到。
scale.genes <- rownames(scRNA)
#如果自己数据集太大,也可到后面只对高变基因中心化
scRNA <- ScaleData(scRNA, features = scale.genes)
- 经过标准化和中心化之后,加上原始表达矩阵已经产生了三套基因表达数据
#原始表达矩阵
GetAssayData(scRNA,slot="counts",assay="RNA")
#标准化之后的表达矩阵
GetAssayData(scRNA,slot="data",assay="RNA")
#中心化之后的表达矩阵
GetAssayData(scRNA,slot="scale.data",assay="RNA")
最后保存结果,供后面的分析
saveRDS(scRNA, file="scRNA.rds")