背景

我們經(jīng)常在探索性分析中使用scRNA-seq數(shù)據(jù)來表征細胞間的異質(zhì)性。諸如聚類和降維之類的過程會根據(jù)細胞的基因表達譜對細胞進行比較，這涉及將每個基因的差異匯總為一對細胞之間的單個（不相似）或相似的度量。用于此計算的基因的選擇對度量的行為和下游分析方法的性能有重大影響。我們希望選擇包含與系統(tǒng)生物學相關且有用信息的基因，同時刪除包含隨機噪聲的基因。這旨在保留有趣的生物學過程的結(jié)構而不會使該結(jié)構模糊不清，并減小數(shù)據(jù)大小以提高后續(xù)步驟的計算效率。

特征選擇的最簡單方法是根據(jù)基因在整個細胞群之間的表達量來選擇變化最大的基因。假定與僅受技術噪聲或基線水平的“無用”生物學變異（例如，來自轉(zhuǎn)錄爆發(fā)）影響的其他基因相比，真正的生物學差異將表現(xiàn)為受影響基因的變異增加。有幾種方法可用于量化每個基因的變異并選擇適當?shù)囊唤M高度可變的基因（HVG）。我們將在下面使用10X PBMC數(shù)據(jù)集進行討論，以進行演示：

#--- 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ù)的變化程度

量化每個基因變異的最簡單方法是簡單地計算細胞群間所有細胞的每個基因的對數(shù)歸一化后的表達值（為簡單起見，稱為“對數(shù)值”）的方差。這具有一個優(yōu)點，即特征基因的選擇基于用于后續(xù)下游步驟的相同對數(shù)值。特別是，對數(shù)值的方差最大的基因?qū)毎g的歐幾里得距離貢獻最大。通過在此處使用對數(shù)值，我們可以確保在整個分析過程中對異質(zhì)性的定量定義是一致的。

每個基因的方差的計算很簡單，但是特征選擇需要對均方差關系進行建模。正如在數(shù)據(jù)均一化中簡要討論的那樣，對數(shù)轉(zhuǎn)換無法實現(xiàn)完美的方差穩(wěn)定化，這意味著基因的方差更多地是由其豐度驅(qū)動的，而不是其潛在的生物異質(zhì)性。為了說明這種影響，我們使用modelGeneVar（）函數(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)

PBMC數(shù)據(jù)集中的差異是平均值的函數(shù)。 每個點代表一個基因，而藍線代表適合所有基因變化度的趨勢。

在任何給定的豐度下，我們都假定大多數(shù)基因的表達譜受隨機技術噪聲的影響。 在此假設下，我們的趨勢代表了技術噪聲作為豐度函數(shù)的估計。 然后，我們將每個基因的總方差分解為技術成分，即該基因豐度時趨勢的擬合值； 生物成分，定義為總差異與技術成分之間的差。 此生物學成分代表每個基因的“有趣”變異，可用作HVG選擇的指標。

# 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