對應(yīng)原版教程第11章http://bioconductor.org/books/release/OSCA/overview.html
從基因表達水平來看，主要是哪些基因造成了每個cluster得以區(qū)分的結(jié)果是這一步的主要目的梢什。實現(xiàn)這一步后默色，對于cluster的生物學(xué)意義解讀將有很大幫助。

源網(wǎng)凛膏，侵刪~

筆記要點

1、cluster marker gene的意義
2脏榆、基于t檢驗
3猖毫、其它檢驗方法

1、cluster marker gene的意義

• 物以類聚须喂，人以群分吁断。細胞因gene表達相似而歸為一類；因不同的表達特征而形成多個cluster坞生。因此仔役，這一步的目的就是找出造成每個cluster與其它clusters分離的marker gene---每個cluster的特征基因集(高表達/低表達)∈羌海可以為下一步的cluster細胞類型注釋提供主要參考信息又兵。
• cluster marker gene的鑒定思路類似于 Bulk RNA-seq的差異分析。不同之處在于往往會有十多個cluster卒废，因此定義cluster的marker gene有蠻多細節(jié)值得考慮沛厨。
• 并且由于scRNA-seq與Bulk RNA-seq的不同，一般不推薦使用諸如DESeq2摔认、edgeR等差異分析思路（原因在后面會提到）逆皮。本篇筆記主要學(xué)習(xí)了使用scran包的findMarkers()函數(shù)提供的基于不同method的鑒定思路。

#示例數(shù)據(jù)
sce.pbmc #獲取方式見原教程
#class: SingleCellExperiment 
#dim: 33694 3985

table(sce.pbmc$label)
#  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
#205 508 541  56 374 125  46 432 302 867  47 155 166  61  84  16

2参袱、基于t檢驗

• Welch t-test是基于兩組的均值是否相同的假設(shè)檢驗思路电谣；為findMarkers()函數(shù)的默認方法
• 比較思路是首先是選擇其中的一個cluster與其余所有的cluster分別兩兩比較秽梅。就比較結(jié)果而言，定義cluster的marker gene有如下三個角度可供選擇（對應(yīng)的零假設(shè)也稍有不同）

2.1 any differences(generous)

• 零假設(shè)：對于基因A辰企，cluster X與其余全部cluster的表達均值均相同
• 換句話說：只要cluster X的基因A的表達水平與其余任意一個cluster的t檢驗結(jié)果顯著风纠，就認為是該cluster的marker gene。

library(scran)
markers.pbmc <- findMarkers(sce.pbmc)
markers.pbmc #分別儲存每個cluster的分析結(jié)果
#List of length 16
#names(16): 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

• 以cluster 7的分析結(jié)果為例

chosen <- "7"
interesting <- markers.pbmc[[chosen]]
colnames(interesting)
# [1] "Top"           "p.value"       "FDR"           "summary.logFC" "logFC.1"      
# [6] "logFC.2"       "logFC.3"       "logFC.4"       "logFC.5"       "logFC.6"      
# [11] "logFC.8"       "logFC.9"       "logFC.10"      "logFC.11"      "logFC.12"     
# [16] "logFC.13"      "logFC.14"      "logFC.15"      "logFC.16"

head(interesting[,1:5],15)
# DataFrame with 15 rows and 5 columns
# Top      p.value          FDR summary.logFC    logFC.1
# <integer>    <numeric>    <numeric>     <numeric>  <numeric>
#   S100A4             1  2.59737e-38  1.27018e-36      -4.27560 -2.1597143
# TAGLN2             1  8.65033e-28  2.44722e-26       5.07327  4.7714408
# FCGR3A             1  8.84356e-63  1.15048e-60      -3.07121 -1.2852520
# GZMA               1 1.15392e-120 7.20000e-118      -1.92877 -2.5023377
# HLA-DQA1           1  3.43640e-83  8.90663e-81      -3.54890 -0.0220385
# ...              ...          ...          ...           ...        ...
# CTSS               2  1.58027e-33  6.15557e-32      -3.51820  -0.242616
# HOPX               2  4.60496e-78  1.05551e-75      -1.99060  -1.990602
# PF4                2  9.57857e-35  3.97464e-33       6.71811   6.862880
# PRELID1            2 1.10561e-107 5.47832e-105      -1.61240  -0.761206
# AC090498.1         2 1.19635e-156 1.55038e-153      -1.93799  -0.872609

a7=as.data.frame(interesting)
dim(a7)  #[1] 33694    19
View(a7)

如下圖牢贸，不同列的釋義

• Top：cluster 7 分別與其它所有cluster的t-test結(jié)果里，P值最顯著的基因镐捧。按理潜索，應(yīng)當有15個Top1 gene。但如圖可以看到只有10個Top1懂酱，我認為是由于存在重復(fù)的gene結(jié)果導(dǎo)致竹习。并且在前面出現(xiàn)過的基因，再后面的TopX就不再出現(xiàn)了列牺。
• p.value：偏離零假設(shè)的程度整陌，具體計算是結(jié)合特定cluster與其余所有cluster的兩兩t檢驗p值的combined p value（Simon method 多重檢驗）
• FDR：p.value的FDR值
• summary.logFC：即最顯著的那次cluster兩兩比較的logFoldChange
• logFC.X：該基因分別與其它所有cluster的logFC結(jié)果
  

這種思路定義的cluster marker gene標準是比較寬松的。因為只有一個基因在任意一次cluster間比較p值顯著瞎领，就會認為是marker gene

2.2 all difference(stringent)

• 零假設(shè)：對于基因A泌辫，cluster X與其余任一cluster的表達均值均相同
• 換句話說：cluster X的基因A的表達水平與其余所有cluster的t檢驗結(jié)果均顯著，才能認為是該cluster的marker gene九默≌鸱牛可以理解為cluster-specific gene。

markers.pbmc2 <- findMarkers(sce.pbmc, pval.type="all")
interesting2 <- markers.pbmc2[[chosen]]
colnames(interesting2)
# [1] "p.value"       "FDR"           "summary.logFC" "logFC.1"       "logFC.2"      
# [6] "logFC.3"       "logFC.4"       "logFC.5"       "logFC.6"       "logFC.8"      
# [11] "logFC.9"       "logFC.10"      "logFC.11"      "logFC.12"      "logFC.13"     
# [16] "logFC.14"      "logFC.15"      "logFC.16

a7_2=as.data.frame(interesting2)
View(a7_2)

• 如下圖結(jié)果驼修，少了top列殿遂，原因很容易理解。
• p.value為該基因的15次t檢驗的p值結(jié)果中的最大值乙各；summary.logFC同樣與之對應(yīng)墨礁；其余列含義可參考2.1
  

這種思路定義的cluster marker gene標準是比較嚴苛的，即當某個cluster的一個基因表達水平與其它所有cluster都具有顯著差異時耳峦，才會被認為是marker gene恩静。

2.3 middle difference(middle)

• 零假設(shè)：對于基因A，cluster X與最多其余一半的cluster表達均值均相同
• 換句話說：cluster X的基因A的表達水平與至少其余一般的cluster的t檢驗結(jié)果都顯著妇萄，才能認為是該cluster的marker gene蜕企。

markers.pbmc3 <- findMarkers(sce.pbmc, pval.type="some")
interesting3 <- markers.pbmc3[[chosen]]
colnames(interesting3)
a7_3=as.data.frame(interesting3)
View(a7_3)

• p.value為該基因的15次t檢驗的p值結(jié)果排名中間的結(jié)果；summary.logFC同樣與之對應(yīng)冠句；其余列含義可參考2.1
  

相較于2.1過于寬松和2.2過于嚴苛的方法轻掩，這種思路則采取了折中的標準進行篩選。

2.4 其它參數(shù)

• direction上調(diào)懦底？下調(diào)唇牧？
  對于marker gene的生物意義來說罕扎，表達上調(diào)的marker gene更具有可解釋性（細胞類型注釋）以及可實驗驗證性。
  可通過設(shè)置direction參數(shù)丐重，只篩選出上調(diào)的marker gene

markers.pbmc.up <- findMarkers(sce.pbmc, direction="up")

但這種篩選方法對于那些由相對下調(diào)的特征基因集特征細胞組成的cluster可能就沒有預(yù)期的結(jié)果腔召。

• lfc差異倍數(shù)
  一般t-test的零假設(shè)為均值相同，也可以設(shè)置參數(shù)lfc設(shè)定差異倍數(shù)為零假設(shè)閾值扮惦，并配合direction參數(shù)可篩選出差異倍數(shù)大于指定閾值的高表達基因臀蛛。

markers.pbmc.up2 <- findMarkers(sce.pbmc, direction="up", lfc=1)

• 批次效應(yīng)block/design
  對于存在batch的測序情況，可以通過block/design指定batch崖蜜，進行分析(每個batch單單獨差異分析浊仆，然后計算合并 combined pvalue)
  以sce.416b為例

#方式一：block參數(shù)
m.out <- findMarkers(sce.416b, block=sce.416b$block, direction="up") 

#方式二：design參數(shù)
design <- model.matrix(~sce.416b$block)
design <- design[,-1,drop=FALSE]

m.alt <- findMarkers(sce.416b, design=design, direction="up")

3、其它檢驗方法

• 之間介紹的pval.type豫领、direction抡柿、lfc等參數(shù)也都適用于下面的方法。

3.1 Wilcoxon rank sum test

• 又稱為 Wilcoxon-Mann-Whitney test等恐、WMW test洲劣，用于評價同一指標在兩組中的分離度(separation)。
• 計算思路簡單來說就是：對于基因A在cluster1與cluster2的表達情況课蔬，任意抽取cluster1的一個細胞的基因A表達量是否高于cluster2的任意一個細胞的基因A表達量囱稽。反映在結(jié)果中，用AUC(area-under-the-curve)指標反應(yīng)购笆。該值越接近1粗悯，則表示cluster1的基因A表達量顯著高于cluster2；越接近0，則表示cluster1的基因A表達量顯著低于cluster2。

markers.pbmc.wmw <- findMarkers(sce.pbmc, test="wilcox", direction="up")
names(markers.pbmc.wmw)
# [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15" "16"

interesting.wmw <- markers.pbmc.wmw[[chosen]]
interesting.wmw[1:10,1:4]

WMW方法區(qū)別于t-test的優(yōu)勢在于不用考慮兩個cluster的細胞數(shù)量差異咐汞，但同時缺點是計算量比較大羹令。

3.2 binomial test

• 這種分析方法重新設(shè)置了表達矩陣，凡是表達量大于0的都視為1。即把count表達矩陣轉(zhuǎn)化為0/1矩陣（0表示不表達、1代表表達）。
• 此時的零假設(shè)就是基因A在兩個cluster的active(有表達)的程度相同撤逢。（至于表達量的高低就不管了）
• foldchange的含義：對于基因A，cluster1的1的比例與cluster2的1的比例的foldchange

markers.pbmc.binom <- findMarkers(sce.pbmc, test="binom", direction="up")
names(markers.pbmc.binom)

interesting.binom <- markers.pbmc.binom[[chosen]]
colnames(interesting.binom)
interesting.binom[1:10,1:4]

這種定義marker gene的標準是十分嚴苛的粮坞。篩選出來的marker基因直接反映了cluster對于該基因的表達有無情況蚊荣。

3.3 關(guān)于傳統(tǒng)Bulk RNA-seq方法

• 對于Bulk RNA-seq的差異分析方法(DESeq2、voom-limma pipeline)莫杈，還是不太適合適用于scRNA-seq的clust間的兩兩比較互例。因為一般將cluster的一個細胞視為1 replication的基礎(chǔ)上，每個cluster往往會有數(shù)百次重復(fù)筝闹，不符合DESeq2媳叨、voom-limma pipeline等適合有限重復(fù)的情況腥光；另一方面DESeq2、voom-limma pipeline等假設(shè)每次重復(fù)的基因表達情況相似糊秆，而scRNA-seq測序過程中的每個細胞的表達水平分布可能受各種因素不相一致武福，即使在同一個cluster里。
• 即使不做推薦痘番，教程里還是演示了采用voom-limma pipeline循環(huán)比較所有cluster間兩兩差異分析的代碼捉片，這里就不展示了。
• 換一個角度思考夫偶，findMarkers()函數(shù)提供的分析方法可以適用于涉及多組涉及的Bulk RNA-seq差異分析方法界睁。

關(guān)于replication這個問題，在上述所有方法中是把每個cluster的cell視為一次重復(fù)兵拢，但這還是一個值得細思的問題。即使這些細胞來自同一個取樣組織逾礁，但本質(zhì)上還是不符合生物學(xué)重復(fù)的概念说铃。換句話說，一次測序的結(jié)果應(yīng)當視為一次重復(fù)嘹履，對于Bulk RNA-seq比較容易理解腻扇，對于scRNA-seq往往會忽略這個因素。