首先我們來看看對AUcell的介紹
AUCell可以識別單細胞RNA序列數(shù)據(jù)中具有活躍基因集(例如signatures鸟廓,基因模塊...)的細胞谚中。 AUCell使用“曲線下面積”(AUC)來計算輸入基因集的關(guān)鍵子集是否在每個細胞的表達基因中富集院刁。 AUC分數(shù)在所有細胞中的分布允許探索特征的相對表達缤言。 由于計分方法是基于排名的滥嘴,因此AUCell不受基因表達單位和標(biāo)準化程序的影響县耽。 此外动羽,由于對cell進行了單獨評估包帚,因此可以輕松地將其應(yīng)用于更大的數(shù)據(jù)集,并根據(jù)需要對表達式矩陣進行分組运吓。
那也就是說渴邦,AUcell是分析感興趣的基因集在所有細胞是否存在富集,道理很簡單拘哨,我們來看看主要分析內(nèi)容
用AUcell進行分析分三步:
1谋梭、 Build the rankings
2、Calculate the Area Under the Curve (AUC)
3倦青、Set the assignment thresholds
我們首先來分析第一步:AUcell首先是對cell怎么排序的瓮床。
For each cell, the genes are ranked from highest to lowest value. The genes with same expression value are shuffled. Therefore, genes with expression ‘0’ are randomly sorted at the end of the ranking. It is important to check that most cells have at least the number of expressed/detected genes that are going to be used to calculate the AUC (aucMaxRank in calcAUC()). The histogram provided by AUCell_buildRankings() allows to quickly check this distribution. plotGeneCount(exprMatrix) allows to obtain only the plot before building the rankings.
這個地方我們可以看出,對于每個基因产镐,從高到底進行排名纤垂,也就是說每個細胞的每個基因都有一個排名,得到一個排名的矩陣磷账。
接下來第二步:
- Calculate enrichment for the gene signatures (AUC)
首先我們要知道什么是AUC曲線:
AUC(Area Under Curve)被定義為ROC下與坐標(biāo)軸圍成的面積峭沦,顯然這個面積的數(shù)值不會大于1。又由于ROC曲線一般都處于y=x這條直線的上方逃糟,所以AUC的取值范圍在0.5和1之間吼鱼。AUC越接近1.0,檢測方法真實性越高;等于0.5時绰咽,則真實性最低菇肃,無應(yīng)用價值。
那么什么是ROC曲線呢取募?
一琐谤、 ROC曲線的由來
很多學(xué)習(xí)器是為測試樣本產(chǎn)生一個實值或概率預(yù)測,然后將這個預(yù)測值與一個分類閾值進行比較玩敏,若大于閾值則分為正類斗忌,否則為反類质礼。例如,神經(jīng)網(wǎng)絡(luò)在一般情形下是對每個測試樣本預(yù)測出一個[0.0织阳,1.0]之間的實值眶蕉,然后將這個值與閾值0.5進行比較,大于0.5則判為正例唧躲,否則為反例造挽。這個閾值設(shè)置的好壞,直接決定了學(xué)習(xí)器的泛化能力弄痹。
在不同的應(yīng)用任務(wù)中饭入,我們可根據(jù)任務(wù)需求來采用不同的閾值。例如肛真,若我們更重視“查準率”谐丢,則可以把閾值設(shè)置的大一些,讓分類器的預(yù)測結(jié)果更有把握毁欣;若我們更重視“查全率”庇谆,則可以把閾值設(shè)置的小一些岳掐,讓分類器預(yù)測出更多的正例凭疮。因此,閾值設(shè)置的好壞串述,體現(xiàn)了綜合考慮學(xué)習(xí)器在不同任務(wù)下的泛化性能的好壞执解。為了形象的描述這一變化,在此引入ROC曲線纲酗,ROC曲線則是從閾值選取角度出發(fā)來研究學(xué)習(xí)器泛化性能的有力工具衰腌。
如果你還對“查準率”和“查全率”不了解,看我之前的文章【錯誤率觅赊、精度右蕊、查準率、查全率和F1度量】詳細介紹
二吮螺、 什么是ROC曲線
ROC全稱是“受試者工作特征”(Receiver OperatingCharacteristic)曲線饶囚。我們根據(jù)學(xué)習(xí)器的預(yù)測結(jié)果,把閾值從0變到最大鸠补,即剛開始是把每個樣本作為正例進行預(yù)測萝风,隨著閾值的增大,學(xué)習(xí)器預(yù)測正樣例數(shù)越來越少紫岩,直到最后沒有一個樣本是正樣例规惰。在這一過程中,每次計算出兩個重要量的值泉蝌,分別以它們?yōu)闄M歇万、縱坐標(biāo)作圖揩晴,就得到了“ROC曲線”。
** ROC曲線的縱軸是“真正例率”(True Positive Rate, 簡稱TPR)堕花,橫軸是“假正例率”(False Positive Rate,簡稱FPR)文狱,**基于上篇文章《錯誤率、精度缘挽、查準率瞄崇、查全率和F1度量》的表1中符號,兩者分別定義為:
顯示ROC曲線的圖稱為“ROC圖”壕曼。圖1給出了一個示意圖苏研,顯然,對角線對應(yīng)于“隨機猜測”模型腮郊,而點(0,1)則對應(yīng)于將所有正例預(yù)測為真正例摹蘑、所有反例預(yù)測為真反例的“理想模型”。
圖1:ROC曲線與AUC面積
現(xiàn)實任務(wù)中通常是利用有限個測試樣例來繪制ROC圖轧飞,此時僅能獲得有限個(真正例率衅鹿,假正例率)坐標(biāo)對,無法產(chǎn)生圖1中的光滑ROC曲線过咬,只能繪制出圖2所示的近似ROC曲線大渤。繪制過程很簡單:給定m+ 個正例和m-個反例,根據(jù)學(xué)習(xí)器預(yù)測結(jié)果對樣例進行排序掸绞,然后把分類閾值設(shè)置為最大泵三,即把所有樣例均預(yù)測為反例,此時真正例率和假正例率均為0衔掸,在坐標(biāo)(0,0)處標(biāo)記一個點烫幕。然后,將分類閾值依次設(shè)為每個樣例的預(yù)測值敞映,即依次將每個樣例劃分為正例较曼。設(shè)前一個標(biāo)記點坐標(biāo)為
三、 ROC曲線的意義
(1)主要作用
1. ROC曲線能很容易的查出任意閾值對學(xué)習(xí)器的泛化性能影響栓霜。
2.有助于選擇最佳的閾值翠桦。ROC曲線越靠近左上角,模型的查全率就越高。最靠近左上角的ROC曲線上的點是分類錯誤最少的最好閾值销凑,其假正例和假反例總數(shù)最少丛晌。
3.可以對不同的學(xué)習(xí)器比較性能。將各個學(xué)習(xí)器的ROC曲線繪制到同一坐標(biāo)中斗幼,直觀地鑒別優(yōu)劣澎蛛,靠近左上角的ROC曲所代表的學(xué)習(xí)器準確性最高。
(2)優(yōu)點
1. 該方法簡單蜕窿、直觀谋逻、通過圖示可觀察分析方法的準確性,并可用肉眼作出判斷桐经。ROC曲線將真正例率和假正例率以圖示方法結(jié)合在一起毁兆,可準確反映某種學(xué)習(xí)器真正例率和假正例率的關(guān)系,是檢測準確性的綜合代表阴挣。
2. 在生物信息學(xué)上的優(yōu)點:ROC曲線不固定閾值悯嗓,允許中間狀態(tài)的存在铆惑,利于使用者結(jié)合專業(yè)知識薄霜,權(quán)衡漏診與誤診的影響牧牢,選擇一個更加的閾值作為診斷參考值。
四誓沸、 AUC面積的由來
如果兩條ROC曲線沒有相交梅桩,我們可以根據(jù)哪條曲線最靠近左上角哪條曲線代表的學(xué)習(xí)器性能就最好。但是蔽介,實際任務(wù)中摘投,情況很復(fù)雜煮寡,如果兩條ROC曲線發(fā)生了交叉虹蓄,則很難一般性地斷言誰優(yōu)誰劣。在很多實際應(yīng)用中幸撕,我們往往希望把學(xué)習(xí)器性能分出個高低來薇组。在此引入AUC面積。
在進行學(xué)習(xí)器的比較時坐儿,若一個學(xué)習(xí)器的ROC曲線被另一個學(xué)習(xí)器的曲線完全“包住”律胀,則可斷言后者的性能優(yōu)于前者;若兩個學(xué)習(xí)器的ROC曲線發(fā)生交叉貌矿,則難以一般性的斷言兩者孰優(yōu)孰劣炭菌。此時如果一定要進行比較,則比較合理的判斷依據(jù)是比較ROC曲線下的面積逛漫,即AUC(Area Under ROC Curve)黑低,如圖1圖2所示。
五、 什么是AUC面積
** AUC就是ROC曲線下的面積克握,衡量學(xué)習(xí)器優(yōu)劣的一種性能指標(biāo)蕾管。**從定義可知,AUC可通過對ROC曲線下各部分的面積求和而得菩暗。假定ROC曲線是由坐標(biāo)為的點按序連接而形成掰曾,參見圖2,則AUC可估算為公式3停团。
六旷坦、 AUC面積的意義
AUC是衡量二分類模型優(yōu)劣的一種評價指標(biāo),表示預(yù)測的正例排在負例前面的概率佑稠。
看到這里塞蹭,是不是很疑惑,根據(jù)AUC定義和計算方法讶坯,怎么和預(yù)測的正例排在負例前面的概率扯上聯(lián)系呢番电?如果從定義和計算方法來理解AUC的含義,比較困難辆琅,實際上AUC和Mann-WhitneyU test(曼-慧特尼U檢驗)有密切的聯(lián)系漱办。從Mann-Whitney U statistic的角度來解釋,AUC就是從所有正樣本中隨機選擇一個樣本婉烟,從所有負樣本中隨機選擇一個樣本娩井,然后根據(jù)你的學(xué)習(xí)器對兩個隨機樣本進行預(yù)測,把正樣本預(yù)測為正例的概率p1似袁,把負樣本預(yù)測為正例的概率p2,p1 > p2的概率就等于AUC洞辣。所以AUC反映的是分類器對樣本的排序能力。根據(jù)這個解釋昙衅,如果我們完全隨機的對樣本分類扬霜,那么AUC應(yīng)該接近0.5。
另外值得注意的是而涉,AUC的計算方法同時考慮了學(xué)習(xí)器對于正例和負例的分類能力著瓶,在樣本不平衡的情況下,依然能夠?qū)Ψ诸惼髯龀龊侠淼脑u價啼县。AUC對樣本類別是否均衡并不敏感材原,這也是不均衡樣本通常用AUC評價學(xué)習(xí)器性能的一個原因。例如在癌癥預(yù)測的場景中季眷,假設(shè)沒有患癌癥的樣本為正例余蟹,患癌癥樣本為負例,負例占比很少(大概0.1%)子刮,如果使用準確率評估威酒,把所有的樣本預(yù)測為正例便可以獲得99.9%的準確率。但是如果使用AUC,把所有樣本預(yù)測為正例兼搏,TPR為1卵慰,F(xiàn)PR為1。這種情況下學(xué)習(xí)器的AUC值將等于0.5佛呻,成功規(guī)避了樣本不均衡帶來的問題裳朋。
最后,我們在討論一下:在多分類問題下能不能使用ROC曲線來衡量模型性能吓著?
我的理解:ROC曲線用在多分類中是沒有意義的鲤嫡。只有在二分類中Positive和Negative同等重要時候,適合用ROC曲線評價绑莺。如果確實需要在多分類問題中用ROC曲線的話暖眼,可以轉(zhuǎn)化為多個“一對多”的問題。即把其中一個當(dāng)作正例纺裁,其余當(dāng)作負例來看待诫肠,畫出多個ROC曲線。
回到我們生信分析的第二步欺缘,對于AUC的計算
To determine whether the gene set is enriched at the top of the gene-ranking for each cell, AUCell uses the “Area Under the Curve” (AUC) of the recovery curve.
In order to calculate the AUC, by default only the top 5% of the genes in the ranking are used (i.e. checks whether the genes in the gene-set or signature are within the top 5%). This allows faster execution on bigger datasets, and reduce the effect of the noise at the bottom of the ranking (e.g. where many genes might be tied at 0 counts). The percentage to take into account can be modified with the argument aucMaxRank. For datasets where most cells express many genes (e.g. a filtered dataset), or these have high expression values, it might be good to increase this threshold. Check the histogram provided by AUCell_buildRankings to get an estimation on where this threshold lies within the dataset.
這里我們深入理解一下栋豫,對于一個細胞進行基因排序之后,得到下面的圖:
也就是說谚殊,看看排名前250位的基因在排名慢慢變大的過程中丧鸯,包含基因集數(shù)量的變化,而形成了一條曲線嫩絮,曲線下面的面積丛肢,就是我們得到的富集分數(shù)。(似乎和剛才講的ROC沒什么關(guān)系~~)
The AUC estimates the proportion of genes in the gene-set that are highly expressed in each cell. Cells expressing many genes from the gene-set will have higher AUC values than cells expressing fewer (i.e. compensating for housekeeping genes, or genes that are highly expressed in all the cells in the dataset). Since the AUC represents the proportion of expressed genes in the signature, we can use the relative AUCs across the cells to explore the population of cells that are present in the dataset according to the expression of the gene-set.
However, determining whether the signature is active (or not) in a given cell is not always trivial. The AUC is not an absolute value, but it depends on the the cell type (i.e. sell size, amount of transcripts), the specific dataset (i.e. sensitivity of the measures) and the gene-set. It is often not straight forward to obtain a pruned signature of clear marker genes that are completely “on” in the cell type of interest and off" in every other cell. In addition, at single-cell level, most genes are not expressed or detected at a constant level.
The ideal situation will be a bi-modal distribution, in which most cells in the dataset have a low “AUC” compared to a population of cells with a clearly higher value (i.e. see “Oligodendrocites” in the next figure). This is normally the case on gene sets that are active mostly in a population of cells with a good representation in the dataset (e.g. ~ 5-30% of cells in the dataset). Similar cases of “marker” gene sets but with different proportions of cells in the datasets are the “neurons” and “microglia” (see figure). When there are very few cells within the dataset, the distribution might look normal-like, but with some outliers to the higher end (e.g. microglia). While if the gene set is marker of a high percentage of cells in the dataset (i.e. neurons), the distribution might start approaching the look of a gene-set of housekeeping genes. As example, the ‘housekeeping’ gene-set in the figure includes genes that are detected in most cells in the dataset.
Note that the size of the gene-set will also affect the results. With smaller gene-genes (fewer genes), it is more likely to get cells with AUC = 0. While this is the case of the “perfect markers” it is also easier to get it by chance with smal datasets. (i.e. Random gene set with 50 genes in the figure). Bigger gene-sets (100-2k) can be more stable and easier to evaluate, as big random gene sets will approach the normal distibution.
To ease the exploration of the distributions, the function AUCell_exploreThresholds() automatically plots all the histograms and calculates several thresholds that could be used to consider a gene-set ‘a(chǎn)ctive’ (returned in aucThr). The distributions are plotted as dotted lines over the histogram and the corresponding thresholds as vertical bars in the matching color. The thicker vertical line indicates the threshold selected by default (selected): the highest value to reduce the false positives.
而將具有較高AUC值反映到細胞中剿干,可以得到
至于代碼很簡單:
library(AUCell)
cells_rankings <- AUCell_buildRankings(exprMatrix)##基因排序
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings, aucMaxRank=nrow(cells_rankings)*0.05) ##計算AUC值蜂怎。
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, nCores=1,assign=TRUE)##挑選閾值
至于基因集的選擇,可以借助于hallmark怨愤,用于研究腫瘤特征分析派敷。