寫(xiě)在前面
TBtools已經(jīng)拆除了絕大部分生信分析的門(mén)檻，讓眾多濕實(shí)驗(yàn)工作者們（無(wú)代碼基礎(chǔ)）可以暢快地分析他們的數(shù)據(jù)。
不求人非常爽！
不用跟公司低效溝通分析細(xì)節(jié)非常爽！
很早之前陳師兄就跟我提過(guò)是否有興趣整一個(gè)單細(xì)胞分析方面的插件，一方面由于我R用的非常不6，另一方面確實(shí)是懶。。。所以擱置了。。 近期又提了一下，于是菜菜的我終于踏上了一段辛酸的寫(xiě)bug之旅。。
總的來(lái)說(shuō)，第一次接觸shiny開(kāi)發(fā)體系，確實(shí)非常不習(xí)慣，跟之前數(shù)據(jù)庫(kù)整的html+css+JavaScript+php完全不一樣，特別是在debug上，菜狗如我，太難了。。。 不過(guò)shiny確實(shí)在實(shí)現(xiàn)可交互功能的時(shí)候非常強(qiáng)大，能夠在ui端可視化地實(shí)時(shí)控制參數(shù)進(jìn)行分析或者自定義繪圖。
磕磕絆絆，盡管質(zhì)量有待提高，最終還是完成了這個(gè)插件，Seurat Plugin for TBtools。將Seurat的基本分析流程進(jìn)行了打包與可視化，能夠在ui界面上直接完成從讀入標(biāo)準(zhǔn)表達(dá)矩陣到數(shù)據(jù)質(zhì)控、降維，然后進(jìn)行細(xì)胞分群、marker基因鑒定與可視化等基本的single cell分析流程。
Seurat Plugin for TBtools
插件安裝方式
注意到這是一個(gè) TBtools R-plugin。與其他R插件類(lèi)似，如果沒(méi)有安裝 Rserver插件，那么需要先安裝 Rserver 插件（可以直接在 TBtools Plugin Store 找到并安裝）。
安裝完成后，打開(kāi)插件。進(jìn)入方式非常簡(jiǎn)單，你只需要點(diǎn)一下 Start 。
1. 上傳數(shù)據(jù)
輸入表達(dá)量數(shù)據(jù)，目前支持兩種格式的表達(dá)量文件
- 基于Cell Ranger輸出的稀疏矩陣，包括三個(gè)文件：barcodes.tsv，genes.tsv，matrix.mtx。
- count矩陣。一般在GEO數(shù)據(jù)庫(kù)或者其他公共數(shù)據(jù)庫(kù)存儲(chǔ)平臺(tái)中（比如浙大樊龍江老師團(tuán)隊(duì)開(kāi)發(fā)的PlantscRNAdb）下載單細(xì)胞公共數(shù)據(jù)，大部分都是count矩陣，就是這個(gè)文件可能會(huì)比較大。。。同樣，我們選擇count data選項(xiàng)，然后上傳相應(yīng)的文件，注意需要是csv格式，即逗號(hào)分隔符，如果是tab分割的，可以TBtools或者notepad++做一下轉(zhuǎn)換。
- 我們也準(zhǔn)備了一個(gè)demo data，使用的是PRJNA497883這套公共數(shù)據(jù)（擬南芥根系）。
2. 數(shù)據(jù)預(yù)處理以及創(chuàng)建Seurat對(duì)象
上傳完數(shù)據(jù)之后，在右邊的“Data Preview”中會(huì)展示count矩陣，用戶(hù)可以自行瀏覽或者檢索。由于單細(xì)胞數(shù)據(jù)量（細(xì)胞數(shù)目）一般非常大，因此我們只展示了前30個(gè)細(xì)胞的數(shù)據(jù)。
接下來(lái)對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，過(guò)濾掉一些不表達(dá)的基因以及基本沒(méi)有表達(dá)基因的細(xì)胞，一定程度可以去除一些異常值，加速下游的計(jì)算。通常這一步我們過(guò)濾兩個(gè)指標(biāo)：
min.cells，保留至少在幾個(gè)細(xì)胞中檢測(cè)到表達(dá)的基因，默認(rèn)為3個(gè)細(xì)胞，可以根據(jù)自己樣本測(cè)序的實(shí)際細(xì)胞數(shù)量來(lái)進(jìn)行調(diào)整，目的是為了去除一些幾乎沒(méi)有表達(dá)的基因。
min.features，保留至少檢測(cè)到有多少個(gè)基因表達(dá)的細(xì)胞，默認(rèn)為200個(gè)基因。如果表達(dá)細(xì)胞表達(dá)的基因過(guò)少，那么很有可能是一個(gè)異常細(xì)胞，將其過(guò)濾。
預(yù)處理完數(shù)據(jù)之后，后臺(tái)自動(dòng)對(duì)其構(gòu)建Seurat對(duì)象。
質(zhì)控
質(zhì)量控制。這一步我們對(duì)數(shù)據(jù)的一些QC指標(biāo)進(jìn)行可視化，用戶(hù)可以根據(jù)可視化圖片判斷自己的數(shù)據(jù)是否正常。此外，我們提供了一個(gè)接口，用戶(hù)可以在QC step1中計(jì)算一些基因在所有細(xì)胞中表達(dá)的比例，例如計(jì)算線(xiàn)粒體基因表達(dá)比例，其是過(guò)濾異常細(xì)胞的一個(gè)常用指標(biāo)。一般我們認(rèn)為線(xiàn)粒體基因表達(dá)過(guò)高的細(xì)胞，是一個(gè)處于應(yīng)激或者其他不好狀態(tài)的細(xì)胞，需要將其過(guò)濾。也可以計(jì)算核糖體基因等的表達(dá)情況。
Step1：計(jì)算線(xiàn)粒體（或者其他基因）的表達(dá)比例
這里我們以線(xiàn)粒體基因?yàn)槔M南芥基因組中（TAIR10）的線(xiàn)粒體基因?yàn)锳TMG開(kāi)頭，因此我們可以用正則表達(dá)式“^ATMG”來(lái)匹配擬南芥中所有線(xiàn)粒體基因。如果是人類(lèi)，線(xiàn)粒體基因則以MT-開(kāi)頭，小鼠以mt-開(kāi)頭。我們提供了兩種基因檢索方式：
正則匹配，匹配線(xiàn)粒體基因等具有特定模式的基因集，例如擬南芥中：^ATMG
用戶(hù)自定義基因集，用戶(hù)輸入自定義的基因集進(jìn)行計(jì)算，每行一個(gè)基因。
Tips：如果你研究的物種并不是擬南芥等模式物種，我們?cè)撛趺创_定它的線(xiàn)粒體基因呢？
- 到Ensembl數(shù)據(jù)庫(kù)或者其他的物種數(shù)據(jù)庫(kù)中找到該物種的注釋文件，在Ensembl中通常會(huì)有一個(gè)線(xiàn)粒體基因注釋文件，下載到本地，打開(kāi)你就可以發(fā)現(xiàn)這個(gè)物種的線(xiàn)粒體基因ID是否有規(guī)律，如果有固定的開(kāi)頭，即可以用正則表達(dá)式去匹配，如^ATMG，如果沒(méi)有，那么就TBtools提取所有的線(xiàn)粒體基因ID，然后使用mode2，輸入線(xiàn)粒體基因ID集，進(jìn)行計(jì)算。
Step2：QC指標(biāo)的可視化
根據(jù)QC指標(biāo)過(guò)濾異常細(xì)胞
在Step2中，我們對(duì)幾個(gè)QC指標(biāo)進(jìn)行了可視化，其中包括Step1中用戶(hù)指定計(jì)算的基因表達(dá)量（如線(xiàn)粒體）?？梢暬闹笜?biāo)如下：
每個(gè)細(xì)胞的線(xiàn)粒體基因比例（如果用戶(hù)有在Step1中計(jì)算）
每個(gè)細(xì)胞檢測(cè)到的基因數(shù)目（nFeature_RNA）
每個(gè)細(xì)胞測(cè)得的UMI總數(shù)（nCount_RNA）
接著根據(jù)線(xiàn)粒體基因比例和每個(gè)細(xì)胞檢測(cè)到的基因數(shù)目進(jìn)行細(xì)胞過(guò)濾。線(xiàn)粒體基因表達(dá)過(guò)高的細(xì)胞一般是異常細(xì)胞，而檢測(cè)到基因數(shù)目異常過(guò)高的細(xì)胞，大概率是雙胞（doublets，一個(gè)文庫(kù)中存在兩個(gè)細(xì)胞的情況，會(huì)對(duì)后續(xù)的分群和細(xì)胞鑒定造成影響和誤導(dǎo)，需要剔除。當(dāng)然有專(zhuān)門(mén)的doublets預(yù)測(cè)軟件，如DoubletFinder等）
注意：示例數(shù)據(jù)使用的是擬南芥根尖的公共數(shù)據(jù)，該套數(shù)據(jù)已經(jīng)經(jīng)過(guò)質(zhì)控，因此QC指標(biāo)基本正常，并不存在異常值。
查看單細(xì)胞數(shù)據(jù)的質(zhì)量
我們也對(duì)特征值進(jìn)行了相關(guān)性分析，并可視化。用戶(hù)可以根據(jù)nFeature_RNA和nCount_RNA的相關(guān)性進(jìn)行大致判斷自己的數(shù)據(jù)是否正常，理論上一個(gè)細(xì)胞測(cè)得的UMI越多，其測(cè)得的基因應(yīng)該也越多，因此這兩個(gè)指標(biāo)應(yīng)該是強(qiáng)正相關(guān)。
數(shù)據(jù)標(biāo)準(zhǔn)化與歸一化
在對(duì)數(shù)據(jù)進(jìn)行降維之前，需要對(duì)數(shù)據(jù)進(jìn)行處理，即三步法：Normalization，F(xiàn)ind Variable Genes以及Scale。通常情況下，默認(rèn)只是標(biāo)準(zhǔn)化2000個(gè)高變基因，運(yùn)算速度更快，不影響 PCA降維和細(xì)胞分群。對(duì)數(shù)據(jù)進(jìn)行scale，可以消除極高表達(dá)的基因的影響。在Seurat plugin中，用戶(hù)可以直接點(diǎn)擊start按鈕進(jìn)行一鍵化跑三步法。
PCA線(xiàn)性降維
對(duì)處理過(guò)后的數(shù)據(jù)進(jìn)行pca降維，默認(rèn)計(jì)算前50個(gè)主成分。用戶(hù)點(diǎn)擊PCA降維按鈕之后，會(huì)實(shí)時(shí)可視化展示PCA結(jié)果。
細(xì)胞聚類(lèi)
在進(jìn)行細(xì)胞聚類(lèi)之前，我們首先要確定好兩個(gè)參數(shù)：1，使用的PC個(gè)數(shù)；2，Resolution（分辨率）。
對(duì)于PC數(shù)目來(lái)說(shuō)，我們需要使用合適的PC數(shù)目來(lái)進(jìn)行聚類(lèi)，太多的PC可能會(huì)對(duì)聚類(lèi)造成一定的干擾，太少的PC又不能很好的解釋數(shù)據(jù)集。因此我們推薦使用Elbow Plot來(lái)確定使用的PC數(shù)目。在Elbow Plot選擇處于拐點(diǎn)處，即趨近于水平處的PC個(gè)數(shù)，來(lái)進(jìn)行聚類(lèi)，這些PC能夠很大程度上解釋整個(gè)數(shù)據(jù)集。
對(duì)于Resolution。Resolution參數(shù)決定下游聚類(lèi)所得到的分群數(shù)，對(duì)于3千左右的細(xì)胞量，seurat官方說(shuō)0.4-1.2能得到較好的結(jié)果。如果數(shù)據(jù)量增大，該參數(shù)也應(yīng)該適當(dāng)增大。 那么，到底選擇多少的Resolution最為合適呢？ 這里我們推薦使用clustree來(lái)確定Resolution。clustree能夠進(jìn)行多次不同的Resolution分群，每個(gè)不同的Resolution的分群結(jié)果進(jìn)行可視化，用戶(hù)可以根據(jù)Resolution的變化而得到的不同細(xì)胞群數(shù)，來(lái)確定適合你數(shù)據(jù)的最佳Resolution。
個(gè)人項(xiàng)目經(jīng)驗(yàn)，如果Resolution選取較小，則分群較少，會(huì)將一些相似的細(xì)胞類(lèi)群合并成一個(gè)大的cluster，從而對(duì)后續(xù)的細(xì)胞類(lèi)型鑒定造成困擾，丟失一些細(xì)胞類(lèi)型信息。而Resolution過(guò)大，則會(huì)將cluster分的過(guò)多，加大后續(xù)細(xì)胞類(lèi)型鑒定的難度。因此我們一般會(huì)將Resolution選擇在中等稍微偏大一點(diǎn)的值，鑒定細(xì)胞類(lèi)型的時(shí)候，可以將具有相同marker基因，且在umap/tSNE上距離相近的cluster歸為同一個(gè)細(xì)胞類(lèi)型。
確定好PC數(shù)目和Resolution之后，則可以點(diǎn)擊按鈕進(jìn)行細(xì)胞聚類(lèi)。
UMAP/tSNE非線(xiàn)性降維
細(xì)胞聚類(lèi)完成之后，用戶(hù)可以點(diǎn)擊按鈕進(jìn)行UMAP/tSNE非線(xiàn)性降維，并對(duì)其進(jìn)行可視化。
marker基因鑒定/差異表達(dá)分析
marker基因鑒定，本質(zhì)上其實(shí)是進(jìn)行差異分析，鑒定每個(gè)cluster的差異表達(dá)基因（差異上調(diào)，往往cluster的marker基因是只在該cluster中特異表達(dá)，而在其他cluster中不表達(dá)。因此鑒定marker的時(shí)候，我們通常會(huì)設(shè)置only.pos=TRUE，即只保留差異上調(diào)基因）。在該模塊中，提供了三種鑒定marker基因的模式。
一鍵化鑒定所有cluster的marker基因（比較慢，運(yùn)行時(shí)間與細(xì)胞數(shù)和cluster數(shù)成正比），即對(duì)每一個(gè)cluster都與剩余其他cluster進(jìn)行差異分析；
用戶(hù)指定鑒定某個(gè)cluster的基因（快速），即對(duì)選定的cluster與其他所有cluster進(jìn)行差異分析；
用戶(hù)指定鑒定某兩個(gè)cluster之間的marker基因，實(shí)際上就是在對(duì)這兩個(gè)選定cluster之間進(jìn)行差異分析，選用這種模式我們就沒(méi)有必要只保留差異上調(diào)基因。
marker基因可視化
鑒定好marker基因之后，用戶(hù)可以在該模塊中，對(duì)感興趣的marker基因進(jìn)行各種可視化，分析其是否確實(shí)在某些cluster中特異表達(dá)，而在其他cluster中不表達(dá)。
- 對(duì)每個(gè)cluster中的Top5 marker基因進(jìn)行熱圖展示，用戶(hù)可以選定需要展示的marker基因數(shù)目
- 對(duì)每個(gè)cluster中的Top5 marker基因進(jìn)行Dot plot展示，用戶(hù)可以選定需要展示的marker基因數(shù)目
- 用戶(hù)輸入基因或基因集，對(duì)其進(jìn)行Violin plot可視化
- 用戶(hù)輸入基因或基因集，對(duì)其進(jìn)行Feature Plot可視化，將marker基因的表達(dá)量映射到UMAP或者tSNE上，可以非常清晰地看到marker基因在cluster中的表達(dá)分布情況
- Ridge plot。
寫(xiě)在最后
磕磕絆絆，最終還是寫(xiě)完了這個(gè)插件。雖然寫(xiě)的確實(shí)很爛，功能也非常簡(jiǎn)單，但在寫(xiě)這些代碼的折騰中，確確實(shí)實(shí)能感受到收獲了新的東西，接觸到了新的網(wǎng)頁(yè)工具開(kāi)發(fā)方式（shiny跟html和php等網(wǎng)站開(kāi)發(fā)確實(shí)不一樣。。。） 。Seurat Plugin for TBtools，確實(shí)還有許多需要改進(jìn)的地方以及新增的功能，后續(xù)也會(huì)抽空慢慢將其完善，副教授說(shuō)“優(yōu)秀的產(chǎn)品都是一步一步優(yōu)化的”哈哈，希望自己也能像產(chǎn)品一樣，”一步一步優(yōu)化“！
接觸了shiny才發(fā)現(xiàn)，用其與R腳本結(jié)合，搭建實(shí)時(shí)可交互的網(wǎng)頁(yè)/插件工具確實(shí)非常強(qiáng)大，后續(xù)也會(huì)慢慢摸索，嘗試將單細(xì)胞分析中常用的幾個(gè)包寫(xiě)成插件，順便繼續(xù)磨煉一下自己的能力。。 我果然還是太菜了。。
博士生涯轉(zhuǎn)眼已過(guò)一半，臨近過(guò)年，希望各位老板舒舒服服過(guò)個(gè)好年，來(lái)年再戰(zhàn)！
