😽

【シングルセル解析】シングルセル単位でgene set scoring (3) ~ AUCell ~

2023/09/06に公開

bioinformatics

scrnaseq

aucell

tech

特定のpathway geneや細胞種markerなどの遺伝子群 (gene set)を使って、gene set scoringを行うことでactivity単位の評価ができる。
これを1細胞単位で行う手法について紹介する。

AddModuleScore ※ 別記事
UCell ※ 別記事
scGSEA ※ 別記事
scGSVA ※ 別記事
スコアの有意差検定 ※ 別記事

AUCell

各細胞の遺伝子ランキング（発現量が多い順）の上位にgeneset遺伝子がどれだけ含まれるかをスコアリングする方法。
他の遺伝子ランキングを使う手法に比べて、発現ランキング上位しか計算に使用しない点が特徴。

シングルセルレベルで転写因子制御をスコアリングするSCINEではこのパッケージ/アルゴリズムを用いている。SCINEの出力を理解するにもAUCellのアルゴリズムを知っておくとよい。

AUCは、1細胞中の発現量上位遺伝子中にどれだけgeneset遺伝子が含まれるかの割合を表現している。
AUCはgenesetの登録遺伝子数や遺伝子の種類、対象細胞の種類などに依存するため、どの程度であればスコアが高いのか低いのかは決められない。そこで、AUCスコアを絶対値として扱うだけではなく、single cell dateset内での相対的な評価も行い、各genesetの陽性細胞/陰性細胞の2値化も試みる。相対評価であるためdateset内にはスコアが低い細胞と高い細胞が十分含まれるシチュエーションが望ましい。

理想的なケースでは、低いAUCを持つ細胞と高いAUCを持つ細胞で二峰性の分布が得られる。ただし、高いAUCを持つ細胞がレアな割合の場合は、正規分布のような見た目な分布で、高いAUCを持つ細胞は外れ値のような場所に位置する。逆に高いAUCを持つ細胞が大半を占める場合は、正規分布の大きな山がメインになる。
もし生物学的にはgenesetがnegativeになるはずであっても、スコアが低い範囲内で相対評価が行われ、geneset positiveと判定されることも有り得る。またhouse keeping遺伝子のようにどの細胞でも高く発現している遺伝子を含む特異性の低いgenesetでは全細胞でスコアが高く出てしまい、相対的評価ではgeneset negativeとカウントされる。

使い方

デフォルトでは、AUCの計算には発現量が多い遺伝子順のランキングでの上位 5% のみが使用される。この上位5%にgeneset遺伝子が含まれる割合が結果に直接影響する (含まれるだけでなく、より上位にgeneset遺伝子があるほどスコアは高くなる)。
つまり、ランキング上位にgeneset遺伝子を多く持つ細胞ほどAUCが高くなり、ランキング上位にgeneset遺伝子を含まない細胞はAUCが低くなる。

ランキング上位を使う利点は、① 計算量軽減/処理高速化 ② 発現量が0ばかりのランキング下位遺伝子のノイズ軽減などが挙げられる。

インストール

Bioconductorから取得可能。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("AUCell")

併せて並列実行やplotに関するパッケージのインストールもvignetteでは推奨されている。

# To support paralell execution:
BiocManager::install(c("doMC", "doRNG","doSNOW"))
# For the main example:
BiocManager::install(c("mixtools", "SummarizedExperiment"))
# For the examples in the follow-up section of the tutorial:
BiocManager::install(c("DT", "plotly", "NMF", "d3heatmap", "shiny", "rbokeh",
                       "dynamicTreeCut","R2HTML","Rtsne", "zoo"))

入力データ - 発現マトリクス

行に遺伝子、列に細胞バーコードの発現マトリクスを用意する。単なるMatrix型ではなく、dgCMatrix型に変換しておく。

matrixをdgCMatrixに変換

library(Matrix)
exprMatrix <- as(exprMatrix, "dgCMatrix")

本記事ではAUCellパッケージvignetteで紹介されているGSE60361のシングルセルデータを使用する。

デモデータの用意

GEOから生データをダウンロードし、データ成型した発現マトリクスをローカルに保存する。

# (Downloads the data)
if (!requireNamespace("GEOquery", quietly = TRUE)) BiocManager::install("GEOquery")
library(GEOquery)
# gse <- getGEO('GSE60361') # does not work, the matrix is in a suppl file
geoFile <- "GSE60361_C1-3005-Expression.txt.gz"
download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60361/suppl/GSE60361_C1-3005-Expression.txt.gz", destfile = geoFile)

library(data.table)
exprMatrix <- fread(geoFile, sep="\t")
geneNames <- unname(unlist(exprMatrix[,1, with=FALSE]))
exprMatrix <- as.matrix(exprMatrix[,-1, with=FALSE])
rownames(exprMatrix) <- geneNames
exprMatrix <- exprMatrix[unique(rownames(exprMatrix)),]

# Remove file(s) downloaded:
file.remove(geoFile)

# Convert to sparse
library(Matrix)
exprMatrix <- as(exprMatrix, "dgCMatrix")

# Save for future use
mouseBrainExprMatrix <- exprMatrix
save(mouseBrainExprMatrix, file="exprMatrix_MouseBrain.RData")

デモ用に細胞数を減らした発現マトリクスで進める。

## すでにローカルに保存していれば次回からはload()で読み込める。
# load("exprMatrix_MouseBrain.RData")

set.seed(333)
exprMatrix <- mouseBrainExprMatrix[sample(rownames(mouseBrainExprMatrix), 5000),]

Seuratオブジェクトからなら次のようにして発現マトリクスを取得して進めるとよい。細胞ごとの発現ランキングを使うので、countsスロットでもdataスロットでも変わらないはず。

exprMatrix <- GetAssayData(seuratObject, slot = "counts", assay = "RNA")
exprMatrix <- seuratObject@assays$RNA@counts

入力データ - geneset

genesetに登録されている遺伝子をベクトルとしてまとめて、さらに複数のgenesetベクトルをリストにまとめる。リストの名前属性にはgeneset名をつけておく。
遺伝子表記（IDやSYMBOL）は入力に使用する発現マトリクスのものを使用する。

genesetのリストの例

本記事ではデモ用に用意されたgenesetを使用する。

デモ用genesetの用意

AUCellの機能を使ってデータ処理するが、実際にはあまり使用しない機能なので覚えなくてよい。

AUCellパッケージに保存されているgmtファイルからgenesetを読み込む。

library(AUCell)
library(GSEABase)
gmtFile <- paste(file.path(system.file('examples', package='AUCell')), "geneSignatures.gmt", sep="/")
geneSets <- getGmt(gmtFile)

subsetGeneSets()機能で、発現マトリクスと共通する遺伝子数を数え上げる。

geneSets <- subsetGeneSets(geneSets, rownames(exprMatrix)) 
cbind(nGenes(geneSets))

重複遺伝子数をgeneset名につけておく。「geneset名 (遺伝子数g)」という表記になる。

geneSets <- setGeneSetNames(geneSets, newNames=paste(names(geneSets), " (", nGenes(geneSets) ,"g)", sep=""))

発現マトリクスから遺伝子をランダムにサンプリングしたgenesetを作成

# Random geneset
set.seed(321)
extraGeneSets <- c(
  GeneSet(sample(rownames(exprMatrix), 50), setName="Random (50g)"),
  GeneSet(sample(rownames(exprMatrix), 500), setName="Random (500g)"))

geneSets <- GeneSetCollection(c(geneSets,extraGeneSets))

ハウスキーピング遺伝子のgenesetを作成

# 遺伝子ごとに発現が0より大きい細胞数を集計
countsPerGene <- apply(exprMatrix, 1, function(x) sum(x>0))

# Housekeeping-like geneset
# 分位数95%より大きい遺伝子から100遺伝子をサンプリング
extraGeneSets <- GeneSet(sample(names(countsPerGene)[which(countsPerGene>quantile(countsPerGene, probs=.95))], 100), setName="HK-like (100g)")

geneSets <- GeneSetCollection(c(geneSets,extraGeneSets))

1. Rank matrixの作成 AUCell_buildRankings()

まずは発現マトリクスから発現ランクマトリクスに変換する。
AUCellでは同ランク値にはrandomな順位が与えられる。低発現帯ほど同値の遺伝子が増えるので、検出遺伝子数が少ない細胞ほどAUCスコアが不安定になる。(もちろん高発現帯であっても同値の遺伝子はあるので計算結果は若干変動する。)
ランダムな処理があるので、再現性が必要であればシード値を固定しておく。(vignetteでは固定していないので若干値が変わる)

set.seed(333)

cells_rankings <- AUCell_buildRankings(exprMatrix, plotStats=FALSE)

cells_rankingsの中身

cells_rankings@assays@data@listData[["ranking"]] にランクマトリクスが保存されている

AUCell_buildRankings()機能のplotStats=TRUEオプションを使うと、細胞ごとの検出細胞数のヒストグラムが表示される。デモデータでは1000遺伝子程度検出されている細胞が最も多く、それ以上の遺伝子数が検出できる細胞は減っている。なので計算対象の上位遺伝子数を1500や2000とするのは不適切なのが分かる。（検出遺伝子数が少ない細胞を考えると、250遺伝子しか検出できない細胞で、上位1000遺伝子を計算対象とすると、251~1000番目の遺伝子は発現値0である。）

plotStats=TRUE オプション使用時に表示されるヒストグラム

2. AUCスコアの計算 AUCell_calcAUC()

各genesetと各細胞のAUCスコアを計算する。引数にはgenesetのリストと発現ランクマトリクスを指定。

cells_AUC <- AUCell_calcAUC(geneSets = geneSets, 
                            rankings = cells_rankings)

上位どれくらいの遺伝子をAUC算出に使用するかはaucMaxRank=で変更可能。デフォルトはaucMaxRank=ceiling(0.05*nrow(cells_rankings))となっている。0.05の箇所を変更して上位何パーセントを指定するかを変える。

3. AUC閾値の探索 AUCell_exploreThresholds()

geneset陽性か陰性かを2値化するための閾値をgenesetごとに設定する。
AUCスコアの閾値となる候補は複数手法で用意し、その中で最も厳しい閾値 (最も高いAUCスコア)を採用する。

引数にはAUCell_calcAUC()の出力を使用する。その他引数は次の通り。

plotHist: デフォルトでTRUEでAUCスコアのヒストグラムとAUC閾値の線、AUC閾値を満たす細胞数などが表示される。
assignCells: デフォルトでFALSEだが、TRUEにするとAUC閾値を満たした陽性細胞を返す。(デモ用にplotする際に使用するのでここではTRUEにしている)
nCores: 並列実行数。デフォルトは1。
thrP: 正規分布で外れ値とする確率。デフォルトは0.01で1%の細胞が外れ値になるようなAUC閾値が計算される。
smallestPopPercent: 予測される最小の陽性細胞割合。AUC閾値候補の計算の際に使用される。デフォルトは0.25。もしAUCスコアの分布が正規分布していれば25%の細胞が陽性、75%の細胞が陰性となるようなAUC閾値が計算される。
nBreaks: 表示されるヒストグラムのビン数。デフォルトは100。
densAdjust: 内部で使われるdensity()関数のadjust=引数に渡される。デフォルトは2。

ランダムな処理があるので再現性が必要であればシード値を固定しておく。

# plotを3行3列に並べて表示する設定。
par(mfrow=c(3,3)) # 元に戻す時はdev.off()

set.seed(333)
cells_assignment <- AUCell_exploreThresholds(cellsAUC = cells_AUC,
                                             plotHist = TRUE, 
                                             assignCells = TRUE
                                             )

ヒストグラムの濃い青色はAUC閾値を満たした陽性細胞を示している。縦線に青/赤/ピンク/緑/灰色などがあるが、それぞれがAUC閾値候補を示しており、実線が採用されたAUC閾値である。

AUCスコアのヒストグラムとAUC閾値の情報

cells_assignmentの中身

AUCell_calcAUC()の返り値は各genesetでの結果が入ったリストである。

AUCell_calcAUC()の返り値

cells_assignment$geneset名$aucThr$thresholdsにはAUC閾値候補、cells_assignment$geneset名$aucThr$selectedには採用されたAUC閾値が記録されている。以下は3つの候補のうちL_k2が採用されたとのこと。

AUC閾値

AUC閾値候補

AUC閾値の候補は最大6種類計算されている。AUCell_calcAUC()のデフォルト引数でthrP = 0.01、smallestPopPercent = 0.25が設定されており、AUC閾値に影響する。特にsmallestPopPercent = 0.25は重要。

AUC閾値候補たち

outlierOfGlobal　濃い緑線
AUCスコアが正規分布している場合に、外れ値のような細胞だけを採用する。
デモデータだと、「Random (500g)」で正規分布の大きな山の右端のものが採用されている。
そのgenesetのAUCスコアの平均、AUCスコアのSDに従う正規分布から累積確率（1-(thrP/nCells)）を満たす確率変数X (AUCスコア)を算出し閾値とする。閾値以上の細胞が陽性細胞である。※ 1-(thrP/nCells)はほぼ1に近いので、thrPを上げなければ陽性細胞は非常に少なく見積もられる。
通常、genesetのAUCスコアに高い低いがあるはずだが、randomなgenesetのように細胞特異性が低い場合にこの閾値が採用されることがある。

tenPercentOfMax　濃い緑線
AUCスコアが0か低い値を持つ細胞が多い時は、左側に高いピークがあり、右裾が長い分布になる。その場合は、AUCの最大値 * 0.1をAUC閾値候補として記録する。
左に偏っていても、右側に高いAUC集団の分布があれば、この閾値が採用されることは無い。
デモデータだと「Random (50g)」でAUC候補として計算されている。

minimumDens　青線
AUCスコアが高い細胞集団、低い細胞集団が得られた場合、二峰性（もしくは三峰性）のピークをもつ分布が得られる。density()関数を使って確率密度関数を推定し、山と山の間の谷の部分にあたるAUCスコアを閾値候補とする。
デモデータだと、4つのgenesetでAUC閾値として採用されている。しかし、「Astrocyte_Lein (7g)」では見た目の分布からは最適だとは言えない。
Global_k1　灰色線
AUCスコア平均、AUCスコアSDに従う正規分布で、累積確率（1 - (thrP/nCells + smallestPopPercent)）を満たす確率変数 (AUCスコア)の値をAUC閾値候補とする。
smallestPopPercent=引数の値が結果に大きく影響する。デフォルトはsmallestPopPercent=0.25なので、「1 - smallestPopPercent = 0.75」を満たすAUC閾値は大きな値になりやすい。(thrP=のデフォルトは0.01であり、thrP/nCellsはとても小さい値なので大きく影響はしない。)
デモデータだと、「Random (50g)」でAUC閾値として採用されている。実際に正規分布していなくても正規分布を仮定した時の値が決まってしまう。

L_k2　赤線
AUCスコアが二峰性の分布を持つ場合、2つの正規分布の混合分布として考えることができる。
平均値が低い山の右にAUC閾値を設定すれば、2つの混合モデルを分ける閾値となるはずである。そこで、左側の山の平均値と標準偏差に従う正規分布から累積確率（1 - (thrP/nCells)）を満たす確率変数（AUCスコア）を計算して、AUC閾値候補とする。「1 - (thrP/nCells)」は1に近いので左側の山をほぼ満たす箇所に閾値が設定できる。
デモデータだと、3つのgenesetでAUC閾値として採用されている。上記のminimumDensは推定された確率密度関数を使って二峰性の分布を探しているが、こちらは混合分布モデルを使って計算している。似た結果になることもあるが、density()では見つけられなかった分布が見つかることがある。しかし、「HK-like (100g)」のように外れ値のような山を見つけて二峰性と判断して線引きしてしまうこともある。

R_k3　ピンク線
場合によってはAUCスコアのピークが3か所ある分布を取ることもあり得る。例えば、genesetに免疫細胞の汎用性マーカーとT細胞マーカーを含んでいると、T細胞集団はAUCスコア強、T細胞以外の免疫細胞はAUCスコア弱、免疫細胞以外の細胞種ではAUCスコアゼロといった分布が起こり得る。
3つの正規分布が混合したモデルでは、最も平均値が高い右側の山の左にAUC閾値を設定する。右側の山の平均値と標準偏差に従う正規分布から累積確率（thrP）を満たす確率変数（AUCスコア）を計算して、AUC閾値候補とする。thrPはデフォルトでは0.01なので、右側の山を殆ど満たさない位置に閾値を設定できる。
デモデータだとどのgenesetにも採用されていないが、「Astrocyte_Lein (7g)」ではR_k3が尤もらしい。

4. AUC閾値の修正

採用されたAUC閾値が相応しい場合には、自ら閾値を再設定する。
AUCell_calcAUC()の処理の中で警告すべきことが結果のリストに保存されているのでどのgenesetを見直すべきかヒントにするとよい。

# 警告メッセージを取り出す
warningMsg <- as.matrix(sapply(cells_assignment, function(x) x$aucThr$comment, simplify = F))

# 警告メッセージがあるものだけを表示
warningMsg[ warningMsg[,1] != "", ]

警告メッセージ

※ ただし、「Astrocyte_Lein (7g)」のように不適切だと思っても、警告メッセージは何も無いということがあるので、なるべく目で閾値を確認した方が良い。

AUCell_plotHist()機能で自身が設定したAUC閾値をヒストグラムで確認することもできる。

thr <- cells_assignment$`Astrocyte_Lein (7g)`$aucThr$thresholds["R_k3",1]
AUCell_plotHist(cells_AUC["Astrocyte_Lein (7g)",], aucThr=thr)

必要であれば、cells_assignmentのAUC陽性細胞情報を更新してもよい。

genesetName <- "Astrocyte_Lein (7g)"
newSelectedCells <- names(which(getAUC(cells_AUC)[genesetName, ] > thr))
cells_assignment[[genesetName]]$assignment <- newSelectedCells

5. AUCスコア/陽性細胞の可視化

AUCellのvignetteに可視化方法が紹介されているのでそちらを確認して欲しい。
この記事ではSeuratオブジェクトにAUCellの結果を反映させる例を紹介する。(Seuratではmeta.dataの列名に半角スペースや括弧を使えないのでピリオドに変換されてしまうが)

デモデータからSeuratオブジェクトを作成

library(Seurat)
library(dplyr)
sc <- CreateSeuratObject(counts = exprMatrix)
sc <- sc %>% 
  NormalizeData() %>% 
  ScaleData() %>% 
  FindVariableFeatures() %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:30)

getAUC()機能ではAUCell_calcAUC()の出力をmatrixに変換することができる。SeuratのAddMetaData()に渡すために行に細胞バーコード、列にgenesetとなるように転置する。

# AUCのマトリクスを転置して取り出し
aucMatrix <- data.frame(t(getAUC(cells_AUC)))

# Seuratオブジェクトに登録
sc <- AddMetaData(object = sc, metadata = aucMatrix)

可視化

FeaturePlot(sc, features = colnames(aucMatrix), ncol = 3)

AUC閾値で二値化したデータの場合。

0,1のマトリクスを作成

cellsAssigned <- lapply(cells_assignment, function(x) x$assignment)
assignmentTable <- reshape2::melt(cellsAssigned, value.name="cell")
colnames(assignmentTable)[2] <- "geneSet"
assignmentMat <- table(assignmentTable[,"geneSet"], assignmentTable[,"cell"])
assignmentMat <- data.frame(t(as.data.frame.matrix(assignmentMat)))

可視化

sc <- AddMetaData(object = sc, metadata = assignmentMat)

FeaturePlot(sc, features = colnames(assignmentMat), ncol = 3, cols = c("lightgray","black"))

【アルゴリズム】

以下では計算の再現を行う。興味がある方は見てみてほしい。

genesetの用意

遺伝子ベクトルをまとめたリストを用意。リストの名前属性に各genesetの名前を付けておく。遺伝子表記は発現マトリクスのものと揃えておく。

Rank matrix

行に遺伝子、列に細胞の発現マトリクスを使用し、apply文で細胞ごとに発現ランク値に変換する。(1細胞内の発現値ランキング)

発現量が多い順にrank値をつけるため、正負を逆転してからrank()関数に渡している。ties.method="randomで"同じ発現量の遺伝子にはrandomにrank値がつくため必要であればシード値を固定しておく。

rankMat <- apply( -exprMatrix, 2, rank,ties.method="random")

AUCスコアの算出

パラメーター設定。上位どれくらいの遺伝子を使用するか。デフォルトが5%。

aucMaxRank <- ceiling(0.05*nrow(rankMat))

genesetごと、細胞ごとに繰り返し処理を行う。

auc_list <- 
  lapply(X = geneSets, FUN = function(x){
    # geneSetのリストからgenesetを1つ取り出す
    geneSet <- x
    
    # 細胞ごとに処理
    apply(rankMat, MARGIN = 2, FUN = function(y){
      # 1細胞の発現ランクベクトルを並び替え
      y <- sort(unlist(y))
      
      # genesetに含まれているかいないかで0,1に変換
      y[!names(y) %in% geneSet] <- 0 # genesetに含まれていないと0
      y[y > 0] <- 1 # それ以外は1
      
      ## AUC計算の分子
      # 0,1ベクトルの累積和ベクトルの総和
      numerator <- sum(cumsum(y[1:(aucMaxRank-1)])) 
            
      ## geneset AUC計算の分母の部分
      # genesetの登録遺伝子数がaucMaxRankより多ければ、aucMaxRankの総和。
      if(length(geneSet) >= aucMaxRank){
        denominator <- sum(seq(aucMaxRank-1))
      } else {
      # aucMaxRankよりgenesetの登録遺伝子数が少なければ,,,
        denominator <- sum(
          c(seq(length(geneSet)),
            rep(length(geneSet), aucMaxRank - length(geneSet) - 1))
          )
      }
      
      # AUCスコア
      auc <- numerator / denominator
      return(auc) 
    })
  })

# listをdata.frameにまとめる。
aucMatrix <- do.call(rbind, auc_list)

発現量上位から遺伝子を見ていき、genesetに登録された遺伝子であれば1を加算するという累積和を行う。ある細胞での例では次のような折れ線グラフが作れる。この線の下の面積を計算の分子に使用する。

上位250遺伝子が計算対象の場合の、ある細胞のAUC

分母は計算対象の全遺伝子が全てgeneset遺伝子だったという超理想的な場合のAUCが分母となる面積である。

上位250遺伝子が計算対象の場合の、分母となるAUC

code

y <- rankMat[,1]
y <- sort(unlist(y))
y[!names(y) %in% geneSet]
y[y > 0] 
plot(cumsum(y[1:(aucMaxRank-1)]), ylim = c(0,250), col = "blue", type = "l")

plot((1:(aucMaxRank-1)), ylim = c(0,250), col = "red", type = "l")

ちなみに、上位250個の遺伝子を対象にする場合、geneset遺伝子が250個以上あればよいが、250個より少ない場合は250個のなるべく上位にgeneset遺伝子が位置する場合のAUCが分母となる。

aucMaxRankよりgenesetの登録遺伝子数が少ない場合

genesetのフィルタリング

全ての細胞でAUCスコアがゼロのようなgenesetは除外する。

# geneSetごとの合計AUCを計算
rowSumAUC <- rowSums(aucMatrix)
# 合計AUCが0より大きいgenesetだけを残す
aucMatrix <- aucMatrix[rowSumAUC > 0, , drop = FALSE]

本来であれば、以降はgenesetごとに繰り返し処理を行うが、解説のため1つのgeneset * 細胞数分のAUCスコアのベクトルを使用する。

i番目のgenesetのAUCスコアのベクトルを使用

gSetNames <- rownames(aucMatrix)
gSetName <- gSetNames[i]

# 1つのgenesetのベクトル取得 (要素にAUCスコア、名前属性に細胞バーコード)
auc <- aucMatrix[gSetName ,]
# AUCスコアが小さい方から昇順に並び替え
auc <- sort(auc)

AUC閾値の探索

geneset陽性か陰性かを2値化するための閾値をgenesetごとに設定する。
AUCスコアの分布に応じて、AUC閾値の計算、取捨選択を行っていく。

パラメーター設定

以下はAUCell_exploreThresholds()機能のデフォルト値である。

smallestPopPercent <- 0.25 # <--想定される最小の細胞集団比率。デフォルトが0.25
thrP <- 0.01 # <-- 分布の中の外れ値確率
densAdjust <- 2 # <-- 確率密度関数を書くときの滑らかさ。いじらなくて良さそう。

中でもsmallestPopPercentは注意。あるgenesetの陽性細胞割合の予測最小値であるが、本来genesetごとに変わるはずだ。綺麗にAUCスコアの二峰性の分布が得られるgenesetの場合は、この引数は重要ではないが、それ以外のAUCスコア分布の場合はこの予測割合が重要となる。

AUCスコアの平均とSDに従う正規分布から 1 - smallestPopPercent (デフォルトだと 0.75) の累積確率を満たすAUC閾値を設定する。もしAUCスコアが正規分布に従っているなら、75%の細胞が陰性、25%の細胞が陽性になるように閾値が設定される。

AUCell_exploreThresholds()結果のオブジェクトでcells_assignment[["geneset名"]][["aucThr"]][["selected"]]が"Global_k1"となっている採用している場合や、表示されるヒストグラム内の閾値線が灰色の場合は、このsmallestPopPercentを使って計算されたAUC閾値が採用されている。
genesetごとのヒストグラムをしっかりと確認して、デフォルトのままで良いのかは判断する。

以下は引数からは指定できないが、内部で使用する変数である。
AUCスコアの分布を元にuseGlobalとright_skewed_distributionのTRUE/FALSEを切り替えていく。

処理内部で扱う設定変数

useGlobal <- FALSE # <-- 1 componentの正規分布でのAUC閾値を採用するかどうか。
right_skewed_distribution <- FALSE # 分布が左に偏っているかどうか

nCells <- length(auc)
commentMsg <- "" # <- 警告を記録するための空ベクトル
aucThrs <- c() # AUC閾値候補を記録するための空ベクトル

AUCスコアが低い細胞ばかりの場合

AUCスコアが低い集団が非常に多い場合は警告が残る。解析が進まないわけではないが、このようなgenesetのAUCスコアは要注意。

- AUCスコアが0.2にも満たない細胞だらけの場合

# AUCスコアが0.2より大きい細胞の割合。signifで2桁に丸める。(小数点以下1桁まで)
highValCells <- signif(sum(auc > 0.2)/length(auc), 2)

# もしAUCが0.2より大きい細胞の割合が5%未満であれば、warningをつける
if (highValCells < 0.05) {
  commentMsg <- paste0(commentMsg, 
                       "Few cells have high AUC values (", 
                       highValCells, 
                       "% cells with AUC> 0.20). ",
                       sep = "")
  cat(commentMsg)
}

- AUCスコア0の細胞数が予測陰性割合より多い場合

# 予測される陰性割合
notPopPercent <- 1 - smallestPopPercent 

if (sum(auc == 0) > (nCells * notPopPercent)) {
  commentMsg <- paste0(commentMsg, 
                       round((sum(auc == 0)/nCells) * 100),
                       "% (more than ", notPopPercent, "%) of AUC are zero. ",
                       sep = "")
  cat(commentMsg)
  
  useGlobal <- TRUE  # Globalな分布でのAUC閾値を候補に残す。
  }

AUCスコア分布

AUCスコアの全体的な分布を確認して、設定用変数を切り替える。

- AUCスコアが正規分布に従うか

まずAUCスコアが正規分布に従っているのか確認する。本来、genesetスコアが低い細胞集団/高い集団のように二峰性に分かれるのであれば、正規分布には従わないはずである。

Kolmogorov-Smirnov Testsで確率分布を検定。
AUCのスコア分布と、AUC平均/AUC標準偏差に従う正規分布とで確率分布が同じかどうかを確かめる。

# AUCスコアの平均とSD取得
meanAUC <- mean(auc)
sdAUC <- sd(auc)

# AUCの平均とSDに従う正規分布からデータを100個取得したシミュレーションデータ
y <- rnorm(max(100, length(auc)), mean = meanAUC, sd = sdAUC)

# AUCスコアの分布と、正規分布のシミュレーションデータとで分布が異なるか検定
res_ks_test <- ks.test(x = auc, 
                       y = y,
                       alternative = "less"
                       )

有意水準以下なら2つの分布が異なる。棄却できなければAUCスコアの分布は正規分布とする。
AUCスコアが正規分布する場合は、randomなgenesetのように細胞特異性が低いことが考えられる。

もし正規分布しているなら、外れ値を分けるAUC閾値を算出する。具体的には累積確率が1-(thrP/nCells)になる確率変数(x軸)の値をAUC閾値候補として保存。(※ 1-(thrP/nCells)は殆ど1に近い。)

if (res_ks_test$p.value > 0.01) {
  # AUC閾値候補計算
  aucThrs["outlierOfGlobal"] <- qnorm(p = 1-(thrP/nCells),
                                      mean = meanAUC,
                                      sd = sdAUC)
　
  commentMsg <- paste0(commentMsg, 
                       "The AUC might follow a normal distribution (random gene-set?). "
                       )
  cat(commentMsg)
  
  # Globalな分布でのAUC閾値候補を残す
  useGlobal <- TRUE  
  }

- AUCスコアが左に偏っているか

AUCスコアの分布が大きく左に偏っているかを確認する。例えば大半の細胞がAUCスコアがゼロの場合はヒストグラムが左に偏り、右に長い裾を持つ。

1. ヒストグラムデータを取得。分割数は100個。

hist()関数にplot=FALSEオプションを入れると、plotされずにデータが返ってくる。breaks=100でbin数を100個に指定。

# 最大値を1にしたAUCスコアのヒストグラム。
# countは各binの頻度。breaksが100なので、0~1の間の100bin。
histogram <- hist(c(0, auc/max(auc)), breaks = 100, plot = FALSE)$count

2. ヒストグラムの最初の5 binに予測される陰性割合の75%以上が含まれる場合

大半の細胞のAUCスコアが0ならヒストグラムが強く左に偏っている。
right_skewed_distributionとuseGlobalをTRUEにして、全体分布から算出されるAUC閾値を候補として残すように設定する。

# 予測される陰性割合
notPopPercent <- 1 - smallestPopPercent 

if ((sum(histogram[1:5])/sum(histogram)) >= notPopPercent * 0.75) {
  right_skewed_distribution <- TRUE
  useGlobal <- TRUE
}

3. ヒストグラムの最初の10 binに、予測される陰性割合の50%以上含まれる場合

同様にright_skewed_distributionとuseGlobalをTRUEにしつつ、AUC閾値候補として「AUCスコアの最大値 * 0.1」を記録する。

if ((sum(histogram[1:10])/sum(histogram)) >= notPopPercent * 0.5) {
  right_skewed_distribution <- TRUE
  useGlobal <- TRUE
  
  aucThrs["tenPercentOfMax"] <- max(auc) * 0.1
}

上記の75%、50%は何を元に決定されたかはわからないが、経験的なものだろう。

AUCの確率密度関数から分布の谷を探す

AUCスコアが高い細胞集団、低い細胞集団が得られた場合、二峰性（もしくは三峰性）のピークをもつ分布が得られる。ピークとピークの間の谷の部分にあたるAUCスコアを探索して閾値候補とする。

具体的には、density()関数でAUCスコアの確率密度関数を推定し、その分布の確率密度(y軸)が「上昇->減少」「減少->上昇」に変わる確率変数(x軸)を探索する。
確率密度関数の滑らかさは adjust=引数によって変わる。(低いほどギザギザするのでターニングポイントとなる座標候補が増える)。AUCellのdensAdjust=引数で指定できる。基本はデフォルトのdensAdjust=2で良さそう。

1. aucの確率密度関数を推定

densCurveの中には、x軸512要素、y軸512要素の結果などが入っている。

density_Curve <- density(auc, adjust = densAdjust, cut = 0)

density()の返り値

2. 確率密度関数のturning pointを探索

次のコマンドで確率密度関数のyが「上昇->減少」「減少->上昇」に変化するターニングポイントを調べることができる。

turningPoints <- diff(sign(diff(density_Curve$y)))

上昇から減少に変わるxをpeak、減少から上昇に変わるxを谷とする。

# 上昇から減少に変わるところ -> peak
peaks_x <- which(turningPoints == -2)

# 減少から上昇に変わるところ。谷の部分。
valleys_x <- which(turningPoints == 2)

場合によっては、谷となる部分が無いこともある、その場合はこのAUC閾値探索はスキップする。

3. AUC閾値候補として相応しいか判断

AUCスコア分布が左に偏っているかどうかでも閾値設定を変える必要がある。もし左側に偏っているなら、左側の山の次の谷がAUCスコア閾値となるし、左側に偏っていなければ一番右側の山の左の谷がAUCスコア閾値となる。

## もし谷の部分があれば、AUC閾値を探索
if( length(valleys_x) > 0){

  # peakの確率密度
  peaks_y <- density_Curve$y[peaks_x]
  # 確率密度の高いピーク（山）を記録
  max_peak_x <- peaks_x[which.max(peaks_y)]
    
  ## AUC分布が左側に偏っているなら、左の大きい山が最も大きなピークとなりがち。その右側の最初の谷をAUC閾値とする
  if (right_skewed_distribution == TRUE){
    # ピークが高い山の次の谷の位置を取得
    tmp <- min(valleys_x[ valleys_x > max_peak_x ])
    # AUCスコアに直す
    density_Threshold <- density_Curve$x[tmp]
  } else {
    ## 分布が左に偏っていなければ、最も大きな山の直前の谷or直後の谷で左にある方をAUC閾値とする。
    tmp1 <- valleys_x[ valleys_x < max_peak_x ] # 最も大きな山以前の谷の位置を取得
    tmp1 <- tmp1[length(tmp1)] # 最後の谷を取得
    tmp2 <- valleys_x[ valleys_x > max_peak_x ] # 最も大きな山以降の谷の位置を取得
    # 一番左側にある谷の位置を採用
    tmp <- min(c(tmp1,tmp2))
    density_Threshold <- density_Curve$x[tmp] # AUCスコアに直す
  }
  
  ## AUC閾値より右側のAUC peakの中で最も大きいpeakの確率密度を確認
  ## AUC閾値より右側のピークが小さければ二峰性ではないのでこのAUC閾値は棄却
  
  # 右側のピークのxの値を取得
  peaks_next <- peaks_x[which(density_Curve$x[peaks_x] > density_Threshold)]
  # 右側のピークのyの値を取得
  peaks_next_y <- density_Curve$y[peaks_next]
  
  # AUC閾値の右側のピークが、全体の最も大きなピークの5%未満ならAUC閾値を白紙に
  if (( max(peaks_next_y ) / max(peaks_y) ) < 0.05) {
    density_Threshold <- NULL
  } else {
  # AUC閾値の右側のピークが5%以上の大きさであれば、AUC閾値候補として残す。
    aucThrs["valley"] <- density_Threshold
  }  
  
}

AUCスコア平均/標準偏差に従う正規分布でAUC閾値設定

AUCスコアが正規分布していなくとも、AUCスコアの平均/標準偏差に従う正規分布を考え、累積確率（1 - (thrP/nCells + smallestPopPercent)）を満たす確率変数 (AUCスコア)の値をAUC閾値候補とする。

smallestPopPercent=引数の値が結果に大きく影響する。デフォルトは"smallestPopPercent=0.25"なので、「1 - smallestPopPercent = 0.75」を満たすAUC閾値は大きな値になりやすい。(thrP=のデフォルトは0.01であり、thrP/nCellsはとても小さい値なので大きく影響はしない。)

global_Probability <- 1 - (thrP/nCells + smallestPopPercent)
aucThrs["Global_k1"] <- qnorm(p = global_Probability,
                              mean = meanAUC, 
                              sd = sdAUC
                              )

このAUC閾値候補は大きな値を取りやすいので、useGlobal変数のTRUE/FALSEに応じて不使用になる。

混合モデルを使った確率分布でAUC閾値設定

AUCスコアが二峰性や三峰性の分布を持つ場合、2つ3つの正規分布の混合分布として考えることができる。
mixtoolsパッケージのnormalmixEM()機能を使って正規分布の混合モデルを推定する。k=で何個のcomponentで分布が作られているか指定する。k=2とすると2つの正規分布が混合したモデルを想定して推定を行う。
tryCatch()を組んでいるように必ずしも推定がうまくいくわけではない。

# componentが2つの混合モデルを推定
MixModel_k2 <- tryCatch(
  mixtools::normalmixEM(x = auc, k = 2),
  error = function(e) {
    return(NULL)
    })
    
# componentが3つの混合モデルを推定
MixModel_k3 <- tryCatch(
  mixtools::normalmixEM(x = auc, k = 3),
  error = function(e) {
    return(NULL)
    })

推定の結果、各componentの分布が従う平均やSDが得られる。lambdaは混合モデルでの分布割合で合計が1になる。

k=2の例

k=3の例

2つの正規分布の混合モデルの場合、平均値が低いcomponent（左側の山）の右にAUC閾値を設定すれば、2つの混合モデルを分ける閾値となるはずである。そこで、左側の山の平均値と標準偏差に従う正規分布から累積確率（1 - (thrP/nCells)）を満たす確率変数（AUCスコア）を計算して、AUC閾値候補とする。「1 - (thrP/nCells)」は1に近いので左側の山をほぼ満たす箇所に閾値が設定できる。

componentが2つの混合モデルでAUC閾値候補計算

if (!is.null(MixModel_k2)) {
  # 左の山(2つの分布のうち、平均muが低い方)
  k2_Lt <- which.min(MixModel_k2$mu)
  # 左の山の平均とSDの分布で累積確率(1-(thrP/nCells))を満たす確率変数(x軸)の値
  aucThrs["MixModel_k2_Lt"] <- qnorm(p = (1 - (thrP/nCells)),
                                     mean = MixModel_k2$mu[k2_Lt],
                                     sd = MixModel_k2$sigma[k2_Lt]
                                     )
}

3つの正規分布が混合したモデルでは、最も平均値が高いcomponent（右側の山）の左にAUC閾値を設定する。右側の山の平均値と標準偏差に従う正規分布から累積確率（thrP）を満たす確率変数（AUCスコア）を計算して、AUC閾値候補とする。thrPはデフォルトでは0.01なので、右側の山を殆ど満たさない位置に閾値を設定できる。

componentが3つの混合モデルでAUC閾値候補計算

if (!is.null(MixModel_k3)) {
  # 一番右の山を取得 (3つの分布のうち、平均muが高い方)
  k3_Rt <- which.max(MixModel_k3$mu)
  # 一番右の山の分布で累積確率thrP (0.01)が得られる確率変数(x軸)の値
  # --> 左の山2つを含んで右の山を含まないような閾値設定
  k3_Rt_threshold <- qnorm(p = thrP,
                           mean = MixModel_k3$mu[k3_Rt], 
                           sd = MixModel_k3$sigma[k3_Rt]
                           )
  
  # 正の値であればR_k3として閾値に記録
  if (k3_Rt_threshold > 0) {
    aucThrs["MixModel_k3_Rt"] <- k3_Rt_threshold
  }
  }

【参考】
https://kusanagi.hatenablog.jp/entry/2017/01/24/193909
https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb
http://tnstdyrhacker.blogspot.com/2015/05/r.html

AUC閾値のフィルタリング

最終的なAUC閾値は、閾値候補の中で最も厳しい閾値（高い値）を採用する。(谷の閾値がある場合は、谷の閾値を使用する)

しっかりと二峰性の分布が確認できる場合の閾値と比べ、全体で正規分布を仮定した閾値（"Global_k1"）は非常に厳しい設定となってしまう。そこで、二峰性の分布と全体分布の位置をヒントにして、AUC閾値候補を残すかどうか決定する。

1. 二峰性の分布の右の山が十分大きい場合

右の山が大きい（AUCスコアが高い細胞が多い）場合、AUC候補として残し、左側の山の方が高い場合に、AUC候補から除く。(理由はよくわからなかった。。。)

if(!is.null(MixModel_k2)){
  component_Lt <- which.min(MixModel_k2$mu)
  component_Rt <- which.max(MixModel_k2$mu)
  
  # 0.4/標準偏差 * labmda
  # 標準偏差が小さい方がより大きな値を取る。ピークがsharpな方が高く出がち。
  # lambaは混合モデルでの分布割合。合計が1になる。lambdaが大きい方が高くでる。
  height1 <- 0.4/MixModel_k2$sigma[component_Lt] * MixModel_k2$lambda[component_Lt]
  height2 <- 0.4/MixModel_k2$sigma[component_Rt] * MixModel_k2$lambda[component_Rt]
    # --> つまりは分布の標準偏差と分布割合のどちらも考慮しての値
  
  # 上記計算結果が右の山の方が高ければ、Global AUC閾値も候補として残す
  if(height1 < height2){
    useGlobal <- TRUE
  }

2. 二峰性の分布と、混合モデルを想定しない正規分布の包含関係

二峰性の分布が得られていても2つの部分分布が均等（同じ高さ、同じ広がり）であれば、混合モデルを想定しない全体分布の中に含まれることになる。その場合は、全体分布でのAUC閾値を候補として残す。

if(!is.null(MixModel_k2)){
  component_Lt <- which.min(MixModel_k2$mu)
  component_Rt <- which.max(MixModel_k2$mu)
  
  ## Global分布と二峰性の分布の包含関係
  # 左側の分布がGlobal分布に含まれているかどうか。
  # 左側の分布の平均+1.5SDの値がGlobal分布の平均より大きければ、左側分布はGlobal分布に含まれている
  Global_include_LtComponent <- (
    Global_k1$mu < (MixModel_k2$mu[component_Lt] + (1.5 * MixModel_k2$sigma[component_Lt])))

  # Global分布に混合モデルの2つの分布が含まれるかどうか。
  # 混合モデルの左の山の平均がGlobal分布の-SD内 かつ 右の山の平均がGlobal分布の+SD内
  Global_include_2component <- (
      # Global平均-SDよりも左山平均が高いか
      ( (Global_k1$mu - Global_k1$sigma) < MixModel_k2$mu[component_Lt] ) &&
      # Global平均+SDよりも右山平均が低いか
      ( (Global_k1$mu + Global_k1$sigma) > MixModel_k2$mu[component_Rt] )
      )
  
  # Global分布が2つの分布を含んでいるなら、Global分布のAUC閾値を候補として残す
  if (Global_include_LtComponent && Global_include_2component) {
    useGlobal <- TRUE
    }
}

分布の位置関係に応じたメッセージが残るようになっている。

if(!is.null(MixModel_k2)){
  if (Global_include_LtComponent && Global_include_2component ) {
    commentMsg <- paste0(commentMsg,
                         "The global distribution overlaps the partial distributions. ")
    }
  if ( (height1 < height2) && !Global_include_2component ) {
    commentMsg <- paste0(commentMsg, 
                         "The right distribution is taller. ")
  }
  }

AUC閾値の最終決定

確率密度関数で谷の閾値が見つかっていればそれが優先して採用される。ただし、谷と採用したAUC閾値より右側にも起伏があれば、分布が複数あることを警告するメッセージが出る。

if ("valley" %in% names(aucThrs)){
  aucThr <- aucThrs["valley"]
}

# もしvalley (最も左側の谷のx軸の箇所)が採用された場合、
# valleyが閾値なのに、確率密度が1を超えるような山が3つ以上あって、さらに起伏がちゃんとあるなら警告を出す。
if ((length(aucThr) > 0) && (names(aucThr) == "valley")) {
  
  # y軸の値が1より大きい山のみを取得
  peaks_x <- peaks_x[which(density_Curve$y[peaks_x] > 1)]
  
  # ピークが2つより多いなら
  if (length(peaks_x) > 2) {
    
    # 谷とその前の山のピークの対応表
    tmp <- cbind(
      valleys_x[seq(length(peaks_x) - 1)], # peaks_xの要素数-1のvalley_x
      peaks_x[-1] # peaks_xの最後の1つは除いたもの
      )
    
    # 谷の確率密度y と その前の山のピークの確率密度y の比が1.5倍より大きければwarning
    FCs <- density_Curve$y[tmp[, 2]]/density_Curve$y[tmp[, 1]]
    if (any(FCs > 1.5)) 
      warning(gSetName, ":\tCheck the AUC histogram. ", 
        "'valley' was selected as the best threshold, ", 
        "but there might be several distributions in the AUC.")
  }
}

谷の閾値が無い場合は、閾値候補のうち最も値が大きいものを使用する。

if (length(aucThr) == 0) {
  aucThr <- aucThrs[which.max(aucThrs)]
}

# それでもまだ未決定なら、AUC閾値は1
if (length(aucThr) == 0) {
  aucThr <- 1
}

複数genesetに繰り返し処理

一応、上記を複数geneset版に繰り返し処理を組んだものを記載しておく。

gSetNames <- rownames(aucMatrix)
res_list <- lapply(gSetNames, function(gSetName){
  
  useGlobal <- FALSE # <-- 1 componentの正規分布でのAUC閾値を採用するかどうか。
  right_skewed_distribution <- FALSE # 分布が左に偏っているかどうか
  commentMsg <- "" # <- 警告を記録するための空ベクトル
  aucThrs <- c() # AUC閾値候補を記録するための空ベクトル
  aucThr <- c() # 採用したAUC閾値を入れる用
  
  # 1つのgenesetのベクトル取得 (要素にAUCスコア、名前属性に細胞バーコード)
  auc <- aucMatrix[gSetName ,]
  
  # AUCスコアが小さい方から昇順に並び替え
  auc <- sort(auc)

  nCells <- length(auc)
  
  # AUCスコアが0.2より大きい細胞の割合。signifで2桁に丸める。(小数点以下1桁まで)
  highValCells <- signif(sum(auc > 0.2)/length(auc), 2)

  # もしAUCが0.2より大きい細胞の割合が5%未満であれば、warningをつける
  if (highValCells < 0.05) {
    commentMsg <- paste0(commentMsg, 
                         "Few cells have high AUC values (", 
                         highValCells, 
                         "% cells with AUC> 0.20). ",
                         sep = "")
    cat(commentMsg)
  }
  
    
  # 予測される陰性割合
  notPopPercent <- 1 - smallestPopPercent 
  
  if (sum(auc == 0) > (nCells * notPopPercent)) {
    commentMsg <- paste0(commentMsg, 
                         round((sum(auc == 0)/nCells) * 100),
                         "% (more than ", notPopPercent, "%) of AUC are zero. ",
                         sep = "")
    cat(commentMsg)
    
    useGlobal <- TRUE  # Globalな分布でのAUC閾値を候補に残す。
  }
  
  # AUCスコアの平均とSD取得
  meanAUC <- mean(auc)
  sdAUC <- sd(auc)
  
  # AUCの平均とSDに従う正規分布からデータを100個取得したシミュレーションデータ
  y <- rnorm(max(100, length(auc)), mean = meanAUC, sd = sdAUC)
  
  # AUCスコアの分布と、正規分布のシミュレーションデータとで分布が異なるか検定
  res_ks_test <- ks.test(x = auc, 
                         y = y,
                         alternative = "less"
                         )
  
  if (res_ks_test$p.value > 0.01) {
    # AUC閾値候補計算
    aucThrs["outlierOfGlobal"] <- qnorm(p = 1-(thrP/nCells),
                                        mean = meanAUC,
                                        sd = sdAUC)
  　
    commentMsg <- paste0(commentMsg, 
                         "The AUC might follow a normal distribution (random gene-set?). "
                         )
    cat(commentMsg)
    
    # Globalな分布でのAUC閾値候補を残す
    useGlobal <- TRUE  
    }
  
  # 最大値を1にしたAUCスコアのヒストグラム。plot=FALSEにするとplotされずにデータが返ってくる。
  # countは各binの頻度。breaksが100なので、0~1の間の100bin。
  histogram <- hist(c(0, auc/max(auc)), breaks = 100, plot = FALSE)$count
  
  # 予測される陰性割合
  notPopPercent <- 1 - smallestPopPercent 
  
  if ((sum(histogram[1:5])/sum(histogram)) >= notPopPercent * 0.75) {
    right_skewed_distribution <- TRUE
    useGlobal <- TRUE
  }
  
  if ((sum(histogram[1:10])/sum(histogram)) >= notPopPercent * 0.5) {
    right_skewed_distribution <- TRUE
    useGlobal <- TRUE
    
    aucThrs["tenPercentOfMax"] <- max(auc) * 0.1
  }
  
  density_Curve <- density(auc, adjust = densAdjust, cut = 0)
  turningPoints <- diff(sign(diff(density_Curve$y)))
  # 上昇から減少に変わるところ -> peak
  peaks_x <- which(turningPoints == -2)
  
  # 減少から上昇に変わるところ。谷の部分。
  valleys_x <- which(turningPoints == 2)
  
  ## もし谷の部分があれば、AUC閾値を探索
  if( length(valleys_x) > 0){
  
    # peakの確率密度
    peaks_y <- density_Curve$y[peaks_x]
    # 確率密度の高いピーク（山）を記録
    max_peak_x <- peaks_x[which.max(peaks_y)]
      
    ## AUC分布が左側に偏っているなら、左の大きい山が最も大きなピークとなりがち。その右側の最初の谷をAUC閾値とする
    if (right_skewed_distribution == TRUE){
      # ピークが高い山の次の谷の位置を取得
      tmp <- min(valleys_x[ valleys_x > max_peak_x ])
      # AUCスコアに直す
      density_Threshold <- density_Curve$x[tmp]
    } else {
      ## 分布が左に偏っていなければ、最も大きな山の直前の谷or直後の谷で左にある方をAUC閾値とする。
      tmp1 <- valleys_x[ valleys_x < max_peak_x ] # 最も大きな山以前の谷の位置を取得
      tmp1 <- tmp1[length(tmp1)] # 最後の谷を取得
      tmp2 <- valleys_x[ valleys_x > max_peak_x ] # 最も大きな山以降の谷の位置を取得
      # 一番左側にある谷の位置を採用
      tmp <- min(c(tmp1,tmp2))
      density_Threshold <- density_Curve$x[tmp] # AUCスコアに直す
    }
    
    ## AUC閾値より右側のAUC peakの中で最も大きいpeakの確率密度を確認
    ## AUC閾値より右側のピークが小さければ二峰性ではないのでこのAUC閾値は棄却
    
    # 右側のピークのxの値を取得
    peaks_next <- peaks_x[which(density_Curve$x[peaks_x] > density_Threshold)]
    # 右側のピークのyの値を取得
    peaks_next_y <- density_Curve$y[peaks_next]
    
    # AUC閾値の右側のピークが、全体の最も大きなピークの5%未満ならAUC閾値を白紙に
    if (( max(peaks_next_y ) / max(peaks_y) ) < 0.05) {
      density_Threshold <- NULL
    } else {
    # AUC閾値の右側のピークが5%以上の大きさであれば、AUC閾値候補として残す。
      aucThrs["valley"] <- density_Threshold
    }  
    
  }
  
  global_Probability <- 1 - (thrP/nCells + smallestPopPercent)
  aucThrs["Global_k1"] <- qnorm(p = global_Probability,
                                mean = meanAUC, 
                                sd = sdAUC
                                )
  
  
  # componentが2つの混合モデルを推定
  MixModel_k2 <- tryCatch(
    mixtools::normalmixEM(x = auc, k = 2),
    error = function(e) {
      return(NULL)
      })
      
  # componentが3つの混合モデルを推定
  MixModel_k3 <- tryCatch(
    mixtools::normalmixEM(x = auc, k = 3),
    error = function(e) {
      return(NULL)
      })
  
  
  if (!is.null(MixModel_k2)) {
    # 左の山(2つの分布のうち、平均muが低い方)
    k2_Lt <- which.min(MixModel_k2$mu)
    # 左の山の平均とSDの分布で累積確率(1-(thrP/nCells))を満たす確率変数(x軸)の値
    aucThrs["MixModel_k2_Lt"] <- qnorm(p = (1 - (thrP/nCells)),
                                       mean = MixModel_k2$mu[k2_Lt],
                                       sd = MixModel_k2$sigma[k2_Lt]
                                       )
  }
  
  if (!is.null(MixModel_k3)) {
    # 一番右の山を取得 (3つの分布のうち、平均muが高い方)
    k3_Rt <- which.max(MixModel_k3$mu)
    # 一番右の山の分布で累積確率thrP (0.01)が得られる確率変数(x軸)の値
    # --> 左の山2つを含んで右の山を含まないような閾値設定
    k3_Rt_threshold <- qnorm(p = thrP,
                             mean = MixModel_k3$mu[k3_Rt], 
                             sd = MixModel_k3$sigma[k3_Rt]
                             )
    
    # 正の値であればR_k3として閾値に記録
    if (k3_Rt_threshold > 0) {
      aucThrs["MixModel_k3_Rt"] <- k3_Rt_threshold
    }
  }
  
  if(!is.null(MixModel_k2)){
    component_Lt <- which.min(MixModel_k2$mu)
    component_Rt <- which.max(MixModel_k2$mu)
    
    # 0.4/標準偏差 * labmda
    # 標準偏差が小さい方がより大きな値を取る。ピークがsharpな方が高く出がち。
    # lambaは混合モデルでの分布割合。合計が1になる。lambdaが大きい方が高くでる。
    height1 <- 0.4/MixModel_k2$sigma[component_Lt] * MixModel_k2$lambda[component_Lt]
    height2 <- 0.4/MixModel_k2$sigma[component_Rt] * MixModel_k2$lambda[component_Rt]
      # --> つまりは分布の標準偏差と分布割合のどちらも考慮しての値
    
    # 上記計算結果が右の山の方が高ければ、Global AUC閾値も候補として残す
    if(height1 < height2){
      useGlobal <- TRUE
    }
    }
    
    if(!is.null(MixModel_k2)){
    component_Lt <- which.min(MixModel_k2$mu)
    component_Rt <- which.max(MixModel_k2$mu)
    
    ## Global分布と二峰性の分布の包含関係
    # 左側の分布がGlobal分布に含まれているかどうか。
    # 左側の分布の平均+1.5SDの値がGlobal分布の平均より大きければ、左側分布はGlobal分布に含まれている
    Global_include_LtComponent <- (
      Global_k1$mu < (MixModel_k2$mu[component_Lt] + (1.5 * MixModel_k2$sigma[component_Lt])))
  
    # Global分布に混合モデルの2つの分布が含まれるかどうか。
    # 混合モデルの左の山の平均がGlobal分布の-SD内 かつ 右の山の平均がGlobal分布の+SD内
    Global_include_2component <- (
        # Global平均-SDよりも左山平均が高いか
        ( (Global_k1$mu - Global_k1$sigma) < MixModel_k2$mu[component_Lt] ) &&
        # Global平均+SDよりも右山平均が低いか
        ( (Global_k1$mu + Global_k1$sigma) > MixModel_k2$mu[component_Rt] )
        )
    
    # Global分布が2つの分布を含んでいるなら、Global分布のAUC閾値を候補として残す
    if (Global_include_LtComponent && Global_include_2component) {
      useGlobal <- TRUE
      }
    }
    
    if(!is.null(MixModel_k2)){
      if (Global_include_LtComponent && Global_include_2component ) {
        commentMsg <- paste0(commentMsg,
                             "The global distribution overlaps the partial distributions. ")
        }
      if ( (height1 < height2) && !Global_include_2component ) {
        commentMsg <- paste0(commentMsg, 
                             "The right distribution is taller. ")
    }
    }
    
    if ("valley" %in% names(aucThrs)){
      aucThr <- aucThrs["valley"]
      }
  
    # もしvalley (最も左側の谷のx軸の箇所)が採用された場合、
    # valleyが閾値なのに、確率密度が1を超えるような山が3つ以上あって、さらに起伏がちゃんとあるなら警告を出す。
    if ((length(aucThr) > 0) && (names(aucThr) == "valley")) {
      
      # y軸の値が1より大きい山のみを取得
      peaks_x <- peaks_x[which(density_Curve$y[peaks_x] > 1)]
      
      # ピークが2つより多いなら
      if (length(peaks_x) > 2) {
        
        # 谷とその前の山のピークの対応表
        tmp <- cbind(
          valleys_x[seq(length(peaks_x) - 1)], # peaks_xの要素数-1のvalley_x
          peaks_x[-1] # peaks_xの最後の1つは除いたもの
          )
        
        # 谷の確率密度y と その前の山のピークの確率密度y の比が1.5倍より大きければwarning
        FCs <- density_Curve$y[tmp[, 2]]/density_Curve$y[tmp[, 1]]
        if (any(FCs > 1.5)) 
          warning(gSetName, ":\tCheck the AUC histogram. ", 
            "'valley' was selected as the best threshold, ", 
            "but there might be several distributions in the AUC.")
      }
    }
  
  
    if (length(aucThr) == 0) {
        aucThr <- aucThrs[which.max(aucThrs)]
      }
      
      # それでもまだ未決定なら、AUC閾値は1
      if (length(aucThr) == 0) {
        aucThr <- 1
      }
    
    positive_cells <- names(auc[auc > aucThr])
      
    return(list(selected = aucThr,
                thresholds = aucThrs, 
                comment = commentMsg,
                assignment = positive_cells))
    }
  )
names(res_list) <- gSetNames