單細胞轉(zhuǎn)錄組(scRNA-seq)分析01 | Scater包的使用

一、介紹

  • 用于單細胞RNA-seq數(shù)據(jù)
  • 提供嚴格的質(zhì)量控制:將原始測序讀數(shù)處理為可用于下游分析的高質(zhì)量表達數(shù)據(jù)集
  • 提供了豐富的繪圖工具套件
  • R包地址:http://bioconductor.org/packages/scater

二钻心、工作流

mark

三苍柏、常用函數(shù)

plotColData

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

# 載入示例數(shù)據(jù)
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 對象中每個特征和細胞的質(zhì)控標準
example_sce <- calculateQCMetrics(example_sce)

# 計算 SingleCellExperiment 對象中read計數(shù)矩陣的歸一化表達值
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/

解釋變量是一種自變量扛邑。這兩個術語通常可互換使用溅潜。但是潦牛,兩者之間的細微差別凹嘲。當一個變量是獨立的宛畦,它不影響在所有的任何其他變量。當變量不是獨立的時候,它是一個解釋變量。 它在臨床研究中非常重要熬的。對于大多數(shù)情況,特別是在統(tǒng)計數(shù)據(jù)中崭闲,這兩個術語基本相同。

假設您有兩個變量來解釋體重增加:快餐和蘇打水。雖然你可能認為吃快餐和喝蘇打水是相互獨立的罗岖,但它們并不是真的炕淮。那是因為快餐店鼓勵你在用餐時買蘇打水。如果你停在某個地方買蘇打水晶疼,那里經(jīng)常會有很多快餐選擇酒贬,比如熱狗。雖然這些變量并非完全相互獨立翠霍,但它們確實會對體重增加產(chǎn)生影響锭吨。它們被稱為解釋變量,因為它們可能為體重增加提供一些解釋寒匙。

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)

[外鏈圖片轉(zhuǎn)存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(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
?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末零如,一起剝皮案震驚了整個濱河市躏将,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌考蕾,老刑警劉巖祸憋,帶你破解...
    沈念sama閱讀 206,126評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異肖卧,居然都是意外死亡蚯窥,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評論 2 382
  • 文/潘曉璐 我一進店門塞帐,熙熙樓的掌柜王于貴愁眉苦臉地迎上來拦赠,“玉大人,你說我怎么就攤上這事葵姥『墒螅” “怎么了?”我有些...
    開封第一講書人閱讀 152,445評論 0 341
  • 文/不壞的土叔 我叫張陵榔幸,是天一觀的道長允乐。 經(jīng)常有香客問我,道長牡辽,這世上最難降的妖魔是什么喳篇? 我笑而不...
    開封第一講書人閱讀 55,185評論 1 278
  • 正文 為了忘掉前任,我火速辦了婚禮态辛,結果婚禮上麸澜,老公的妹妹穿的比我還像新娘。我一直安慰自己奏黑,他們只是感情好炊邦,可當我...
    茶點故事閱讀 64,178評論 5 371
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著熟史,像睡著了一般馁害。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上蹂匹,一...
    開封第一講書人閱讀 48,970評論 1 284
  • 那天碘菜,我揣著相機與錄音,去河邊找鬼限寞。 笑死忍啸,一個胖子當著我的面吹牛,可吹牛的內(nèi)容都是我干的履植。 我是一名探鬼主播计雌,決...
    沈念sama閱讀 38,276評論 3 399
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼玫霎!你這毒婦竟也來了凿滤?” 一聲冷哼從身側(cè)響起妈橄,我...
    開封第一講書人閱讀 36,927評論 0 259
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎翁脆,沒想到半個月后眷蚓,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 43,400評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡鹃祖,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,883評論 2 323
  • 正文 我和宋清朗相戀三年溪椎,在試婚紗的時候發(fā)現(xiàn)自己被綠了普舆。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片恬口。...
    茶點故事閱讀 37,997評論 1 333
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖沼侣,靈堂內(nèi)的尸體忽然破棺而出祖能,到底是詐尸還是另有隱情,我是刑警寧澤蛾洛,帶...
    沈念sama閱讀 33,646評論 4 322
  • 正文 年R本政府宣布养铸,位于F島的核電站,受9級特大地震影響轧膘,放射性物質(zhì)發(fā)生泄漏钞螟。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 39,213評論 3 307
  • 文/蒙蒙 一谎碍、第九天 我趴在偏房一處隱蔽的房頂上張望鳞滨。 院中可真熱鬧,春花似錦蟆淀、人聲如沸拯啦。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽褒链。三九已至,卻和暖如春疑苔,著一層夾襖步出監(jiān)牢的瞬間甫匹,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評論 1 260
  • 我被黑心中介騙來泰國打工惦费, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留兵迅,地道東北人。 一個月前我還...
    沈念sama閱讀 45,423評論 2 352
  • 正文 我出身青樓趁餐,卻偏偏與公主長得像喷兼,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子后雷,可洞房花燭夜當晚...
    茶點故事閱讀 42,722評論 2 345

推薦閱讀更多精彩內(nèi)容