劉小澤寫于19.4.4+4.7
結(jié)合一篇文章+scater說明書，并做了一些個人修改
首先上文章和scater包說明書
來自https://f1000research.com/articles/5-2122/v2轩缤，題目是：A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor 作者來自劍橋癌癥研究所傍念、EMBL中心晚缩、Sanger研究所，發(fā)表于2016年
文章主要借助Bioconductor包進行質(zhì)量控制、數(shù)據(jù)挖掘和標準化等基本步驟拔疚，還有細胞周期判斷、差異基因鑒定既荚、降維聚類識別亞群稚失、marker基因檢測等進階步驟。分析了幾個公開可用的數(shù)據(jù)集恰聘，包括造血干細胞句各、腦源性細胞吸占、輔助性T細胞和小鼠胚胎干細胞
注意：這篇文章用到的包版本比較舊，因此我們需要更新它的函數(shù)
文章中用到的主要是scater包凿宾，那么同時打開scater說明書對照看：https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/vignette-qc.html#cell-level-qc-metrics
分析前須知
- scRNA-seq數(shù)據(jù)比bulk RNA數(shù)據(jù)噪音更大矾屯。單個細胞中RNA含量非常低，更別說還要反轉(zhuǎn)錄建庫測序初厚，因此增加了Dropout(細胞中有基因卻檢測不到表達量)的情況
- scRNA的優(yōu)勢就是研究細胞的異質(zhì)性件蚕。例如，識別新的細胞亞型产禾，表征分化過程排作，將細胞與細胞周期階段對應，或檢測跨人群驅(qū)動差異的HVGs（highly variable genes）
- 流程從計數(shù)矩陣開始（其實上游分析也很重要亚情，目前基本公司直接會把表達矩陣給用戶纽绍，但是當結(jié)果不滿意時，可以考慮上游因素）势似。此工作流包含了：質(zhì)控（去掉有問題的細胞）拌夏；標準化（細胞特異性偏差，有沒有spikein等）履因；表達矩陣探索（根據(jù)基因表達數(shù)據(jù)判斷細胞周期障簿、推斷亞群）。最后栅迄，利用HVG和Marker基因?qū)Ω信d趣的基因進行排序站故。
- 基本流程可能相同，但是不同項目中需要不同的搭配毅舆。文章使用了多個公共scRNA-seq數(shù)據(jù)集西篓，看看如何搭配這些流程
造血干細胞分析
數(shù)據(jù)來自文章：https://f1000research.com/articles/5-2122/v2#ref-48
關鍵詞：小鼠造血干細胞(haematopoietic stem cells ，HSCs ) 憋活、Smart-seq2 岂津、96孔板、ERCC悦即、每個基因表達量由比對到外顯子區(qū)域的reads數(shù)決定吮成、GSE61533
單細胞轉(zhuǎn)錄組和普通bulk轉(zhuǎn)錄組比對、定量流程相似辜梳，只是需要注意：
The only additional consideration is that the spike-in information must be included in the pipeline
一般來說粱甫，spike-in序列在比對之前，在基因組構(gòu)建索引時是作為附加的FASTA文件作瞄；定量之前茶宵，GTF文件中也是要加上spike-in轉(zhuǎn)錄本和內(nèi)源基因信息的
加載數(shù)據(jù)
為了更具有代表性，文章給出了讀取gzip的excel文件方法宗挥。讀取的表達矩陣中每一行代表一個內(nèi)源基因或spike-in轉(zhuǎn)錄本乌庶，每一列表示一個單獨的HSC細胞种蝶，然后將表達矩陣保存在 scater 的SingleCellExperiment對象中
注意：這篇文章中用到的還是scater1.0版本的
newSCESet函數(shù)，目前scater版本已經(jīng)到了1.10.1安拟，函數(shù)也升級為了SingleCellExperiment
library(R.utils)
gunzip("GSE61533_HTSEQ_count_results.xls.gz", remove=FALSE, overwrite=TRUE)
library(gdata)
# 文件24.5M蛤吓，讀取大概用了三分鐘
counts <- read.xls('GSE61533_HTSEQ_count_results.xls', sheet=1, header=TRUE, row.names=1)
# 因為大文件讀進來不容易，所以做個備份還是有幫助的
counts_bk <- counts
dim(counts)
counts[1:3,1:3]
library(scater)
sce <- SingleCellExperiment(assays = list(counts = counts),
rowData = data.frame(gene = rownames(counts)),
colData = data.frame(cell = colnames(counts)))
看看有沒有報錯糠赦？是不是提示：
Error in all_dims[, 1L] : incorrect number of dimensions
想了一會会傲，突然想起來是不是數(shù)據(jù)格式的問題，因為提示dimensions有錯誤拙泽，一般數(shù)據(jù)框和矩陣之間轉(zhuǎn)換會出現(xiàn)這個問題淌山。另外SingleCellExperiment是需要使用矩陣的(這個從它的示例數(shù)據(jù)中就可以看到)
那好，試試吧
counts <- as.matrix(counts)
sce <- SingleCellExperiment(assays = list(counts = counts),
rowData = data.frame(gene = rownames(counts)),
colData = data.frame(cell = colnames(counts)))
sce
# 這次成功了
看看行名中有沒有ERCC spike-in或者MT(線粒體)基因(注意：每個數(shù)據(jù)中的線粒體縮寫不一樣顾瞻，有的全部大寫MT泼疑，有的是小寫mt)
# 檢測ERCC和MT，并直接統(tǒng)計數(shù)量
is.spike <- grepl("^ERCC", rownames(sce));sum(is.spike) # 有92個
is.mito <- grepl("^mt-", rownames(sce));sum(is.mito) # 有37個
質(zhì)控
質(zhì)控可以對colData中的cells和rawData中的features 分別進行統(tǒng)計荷荤，默認根據(jù)sce對象assays中的counts值進行質(zhì)控退渗，當存在其它表達量值比如logCounts等，可以用calculateQCMetrics的exprs_values參數(shù)進行轉(zhuǎn)換
結(jié)果得到一個包含一系列統(tǒng)計值的蕴纳，比如： total number of counts 会油、the proportion of counts in mitochondrial genes、spike-in transcripts等古毛，所有統(tǒng)計值都儲存在colData中
sce <- calculateQCMetrics(sce,
feature_controls=list(ERCC=is.spike, Mt=is.mito))
# 除了指定feature_controls,還可以指定cell_controls(比如blank wells or bulk controls)
colnames(rowData(sce))
head(colnames(colData(sce)))
結(jié)果主要有兩方面：
-
細胞層次QC metrics翻翩，其中比較常用的統(tǒng)計值有：
total_counts： total number of counts for the cell (i.e., the library size)total_features_by_counts: number of features for the cell that have counts above the detection limit (default of zero)pct_counts_X:percentage of all counts that come from the feature control set named X如果我們不用counts矩陣進行計算，那么結(jié)果中所有帶counts的字段都會被替換成為我們使用的矩陣的名稱
-
基因?qū)哟蜵C metrics稻薇，主要包括：
mean_counts： the mean count of the gene/featurepct_dropout_by_counts: the percentage of cells with counts of zero for each genepct_counts_Y: percentage of all counts that come from the cell control set named Y(這里如果不加cell_control就沒有這個結(jié)果)
輸入sce可以看到嫂冻，spikeNames(0)還是空值，但上面我們已經(jīng)檢測到了ERCC存在塞椎，因此這里我們要明確指出ERCC代表的是一個spike-in桨仿，這很重要 ，因為下游分析時需要對spike-in進行一些處理忱屑，比如方差估計和標準化蹬敲。可以利用isSpike指定
# 對照組/控制條件可以根據(jù)基因(features)定莺戒，比如spike-in轉(zhuǎn)錄本, 線粒體基因；也可以根據(jù)細胞定急波，如空的測序孔从铲、損傷的細胞等
library(SingleCellExperiment)
isSpike(sce,"ERCC") <- is.spike
# 如果要加其他的spike-in，比如adam
# isSpike(sce, "Adam") <- grepl("^Adam[0-9]", rownames(sce))
# 如果要取消spike-in
# isSpike(temp,"ERCC") <- NULL
QC圖
目的都是為了幫助判斷數(shù)據(jù)質(zhì)量如何
檢查表達量最高的一些基因(默認前50)
圖中每一行是一個基因澄暮，每一個bar(線條)就是這個基因的單個細胞中的表達量名段，圓點是每個基因的平均表達量
plotHighestExprs(sce, exprs_values = "counts")
這個結(jié)果一般是會包含mitochondrial genes, actin, ribosomal protein, MALAT1 ,spike-in 阱扬，但是表達量前50基因中spike-in如果太多，就說明一開始加的spike-in太多伸辟；另外如果出現(xiàn)太多的"假基因"或者預測的基因麻惶，就說明可能比對過程出了問題。
Frequency of expression against mean
有表達基因的細胞所占的比例（縱坐標）與基因表達量平均值(橫坐標)構(gòu)成信夫，對于大多數(shù)基因窃蹋，X和Y坐標應該是正相關（也就是說，隨著細胞數(shù)量的增加静稻，基因表達量也是不斷增加的）警没。如果出現(xiàn)負相關的離群點(outlier)【少量可以接受，但總趨勢還是要正相關】振湾，具體原因需要深入研究杀迹，
比如：橫坐標低縱坐標高的情況：如果比對過程產(chǎn)生了許多"表達"的假基因，但實際上平均到所有細胞這部分就比較少(因為并不是所有細胞都會有假基因)押搪；(基因表達量的增長跟不上細胞數(shù)量的增長)
橫坐標高縱坐標低的情況： PCR擴增偏差或罕見細胞群的存在树酪，導致少數(shù)細胞中平均表達量很高(細胞數(shù)量的增長跟不上基因表達量的增長)
表達基因與feature controls的對比
利用colData(sce)中的縱坐標是feature control的比例，橫坐標是表達的基因的數(shù)量大州。比較好的數(shù)據(jù)就是有大量表達的基因同時feature controls很少续语；如果feature controls較高而表達基因數(shù)很少，說明細胞測序失敗
plotColData(sce, x = "total_features_by_counts",
y = "pct_counts_feature_control") +
theme(legend.position = "top") +
stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)
# 另外如果有分組信息的話摧茴，還可以按分組進行上色
根據(jù)基因信息（rowData）畫圖
> colnames(rowData(sce))
[1] "gene" "is_feature_control"
[3] "is_feature_control_ERCC" "is_feature_control_Mt"
[5] "mean_counts" "log10_mean_counts"
[7] "n_cells_by_counts" "pct_dropout_by_counts"
[9] "total_counts" "log10_total_counts"
> plotRowData(sce, x = "n_cells_by_counts", y = "mean_counts")
也可以用multiplot同時畫多個圖
p1 <- plotColData(sce, x = "total_counts",
y = "total_features_by_counts")
p2 <- plotColData(sce, x = "pct_counts_feature_control",
y = "total_features_by_counts")
p3 <- plotColData(sce, x = "pct_counts_feature_control",
y = "pct_counts_in_top_50_features")
multiplot(p1, p2, p3, cols = 3)
過濾之按細胞過濾
總目標：去掉不好的細胞(文庫太小的)绵载、去掉不好的基因(低表達的)、去掉線粒體占比太高的苛白、去掉spike-in占比太高的
解決前兩項問題：所有細胞中文庫大小和有表達基因的數(shù)目
標準就是：總表達量(文庫)小的不行娃豹，表達基因太少的也不行
先看看過濾前的原始數(shù)據(jù)
> par(mfrow=c(1,2))
> hist(sce$total_counts/1e6, xlab="Library sizes (millions)", main="",
+ breaks=20, col="grey80", ylab="Number of cells")
> hist(sce$total_features_by_counts, xlab="Number of expressed genes", main="",
+ breaks=20, col="grey80", ylab="Number of cells")
過濾的時候不能直接使用絕對值，因為可能出現(xiàn)這種情況：不管細胞質(zhì)量如何购裙，只要測序深度高懂版，就能得到更多的reads。手動設置閾值往往存在把控不到位的情況躏率，對于沒有分析經(jīng)驗的我們來講躯畴，使用isOutlier進行閾值設定是更好的選擇，也就是先將文庫的中位數(shù)進行l(wèi)og轉(zhuǎn)換(排除太大或太小值的影響)薇芝，然后減去3倍的MAD值(median absolute deviations)蓬抄。
isOutlierdefines the threshold at a certain number of median absolute deviations (MADs) away from the median
為何設定3倍的MAD？ ：通常把偏離中位數(shù)三倍以上的數(shù)據(jù)作為異常值夯到，和均值標準差方法比嚷缭，中位數(shù)和MAD的計算不受極端異常值的影響，結(jié)果更加穩(wěn)健。
libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
解決后兩項問題：線粒體占比和spike-in占比
par(mfrow=c(1,2))
hist(sce$pct_counts_Mt, xlab="Mitochondrial proportion (%)",
ylab="Number of cells", breaks=20, main="", col="grey80")
hist(sce$pct_counts_ERCC, xlab="ERCC proportion (%)",
ylab="Number of cells", breaks=20, main="", col="grey80")
同樣設置偏離中位數(shù)三倍以上的數(shù)據(jù)為離群值
mito.drop <- isOutlier(sce$pct_counts_Mt, nmads=3, type="higher")
spike.drop <- isOutlier(sce$pct_counts_ERCC, nmads=3, type="higher")
最后只篩選去掉以上四種離群值的數(shù)據(jù)：
sce <- sce[,!(libsize.drop | feature.drop | mito.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
ByMito=sum(mito.drop), BySpike=sum(spike.drop), Remaining=ncol(sce))
# ByLibSize ByFeature ByMito BySpike Remaining
1 2 2 6 3 86
## 可以看到去掉上面四項內(nèi)容后還有86個細胞文庫
過濾后再次檢測
因為我們之前做的過濾是假設所有細胞都是高質(zhì)量的阅爽，但是事實上可能會有技術原因?qū)е碌牡唾|(zhì)量細胞路幸，最好的檢測辦法是基于QC metrics進行PCA，檢測離群點
> sce <- runPCA(sce, use_coldata = TRUE,
+ detect_outliers = TRUE)
> plotReducedDim(sce, use_dimred="PCA_coldata")
> summary(sce$outlier) # 發(fā)現(xiàn)還有離群點付翁，需要去掉
Mode FALSE TRUE
logical 85 1
> sce <- sce[,-sce$outlier]
> sce <- runPCA(sce, use_coldata = TRUE,
+ detect_outliers = TRUE)
> summary(sce$outlier)
Mode FALSE
logical 85
細胞周期預測
預測的方法參考： Scialdone et al. (2015) 利用小鼠胚胎干細胞每兩個基因表達量之間進行比較简肴，結(jié)果作為一組，最后得到scran包中的預先訓練數(shù)據(jù)集百侧。然后基于與細胞周期相關的轉(zhuǎn)錄過程的保守性砰识，利用cyclone函數(shù)對其他類型的細胞進行細胞周期預測
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
# 除了mouse還有human的數(shù)據(jù)
library(org.Mm.eg.db)
anno <- select(org.Mm.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL")
ensembl <- anno$ENSEMBL[match(rownames(sce), anno$SYMBOL)]
assignments <- cyclone(sce, mm.pairs, gene.names=ensembl)
plot(assignments$score$G1, assignments$score$G2M, xlab="G1 score", ylab="G2/M score", pch=16)
# 于是結(jié)果又過濾掉3個細胞
sce <- sce[,assignments$phases=="G1"]
得到的結(jié)果中，如果G1 score大于G2/M且大于0.5移层，就是G1期仍翰；如果G2/M大于G1且大于0.5，就是G2/M期观话；如果結(jié)果都不大于0.5予借，就是S期。[補充：關于細胞分期]
過濾之按基因過濾
Bourgon et al., 2010 研究表明频蛔，基因表達量為0或者接近0不能進行有效地統(tǒng)計分析灵迫。這里低豐度基因被定義為平均表達量低于過濾閾值1的基因，去除這些基因可以降低離散性晦溪，減少計算量瀑粥，而不會造成大量信息丟失。
ave.counts <- rowMeans(counts(sce))
keep <- ave.counts >= 1
sum(keep)
單細胞特有的標準化
deconvolution方法計算標準化因子
Stegle et al., 2015 研究表明三圆，read counts 在細胞間因捕獲效率和測序深度而產(chǎn)生不同狞换。在下游分析細胞間的偏差時，需要消除細胞間的這種無關生物因素的偏差（bias）舟肉。一般是先假設大部分基因在細胞間不是差異表達的修噪，那么兩個細胞間的多數(shù)非表達基因產(chǎn)生的count size差異就主要是由于bias產(chǎn)生，這部分就是需要標準化的路媚。
每個文庫中需要被標準化的counts值就是"size factors" 黄琼。計算size factors有多種方法，比如在bulk轉(zhuǎn)錄組中：利用DESeq2的 estimateSizeFactorsFromMatrix整慎，或者edgeR的calcNormFactors脏款。但是單細胞數(shù)據(jù)有個特點，它包含了許多的低表達量或0值裤园，因此可以先將多個細胞的count值加起來撤师，提高count size然后再計算總size factor，再利用deconvolved的算法得到細胞水平的factor
sce <- computeSumFactors(sce, sizes=c(20, 40, 60, 80))
summary(sizeFactors(sce))
plot(sizeFactors(sce), sce$total_counts/1e6, log="xy",
ylab="Library size (millions)", xlab="Size factor")
圖中可以看到拧揽，size factors和 library sizes for all cells強相關丈氓。也就說明了細胞間的不同主要是由于捕獲效率或測序深度的影響。如果是真正的細胞間差異基因强法，它們的total count 和size factor之間應該是一條非線性相關的趨勢線万俗，線周圍是升高的散點。
對spike-in transcripts計算size factors
計算內(nèi)源基因size factor的方法不適用于spike-in transcripts饮怯，試想一下沒有文庫定量的實驗(每個文庫在混合和混養(yǎng)測序之前的cDNA含量都不等)闰歪，細胞中內(nèi)源基因的RNA越多，就越需要更大的size factor去進行標準化蓖墅。
但是在文庫制備時库倘，每個細胞中加入的是等量的spike-in，因此它的含量和不同文庫中RNA含量無關论矾。如果也要按照基因?qū)用嬗嬎愕膕ize factor進行標準化教翩，一般會導致標準化過度，表達定量不準確贪壳。
同樣的思路也適用于有文庫定量的情況饱亿。在cDNA總量不變的情況下，任何內(nèi)源性RNA含量的增加都會抑制spike -in轉(zhuǎn)錄本的數(shù)目闰靴。因此彪笼， spike-in計數(shù)的bias將與基于基因的size factor得到的bias相反。
sce <- computeSpikeFactors(sce, type="ERCC", general.use=FALSE)
# general.use=FALSE保證了得到的size factor只適用于spike-in transcripts 而非內(nèi)源基因
將計算的size factor用于基因表達標準化
原來的count值被用來計算標準化的log count值蚂且，下游分析基于這個配猫。計算過程是這樣的：每個細胞中的每個count值先進行log(count+1)，然后除以特定細胞的size factor
sce <- normalize(sce)
標準化后杏死，就可以探索實驗因子與表達量的關系
看看是否存在技術性因素對基因表達的異質(zhì)性有實質(zhì)性貢獻泵肄。這里看下spike-in的影響
plotExplanatoryVariables(sce,
variables = c("total_counts_ERCC",
"log10_total_counts_ERCC"))
歡迎關注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球淑翼，想讓它成為一個不拽術語腐巢、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
