由于網(wǎng)速的原因，昨天的測(cè)試數(shù)據(jù)沒(méi)有下載成功，在我掛機(jī)下載了一整個(gè)晚上之后顯示
image.png
我也只能作罷。。。
先看看第二步吧
Analyzing single-cell RNA-seq data containing read counts
1 Overview
- In this workflow, we use a relatively simple dataset (Lun et al. 2017) to introduce most of the concepts of scRNA-seq data analysis.
This dataset contains two plates of 416B cells (an immortalized mouse myeloid progenitor cell line), processed using the Smart-seq2 protocol (Picelli et al. 2014). A constant amount of spike-in RNA from the External RNA Controls Consortium (ERCC) was also added to each cell’s lysate prior to library preparation.
- Counts for all genes/transcripts in each cell are available from ArrayExpress using the accession number E-MTAB-5522. We download both the count tables (in the “processed files”) as well as the metadata file using the BiocFileCache package. This saves the files to a local cache (
raw_data) and avoids re-downloading them if they are already present.
BiocManager::install("BiocFileCache")
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
lun.zip <- bfcrpath(bfc,
file.path("https://www.ebi.ac.uk/arrayexpress/files",
"E-MTAB-5522/E-MTAB-5522.processed.1.zip"))
lun.sdrf <- bfcrpath(bfc,
file.path("https://www.ebi.ac.uk/arrayexpress/files",
"E-MTAB-5522/E-MTAB-5522.sdrf.txt"))
unzip(lun.zip, exdir=tempdir())
2 Setting up the data
2.1 Loading in the count matrix
- 第一步 將計(jì)數(shù)矩陣加載到內(nèi)存中。分別為研究中使用的每個(gè)細(xì)胞板生成了一個(gè)矩陣
plate1、plate2。在每個(gè)矩陣中，每行代表基因或轉(zhuǎn)錄本，每列表示一個(gè)細(xì)胞（跟RNAseq表達(dá)矩陣類(lèi)似，行名是基因名ensemble/entrizID，列名是樣本名GSM）。隨后，矩陣的每個(gè)條目中的計(jì)數(shù)表示映射到特定細(xì)胞中特定基因/轉(zhuǎn)錄的讀取次數(shù)。
plate1 <- read.delim(file.path(tempdir(), "counts_Calero_20160113.tsv"),
header=TRUE, row.names=1, check.names=FALSE)
plate2 <- read.delim(file.path(tempdir(), "counts_Calero_20160325.tsv"),
header=TRUE, row.names=1, check.names=FALSE)
gene.lengths <- plate1$Length # First column is the gene length.
plate1 <- as.matrix(plate1[,-1]) # Discarding gene length (as it is not a cell).
plate2 <- as.matrix(plate2[,-1])
rbind(Plate1=dim(plate1), Plate2=dim(plate2))
- 第二步 將兩個(gè)矩陣合并到一個(gè)對(duì)象中，以便進(jìn)一步處理。ps
(這是在驗(yàn)證兩個(gè)矩陣之間的基因順序相同后完成的。)
stopifnot(identical(rownames(plate1), rownames(plate2)))
all.counts <- cbind(plate1, plate2)
- 為方便起見(jiàn)，計(jì)數(shù)矩陣存儲(chǔ)在
SingleCellExperiment包中的SingleCellExperiment對(duì)象中。這允許將不同類(lèi)型的行級(jí)和列級(jí)元數(shù)據(jù)與計(jì)數(shù)一起存儲(chǔ)，以便在整個(gè)工作流中同步操作。(這一步也和DEGlist的獲得方法類(lèi)似)
library(SingleCellExperiment)
sce <- SingleCellExperiment(list(counts=all.counts))
rowData(sce)$GeneLength <- gene.lengths
sce
- 從行名稱(chēng)中標(biāo)識(shí)與 ERCC spike-in 轉(zhuǎn)錄腳本對(duì)應(yīng)的行。我們將此信息存儲(chǔ)在
SingleCellExperiment對(duì)象中，以便將來(lái)使用。這是必需的，因?yàn)閟pike-in在下游步驟如規(guī)一化過(guò)程中需要特殊處理。
isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))
summary(isSpike(sce, "ERCC"))
- 此數(shù)據(jù)集稍有不尋常，因?yàn)樗瑏?lái)自另一組spike-in轉(zhuǎn)錄本（Spike-In RNA Variants (SIRV) set)的信息。為簡(jiǎn)單起見(jiàn)，我們將僅在此分析中使用
ERCC spike-in。因此，在進(jìn)一步分析之前，我們必須刪除與 SIRV 轉(zhuǎn)錄本對(duì)應(yīng)的行，只需通過(guò)子設(shè)置SingleCellExperiment對(duì)象即可完成。
is.sirv <- grepl("^SIRV", rownames(sce))
sce <- sce[!is.sirv,]
summary(is.sirv)
Comments from Aaron:
- 某些要素計(jì)數(shù)工具將在計(jì)數(shù)矩陣中報(bào)告映射統(tǒng)計(jì)信息（例如，未對(duì)齊或未分配的讀取數(shù)）。雖然這些值對(duì)質(zhì)量控制很有用，但如果將其視為基因表達(dá)值，則它們會(huì)具有誤導(dǎo)性。因此，在進(jìn)一步分析之前，應(yīng)將其移除（或至少移動(dòng)到 colData）。
- 請(qǐng)注意，在計(jì)數(shù)矩陣的行名稱(chēng)為基因符號(hào)的人類(lèi)數(shù)據(jù)中使用
^ERCC正則表達(dá)式。ERCC基因家族實(shí)際上存在于人類(lèi)注釋中，因此這將導(dǎo)致基因作為spike-in轉(zhuǎn)錄本的錯(cuò)誤識(shí)別。通過(guò)發(fā)布具有標(biāo)準(zhǔn)標(biāo)識(shí)符的計(jì)數(shù)矩陣（例如，Ensembl、Entrez），可以避免此問(wèn)題。
2.2 Incorporating cell-based annotation
- 從 sdrf.txt 文件加載每個(gè)庫(kù)/單元格的元數(shù)據(jù)。請(qǐng)務(wù)必檢查元數(shù)據(jù)表的行與計(jì)數(shù)矩陣的列的順序相同。否則，不正確的元數(shù)據(jù)將分配給每個(gè)單元格。
metadata <- read.delim(lun.sdrf, check.names=FALSE, header=TRUE)
m <- match(colnames(sce), metadata[["Source Name"]]) # Enforcing identical order.
stopifnot(all(!is.na(m))) # Checking that nothing's missing.
metadata <- metadata[m,]
head(colnames(metadata))
- 只保留相關(guān)的元數(shù)據(jù)字段，以避免在
SingleCellExperiment對(duì)象的 colData 中存儲(chǔ)不必要的信息。特別是，此處保留每個(gè)細(xì)胞的原產(chǎn)板 (i.e., block)和表型。第二個(gè)字段是有關(guān)聯(lián)的，因?yàn)樗屑?xì)胞都表達(dá)oncogene，CBFB-MYH11基因，但是這種oncogene基因只在細(xì)胞的一個(gè)亞型中誘導(dǎo)表達(dá)。
colData(sce)$Plate <- factor(metadata[["Factor Value[block]"]])
pheno <- metadata[["Factor Value[phenotype]"]]
levels(pheno) <- c("induced", "control")
colData(sce)$Oncogene <- pheno
table(colData(sce)$Oncogene, colData(sce)$Plate)
2.3 Incorporating gene-based annotation
- 特征計(jì)數(shù)工具(Feature-counting)通常根據(jù) Ensembl 或 Entrez 的標(biāo)準(zhǔn)標(biāo)識(shí)符報(bào)告基因。這些標(biāo)識(shí)符的使用，因?yàn)樗鼈兪敲鞔_的，高度穩(wěn)定。然而，與文獻(xiàn)中較常用的基因符號(hào)( gene symbols )相比，它們很難理解。但是給定Ensembl標(biāo)識(shí)符，方便我們使用注釋包annotation packages（例如org.Mm.eg.db.）獲得相應(yīng)的基因符號(hào)。
library(org.Mm.eg.db)
symb <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="ENSEMBL", column="SYMBOL")
rowData(sce)$ENSEMBL <- rownames(sce)
rowData(sce)$SYMBOL <- symb
head(rowData(sce))
- 下一步就是過(guò)濾去重
- 通常需要將
sce的行名稱(chēng)重命名為基因符號(hào)，因?yàn)檫@些符號(hào)更易于解釋。但是，這需要一些工作來(lái)注釋缺省和重復(fù)的gene symbols。下面的代碼將用 Ensembl ID替換NA的gene symbols，并將重復(fù)的gene symbols替換為 Ensembl ID。
library(scater)
rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ENSEMBL, rowData(sce)$SYMBOL)
head(rownames(sce))
- 我們還使用
TxDb.Mmusculus.UCSC.mm10.ensGene包來(lái)確定每個(gè)基因的染色體位置。這在以后將非常有用，因?yàn)閹讉€(gè)質(zhì)量控制指標(biāo)將從與線粒體基因?qū)?yīng)的行中計(jì)算。
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.ensGene")
library(TxDb.Mmusculus.UCSC.mm10.ensGene)
location <- mapIds(TxDb.Mmusculus.UCSC.mm10.ensGene, keys=rowData(sce)$ENSEMBL,
column="CDSCHROM", keytype="GENEID")
rowData(sce)$CHR <- location
summary(location=="chrM")
或者，可以使用
scater中的getBMFeatureAnnos函數(shù)將來(lái)自BioMart資源的注釋直接添加到對(duì)象中。這可能比上面所示的方法更方便，但取決于 BioMart 數(shù)據(jù)庫(kù)的可用性與互聯(lián)網(wǎng)連接情況，這個(gè)網(wǎng)速就很玄妙了。
3 Quality control on the cells
3.1 Defining the quality control metrics
- 需要去除低質(zhì)量細(xì)胞，以確保技術(shù)差異不會(huì)干擾下游分析結(jié)果。我們使用多個(gè)質(zhì)量控制 （QC） 指標(biāo)：
- 庫(kù)的大小定義為所有要素（即基因和spike-in 轉(zhuǎn)錄本）的計(jì)數(shù)總和。小庫(kù)的細(xì)胞質(zhì)量低，一般是因?yàn)镽NA在制備過(guò)程中沒(méi)有被有效捕獲（即未能轉(zhuǎn)化為cDNA和擴(kuò)增）。
- 每個(gè)細(xì)胞中的表征數(shù)定義為該細(xì)胞具有非零計(jì)數(shù)的表征數(shù)。任何表達(dá)基因很少的細(xì)胞都可能質(zhì)量較差，因?yàn)椴煌霓D(zhuǎn)錄組群尚未被成功捕獲。
- 映射到spike-in轉(zhuǎn)錄本的讀取比例是相對(duì)于每個(gè)單元格的庫(kù)大小計(jì)算的。高比例表明細(xì)胞質(zhì)量低劣，其中內(nèi)源RNA在加工過(guò)程中丟失（例如，由于細(xì)胞裂解或RNA降解）。每個(gè)細(xì)胞的spike-in RNA量相同，因此spike-in 計(jì)數(shù)中的富集是內(nèi)源RNA損失的表現(xiàn)。
- 在沒(méi)有spike-in轉(zhuǎn)錄本的情況下，也可以使用與線粒體基因組中基因映射的讀取比例。高比例表明細(xì)胞質(zhì)量低(Islam et al. 2014; Ilicic et al. 2016)，可能是因?yàn)榇┛准?xì)胞的細(xì)胞質(zhì)RNA的喪失。其理由是線粒體比單個(gè)轉(zhuǎn)錄分子大，不太可能通過(guò)細(xì)胞膜中的撕裂而逃脫。
- 對(duì)于每個(gè)細(xì)胞，我們使用
scater包（McCarthy 等人 2017 年）中的QCMetrics函數(shù)計(jì)算這些質(zhì)量控制指標(biāo)。這些質(zhì)量控制指標(biāo)存儲(chǔ)在SingleCellExperiment的行和列元數(shù)據(jù)中，以供將來(lái)參考。 - 這些指標(biāo)的分布如圖1所示，按基因誘導(dǎo)狀態(tài)和原板分層。目的是去除具有小庫(kù)、低提取表征數(shù)量和高spike-in（或線粒體）比例的低質(zhì)量細(xì)胞。這種細(xì)胞會(huì)干擾下游分析，例如，形成不同的簇，使結(jié)果解釋復(fù)雜化。
sce$PlateOnco <- paste0(sce$Oncogene, ".", sce$Plate)
multiplot(
plotColData(sce, y="total_counts", x="PlateOnco"),
plotColData(sce, y="total_features_by_counts", x="PlateOnco"),
plotColData(sce, y="pct_counts_ERCC", x="PlateOnco"),
plotColData(sce, y="pct_counts_Mt", x="PlateOnco"),
cols=2)
Figure 1Distributions of various QC metrics for all cells in the 416B dataset. This includes the library sizes, number of expressed genes, and proportion of reads mapped to spike-in transcripts or mitochondrial genes.
- 研究質(zhì)量控制指標(biāo)如何相互配合也是很有價(jià)值的（圖2）。通常，它們將大致一致，即總計(jì)數(shù)低的細(xì)胞也將具有低表示特征的數(shù)量和高ERCC/線粒體比例。明顯的差異可能對(duì)應(yīng)于細(xì)胞批次之間的技術(shù)差異（見(jiàn)下文）或RNA含量的真正生物學(xué)差異。
par(mfrow=c(1,3))
plot(sce$total_features_by_counts, sce$total_counts/1e6, xlab="Number of expressed genes",
ylab="Library size (millions)")
plot(sce$total_features_by_counts, sce$pct_counts_ERCC, xlab="Number of expressed genes",
ylab="ERCC proportion (%)")
plot(sce$total_features_by_counts, sce$pct_counts_Mt, xlab="Number of expressed genes",
ylab="Mitochondrial proportion (%)")
Figure 2: Behaviour of each QC metric compared to the total number of expressed features. Each point represents a cell in the 416B dataset.
3.2 Identifying outliers for each metric
- 為這些指標(biāo)選擇閾值并不簡(jiǎn)單，因?yàn)槠浣^對(duì)值取決于實(shí)驗(yàn)流程。例如，無(wú)論細(xì)胞的質(zhì)量如何，到更深的測(cè)序深度都會(huì)導(dǎo)致更多的reads和更多expressed features。同樣，在實(shí)驗(yàn)中使用更多的spike-in RNA將導(dǎo)致更高的spike-in比例。為了獲得自適應(yīng)閾值，我們假設(shè)大多數(shù)數(shù)據(jù)集由高質(zhì)量細(xì)胞組成，并標(biāo)識(shí)出各種 QC 指標(biāo)異常的細(xì)胞。
- 異常值是根據(jù)所有細(xì)胞中每個(gè)指標(biāo)的中值的絕對(duì)偏差 （MAD） 定義的。我們刪除日志庫(kù)大小低于日志庫(kù)中位數(shù) 3 個(gè)以上的細(xì)胞。對(duì)數(shù)變換可提高較小值的分辨率，尤其是在原始值的 MAD 與中位數(shù)相當(dāng)或大于中位數(shù)時(shí)。我們還去除了對(duì)數(shù)轉(zhuǎn)換的表達(dá)基因數(shù)量低于中值3 MAD的細(xì)胞。
libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower",
log=TRUE, batch=sce$PlateOnco)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower",
log=TRUE, batch=sce$PlateOnco)
-
batch=參數(shù)確保在指定板/基因的每個(gè)級(jí)別內(nèi)異常值。這允許isoutlier()適應(yīng)不同板的 QC 指標(biāo)的系統(tǒng)差異（圖 1），這可能是由于處理技術(shù)差異（例如，測(cè)序深度的差異）而不是質(zhì)量的任何變化而產(chǎn)生的。同樣的推理也適用于腫瘤基因誘導(dǎo)狀態(tài)，其中誘導(dǎo)細(xì)胞可能由于生物原因自然表達(dá)的基因較少。如果無(wú)法解釋這些系統(tǒng)差異，就會(huì)使MAD估計(jì)值膨脹，并損害低質(zhì)量細(xì)胞的去除。 - 我們以類(lèi)似的方式識(shí)別基于比例的指標(biāo)的異常值。在這里，不需要轉(zhuǎn)換，因?yàn)槲覀冋谧R(shí)別大型異常值，在原始尺度上，區(qū)別應(yīng)該相當(dāng)清楚。我們不需要使用線粒體比例，因?yàn)槲覀円呀?jīng)具有此數(shù)據(jù)集的spike-in比例（用于類(lèi)似目的）。這避免了由于細(xì)胞類(lèi)型之間的線粒體含量的真正差異而產(chǎn)生的潛在問(wèn)題，這些細(xì)胞類(lèi)型可能會(huì)混淆異常值識(shí)別。
spike.drop <- isOutlier(sce$pct_counts_ERCC, nmads=3, type="higher",
batch=sce$PlateOnco)
-按列進(jìn)行設(shè)置將僅保留通過(guò)上述每個(gè)篩選器的高質(zhì)量細(xì)胞。我們檢查每個(gè)篩選器移除的細(xì)胞數(shù)以及保留的細(xì)胞總數(shù)。去除相當(dāng)比例的細(xì)胞（> 10%）可能表示數(shù)據(jù)質(zhì)量存在整體問(wèn)題。
keep <- !(libsize.drop | feature.drop | spike.drop)
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
BySpike=sum(spike.drop), Remaining=sum(keep))
結(jié)果是
ByLibSize ByFeature BySpike Remaining
1 5 4 6 183
- 然后，我們對(duì)SingleCellExperiment對(duì)象進(jìn)行子集化，以?xún)H保留假定的高質(zhì)量細(xì)胞。我們還將原始對(duì)象保存到文件中以供以后使用。
sce$PassQC <- keep
saveRDS(sce, file="416B_preQC.rds")
sce <- sce[,keep]
dim(sce)
Comments from Aaron:
- 有關(guān)基于異常值的低質(zhì)量單元檢測(cè)所依據(jù)的假設(shè)的更詳細(xì)討論，請(qǐng)參閱本節(jié)。
isOutlier()還將返回每個(gè)指標(biāo)的確切篩選器閾值（如果指定了batch=，則在每個(gè)批次中）。這些對(duì)于檢查自動(dòng)選擇的閾值是否合適可能很有用。
attr(libsize.drop, "thresholds")
attr(spike.drop, "thresholds")
4 Classification of cell cycle phase
- 我們使用Scialdone et al. (2015)描述的預(yù)測(cè)方法，根據(jù)基因表達(dá)數(shù)據(jù)將細(xì)胞分類(lèi)為細(xì)胞周期階段。使用訓(xùn)練數(shù)據(jù)集，計(jì)算了每對(duì)基因之間兩個(gè)基因表達(dá)差異的符號(hào)。選擇在細(xì)胞周期階段與符號(hào)變化成對(duì)作為標(biāo)記。然后，測(cè)試數(shù)據(jù)集中的細(xì)胞可以分類(lèi)到適當(dāng)?shù)碾A段，具體取決于每個(gè)標(biāo)記對(duì)的觀測(cè)符號(hào)是否與一個(gè)相位或另一個(gè)的階段一致。
- 此方法用
scran包中cyclone函數(shù)實(shí)現(xiàn)。該包包含一組預(yù)先訓(xùn)練的鼠的標(biāo)準(zhǔn)數(shù)據(jù)，我們可以在readRDS函數(shù)中加載這些數(shù)據(jù)。我們使用數(shù)據(jù)集中每個(gè)基因的 Ensembl 標(biāo)識(shí)符來(lái)匹配標(biāo)準(zhǔn)數(shù)據(jù)的基因?qū)械幕颉?/li>
set.seed(100)
BiocManager::install("scran")
library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
assignments <- cyclone(sce, mm.pairs, gene.names=rowData(sce)$ENSEMBL)
- HSC 數(shù)據(jù)集中每個(gè)細(xì)胞的
cyclone結(jié)果如圖 3 所示。為每個(gè)細(xì)胞分配每個(gè)階段的分?jǐn)?shù)，分?jǐn)?shù)越高，與細(xì)胞處于該階段的概率較高對(duì)應(yīng)。我們專(zhuān)注于 G1 和 G2/M 指標(biāo)，因?yàn)檫@些是最具信息性的分類(lèi)。
plot(assignments$score$G1, assignments$score$G2M,
xlab="G1 score", ylab="G2/M score", pch=16)
Figure 3: Cell cycle phase scores from applying the pair-based classifier on the 416B dataset. Each point represents a cell, plotted according to its scores for G1 and G2/M phases.
- 如果 G1 分?jǐn)?shù)高于 0.5 且大于 G2/M 分?jǐn)?shù)，則細(xì)胞被歸類(lèi)為 G1 階段;G2/M 階段如果 G2/M 分?jǐn)?shù)高于 0.5 且大于 G1 分?jǐn)?shù);在S階段，如果兩個(gè)分?jǐn)?shù)都高于0.5。在這里，絕大多數(shù)細(xì)胞被歸類(lèi)為G1階段。我們將這些歸類(lèi)保存到
SingleCellExperiment對(duì)象中以供以后使用。
sce$phases <- assignments$phases
table(sce$phases)
結(jié)果如下
> table(sce$phases)
G1 G2M S
98 62 23
-
scran包中提供預(yù)先訓(xùn)練的數(shù)據(jù)集（Pre-trained classifiers），包括人類(lèi)和小鼠數(shù)據(jù)。雖然這里使用的小鼠預(yù)先訓(xùn)練的數(shù)據(jù)集是根據(jù)胚胎干細(xì)胞的數(shù)據(jù)訓(xùn)練得來(lái)的，但對(duì)于其他細(xì)胞類(lèi)型來(lái)說(shuō)仍然準(zhǔn)確（Scialdone et al. 2015）。這可能是由于與細(xì)胞周期相關(guān)的轉(zhuǎn)錄程序（Bertoli, Skotheim, and Bruin 2013; Conboy et al. 2007）?；谂鋵?duì)的方法也是一個(gè)非參數(shù)過(guò)程，對(duì)數(shù)據(jù)集之間的大多數(shù)技術(shù)差異是可靠的。
Comments from Aaron:
- 為了消除由于細(xì)胞周期階段造成的混淆效應(yīng)，我們可以篩選細(xì)胞，只保留特定階段（通常為 G1）的細(xì)胞，以便進(jìn)行下游分析?；蛘撸绻豢珊雎缘募?xì)胞數(shù)位于其他階段，我們可以使用分配的階段作為 blocking factor。這可以防止細(xì)胞循環(huán)效應(yīng)，而不會(huì)丟棄信息，稍后將詳細(xì)討論。
- classifiers對(duì)于與訓(xùn)練集中使用的數(shù)據(jù)有很大不同的數(shù)據(jù)可能不準(zhǔn)確，例如，由于使用不同的protocol。在這種情況下，用戶(hù)可以使用
sandbag函數(shù)從自己的訓(xùn)練數(shù)據(jù)構(gòu)造自定義classifiers。對(duì)于沒(méi)有預(yù)先訓(xùn)練的classifiers的其他模型生物體來(lái)說(shuō)，這也是必要的。- 在應(yīng)用
cyclone之前，不要過(guò)濾掉低豐度基因。即使基因未在任何細(xì)胞中表達(dá)，如果基因是相位特異性的，它仍然可用于分類(lèi)。相對(duì)于其他基因，它的表達(dá)水平雖低但仍然會(huì)產(chǎn)生信息，而過(guò)濾掉它們會(huì)減少權(quán)重（reduce power，姑且先這么理解）。
5 Examining gene-level expression metrics
5.1 Inspecting the most highly expressed genes
- 富集分析高表達(dá)的基因（圖4）。這通常應(yīng)該由構(gòu)成表達(dá)的轉(zhuǎn)錄本，如核糖體或線粒體蛋白的轉(zhuǎn)錄本主導(dǎo)。如果其他類(lèi)別的功能與預(yù)期的生物學(xué)不一致，則它們的存在可能會(huì)引起關(guān)注。例如，包含許多spike-in轉(zhuǎn)錄本的數(shù)據(jù)集表明，在庫(kù)制備過(guò)程中添加了過(guò)多的spike-in RNA，而沒(méi)有核糖體蛋白和/或其假基因的存在則表明其排列不理想。
fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotHighestExprs(sce, n=50) + fontsize
圖 4：分配給 416B 數(shù)據(jù)集中前 50 個(gè)高表達(dá)features的總計(jì)數(shù)百分比。對(duì)于每個(gè)features，每個(gè)條形表示分配給單個(gè)細(xì)胞的基因的百分比，而圓表示所有細(xì)胞的平均值。條形由每個(gè)單元格中的features總數(shù)著色，而圓形則根據(jù)features是否標(biāo)記為control feature進(jìn)行著色。
5.2 Filtering out low-abundance genes
- 低豐度基因存在問(wèn)題，因?yàn)榱慊蚪咏阌?jì)數(shù)不包含太多用于可靠統(tǒng)計(jì)推理的信息。在涉及假設(shè)檢驗(yàn)的應(yīng)用中，這些基因通常不能提供足夠的證據(jù)來(lái)拒絕零假設(shè)，但它們?nèi)匀辉黾恿硕嘀貦z驗(yàn)校正的嚴(yán)重程度。計(jì)數(shù)的離散性也可能干擾統(tǒng)計(jì)過(guò)程，例如，通過(guò)影響連續(xù)近似的準(zhǔn)確性。因此，在應(yīng)用下游方法之前，許多RNA-seq分析策略中通常會(huì)去除低豐度基因。
- 篩選策略的"最優(yōu)"選擇取決于下游應(yīng)用程序。與移除動(dòng)力不足的測(cè)試所需的分立性相比，通常需要更激進(jìn)的濾波器來(lái)去除離散性（e.g., for normalization）。對(duì)于假設(shè)檢驗(yàn)，篩選器統(tǒng)計(jì)量還應(yīng)獨(dú)立于零假設(shè)下的檢驗(yàn)統(tǒng)計(jì)量（Bourgon, Gentleman, and Huber 2010）。因此，我們（或相關(guān)函數(shù)）將根據(jù)需要在每個(gè)步驟進(jìn)行篩選，而不是對(duì)整個(gè)分析應(yīng)用單個(gè)篩選器（single filter）。
- 幾個(gè)指標(biāo)可以用來(lái)定義低豐度基因。最明顯的是每個(gè)基因的平均計(jì)數(shù)，該數(shù)在數(shù)據(jù)集中的所有單元格中計(jì)算。我們使用
calcAverage()函數(shù)計(jì)算這一點(diǎn)，該函數(shù)還對(duì)細(xì)胞之間的庫(kù)大小差異執(zhí)行一些調(diào)整。我們通常觀察到中度表達(dá)基因的峰值，遵循低表達(dá)基因的plateau （從圖上看像是一個(gè)平原）（圖5）。
ave.counts <- calcAverage(sce, use_size_factors=FALSE)
hist(log10(ave.counts), breaks=100, main="", col="grey80",
xlab=expression(Log[10]~"average count"))
Figure 5: Histogram of log-average counts for all genes in the 416B dataset.
- 可以對(duì)此值應(yīng)用最小閾值來(lái)篩選出低表示的基因。下面的示例演示了我們?nèi)绾蝿h除平均計(jì)數(shù)小于 1 的基因。
demo.keep中的TRUE值數(shù)對(duì)應(yīng)于篩選后保留的行/基因數(shù)。
demo.keep <- ave.counts >= 1
filtered.sce <- sce[demo.keep,]
summary(demo.keep)
- 我們還檢查表達(dá)每個(gè)基因的細(xì)胞數(shù)量。這與大多數(shù)基因的平均計(jì)數(shù)密切相關(guān)，因?yàn)榛蛟谠S多細(xì)胞中的表達(dá)將導(dǎo)致更高的平均值（圖6）。在極少數(shù)細(xì)胞中表達(dá)的基因往往會(huì)不被關(guān)注，因?yàn)樗鼈兪怯蓴U(kuò)增偽影驅(qū)動(dòng)的（盡管它們也可能來(lái)自稀有種群）。然后，我們可以去除在少于n個(gè)細(xì)胞中表達(dá)的基因。
num.cells <- nexprs(sce, byrow=TRUE)
smoothScatter(log10(ave.counts), num.cells, ylab="Number of cells",
xlab=expression(Log[10]~"average count"))
Figure 6: The number of cells expressing each gene in the 416B dataset, plotted against the log-average count. Intensity of colour corresponds to the number of genes at any given location.
- 如上所述，我們將在每個(gè)步驟中應(yīng)用這些篩選器，而不是通過(guò)子設(shè)置
sce全局應(yīng)用它們。這可確保在每個(gè)應(yīng)用程序中使用最合適的篩選器。盡管如此，我們刪除任何細(xì)胞中不表達(dá)的基因，以減少下游步驟的計(jì)算工作。這些基因不提供信息，任何過(guò)濾策略都會(huì)被移除。
to.keep <- num.cells > 0
sce <- sce[to.keep,]
summary(to.keep)
6 Normalization of cell-specific biases
6.1 Using the deconvolution method to deal with zero counts
- 讀取計(jì)數(shù)取決于細(xì)胞之間的捕獲效率和測(cè)序深度的差異(Stegle, Teichmann, and Marioni 2015)。在下游定量分析之前，需要標(biāo)準(zhǔn)化來(lái)消除這些細(xì)胞特異性偏差。這通常是通過(guò)假設(shè)大多數(shù)基因在細(xì)胞之間沒(méi)有差異表達(dá)（DE）來(lái)實(shí)現(xiàn)的。假定兩個(gè)細(xì)胞之間非DE大多數(shù)基因之間的計(jì)數(shù)大小的任何系統(tǒng)差異都代表偏差，并通過(guò)縮放去除。更具體地說(shuō)，"大小因子"的計(jì)算表示每個(gè)庫(kù)中應(yīng)縮放計(jì)數(shù)的程度。
- 大小因子可以用幾種不同的方法計(jì)算，例如，使用
DESeq2包中的estimateSizeFactorsFromMatrix函數(shù) (Anders and Huber 2010; Love, Huber, and Anders 2014)，或用edgeR包的calcNormFactors函數(shù)。但是，單細(xì)胞數(shù)據(jù)對(duì)于這些基于批量數(shù)據(jù)的方法來(lái)說(shuō)可能存在問(wèn)題，因?yàn)榈陀?jì)數(shù)和零計(jì)數(shù)占主導(dǎo)地位。為了克服這一點(diǎn)，我們匯集了來(lái)自多個(gè)單細(xì)胞的計(jì)數(shù)，以增加計(jì)數(shù)的大小，以便精確估計(jì)大小因子(Lun, Bach, and Marioni 2016)。然后，基于池的大小因子被“deconvolved”到基于細(xì)胞的因子中，以規(guī)范化每個(gè)細(xì)胞的表達(dá)式配置文件。
sce <- computeSumFactors(sce)
summary(sizeFactors(sce))
- 大小因子與所有細(xì)胞的庫(kù)大小密切相關(guān)（圖 7）。這表明，細(xì)胞之間的大多數(shù)系統(tǒng)差異是由捕獲效率或測(cè)序深度的差異引起的。細(xì)胞之間的任何 DE 都會(huì)在總計(jì)數(shù)和大小因子之間產(chǎn)生非線性趨勢(shì)，并且/或增加圍繞趨勢(shì)的散射。我們?cè)诨蛘T導(dǎo)后觀察到一些證據(jù)，其中誘導(dǎo)后的大小因子有計(jì)劃地降低。這與誘導(dǎo)后基因向上調(diào)節(jié)引入的成分偏差(Robinson and Oshlack 2010)是一致的。
plot(sce$total_counts/1e6, sizeFactors(sce), log="xy",
xlab="Library size (millions)", ylab="Size factor",
col=c("red", "black")[sce$Oncogene], pch=16)
legend("bottomright", col=c("red", "black"), pch=16, cex=1.2,
legend=levels(sce$Oncogene))
到這一步的時(shí)候開(kāi)始出錯(cuò)了，提示legend 為0，我單獨(dú)運(yùn)行了等號(hào)后面的是有東西的，不知道是不是前面哪里出了問(wèn)題，我從新從頭運(yùn)行了一次然后出現(xiàn)前面3個(gè)圖合在了一起，估計(jì)是這里的設(shè)置問(wèn)題
test
但是仍然出錯(cuò)，
image.png
Error in legend("bottomright", col = c("red", "black"), pch = 16, cex = 1.2, :
'legend' is of length 0
回頭去查一下看看是哪里的問(wèn)題
Comments from Aaron:
- 雖然deconvolution approach對(duì)scRNA-seq數(shù)據(jù)中的高頻為零是可靠的，但如果計(jì)數(shù)過(guò)多為零，它最終將失敗。這表現(xiàn)為manifests as negative size factors that are obviously nonsensical。為了避免這種情況，
computeSumFactors函數(shù)在計(jì)算size factors 之前將自動(dòng)刪除低豐度基因。library size-adjusted的平均計(jì)數(shù)低于指定閾值（最小值）的基因?qū)⒈缓雎浴?duì)于讀取計(jì)數(shù)數(shù)據(jù)，默認(rèn)值 1 通常令人滿(mǎn)意。- 基于細(xì)胞的QC應(yīng)始終在normalization前執(zhí)行，以去除表達(dá)基因數(shù)量極低的細(xì)胞。否則，
computeSumFactors()可能會(huì)為低質(zhì)量細(xì)胞生成size factors。這是因?yàn)檫@些細(xì)胞中存在太多的零，從而降低了池以消除零的有效性。有關(guān)函數(shù)的更多詳細(xì)信息，請(qǐng)參閱?computeSumFactors。sizes參數(shù)可用于指定用于計(jì)算sizes因子的池sizes數(shù)。sizes越多，估計(jì)量就越精確，但需要花費(fèi)一些計(jì)算時(shí)間和內(nèi)存。通常，所有sizes應(yīng)大于 20 個(gè)細(xì)胞，以確保每個(gè)池中有足夠的非零表達(dá)式值。用于有效池化的細(xì)胞總數(shù)也應(yīng)至少為 100。- 對(duì)于高度異構(gòu)的數(shù)據(jù)集，建議對(duì)細(xì)胞執(zhí)行粗略聚類(lèi)以削弱非 DE 假設(shè)。這可以通過(guò)
quickCluster()函數(shù)完成，并且通過(guò)clusters參數(shù)傳遞給computeSumFactors()函數(shù)的結(jié)果。我們?cè)谙乱粋€(gè)工作流中使用更大的數(shù)據(jù)集演示此方法。
6.2 Computing separate size factors for spike-in transcripts
- 從內(nèi)源基因計(jì)數(shù)計(jì)算的Size factors通常不適合normalizing spike-in轉(zhuǎn)錄本的計(jì)數(shù)?？紤]一個(gè)without library quantification，即，在池化和multiplexed sequencing之前，每個(gè)庫(kù)中的 cDNA 量未相等。在這里，含有更多RNA的細(xì)胞對(duì)內(nèi)源性基因的計(jì)數(shù)較大，因此， larger size factors to scale down those counts。然而，在庫(kù)制備期間，每個(gè)細(xì)胞中都會(huì)添加相同數(shù)量的spike-in RNA。這意味著spike-in轉(zhuǎn)錄本的計(jì)數(shù)不受RNA含量的影響。嘗試用基于基因的size factors使spike-in計(jì)數(shù)normalization將導(dǎo)致表達(dá)的過(guò)度normalization和不正確的量化。在進(jìn)行庫(kù)量化的情況下，也適用類(lèi)似的推理。對(duì)于恒定的總量cDNA，內(nèi)源RNA含量的任何增加都會(huì)抑制spike-in轉(zhuǎn)錄本的覆蓋率。因此，spike-in計(jì)數(shù)中的偏差將與基于基因的size factors捕獲的偏差相反。
- 為了確保正確執(zhí)行normalization，我們?yōu)閟pike-in set 計(jì)算一組單獨(dú)的size factors。對(duì)于每個(gè)細(xì)胞，特定于峰值的size factors定義為spike-in set中所有轉(zhuǎn)錄本的總數(shù)。這假定沒(méi)有一個(gè)spike-in轉(zhuǎn)錄本本身差異地表示，這是合理的，因?yàn)槊總€(gè)細(xì)胞中應(yīng)添加相同量和成分的spike-in RNA(Lun et al. 2017)。（See below for a more detailed discussion on spike-in normalization。）通過(guò)設(shè)置
computeSpikeFactors中的general.use=FALSE，這些size factors存儲(chǔ)在SingleCellExperiment對(duì)象的單獨(dú)字段中。這確保了它們將只與spike-in轉(zhuǎn)錄本一起使用，而不是與內(nèi)源基因一起使用。
6.3 Applying the size factors to normalize gene expression
- count data用于計(jì)算normalized的log-expression值，用于下游分析。每個(gè)值定義為每個(gè)計(jì)數(shù)與相應(yīng)細(xì)胞的size factors的
，在計(jì)算前
+1以避免零計(jì)數(shù)時(shí)值為空。按size factors劃分可確保消除任何特定于細(xì)胞的偏差。如果sce中存在spike-in-specific size factors，則它們將自動(dòng)應(yīng)用于將spike-in轉(zhuǎn)錄本與內(nèi)源基因分開(kāi)normalize。
sce <- normalize(sce)
- log-transformation很有用，因?yàn)樗馕吨抵械娜魏尾町惗贾苯颖硎炯?xì)胞之間表達(dá)差異的
值的changes。這通常比絕對(duì)的覆蓋范圍差異更相關(guān)，需要根據(jù)總體豐度來(lái)解釋這些差異。對(duì)數(shù)變換還提供一些方差穩(wěn)定性度量 (Law et al. 2014)，因此具有較大方差的高豐度基因不會(huì)主導(dǎo)下游分析。計(jì)算值還存儲(chǔ)為
"logcounts"矩陣還可以成為用于其他分析的元素。
7 Modelling the technical noise in gene expression
- 基因間觀察到的表達(dá)值的差異可能由真正的生物學(xué)差異或無(wú)意義的技術(shù)性噪聲。為了區(qū)分這兩個(gè)可能性，我們需要對(duì)每個(gè)基因的表達(dá)值方差的技術(shù)成分進(jìn)行建模。我們使用一組spike-in轉(zhuǎn)錄本來(lái)執(zhí)行此操作，這些spike-in轉(zhuǎn)錄本以相同的數(shù)量添加到每個(gè)細(xì)胞樣品中。因此，spike-in轉(zhuǎn)錄本應(yīng)沒(méi)有生物變異性，即其計(jì)數(shù)的任何差異都應(yīng)是技術(shù)來(lái)源。
- 我們使用
trendVar()函數(shù)將均值相關(guān)趨勢(shì)與spike-in轉(zhuǎn)錄本的對(duì)數(shù)表達(dá)式值的方差擬合。我們?cè)O(shè)置block=，以確保板之間的技術(shù)差異不會(huì)膨脹方差。鑒于一個(gè)基因的平均豐度，該趨勢(shì)的擬合值將用作該基因技術(shù)成分的估計(jì)值。方差的生物分量最終通過(guò)decomposeVar函數(shù)的每個(gè)基因的總方差中減去技術(shù)分量來(lái)計(jì)算。
var.fit <- trendVar(sce, parametric=TRUE, block=sce$Plate,
loess.args=list(span=0.3))
var.out <- decomposeVar(sce, var.fit)
head(var.out)
- 我們通過(guò)可視化檢查趨勢(shì)，確認(rèn)它對(duì)應(yīng)于峰值方差（圖 8）。波浪狀形狀是log-expression值的平均方差趨勢(shì)的典型。當(dāng)平均值從零增加時(shí)，方差的線性增加可以被觀察到，因?yàn)楫?dāng)計(jì)數(shù)增加時(shí)，大方差是可能的。在非常高的豐度下，采樣噪聲的影響會(huì)由于大量數(shù)定律而減小，從而導(dǎo)致方差減小。
plot(var.out$mean, var.out$total, pch=16, cex=0.6, xlab="Mean log-expression",
ylab="Variance of log-expression")
curve(var.fit$trend(x), col="dodgerblue", lwd=2, add=TRUE)
cur.spike <- isSpike(sce)
points(var.out$mean[cur.spike], var.out$total[cur.spike], col="red", pch=16)
Figure 8: Variance of normalized log-expression values for each gene in the 416B dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red).
- 我們檢查具有最大生物成分的基因的表達(dá)值分布。這可確保方差估計(jì)不由一兩個(gè)異常值的細(xì)胞驅(qū)動(dòng)（圖 9）。(這個(gè)小提琴圖還是挺漂亮的像幾個(gè)穿旗袍的美女，哈哈)
chosen.genes <- order(var.out$bio, decreasing=TRUE)[1:10]
plotExpression(sce, features=rownames(var.out)[chosen.genes]) + fontsize
Figure 9: Violin plots of normalized log-expression values for the top 10 genes with the largest biological components in the 416B dataset. Each point represents the log-expression value in a single cell.
Comments from Aaron:
- 在實(shí)踐中，趨勢(shì)擬合由于記錄量少和豐度分布不均而變得復(fù)雜。有關(guān)如何優(yōu)化擬合的更多詳細(xì)信息，請(qǐng)參閱此處。
- 在沒(méi)有spike-in的情況下，用戶(hù)可以設(shè)置
use.spikes=FALSE，以適合內(nèi)源基因方差的趨勢(shì)（參見(jiàn)此處）?；蛘?，我們可以基于泊森技術(shù)噪聲的假設(shè)創(chuàng)建趨勢(shì)，如此處所述。- Negative biological components通常來(lái)自
decomposeVar。這些在直覺(jué)上是毫無(wú)意義的，因?yàn)榛虿豢赡茉诩夹g(shù)噪聲下具有總變異。盡管如此，這些值的發(fā)生是由于對(duì)總方差的估算不精確，尤其是對(duì)于數(shù)量少的細(xì)胞。decomposeVar還生成 p 值，可用于在錯(cuò)誤發(fā)現(xiàn)率 （FDR） 的特定閾值下定義 HVG。稍后我們將更詳細(xì)地討論這一點(diǎn)，因?yàn)樵跀?shù)據(jù)探索期間，特征選擇不需要 HVG 的正式檢測(cè)。
8 Removing the batch effect
- 如前所述，數(shù)據(jù)是在兩個(gè)板塊上收集的。板間處理的微小無(wú)法控制的差異可能導(dǎo)致批次效應(yīng)，即不同板上細(xì)胞之間的表達(dá)系統(tǒng)差異。這種差異并不有趣，可以通過(guò)從 limma 包中應(yīng)用
removeBatchEffect()函數(shù)(Ritchie et al. 2015)來(lái)刪除。這消除了原板的影響，同時(shí)考慮了對(duì)子基因誘導(dǎo)的（感興趣的）影響。
library(limma)
assay(sce, "corrected") <- removeBatchEffect(logcounts(sce),
design=model.matrix(~sce$Oncogene), batch=sce$Plate)
assayNames(sce)
- 對(duì)于非基于模型的下游過(guò)程（例如聚類(lèi)和大多數(shù)尺寸減少形式），需要手動(dòng)批處理校正。但是，如果分析方法可以接受設(shè)計(jì)矩陣，則最好在設(shè)計(jì)矩陣中阻止滋擾因子，而使用
removeBatchEffect()。這是因?yàn)楹笳邲](méi)有考慮到殘余自由度的喪失，也沒(méi)有考慮到估計(jì)blocking factor terms的不確定性。
Comments from Aaron:
removeBatchEffect()執(zhí)行線性回歸，并將與blocking factors對(duì)應(yīng)的系數(shù)設(shè)置為零。如果每個(gè)批次中的人口組成已知（設(shè)置design=）或跨批次相同，則此功能有效。此數(shù)據(jù)集的這種假設(shè)是合理的，涉及兩個(gè)板上的均質(zhì)細(xì)胞系總體。然而，在大多數(shù)scRNA-seq應(yīng)用中，不同批次的變化因素并不相同，而且事先不知道。這促使使用更復(fù)雜的批處理校正方法，如mnnCorrect()。
9 Denoising expression values using PCA
- 一旦對(duì)技術(shù)噪聲進(jìn)行建模，我們就可以使用主要元件分析 （PCA） 來(lái)消除隨機(jī)技術(shù)噪聲。假設(shè)每個(gè)單元格表示高維表達(dá)式空間中的點(diǎn)，其中點(diǎn)的分布表示總方差。PCA 標(biāo)識(shí)此空間中的坐標(biāo)軸，這些軸可捕獲盡可能多的此方差。每個(gè)軸都是一個(gè)主組件 （PC），其中任何早期的 PC 都會(huì)比后來(lái)的 PC 解釋更多的差異。
- 我們假設(shè)，涉及共同調(diào)控的基因組的生物過(guò)程將占數(shù)據(jù)差異最大的部分。如果是這種情況，此過(guò)程應(yīng)由一個(gè)或多個(gè)較早的 PC 表示。相反，隨機(jī)技術(shù)噪聲影響每個(gè)基因獨(dú)立，將由以后的PC表示。
denoisePCA()函數(shù)將刪除后來(lái)的PC，直到丟棄的總方差等于PCA中使用的所有基因的技術(shù)成分之和。 - 該函數(shù)返回一個(gè)
SingleCellExperiment對(duì)象，其中包含reducedDims插槽中每個(gè)單元格的 PC 分?jǐn)?shù)。其目的是消除技術(shù)噪音，并豐富保留 PC 中的生物信號(hào)。這提高了下游過(guò)程（如聚類(lèi)）中基礎(chǔ)生物學(xué)的分辨率。
Comments from Aaron:
denoisePCA()將只使用具有正生物成分的基因，即方差大于擬合趨勢(shì)。這可確保要丟棄的總技術(shù)差異不會(huì)大于數(shù)據(jù)中的總方差。- 對(duì)于
technical=參數(shù)，該函數(shù)還將直接接受趨勢(shì)函數(shù)（即var.fit$trend）或每個(gè)基因的技術(shù)組件的矢量。在這里，我們提供來(lái)自decomposeVar()的DataFrame，以允許函數(shù)在批處理校正后調(diào)整剩余自由度的損失。具體來(lái)說(shuō)，批處理矩陣中的方差略低，需要對(duì)技術(shù)組件進(jìn)行一些重新縮放以進(jìn)行補(bǔ)償。- 此處不對(duì)豐度執(zhí)行濾波，這可確保仍可檢測(cè)到與罕見(jiàn)亞種群對(duì)應(yīng)的PC。離散性問(wèn)題較少，因?yàn)榈拓S度基因的方差也較低，從而降低了它們對(duì) PCA 的貢獻(xiàn)。
- 還可以獲取原始表達(dá)式矩陣的低等級(jí)近似值，捕獲與保留的 PC 等效的方差。這對(duì)于在需要基因表達(dá)值的下游過(guò)程之前進(jìn)行去角非常有用。
10 Visualizing data in low-dimensional space
10.1 With PCA
- 我們通過(guò)構(gòu)建前三個(gè)主成分的 PCA 圖來(lái)可視化細(xì)胞之間的關(guān)系（圖 10）。具有相似表達(dá)式配置文件的細(xì)胞應(yīng)位于圖中，而不同的細(xì)胞應(yīng)相距很遠(yuǎn)。在這種情況下，我們觀察到基于基因誘導(dǎo)狀態(tài)的細(xì)胞的清晰分離，與轉(zhuǎn)錄組的預(yù)期效果一致。
plotReducedDim(sce, use_dimred="PCA", ncomponents=3,
colour_by="Oncogene") + fontsize
Figure 10: Pairwise PCA plots of the first three PCs in the 416B dataset, constructed from normalized log-expression values of genes with positive biological components. Each point represents a cell, coloured according to oncogene induction status.
-相比之下，我們觀察到細(xì)胞沒(méi)有按批次明確分離（圖11）。這表明使用removeBatchEffect()的批處理更正步驟已成功。
plotReducedDim(sce, use_dimred="PCA", ncomponents=3,
colour_by="Plate") + fontsize
Figure 11: Pairwise PCA plots of the first three PCs in the 416B dataset, constructed from normalized log-expression values of genes with positive biological components. Each point represents a cell, coloured according to the plate of origin.
- 請(qǐng)注意，
plotReducedDim()將使用已通過(guò)denoisePCA()存儲(chǔ)在sce中的PCA結(jié)果。這使我們能夠快速生成具有不同美學(xué)的新繪圖，而無(wú)需重復(fù)整個(gè) PCA 計(jì)算。同樣，如果plotPCA()在SingleCellExperiment中可用，將使用現(xiàn)有結(jié)果，否則將重新計(jì)算它們。用戶(hù)應(yīng)設(shè)置rerun=TRUE，以在現(xiàn)有結(jié)果存在的情況下強(qiáng)制重新計(jì)算 PCs。
Comments from Aaron:
- 對(duì)于每種可視化方法，可以將其他特定于單元格的信息合并到每個(gè)點(diǎn)的大小或形狀中。這是使用大多數(shù)繪圖函數(shù)中的
size_by=和shape_by=參數(shù)來(lái)完成的。- 可以顯示更多分量，但這些組件的信息量通常較少，因?yàn)樗鼈兘庾x的方差較少。它們通常也更難解讀，因?yàn)樗鼈儽欢x為與早期 PCs 正交（因此取決于在這些 PCs 中檢測(cè)到的內(nèi)容）。
10.2With t-SNE
- 另一種廣泛使用的降維方法是t-stochastic neighbour embedding (t-SNE)方法(Van der Maaten and Hinton 2008)。t-SNE往往比PCA更好地在更多樣化的種群分離細(xì)胞。這是因?yàn)榍罢呖梢灾苯硬东@高維空間中的非線性關(guān)系，而后者必須在線性軸上表示非線性關(guān)系。但是，這種改進(jìn)的代價(jià)是付出更多的計(jì)算工作量，并且要求用戶(hù)考慮random seed and perplexity等參數(shù)（請(qǐng)參閱comments）。
- 我們使用
plotTSNE()函數(shù)演示了圖12中t-SNE圖的生成。我們?cè)O(shè)置use_dimred="PCA"來(lái)對(duì)數(shù)據(jù)的低級(jí)近似執(zhí)行t-SNE，允許算法利用前面的降噪步驟。
set.seed(100)
out5 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=5),
colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 5")
set.seed(100)
out10 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=10),
colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 10")
set.seed(100)
out20 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=20),
colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 20")
multiplot(out5, out10, out20, cols=3)
又出現(xiàn)跟上次一樣點(diǎn)圖沒(méi)有電的情況的figure12
-
t-SNE 是一種隨機(jī)方法，因此用戶(hù)應(yīng)多次運(yùn)行該算法，以確保結(jié)果具有代表性。腳本應(yīng)設(shè)置種子（通過(guò)
set.seed()命令），以確保所選結(jié)果可重現(xiàn)。還建議測(cè)試"perplexity"參數(shù)的不同設(shè)置，因?yàn)檫@將影響低維空間中點(diǎn)的分布。 - 在這里，我們調(diào)用 20 的perplexity的
runTSNE()， 將 t-SNE 結(jié)果存儲(chǔ)在我們的SingleCellExperiment對(duì)象中。這樣可以避免在我們想要使用plotTSNE()創(chuàng)建新繪圖時(shí)重復(fù)計(jì)算，因?yàn)閷⒏挠么鎯?chǔ)的結(jié)果。同樣，用戶(hù)可以將rerun=TRUE設(shè)置為在存在存儲(chǔ)的結(jié)果時(shí)強(qiáng)制重新計(jì)算。 - 這里沒(méi)有考慮許多其他降維策略，但也可以使用，例如多維縮放、擴(kuò)散圖。它們各有優(yōu)缺點(diǎn)，例如，擴(kuò)散圖（參見(jiàn)
plotDiffusionMap）將細(xì)胞沿連續(xù)軌跡放置，并適合可視化分化等分級(jí)過(guò)程(Angerer et al. 2016)。
Comments from Aaron:
- 有關(guān)如何解釋 t-SNE 繪圖的良指南可在http://distill.pub/2016/misread-tsne/中找到。這演示了二維嵌入中的聚類(lèi)之間的距離如何沒(méi)有多大意義，并且群集的明顯"size"（即spread）也沒(méi)有什么意義。
11 Clustering cells into putative subpopulations
11.1Defining cell clusters from expression data
- 去噪聲的log-expression值用于將細(xì)胞集中到假定的亞群中。具體而言，我們使用 Ward’s criterion來(lái)最小化每個(gè)群集內(nèi)的總方差，對(duì)單元之間的歐氏距離執(zhí)行分層聚類(lèi)。這會(huì)產(chǎn)生一個(gè)樹(shù)狀圖，將具有相似表達(dá)模式的細(xì)胞組合在一起，在選定的基因上。
pcs <- reducedDim(sce, "PCA")
my.dist <- dist(pcs)
my.tree <- hclust(my.dist, method="ward.D2")
- 聚類(lèi)是通過(guò)將dynamic tree cut (Langfelder, Zhang, and Horvath 2008) 來(lái)顯式定義的。這利用樹(shù)突圖中的分支形狀來(lái)優(yōu)化群集定義，并且比
cutree更適合復(fù)雜樹(shù)形圖。通過(guò)在cutreeDynamic中手動(dòng)指定cutHeight，可以更好地控制經(jīng)驗(yàn)聚類(lèi)。我們還將minClusterSize設(shè)置為低于默認(rèn)值 20 的值，以避免遠(yuǎn)距離小群集的虛假聚合。
library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist),
minClusterSize=10, verbose=0))
- 我們根據(jù)已知因素檢查每個(gè)簇中的細(xì)胞分布。每個(gè)群集由兩個(gè)批處理的單元組成，表示聚類(lèi)不是由批處理效應(yīng)驅(qū)動(dòng)的。觀察到每個(gè)簇的組成與
Oncogene有關(guān)的差異，與基因誘導(dǎo)的生物效應(yīng)一致。
table(my.clusters, sce$Plate)
table(my.clusters, sce$Oncogene)
- 我們可視化圖 13 中 t-SNE 圖上所有細(xì)胞的聚類(lèi)分配。相鄰細(xì)胞通常分配給同一群集，表示群集過(guò)程已正確應(yīng)用。
sce$cluster <- factor(my.clusters)
plotTSNE(sce, colour_by="cluster") + fontsize
Figure 13: t-SNE plot of the denoised PCs of the 416B dataset. Each point represents a cell and is coloured according to the cluster identity to which it was assigned.
-
我們使用輪廓寬度檢查聚類(lèi)的分離度（圖 14）。具有較大的正輪廓寬度的單元格比不同群集中的單元格更接近同一群集中的其他單元格。相反，寬度為負(fù)的單元格比分配給它的群集中的其他單元格更接近其他群集。理想情況下，每個(gè)群集將包含許多具有較大正寬度的單元格，這表明它與其他群集很好地分離。
Figure 14: Barplot of silhouette widths for cells in each cluster. Each cluster is assigned a colour and cells with positive widths are coloured according to the colour of its assigned cluster. Any cell with a negative width is coloured according to the colour of the cluster that it is closest to. The average width for all cells in each cluster is shown, along with the average width for all cells in the dataset. - 輪廓寬度可用于確定參數(shù)值，以最大化群集之間的分隔。例如，我們可以在
cutreeDynamic中更改剪切高度或分割深度，以最大化所有單元格的平均輪廓寬度。這通常為進(jìn)一步檢查提供令人滿(mǎn)意的初始聚類(lèi)。但是，請(qǐng)記住，聚類(lèi)的粒度與顯微鏡上的放大倍率非常類(lèi)似。不同的數(shù)據(jù)視圖可以使用不同的粒度獲得，其中一些在分離測(cè)量上可能不夠理想。如果群集不能為特定生物問(wèn)題提供所需的粒度，則用戶(hù)不應(yīng)將聚類(lèi)與最大分離性集中在一起。 - 圖 14 中大多數(shù)單元的輪廓正寬度相對(duì)較小，表明簇之間的分離較弱。這可能是過(guò)度聚類(lèi)的癥狀，其中在子基因歸納狀態(tài)上明確定義的聚類(lèi)進(jìn)一步被分割成分離得不太良好的子集。盡管如此，我們將在
my.cluster中繼續(xù)當(dāng)前聚類(lèi)方案，因?yàn)樗鼮檫M(jìn)一步描述異質(zhì)性提供了合理的分區(qū)。
Comments from Aaron:
- 另一種聚類(lèi)策略是使用從相關(guān)關(guān)系派生的距離矩陣（例如，如
quickCluster）。這對(duì)噪聲和規(guī)范化錯(cuò)誤更加可靠，但對(duì)表達(dá)式配置文件中的細(xì)微變化也不太敏感。- 沃德的標(biāo)準(zhǔn)和完整的連桿產(chǎn)生球形，緊湊的集群。特別是，完全聯(lián)系有利于形成直徑相同的簇。在某些情況下，這可能是可取的，但當(dāng)亞種群的方差不同時(shí)，這不太合適。因此，我們通常使用沃德的標(biāo)準(zhǔn)進(jìn)行初始聚類(lèi)。當(dāng)然，嘗試其他方法（建議）很簡(jiǎn)單（建議），前提是執(zhí)行一些評(píng)估，例如，使用輪廓寬度。
11.2Detecting marker genes between clusters
- 一旦通過(guò)聚類(lèi)識(shí)別假定的亞群體，我們就可以使用
findMarkers函數(shù)識(shí)別每個(gè)聚類(lèi)的標(biāo)記基因。這將在每個(gè)基因的對(duì)數(shù)表達(dá)值和每對(duì)簇之間執(zhí)行Welch t-tests(Soneson and Robinson 2018)。目的是測(cè)試每個(gè)群集中的 DE 與其他群集相比，同時(shí)阻止無(wú)意義因素，如原板差異（詳情請(qǐng)參閱此處）。top DE基因可能是很好的候選標(biāo)記，因?yàn)樗鼈兛梢杂行У貐^(qū)分不同簇中的細(xì)胞。
markers <- findMarkers(sce, my.clusters, block=sce$Plate)
- 對(duì)于每個(gè)群集，相關(guān)比較的 DE 結(jié)果合并到單個(gè)輸出表中。這允許通過(guò)從每個(gè)集群之間的成對(duì)比較中獲取Top DE 基因來(lái)輕松定義一組標(biāo)記基因。例如，要從每個(gè)比較的前 10 個(gè)基因中構(gòu)造群集 1 的標(biāo)記集，將篩選
marker.set保留Top小于或等于 10 的行。還報(bào)告了每個(gè)基因的其他統(tǒng)計(jì)信息，包括調(diào)整后的 p 值（見(jiàn)下文）和相對(duì)于其他每個(gè)群集的對(duì)數(shù)折疊變化。
marker.set <- markers[["1"]]
head(marker.set, 10)
- 我們保存候選標(biāo)記基因列表以供進(jìn)一步檢查。
write.table(marker.set, file="416B_marker_1.tsv", sep="\t",
quote=FALSE, col.names=NA)
- 我們可視化top候選項(xiàng)的表達(dá)式配置文件，以驗(yàn)證 DE 簽名是否健壯（圖 15）。與部分或所有其他群集相比，大多數(shù)top 標(biāo)記在cluster1 的單元中具有強(qiáng)且一致的向上或向下調(diào)節(jié)。粗略地檢查熱圖表明，cluster1包含的腫瘤基因誘導(dǎo)細(xì)胞，DNA復(fù)制和細(xì)胞周期基因有很強(qiáng)的低調(diào)節(jié)。這與衰老作為抗腫瘤反應(yīng)的潛在誘導(dǎo)是一致的(Wajapeyee et al. 2010)。通過(guò)基因集濃縮分析，例如，使用來(lái)自
limma的kegga或goana，可以對(duì)這些標(biāo)記的功能進(jìn)行更全面的調(diào)查。
top.markers <- rownames(marker.set)[marker.set$Top <= 10]
plotHeatmap(sce, features=top.markers, columns=order(sce$cluster),
colour_columns_by=c("cluster", "Plate", "Oncogene"),
cluster_cols=FALSE, center=TRUE, symmetric=TRUE, zlim=c(-5, 5), show_colnames = F)
Figure 15: Heatmap of mean-centred and normalized log-expression values for the top set of markers for cluster 1 in the 416B dataset. Column colours represent the cluster to which each cell is assigned, the plate of origin or the oncogene induction status of each cell, as indicated by the legend.
- 圖 15 中的許多標(biāo)記在所選群集中不是唯一向上或向下調(diào)節(jié)的。對(duì)unique DE的測(cè)試往往過(guò)于嚴(yán)格，因?yàn)樗雎粤嗽趦蓚€(gè)或多個(gè)簇中表達(dá)的重要基因。例如，在僅CD4+、僅CD8+、雙陽(yáng)性和雙陰性T細(xì)胞的混合種群中，Cd4或Cd8都不會(huì)被檢測(cè)為亞群特異性標(biāo)記，因?yàn)槊總€(gè)基因在兩個(gè)亞種群中表達(dá)。通過(guò)我們的方法，這兩個(gè)基因?qū)⒈皇叭楹蜻x標(biāo)記，因?yàn)樗鼈儗⒃谥辽僖粚?duì)亞種群之間被認(rèn)為是DE。然后可以選擇標(biāo)記的組合來(lái)描述亞群，這比嘗試找到unique的DE基因更靈活。
- 我們強(qiáng)烈建議選擇一些標(biāo)記用于具有獨(dú)立復(fù)制細(xì)胞群的驗(yàn)證研究。其目的是確定表達(dá)上調(diào)節(jié)標(biāo)記但不表示向下調(diào)節(jié)標(biāo)記的細(xì)胞的相應(yīng)子集。理想情況下，在驗(yàn)證過(guò)程中也會(huì)使用不同的表達(dá)量化技術(shù)，例如熒光原位雜交或定量PCR。這證實(shí)了亞群確實(shí)存在，并且不僅僅
scRNA-seq或計(jì)算分析的結(jié)果。
Comments from Aaron:
- 通過(guò)設(shè)置
direction="up", findMarkers將只返回每個(gè)cluster中與其他集群相比被調(diào)控的基因。這在高度異質(zhì)的人群中非常方便，可以專(zhuān)注于能夠立即識(shí)別每個(gè)cluster的基因。雖然未表達(dá)基因可能也有信息，但它對(duì)positive identification用處不大。findMarker也可以定向查找所選聚類(lèi)和所有其他群集之間的 DE 基因。這應(yīng)該通過(guò)設(shè)置pval.type="all"來(lái)完成，該值將每個(gè)基因的 p 值定義為涉及所選群集的所有對(duì)對(duì)比較中的最大值。結(jié)合direction="up"，這可用于標(biāo)識(shí)每個(gè)群集的唯一標(biāo)記。但是，這對(duì)過(guò)度聚類(lèi)很敏感，因?yàn)槿绻麑⒕垲?lèi)拆分為兩個(gè)較小的子群集，則唯一標(biāo)記基因?qū)⒉辉俅嬖凇?/li>- 必須強(qiáng)調(diào)，此處計(jì)算的（adjusted）p 值不能正確解釋為重要度量。這是因?yàn)槿杭菑臄?shù)據(jù)中根據(jù)經(jīng)驗(yàn)識(shí)別的，請(qǐng)參閱此處。
12 Concluding remarks
- 基本分析完成后，將
SingleCellExperiment對(duì)象保存到本地使用saveRDS函數(shù)存檔通常很有用。然后，可以使用readRDS函數(shù)輕松地將該對(duì)象還原到新的 R 任務(wù)中。這允許進(jìn)行進(jìn)一步的工作，而無(wú)需重復(fù)上述所有處理步驟。
saveRDS(file="416B_data.rds", sce)
- 有多種方法可用于對(duì)已處理的表達(dá)式數(shù)據(jù)執(zhí)行更復(fù)雜的分析。例如，細(xì)胞可以在偽時(shí)間（例如，沿著分化途徑的進(jìn)度）與
monocle(Trapnell et al. 2014)或TSCAN(Ji and Ji 2016)排序;細(xì)胞狀態(tài)層次結(jié)構(gòu)可以用sincellpackage (Julia, Telenti, and Rausell 2015)來(lái)表示;和振蕩行為可以用Oscope(Leng et al. 2015)。HVGs 可用于基因集擴(kuò)充分析，以識(shí)別具有異構(gòu)活性的生物通路和過(guò)程，使用專(zhuān)為 topGO 等批量數(shù)據(jù)設(shè)計(jì)的包，或使用專(zhuān)用的 scde (Fan et al. 2016)。這些分析的完整描述超出了此工作流的范圍，因此建議感興趣的用戶(hù)查閱相關(guān)文檔。
All software packages used in this workflow are publicly available from the Comprehensive R Archive Network (https://cran.r-project.org) or the Bioconductor project (http://bioconductor.org). The specific version numbers of the packages used are shown below, along with the version of the R installation.
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] cluster_2.1.0 dynamicTreeCut_1.63-1
[3] limma_3.40.6 scran_1.12.1
[5] TxDb.Mmusculus.UCSC.mm10.ensGene_3.4.0 GenomicFeatures_1.36.4
[7] scater_1.12.2 ggplot2_3.2.1
[9] org.Mm.eg.db_3.6.0 AnnotationDbi_1.46.0
[11] SingleCellExperiment_1.6.0 SummarizedExperiment_1.14.0
[13] DelayedArray_0.10.0 BiocParallel_1.18.0
[15] matrixStats_0.54.0 Biobase_2.44.0
[17] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0
[19] IRanges_2.18.1 S4Vectors_0.22.0
[21] BiocGenerics_0.30.0 BiocFileCache_1.8.0
[23] dbplyr_1.4.2 Rsubread_1.34.7
loaded via a namespace (and not attached):
[1] bitops_1.0-6 bit64_0.9-7 RColorBrewer_1.1-2 progress_1.2.2
[5] httr_1.4.1 tools_3.6.1 backports_1.1.4 R6_2.4.0
[9] irlba_2.3.3 KernSmooth_2.23-15 vipor_0.4.5 DBI_1.0.0
[13] lazyeval_0.2.2 colorspace_1.4-1 withr_2.1.2 tidyselect_0.2.5
[17] gridExtra_2.3 prettyunits_1.0.2 bit_1.1-14 curl_4.0
[21] compiler_3.6.1 BiocNeighbors_1.2.0 rtracklayer_1.44.4 labeling_0.3
[25] scales_1.0.0 rappdirs_0.3.1 stringr_1.4.0 digest_0.6.20
[29] Rsamtools_2.0.0 XVector_0.24.0 pkgconfig_2.0.2 rlang_0.4.0
[33] rstudioapi_0.10 RSQLite_2.1.2 DelayedMatrixStats_1.6.0 dplyr_0.8.3
[37] RCurl_1.95-4.12 magrittr_1.5 BiocSingular_1.0.0 GenomeInfoDbData_1.2.1
[41] Matrix_1.2-17 Rcpp_1.0.2 ggbeeswarm_0.6.0 munsell_0.5.0
[45] viridis_0.5.1 stringi_1.4.3 edgeR_3.26.8 zlibbioc_1.30.0
[49] Rtsne_0.15 plyr_1.8.4 grid_3.6.1 blob_1.2.0
[53] dqrng_0.2.1 crayon_1.3.4 lattice_0.20-38 Biostrings_2.52.0
[57] cowplot_1.0.0 hms_0.5.1 locfit_1.5-9.1 zeallot_0.1.0
[61] pillar_1.4.2 igraph_1.2.4.1 reshape2_1.4.3 biomaRt_2.40.4
[65] XML_3.98-1.20 glue_1.3.1 packrat_0.5.0 BiocManager_1.30.4
[69] vctrs_0.2.0 gtable_0.3.0 purrr_0.3.2 assertthat_0.2.1
[73] rsvd_1.0.2 viridisLite_0.3.0 pheatmap_1.0.12 tibble_2.1.3
[77] GenomicAlignments_1.20.1 beeswarm_0.2.3 memoise_1.1.0 statmod_1.4.32
