劉小澤寫于19.8.19
主要介紹官方教程PBMC 3K的下游分析（https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html）
上篇：單細胞Seurat包升級之2,700 PBMCs分析（上）
數(shù)據(jù)標準化（scale）
它的目的是實現(xiàn)數(shù)據(jù)的線性轉(zhuǎn)換(scaling)堕绩，這也是降維處理(如PCA)之前的標準預(yù)處理腔剂。
主要利用了ScaleData函數(shù)痴脾，回憶一下，之前歸一化（normalize）這里做的操作是log處理辆毡，它是對所有基因的表達量統(tǒng)一對待的，最后放在了一個起跑線上甜害。但是為了真正找到差異舶掖，我們還要基于這個起跑線，考慮不同樣本對表達量的影響
它做了：
- 將每個基因在所有細胞的表達量進行調(diào)整尔店，使得每個基因在所有細胞的表達量均值為0
- 將每個基因的表達量進行標準化眨攘，讓每個基因在所有細胞中的表達量方差為1
# 先使用全部基因
all.genes <- rownames(pbmc)
> length(all.genes)
[1] 13714
pbmc <- ScaleData(pbmc, features = all.genes)
# 結(jié)果存儲在pbmc[["RNA"]]@scale.data中
感覺這一步太慢，能不能加速一點嚣州？
這一步主要目的就是下一步的降維鲫售，使用全部基因是比較慢。如果想提高速度该肴，可以只對鑒定的HVGs（之前FindVariableFeatures設(shè)置的是2000個）進行操作情竹，其實這個函數(shù)默認就是對部分差異基因進行操作，：pbmc <- ScaleData(pbmc) 匀哄。這樣操作的結(jié)果中降維和聚類不會被影響秦效，但是只對HVGs基因進行標準化，下面畫的熱圖依然會受到全部基因中某些極值的影響涎嚼。所以為了熱圖的結(jié)果阱州，還是對所有基因進行標準化比較好。
ScaleData的操作過程是怎樣的?
看一下幫助文檔就能了解：
ScaleData這個函數(shù)有兩個默認參數(shù)：do.scale = TRUE和do.center = TRUE铸抑，然后需要輸入進行scale/center的基因名（默認是取前面FindVariableFeatures的結(jié)果）贡耽。
scale和center這兩個默認參數(shù)應(yīng)該不陌生了衷模，center就是對每個基因在不同細胞的表達量都減去均值鹊汛；
scale就是對每個進行center后的值再除以標準差（就是進行了一個z-score的操作）
再看一個來自BioStar的討論：https://www.biostars.org/p/349824/
問題1：為什么在NormalizeData()之后，還需要進行ScaleData操作阱冶？
前面NormalizeData進行的歸一化主要考慮了測序深度/文庫大小的影響（reads *10000 divide by total reads, and then log）刁憋，但實際上歸一化的結(jié)果還是存在基因表達量分布區(qū)間不同的問題；而ScaleData做的就是在歸一化的基礎(chǔ)上再添zero-centres （mean/sd =》 z-score）木蹬，將各個表達量放在了同一個范圍中至耻，真正實現(xiàn)了”表達量的同一起跑線“。
另外镊叁，新版的ScaleData函數(shù)還支持設(shè)定回歸參數(shù)尘颓，如（nUMI、percent.mito）晦譬，來校正一些技術(shù)噪音
問題2： FindVariableGenes() or RunPCA() or FindCluster()這些參數(shù)是基于歸一化Normalized_Data還是標準化Scaled_Data 疤苹？
所有操作都是基于標準化的數(shù)據(jù)Scaled_Data ，因為這個數(shù)據(jù)是針對基因間比較的敛腌。例如有兩組基因表達量如下:
g1 10 20 30 40 50
g2 20 40 60 80 100
雖然它們看起來是g2的表達量是g1的兩倍卧土，但真要降維聚類時惫皱，就要看scale的結(jié)果：
> scale(c(10,20,30,40,50))
[,1]
[1,] -1.2649111
[2,] -0.6324555
[3,] 0.0000000
[4,] 0.6324555
[5,] 1.2649111
attr(,"scaled:center")
[1] 30
attr(,"scaled:scale")
[1] 15.81139
> scale(c(20,40,60,80,100))
[,1]
[1,] -1.2649111
[2,] -0.6324555
[3,] 0.0000000
[4,] 0.6324555
[5,] 1.2649111
attr(,"scaled:center")
[1] 60
attr(,"scaled:scale")
[1] 31.62278
二者標準化后的分布形態(tài)是一樣的（只是中位數(shù)不同），而我們后來聚類也是看數(shù)據(jù)的分布尤莺，分布相似聚在一起旅敷。所以歸一化差兩倍的兩個基因，根據(jù)標準化的結(jié)果依然聚在了一起
問題3： ScaleData()在scRNA分析中一定要用嗎颤霎？
實際上標準化這個過程不是單細胞分析固有的媳谁，在很多機器學習、降維算法中比較不同feature的距離時經(jīng)常使用捷绑。當不進行標準化時韩脑，有時feature之間巨大的差異（就像??問題2中差兩倍的基因表達量，實際上可能差的倍數(shù)會更多）會影響分析結(jié)果粹污，讓我們誤以為它們的數(shù)據(jù)分布相差很遠段多。
很多時候，我們會混淆normalization和scale：那么看一句英文解釋：
Normalization "normalizes" within the cell for the difference in sequenicng depth / mRNA throughput. 主要著眼于樣本的文庫大小差異
Scaling "normalizes" across the sample for differences in range of variation of expression of genes . 主要著眼于基因的表達分布差異
另外壮吩，看scale()函數(shù)的幫助文檔就會發(fā)現(xiàn)进苍，它實際上是對列進行操作：
那么具體操作中是要先對表達矩陣轉(zhuǎn)置，然后scale鸭叙，最后轉(zhuǎn)回來
t(scale(t(dat[genes,])))
問題4：RunCCA()函數(shù)需要歸一化數(shù)據(jù)還是標準化數(shù)據(jù)觉啊？
通過看這個函數(shù)的幫助文檔：
RunCCA(object, object2, group1, group2, group.by, num.cc = 20, genes.use,
scale.data = TRUE, rescale.groups = FALSE, ...)
顯示scale.data = TRUE，那么就需要標準化后的數(shù)據(jù)
怎樣移除不想要的差異來源沈贝？
在進行降維聚類時杠人，主要就是考慮數(shù)據(jù)差異對分布產(chǎn)生的影響。由于技術(shù)誤差導致的差異是我們不感興趣的宋下，排除這一部分其實也在上面的問題1中提到了嗡善，如果說不想線粒體污染導致的差異影響整體差異分布：
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
這個函數(shù)僅針對初級用戶，如果想深入探索技術(shù)誤差的排除学歧，官方給出了另一個專業(yè)版函數(shù)：SCTransform罩引，它也包含vars.to.regress這個選項。詳細的說明在：https://satijalab.org/seurat/v3.0/sctransform_vignette.html
進行線性降維
最常用的線性降維就是PCA方法枝笨，使用標準化的數(shù)據(jù)
它默認選擇之前鑒定的差異基因（2000個）作為input袁铐，但可以使用features進行設(shè)置；默認分析的主成分數(shù)量是50個
# 這里就使用默認值
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
# 結(jié)果保存在reductions 這個接口中
檢查下PCA結(jié)果：
> print(pbmc[["pca"]], dims = 1:3, nfeatures = 5)
# 或者用 print(pbmc@reductions$pca, dims = 1:3, nfeatures = 5)
PC_ 1
Positive: CST3, TYROBP, LST1, AIF1, FTL
Negative: MALAT1, LTB, IL32, IL7R, CD2
PC_ 2
Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
Negative: NKG7, PRF1, CST7, GZMB, GZMA
PC_ 3
Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
Negative: PPBP, PF4, SDPR, SPARC, GNG11
用VizDimLoadings對print的結(jié)果可視化：
VizDimLoadings(pbmc, dims = 1:3, reduction = "pca")
用DimPlot進行降維后的可視化
提取兩個主成分（默認前兩個，當然可以修改dims選項）繪制散點圖
DimPlot(pbmc, reduction = "pca")
這個圖中每個點就是一個樣本，它根據(jù)PCA的結(jié)果坐標進行畫圖屠缭，這個坐標保存在了cell.embeddings中：
> head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
PC_1 PC_2
AAACATACAACCAC -4.7296855 -0.5184265
AAACATTGAGCTAC -0.5174029 4.5918957
其實還支持定義分組：根據(jù)pbmc@meta.data中的ident分組：
> table(pbmc@meta.data$orig.ident)
pbmc3k
2638
# 當然這里只有一個分組
探索異質(zhì)性來源—DimHeatmap
每個細胞和基因都根據(jù)PCA結(jié)果得分進行了排序咪惠，默認畫前30個基因（nfeatures設(shè)置），1個主成分(dims設(shè)置)，細胞數(shù)沒有默認值（cells設(shè)置）
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
# 其中balanced表示正負得分的基因各占一半
降維的真實目的是盡可能減少具有相關(guān)性的變量數(shù)目，例如原來有700個樣本叫惊，也就是700個維度/變量帝洪，但實際上根據(jù)樣本中的基因表達量來歸類似舵，它們或多或少都會有一些關(guān)聯(lián)，這樣有關(guān)聯(lián)的一些樣本就會聚成一小撮葱峡，表示它們的”特征距離“相近砚哗。最后直接分析這些”小撮“，就用最少的變量代表了實際最多的特征
既然每個主成分都表示具有相關(guān)性的一小撮變量（稱之為”metafeature“砰奕，元特征集）蛛芥，那么問題來了：
降維后怎么選擇合適數(shù)量的主成分來代表整個數(shù)據(jù)集？
選10個军援？20個仅淑？50個？最簡單的方法就是：
ElbowPlot(pbmc)
它是根據(jù)每個主成分對總體變異水平的貢獻百分比排序得到的圖胸哥，我們主要關(guān)注”肘部“的PC涯竟，它是一個轉(zhuǎn)折點（也即是這里的PC9-10），說明取前10個主成分可以比較好地代表總體變化
但官方依然建議我們空厌，下游分析時多用幾個成分試試（10個庐船、15個甚至50個），如果起初選擇的合適嘲更，結(jié)果不會有太大變化筐钟；另外推薦開始不確定時要多選一些主成分，也不要直接就定下5個成分進行后續(xù)分析赋朦，那樣很有可能會出錯
細胞聚類
Seurat 3版本提供了graph-based 聚類算法篓冲，參考兩篇文章：SNN-Cliq, Xu and Su, Bioinformatics, 2015] 和 PhenoGraph, Levine et al., Cell, 2015]. 簡言之，這些graph的方法都是將細胞嵌入一個算法圖層中宠哄，例如：K-nearest neighbor (KNN) graph 中壹将，每個細胞之間的連線就表示相似的基因表達模式，然后按照這個相似性將圖分隔成一個個的高度相關(guān)的小群體（名叫‘quasi-cliques’ or ‘communities’）琳拨；
又比如瞭恰，在PhenoGraph中屯曹，首先根據(jù)PCA的歐氏距離結(jié)果狱庇，構(gòu)建了KNN圖，找到了各個近鄰組成的小社區(qū)恶耽；然后根據(jù)近鄰中兩兩住戶（細胞）之間的交往次數(shù)（shared overlap）計算每條線的權(quán)重（術(shù)語叫Jaccard similarity） 密任。這些計算都是用包裝好的函數(shù)FindNeighbors() 得到的，它的輸入就是前面降維最終確定的主成分
# 因為我們前面挑選了10個PCs偷俭，所以這里dims定義為10個
pbmc <- FindNeighbors(pbmc, dims = 1:10)
以上是憑個人理解發(fā)揮浪讳，大概就是這么個意思，不涉及算法層面推導
得到權(quán)重后涌萤，為了對細胞進行聚類淹遵，使用了計算模塊的算法（默認使用Louvain口猜，另外還有包括SLM的其他三種），使用FindClusters() 進行聚類透揣。其中包含了一個resolution的選項济炎，它會設(shè)置一個”間隔“值，該值越大辐真，間隔越大须尚，得到的cluster數(shù)量越多
怎么理解這個
resolution參數(shù)？
可以想象成配眼鏡的時候侍咱，鏡片度數(shù)越高耐床，分辨率越高，看的越清楚楔脯，看到的細節(jié)越豐富（cluster越多）撩轰；反之，如果分辨率調(diào)的很低昧廷，結(jié)果就看的模模糊糊一大坨
一般來說钧敞，這個值在細胞數(shù)量為3000左右時設(shè)為0.4-1.2 會有比較好的結(jié)果
pbmc <- FindClusters(pbmc, resolution = 0.5)
# 結(jié)果聚成幾類，用Idents查看
> head(Idents(pbmc), 5)
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG
1 3 1 2
AAACCGTGTATGCG
6
Levels: 0 1 2 3 4 5 6 7 8
# 發(fā)現(xiàn)3000細胞麸粮，設(shè)resolution=0.5時溉苛，得到Number of communities: 9
另一種降維方法—非線性降維（UMAP/tSNE）
線性降維PCA（1933）一個特點就是：它將高維中不相似的數(shù)據(jù)點放在較低維度時，會盡可能多地保留原始數(shù)據(jù)的差異弄诲，因此導致這些不相似的數(shù)據(jù)相距甚遠愚战，大多的數(shù)據(jù)重疊放置
tSNE（2008）兼顧了高維降維+低維展示，降維后的那些不相似的數(shù)據(jù)點依然可以放得靠近一些齐遵，保證了點既在局部分散寂玲，又在全局聚集
后來開發(fā)的UMAP（2018）相比tSNE又能展示更多的維度(參考文獻：https://sci-hub.tw/https://doi.org/10.1038/nbt.4314)
做個對比：
# 使用UMAP
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap") # 還可以設(shè)置label = TRUE讓數(shù)字顯示在每個cluster上
# 對計算過程很費時的結(jié)果保存一下
save(pbmc, file = 'pbmc_umap.Rdata')
需要注意的點：
注意1：tSNE算法具有隨機性，為了保證重復(fù)性梗摇，要固定一個隨機種子
注意2：如何在R中使用UMAP拓哟？
以下測試是在個人電腦上，需要管理員權(quán)限
因為UMAP是基于Python的伶授，如果要用它断序，必須先保證有一個python的安裝環(huán)境，先試驗一下：
> reticulate::py_install(packages ='umap-learn')
# 結(jié)果報錯糜烹，說讓我升級一下virtualenv(mac和linux默認使用virtualenv违诗，而Windows只能使用conda)
Error: Prerequisites for installing Python packages not available.
Execute the following at a terminal to install the prerequisites:
$ sudo /usr/local/bin/pip install --upgrade virtualenv
然后直接在Rstudio中打開Terminal，輸入：sudo /usr/local/bin/pip install --upgrade virtualenv
再運行一遍上面的R代碼（成功與否取決于個人的網(wǎng)絡(luò)疮蹦，因為它要通過外鏈去下載诸迟，wlan有限制的話就用個人熱點吧）：
一般總共二三十兆的文件，需要下載個十幾分鐘，速度確實很慢阵苇，而這個過程很有可能出現(xiàn)斷線壁公，只能重新運行安裝。已下載完的文件會直接識別不會重復(fù)下載绅项。檢查是否下載成功標志就是：
library(umap)能否成功贮尖，這個是Seurat調(diào)用的基礎(chǔ)。最后特別注意趁怔，要重啟一下Rstudio再運行RunUMAP
找差異表達基因（cluster biomarker）
目的是根據(jù)表達量不同這個特征而分出不同細胞群的基因們（就是找具有生物學意義的HVGs=》即biomarkers）
marker可以理解成標記湿硝，就是一看到有這些基因，就說明是這個細胞群體
主要利用FindAllMarkers()函數(shù)润努，它默認會對單獨的一個細胞群與其他幾群細胞進行比較找到正关斜、負表達biomarker（這里的正負有點上調(diào)、下調(diào)基因的意思铺浇；正marker表示在我這個cluster中表達量高痢畜，而在其他的cluster中低）
# 找到cluster1中的marker基因
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
上面這個代碼中min.pct的意思是：一個基因在任何兩群細胞中的占比最低不能低于多少，當然這個可以設(shè)為0鳍侣，但會大大加大計算量
另外丁稀，如果要比較cluster5和cluster0、cluster3的marker基因：
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
一步到位的辦法FindAllMarkers（對所有cluster都比較一下倚聚，并只挑出來正表達marker）：
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# 這一步過濾好好理解(進行了分類线衫、排序、挑前2個)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
Seurat提供了多種差異分析的檢驗方法惑折，在test.use選項中設(shè)置（詳見：DE vignette ）例如：
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
多種可視化方法：
-
VlnPlot：用小提琴圖對某些基因在全部cluster中表達量進行繪制VlnPlot(pbmc, features = c("MS4A1", "CD79A")) -
FeaturePlot：最常用的可視化=》將基因表達量投射到降維聚類結(jié)果中FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
另外還有：RidgePlot, CellScatter, and DotPlot
-
DoHeatmap對給定細胞和基因繪制熱圖授账，例如對每個cluster繪制前20個markertop10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC) DoHeatmap(pbmc, features = top10$gene) + NoLegend()
賦予每個cluster細胞類型
重點在于背景知識的理解，能夠?qū)arker與細胞類型對應(yīng)起來惨驶，例如這里從各個cluster中找到的marker對應(yīng)到了不同的細胞（這一過程是比較考驗課題理解程度的）：
有了這個白热，繪圖就很簡單：
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球粗卜，想讓它成為一個不拽術(shù)語屋确、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
