背景
我們經(jīng)常在探索性分析中使用scRNA-seq數(shù)據(jù)來表征細(xì)胞間的異質(zhì)性。諸如聚類和降維之類的過程會(huì)根據(jù)細(xì)胞的基因表達(dá)譜對(duì)細(xì)胞進(jìn)行比較病游，這涉及將每個(gè)基因的差異匯總為一對(duì)細(xì)胞之間的單個(gè)（不相似）或相似的度量衷戈。用于此計(jì)算的基因的選擇對(duì)度量的行為和下游分析方法的性能有重大影響狭吼。我們希望選擇包含與系統(tǒng)生物學(xué)相關(guān)且有用信息的基因，同時(shí)刪除包含隨機(jī)噪聲的基因殖妇。這旨在保留有趣的生物學(xué)過程的結(jié)構(gòu)而不會(huì)使該結(jié)構(gòu)模糊不清刁笙，并減小數(shù)據(jù)大小以提高后續(xù)步驟的計(jì)算效率。
特征選擇的最簡(jiǎn)單方法是根據(jù)基因在整個(gè)細(xì)胞群之間的表達(dá)量來選擇變化最大的基因拉一。假定與僅受技術(shù)噪聲或基線水平的“無用”生物學(xué)變異（例如采盒，來自轉(zhuǎn)錄爆發(fā)）影響的其他基因相比，真正的生物學(xué)差異將表現(xiàn)為受影響基因的變異增加蔚润。有幾種方法可用于量化每個(gè)基因的變異并選擇適當(dāng)?shù)囊唤M高度可變的基因（HVG）磅氨。我們將在下面使用10X PBMC數(shù)據(jù)集進(jìn)行討論，以進(jìn)行演示：
#--- loading ---#
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)
library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
#--- gene-annotation ---#
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID,
column="SEQNAME", keytype="GENEID")
#--- cell-detection ---#
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
#--- quality-control ---#
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]
#--- normalization ---#
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
sce.pbmc
## class: SingleCellExperiment
## dim: 33694 3985
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3985): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
## TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
## colData names(3): Sample Barcode sizeFactor
## reducedDimNames(0):
## altExpNames(0):
還有416B數(shù)據(jù)集：
#--- loading ---#
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")
sce.416b$block <- factor(sce.416b$block)
#--- gene-annotation ---#
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SYMBOL")
rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SEQNAME")
library(scater)
rownames(sce.416b) <- uniquifyFeatureNames(rowData(sce.416b)$ENSEMBL,
rowData(sce.416b)$SYMBOL)
#--- quality-control ---#
mito <- which(rowData(sce.416b)$SEQNAME=="MT")
stats <- perCellQCMetrics(sce.416b, subsets=list(Mt=mito))
qc <- quickPerCellQC(stats, percent_subsets=c("subsets_Mt_percent",
"altexps_ERCC_percent"), batch=sce.416b$block)
sce.416b <- sce.416b[,!qc$discard]
#--- normalization ---#
library(scran)
sce.416b <- computeSumFactors(sce.416b)
sce.416b <- logNormCounts(sce.416b)
sce.416b
## class: SingleCellExperiment
## dim: 46604 185
## metadata(0):
## assays(2): counts logcounts
## rownames(46604): 4933401J01Rik Gm26206 ... CAAA01147332.1
## CBFB-MYH11-mcherry
## rowData names(4): Length ENSEMBL SYMBOL SEQNAME
## colnames(185): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(10): Source Name cell line ... block sizeFactor
## reducedDimNames(0):
## altExpNames(2): ERCC SIRV
量化基因的變化程度
log轉(zhuǎn)化數(shù)據(jù)的變化程度
量化每個(gè)基因變異的最簡(jiǎn)單方法是簡(jiǎn)單地計(jì)算細(xì)胞群間所有細(xì)胞的每個(gè)基因的對(duì)數(shù)歸一化后的表達(dá)值（為簡(jiǎn)單起見嫡纠，稱為“對(duì)數(shù)值”）的方差烦租。這具有一個(gè)優(yōu)點(diǎn)，即特征基因的選擇基于用于后續(xù)下游步驟的相同對(duì)數(shù)值除盏。特別是叉橱，對(duì)數(shù)值的方差最大的基因?qū)?duì)細(xì)胞之間的歐幾里得距離貢獻(xiàn)最大。通過在此處使用對(duì)數(shù)值者蠕，我們可以確保在整個(gè)分析過程中對(duì)異質(zhì)性的定量定義是一致的窃祝。
每個(gè)基因的方差的計(jì)算很簡(jiǎn)單，但是特征選擇需要對(duì)均方差關(guān)系進(jìn)行建模踱侣。正如在數(shù)據(jù)均一化中簡(jiǎn)要討論的那樣粪小，對(duì)數(shù)轉(zhuǎn)換無法實(shí)現(xiàn)完美的方差穩(wěn)定化，這意味著基因的方差更多地是由其豐度驅(qū)動(dòng)的抡句，而不是其潛在的生物異質(zhì)性探膊。為了說明這種影響，我們使用modelGeneVar（）函數(shù)描繪所有基因的豐度差異的趨勢(shì)（圖1）待榔。
library(scran)
dec.pbmc <- modelGeneVar(sce.pbmc)
# Visualizing the fit:
fit.pbmc <- metadata(dec.pbmc)
plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression",
ylab="Variance of log-expression")
curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)
在任何給定的豐度下腌闯，我們都假定大多數(shù)基因的表達(dá)譜受隨機(jī)技術(shù)噪聲的影響。 在此假設(shè)下雕憔，我們的趨勢(shì)代表了技術(shù)噪聲作為豐度函數(shù)的估計(jì)姿骏。 然后，我們將每個(gè)基因的總方差分解為技術(shù)成分橘茉，即該基因豐度時(shí)趨勢(shì)的擬合值工腋； 生物成分姨丈，定義為總差異與技術(shù)成分之間的差。 此生物學(xué)成分代表每個(gè)基因的“有趣”變異擅腰，可用作HVG選擇的指標(biāo)蟋恬。
# Ordering by most interesting genes for inspection.
dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),]
## DataFrame with 33694 rows and 6 columns
## mean total tech bio p.value FDR
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## LYZ 1.95605 5.05854 0.835343 4.22320 1.10534e-270 2.17409e-266
## S100A9 1.93416 4.53551 0.835439 3.70007 2.71036e-208 7.61572e-205
## S100A8 1.69961 4.41084 0.824342 3.58650 4.31570e-201 9.43173e-198
## HLA-DRA 2.09785 3.75174 0.831239 2.92050 5.93940e-132 4.86758e-129
## CD74 2.90176 3.36879 0.793188 2.57560 4.83929e-113 2.50484e-110
## ... ... ... ... ... ... ...
## TMSB4X 6.08142 0.441718 0.679215 -0.237497 0.992447 1
## PTMA 3.82978 0.486454 0.731275 -0.244821 0.990002 1
## HLA-B 4.50032 0.486130 0.739577 -0.253447 0.991376 1
## EIF1 3.23488 0.482869 0.768946 -0.286078 0.995135 1
## B2M 5.95196 0.314948 0.654228 -0.339280
