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