单细胞转录组(scRNA-seq)分析01 | Scater包的使用

一、介绍

  • 用于单细胞RNA-seq数据
  • 提供严格的质量控制:将原始测序读数处理为可用于下游分析的高质量表达数据集
  • 提供了丰富的绘图工具套件
  • R包地址:http://bioconductor.org/packages/scater

二、工作流

mark

三、常用函数

plotColData

# 导入包
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))

# 载入示例数据
data("sc_example_counts")
data("sc_example_cell_info")

# 构建 SingleCellExperiment 对象
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)

# 计算 SingleCellExperiment 对象中每个特征和细胞的质控标准
example_sce <- calculateQCMetrics(example_sce)

# 计算 SingleCellExperiment 对象中read计数矩阵的归一化表达值
example_sce <- normalize(example_sce)

plotColData(example_sce, y = "total_features_by_counts",
            x = "log10_total_counts", colour_by = "Mutation_Status")
mark
plotColData(example_sce, y = "total_features_by_counts",
            x = "log10_total_counts", colour_by = "Mutation_Status",
            size_by = "Gene_0001", shape_by = "Treatment")
mark
plotColData(example_sce, y = "Treatment",
            x = "log10_total_counts", colour_by = "Mutation_Status")
mark
plotColData(example_sce, y = "total_features_by_counts",
            x = "Cell_Cycle", colour_by = "Mutation_Status")
mark

plotExplanatoryVariables

解释变量(ExplanatoryVariables):https://www.statisticshowto.datasciencecentral.com/explanatory-variable/

解释变量是一种自变量。这两个术语通常可互换使用。但是,两者之间的细微差别。当一个变量是独立的,它不影响在所有的任何其他变量。当变量不是独立的时候,它是一个解释变量。 它在临床研究中非常重要。对于大多数情况,特别是在统计数据中,这两个术语基本相同。

假设您有两个变量来解释体重增加:快餐和苏打水。虽然你可能认为吃快餐和喝苏打水是相互独立的,但它们并不是真的。那是因为快餐店鼓励你在用餐时买苏打水。如果你停在某个地方买苏打水,那里经常会有很多快餐选择,比如热狗。虽然这些变量并非完全相互独立,但它们确实会对体重增加产生影响。它们被称为解释变量,因为它们可能为体重增加提供一些解释。

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
plotExplanatoryVariables(example_sce)
mark

plotExpression

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- calculateQCMetrics(example_sce)
sizeFactors(example_sce) <- colSums(counts(example_sce))
example_sce <- normalize(example_sce)

# 前十五个基因的表达值
plotExpression(example_sce, 1:15)
mark
plotExpression(example_sce, c("Gene_0001", "Gene_0004"), x="Mutation_Status")
mark
plotExpression(example_sce, c("Gene_0001", "Gene_0004"), x="Gene_0002")
mark
plotExpression(example_sce, 1:6, colour_by = "Mutation_Status")
mark
plotExpression(example_sce, 1:6, colour_by = "Mutation_Status",
               shape_by = "Treatment", size_by = "Gene_0010")
mark
plotExpression(example_sce, 1:4, "Gene_0004", show_smooth = TRUE)
mark

plotExprsFreqVsMean

表达频率(即表达细胞的百分比)Vs SingleCellExperiment对象中每个特征的平均表达水平

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
example_sce <- calculateQCMetrics(example_sce,
                                  feature_controls = list(set1 = 1:500))
plotExprsFreqVsMean(example_sce)
mark
plotExprsFreqVsMean(example_sce, size_by = "is_feature_control")
mark

plotExprsVsTxLength

Plot mean expression values for all features in a SingleCellExperiment object against transcript length values.

data("sc_example_counts")
data("sc_example_cell_info")
rd <- DataFrame(gene_id = rownames(sc_example_counts),
                feature_id = paste("feature", rep(1:500, each = 4), sep = "_"),
                median_tx_length = rnorm(2000, mean = 5000, sd = 500),
                other = sample(LETTERS, 2000, replace = TRUE)
)
rownames(rd) <- rownames(sc_example_counts)
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info, rowData = rd
)
example_sce <- normalize(example_sce)
plotExprsVsTxLength(example_sce, "median_tx_length")
mark
plotExprsVsTxLength(example_sce, "median_tx_length", show_smooth = TRUE)
mark
plotExprsVsTxLength(example_sce, "median_tx_length", show_smooth = TRUE,
                    colour_by = "other", show_exprs_sd = TRUE)
mark
## using matrix of tx length values in assays(object)
mat <- matrix(rnorm(ncol(example_sce) * nrow(example_sce), mean = 5000,
                    sd = 500), nrow = nrow(example_sce))
dimnames(mat) <- dimnames(example_sce)
assay(example_sce, "tx_len") <- mat
plotExprsVsTxLength(example_sce, "tx_len", show_smooth = TRUE,
                    length_is_assay = TRUE, show_exprs_sd = TRUE)
mark
## using a vector of tx length values
plotExprsVsTxLength(example_sce,
                    data.frame(rnorm(2000, mean = 5000, sd = 500)))
mark

plotHeatmap

Create a heatmap of expression values for each cell and specified features in a SingleCellExperiment
object.

example(normalizeSCE) # borrowing the example objects in here.
plotHeatmap(example_sce, features=rownames(example_sce)[1:10])
mark
plotHeatmap(example_sce, features=rownames(example_sce)[1:10],
            center=TRUE, symmetric=TRUE)
mark
plotHeatmap(example_sce, features=rownames(example_sce)[1:10],
            colour_columns_by=c("Mutation_Status", "Cell_Cycle"))
mark

plotHighestExprs

Plot the features with the highest average expression across all cells, along with their expression in
each individual cell.

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- calculateQCMetrics(example_sce,
                                  feature_controls = list(set1 = 1:500)
)
plotHighestExprs(example_sce, colour_cells_by ="total_features_by_counts")
mark
plotHighestExprs(example_sce, controls = NULL)
mark
plotHighestExprs(example_sce, colour_cells_by="Mutation_Status")
mark

plotPlatePosition

Plots cells in their position on a plate, coloured by metadata variables or feature expression values
from a SingleCellExperiment object.

## prepare data
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
example_sce <- calculateQCMetrics(example_sce)
## define plate positions
example_sce$plate_position <- paste0(
  rep(LETTERS[1:5], each = 8),
  rep(formatC(1:8, width = 2, flag = "0"), 5)
)
## plot plate positions
plotPlatePosition(example_sce, colour_by = "Mutation_Status")
mark
plotPlatePosition(example_sce, shape_by = "Treatment", colour_by = "Gene_0004")
mark
plotPlatePosition(example_sce, shape_by = "Treatment", size_by = "Gene_0001",
                  colour_by = "Cell_Cycle")
mark

plotQC

Produce QC diagnostic plots

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info)
example_sce <- normalize(example_sce)
example_sce <- calculateQCMetrics(example_sce)
plotQC(example_sce, type="high", colour_cells_by="Mutation_Status")
mark

plotReducedDim

Plot cell-level reduced dimension results stored in a SingleCellExperiment object.

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
example_sce <- runPCA(example_sce, ncomponents=5)
plotReducedDim(example_sce, "PCA")
mark
plotReducedDim(example_sce, "PCA", colour_by="Cell_Cycle")
mark
plotReducedDim(example_sce, "PCA", colour_by="Gene_0001")
mark
plotReducedDim(example_sce, "PCA", ncomponents=5)
mark
plotReducedDim(example_sce, "PCA", ncomponents=5, colour_by="Cell_Cycle",
               shape_by="Treatment")
mark

plotRLE

Produce a relative log expression (RLE) plot of one or more transformations of cell expression values.

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
plotRLE(example_sce, colour_by = "Mutation_Status", style = "minimal")
mark
plotRLE(example_sce, colour_by = "Mutation_Status", style = "full",
        outlier.alpha = 0.1, outlier.shape = 3, outlier.size = 0)
mark

plotRowData

Plot row-level (i.e., gene) metadata from a SingleCellExperiment object.

data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
example_sce <- calculateQCMetrics(example_sce,
                                  feature_controls = list(ERCC=1:40))
example_sce <- normalize(example_sce)
plotRowData(example_sce, y="n_cells_by_counts", x="log10_total_counts")
mark
plotRowData(example_sce, y="n_cells_by_counts",
            size_by ="log10_total_counts",
            colour_by = "is_feature_control")
mark

plotScater

Plot the relative proportion of the library size that is accounted for by the most highly expressed features for each cell in a SingleCellExperiment object.

## Set up an example SingleCellExperiment
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts),
  colData = sc_example_cell_info
)
plotScater(example_sce)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9Y6IsBEG-1571570422795)(http://baimoc.ziptop.top/blog/20190406/smYj6xT8Fyil.png)]

plotScater(example_sce, exprs_values = "counts", colour_by = "Cell_Cycle")
mark
plotScater(example_sce, block1 = "Treatment", colour_by = "Cell_Cycle")
mark
cpm(example_sce) <- calculateCPM(example_sce, use_size_factors = FALSE)
plotScater(example_sce, exprs_values = "cpm", block1 = "Treatment",
           block2 = "Mutation_Status", colour_by = "Cell_Cycle")
mark

Reduced dimension plots

PCA

## Set up an example SingleCellExperiment
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
assays = list(counts = sc_example_counts),
colData = sc_example_cell_info
)
example_sce <- normalize(example_sce)
## Examples plotting PC1 and PC2
plotPCA(example_sce)
mark
plotPCA(example_sce, colour_by = "Cell_Cycle")
mark
plotPCA(example_sce, colour_by = "Cell_Cycle", shape_by = "Treatment")
mark
plotPCA(example_sce, colour_by = "Cell_Cycle", shape_by = "Treatment",
size_by = "Mutation_Status")
mark
## Force legend to appear for shape:
example_subset <- example_sce[, example_sce$Treatment == "treat1"]
plotPCA(example_subset, colour_by = "Cell_Cycle", shape_by = "Treatment",
by_show_single = TRUE)
mark
## Examples plotting more than 2 PCs
plotPCA(example_sce, ncomponents = 4, colour_by = "Treatment",
shape_by = "Mutation_Status")
mark
## Same for TSNE:
plotTSNE(example_sce, run_args=list(perplexity = 10))
mark
## Same for DiffusionMaps:
plotDiffusionMap(example_sce)
mark
## Same for MDS plots:
plotMDS(example_sce)
mark
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,658评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,482评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,213评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,395评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,487评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,523评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,525评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,300评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,753评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,048评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,223评论 1 343
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,905评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,541评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,168评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,417评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,094评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,088评论 2 352

推荐阅读更多精彩内容