1. 演示数据和包的加载
# if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")## BiocManager::install("scater")# BiocManager::install("scRNAseq")library(scater)library(scRNAseq)# 创建SingleCellExperiment对象example_sce <- ZeiselBrainData()
直接创建SingleCellExperiment对象
SingleCellExperiment(
...,
reducedDims = list(),
altExps = list(),
rowPairs = list(),
colPairs = list(),
mainExpName = NULL
)
2.质控作图，细胞与基因的筛选
## Diagnostic plots for quality controlexample_sce <- addPerCellQC(example_sce,subsets=list(Mito=grep("mt-", rownames(example_sce))))plotColData(example_sce, x = "sum", y="detected", colour_by="tissue")plotColData(example_sce, x = "sum", y="subsets_Mito_percent",other_fields="tissue") + facet_wrap(~tissue)plotHighestExprs(example_sce, exprs_values = "counts")## remove damaged cells and poorly sequenced libraries.自己制定筛选条件# 查看是否有spike-insrownames(example_sce)[grep("^ERCC-", rownames(example_sce))]## 选择总counts数大于5,000，表达的基因数大于500的细胞。keep.total <- example_sce$sum > 5e3#table(keep.total)keep.n <- example_sce$`total mRNA mol` > 500# table(keep.n)# 根据设定的条件进行过滤example_sce <- example_sce[,keep.total & keep.n]#dim(example_sce)## 选择至少在三个细胞中表达的基因keep_feature <- nexprs(example_sce, byrow=TRUE) >= 3# table(keep_feature)example_sce <- example_sce[keep_feature,]# log转化并归一化example_sce <- logNormCounts(example_sce)# plotExplanatoryVariables(example_sce)vars <- getVarianceExplained(example_sce,variables=c("tissue", "total mRNA mol", "sex", "age"))head(vars)class(vars)plotExplanatoryVariables(vars)
3. 基因表达值可视化作图
## Visualizing expression valuescolnames(example_sce@colData)## 不同分组中一个或多个基因的表达量plotExpression(example_sce, rownames(example_sce)[1:3])# 分组统计plotExpression(example_sce, rownames(example_sce)[1:3],x = "level1class")# 加颜色plotExpression(example_sce, rownames(example_sce)[1:3],x = "level1class", colour_by="tissue")## 基因表达量之间的相关性plotExpression(example_sce, rownames(example_sce)[1:3],x = rownames(example_sce)[10])
4. 降维并作图
## Dimensionality reduction# calculatePCA(# x,# ncomponents = 50,# ntop = 500,# subset_row = NULL,# scale = FALSE,# transposed = FALSE,# BSPARAM = bsparam(),# BPPARAM = SerialParam()# )example_sce <- runPCA(example_sce)# uses the top 500 genes with the highest variances to compute the first PCsstr(reducedDim(example_sce, "PCA"))plotReducedDim(example_sce, dimred = "PCA",colour_by = "level1class")plotPCA(example_sce,colour_by = "level1class")plotPCA(example_sce, colour_by="Mog")# PCA1-4 两两之间作图plotPCA(example_sce, ncomponents = 4,colour_by = "level1class")# 指定基因，指定名称，指定纬度数据example_sce <- runPCA(example_sce, name="PCA2",subset_row=rownames(example_sce)[1:1000],ncomponents=25)str(reducedDim(example_sce, "PCA2"))set.seed(1000)example_sce <- runTSNE(example_sce, perplexity=10) # 默认ncomponents = 2example_sce <- runTSNE(example_sce)#numeric; Perplexity parameter (should not be bigger than# 3 * perplexity < nrow(X) - 1head(reducedDim(example_sce, "TSNE"))#dev.new()# 颜色显示分组# plotTSNE(example_sce, colour_by = "level1class")plotReducedDim(example_sce, dimred = "TSNE",colour_by = "level1class")# 按基因表达量显示颜色plotTSNE(example_sce, colour_by = "Snap25")# "Snap25" %in% rownames(example_sce)
5.ggcells，ggfeatures作图
ggcells，ggfeatures函数将智能地从colData（）、assays（）、altExps（）或reducedDims（）中检索字段，以创建单个data.frame供立即使用。
# These functions generate a data.frame from the contents of# a SingleCellExperiment and pass it to ggplot.ggcells(example_sce, mapping=aes(x=level1class, y=Snap25)) +geom_boxplot() +facet_wrap(~tissue)# 显示Snap25表达量并按组织来源分别显示的TSNE图。ggcells(example_sce, mapping=aes(x=TSNE.1, y=TSNE.2, colour=Snap25)) +geom_point() +stat_density_2d() +facet_wrap(~tissue) +scale_colour_distiller(direction=1)#examine the relationship between the size factors# on the normalized gene expressionggcells(example_sce, mapping=aes(x=sizeFactor, y=Actb)) +geom_point() +geom_smooth()# 使列名（cell barcodes）有效colnames(example_sce) <- make.names(colnames(example_sce))make.names(".2way") # 输出[1] "X.2way"ggfeatures(example_sce, mapping=aes(x=featureType, y=X1772062111_E06)) +geom_violin()
参考
/packages/release/bioc/vignettes/scater/inst/doc/overview.html
如果觉得《scater分析单细胞转录组数据代码》对你有帮助，请点赞、收藏，并留下你的观点哦！
