seurat 涉及的數(shù)據(jù)分析包括很多步驟。之前只顧著干活兒钾唬，也沒有系統(tǒng)的整理過分析中的具體內(nèi)容。這里就參照網(wǎng)上大神們分享的帖子侠驯，來梳理一下抡秆。
一、讀入數(shù)據(jù)吟策。
Read10X() 和 CreateSeuratObject()函數(shù)就是根據(jù)輸入矩陣/數(shù)據(jù)框儒士，創(chuàng)建Seurat對(duì)象的。重要步驟是 設(shè)置 ident 和添加 meta.data檩坚。
dat <- Read10X(data.dir = "your/work/path")
organoids <- CreateSeuratObject(counts = dat,
project = "Organoids",
min.cells = 3,
min.features = 200)
- min.cells 表示一個(gè)基因至少要在3個(gè)細(xì)胞中被檢測(cè)到着撩，否則不要诅福。
- min.features 參數(shù)指定每個(gè)細(xì)胞需要檢測(cè)的最小基因數(shù)量。此參數(shù)將過濾掉質(zhì)量較差的細(xì)胞睹酌，這些細(xì)胞可能只是封裝了隨機(jī)barcodes权谁，而沒有任何真實(shí)的細(xì)胞。通常憋沿，檢測(cè)到的基因少于100或者200個(gè)的細(xì)胞不會(huì)被考慮進(jìn)行分析旺芽。
這里還是設(shè)計(jì)一個(gè)知識(shí)點(diǎn)就是R里面的S3類和S4類。list一般情況下被認(rèn)為是S3類辐啄，S4類是指使用slots存儲(chǔ)數(shù)據(jù)的格式采章。（如果說的不對(duì)歡迎中糾錯(cuò)）。這里讀進(jìn)去的數(shù)值是三個(gè)文本文件創(chuàng)建的稀疏矩陣壶辜。
什么是稀疏矩陣悯舟？
在[矩陣]中，若數(shù)值為0的元素?cái)?shù)目遠(yuǎn)遠(yuǎn)多于非0元素的數(shù)目砸民，并且非0元素分布沒有規(guī)律時(shí)抵怎，則稱該矩陣為稀疏矩陣；與之相反岭参，若非0元素?cái)?shù)目占大多數(shù)時(shí)反惕，則稱該矩陣為稠密矩陣。
如果自己手上有單細(xì)胞數(shù)據(jù)演侯，那個(gè)matrix文件里面包含很多0姿染。
二、NormalizeData
因?yàn)樵跍y(cè)序之前會(huì)對(duì)抓取到的RNA進(jìn)行PCR擴(kuò)增秒际，所以需要考慮文庫深度的對(duì)測(cè)序的影響悬赏，所以需要對(duì)上一步得到的稀疏矩陣進(jìn)行Normalize。
Normalize的方式：每個(gè)細(xì)胞每個(gè)基因的特征計(jì)數(shù)除以該細(xì)胞的特征總計(jì)數(shù)娄徊，再乘以scale.factor（默認(rèn)10000）闽颇，然后使用log1p進(jìn)行對(duì)數(shù)轉(zhuǎn)換。（log1p=log(n+1)）
Normalize之后的數(shù)據(jù)儲(chǔ)存在seurat[['RNA']]@data這里寄锐。
combined.data <- NormalizeData(combined.data,
normalization.method = "LogNormalize",
scale.factor = 10000)
1进萄、單個(gè)樣本Normalize和多個(gè)樣本合并之后Normalize的結(jié)果是否存在差異。
首先我們合并數(shù)據(jù)的時(shí)候一般直接用的是merge锐峭，所以不同樣本的細(xì)胞的數(shù)值不會(huì)發(fā)生變化。
2可婶、scale factor為什么是10000沿癞？
這里截取whitebird的解釋。
3矛渴、為什么要進(jìn)行l(wèi)og轉(zhuǎn)化椎扬？
做過傳統(tǒng)轉(zhuǎn)錄組分析的家人們都明白惫搏，用轉(zhuǎn)錄組數(shù)據(jù)的FPKM和TPM繪制熱圖等等的時(shí)候，因?yàn)閿?shù)值的變化范圍太過巨大蚕涤，都需要進(jìn)行一個(gè)log轉(zhuǎn)換筐赔，讓數(shù)據(jù)壓縮在一個(gè)區(qū)間。
其次揖铜，也是最重要的改變數(shù)據(jù)分布：測(cè)序數(shù)值本身不符合正態(tài)分布茴丰，log轉(zhuǎn)換能讓數(shù)據(jù)趨近于正態(tài)分布，方便后續(xù)的進(jìn)一步分析天吓。
三贿肩、FindVariableFeatures()
高變異基因 就是highly variable features（HVGs），就是在細(xì)胞與細(xì)胞間進(jìn)行比較龄寞，選擇表達(dá)量差別最大的基因汰规，Seurat使用FindVariableFeatures函數(shù)鑒定高可變基因，這些基因在不同細(xì)胞之間的表達(dá)量差異很大（在一些細(xì)胞中高表達(dá)物邑，在另一些細(xì)胞中低表達(dá)）溜哮。默認(rèn)情況下，會(huì)返回2,000個(gè)高可變基因用于下游的分析色解，如PCA等茂嗓。
1、vst的基本流程
算法實(shí)現(xiàn)在 FindVariableFeatures.default() 中冒签。
目的是在var~mean曲線中在抛，不同mean值區(qū)域都能挑選var較大的基因。
- 使用loess（局部加權(quán)回歸）擬合平滑曲線模型
- 獲取模型計(jì)算的值作為y=var.exp值
- var.standarlized = get variance after feature standardization: (每個(gè)基因 - mean)/sd 后 取var(). 注意sd=sqrt(var.exp)
- 按照 var.standarlized 降序排列萧恕，取前n(比如2000)個(gè)基因作為高變基因刚梭。
combined.data <- FindVariableFeatures(combined.data,
selection.method = "vst",
nfeatures = 2000)
四、ScaleData（）歸一化
單細(xì)胞基因表達(dá)counts矩陣數(shù)據(jù)經(jīng)過NormalizeData()處理后票唆，還需要進(jìn)行scale朴读。
ScaleData()函數(shù)將基因的表達(dá)轉(zhuǎn)換為Z分?jǐn)?shù)（值以 0 為中心，方差為 1）走趋。
它存儲(chǔ)在 seurat_obj[['RNA']]@scale.data衅金，用于下游的PCA降維。
默認(rèn)是僅在高可變基因上運(yùn)行標(biāo)準(zhǔn)化簿煌。
combined.data <- ScaleData(combined.data)
最開始分析單細(xì)胞的時(shí)候氮唯，這里有點(diǎn)疑惑。為什么前面已經(jīng)Normalize姨伟，這里還要scale一下惩琉？
Scaling與Normalization的區(qū)別
- scale改變的是數(shù)據(jù)的范圍，normalize改變的是數(shù)據(jù)的分布夺荒。
- scale是將數(shù)據(jù)的分布限定在一個(gè)范圍內(nèi)瞒渠，這樣子方便比較良蒸。normalize卻是將偏態(tài)分布轉(zhuǎn)換成趨近于正態(tài)分布。
R語言的Z score計(jì)算是通過scale()函數(shù)求得伍玖，Seurat包中ScaleData()函數(shù)也基本參照了scale()函數(shù)的功能嫩痰。
scale方法中的兩個(gè)參數(shù)：center和scale
- center和scale默認(rèn)為真，即T或者TRUE
- center為真表示數(shù)據(jù)中心化：數(shù)據(jù)集中的各項(xiàng)數(shù)據(jù)減去數(shù)據(jù)集的均值
如有數(shù)據(jù)集1, 2, 3, 6, 3窍箍，其均值為3
那么中心化之后的數(shù)據(jù)集為1-3,2-3,3-3,6-3,3-3串纺，即：-2,-1,0,3,0 - scale為真表示數(shù)據(jù)標(biāo)準(zhǔn)化：指中心化之后的數(shù)據(jù)在除以數(shù)據(jù)集的標(biāo)準(zhǔn)差，即數(shù)據(jù)集中的各項(xiàng)數(shù)據(jù)減去數(shù)據(jù)集的均值再除以數(shù)據(jù)集的標(biāo)準(zhǔn)差仔燕。
如有數(shù)據(jù)集1, 2, 3, 6, 3造垛，其均值為3，其標(biāo)準(zhǔn)差為1.87
那么標(biāo)準(zhǔn)化之后的數(shù)據(jù)集為(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87晰搀，即：-1.069,-0.535,0,1.604,0
Z score的概念是指原始數(shù)據(jù)距離均值有多少個(gè)標(biāo)準(zhǔn)差五辽。當(dāng)以標(biāo)準(zhǔn)差為單位進(jìn)行測(cè)量時(shí)，Z score衡量的是一個(gè)數(shù)值偏離總體均值以上或以下多少個(gè)標(biāo)準(zhǔn)差外恕。如果原始數(shù)值高于均值杆逗，則 Z score得分為正，如果低于均值鳞疲，則Z score為負(fù)罪郊。
Z score其實(shí)是標(biāo)準(zhǔn)正態(tài)分布（Standard Normal Distribution），即平均值μ=0尚洽，標(biāo)準(zhǔn)差 σ=1 的正態(tài)分布悔橄。SND標(biāo)準(zhǔn)正態(tài)分布的直方圖如下所示：
五、RunPCA() 降維處理
Seurat使用RunPCA函數(shù)對(duì)標(biāo)準(zhǔn)化后的表達(dá)矩陣進(jìn)行PCA降維處理腺毫。默認(rèn)情況下癣疟，只對(duì)前面選出的2000個(gè)高可變基因進(jìn)行線性降維，也可以通過feature參數(shù)指定想要降維的數(shù)據(jù)集潮酒。
RunPCA之后會(huì)返回5個(gè)PC和對(duì)應(yīng)的一大堆基因睛挚。
每個(gè)PC對(duì)應(yīng)60個(gè)基因，分為Positive和Negative急黎，各30個(gè)扎狱。
Positive和Negative就是PC軸的正負(fù)映射關(guān)系，正值為Positive勃教，負(fù)值為Negative淤击。
返回的是正值和負(fù)值絕對(duì)值最大的top30，可以理解為對(duì)所有細(xì)胞區(qū)分度最大的基因故源。
combined.data <- RunPCA(combined.data)
View(pbmc@reductions[["pca"]]@feature.loadings)
View(pbmc@reductions[["pca"]]@cell.embeddings)
這是PCA之后得到的兩個(gè)結(jié)果遭贸。
第一個(gè)是每個(gè)PC對(duì)應(yīng)的基因。
第二個(gè)是每個(gè)細(xì)胞對(duì)應(yīng)的PC上的坐標(biāo)心软。
得到的PCA的結(jié)果還可以用以下兩種方式來看壕吹。
DimPlot(pbmc,reduction = 'pca')
p1=FeaturePlot(pbmc,features = "PC_1", order = T)
p2=FeaturePlot(pbmc,features = "PC_2", order = T)
p1|p2
六、FindNeighbors()
首先計(jì)算每個(gè)細(xì)胞的KNN删铃，也就是計(jì)算每個(gè)細(xì)胞之間的相互距離耳贬，依據(jù)細(xì)胞之間鄰居的overlap來構(gòu)建snn graph。
計(jì)算給定數(shù)據(jù)集的k.param最近鄰猎唁。也可以選擇(通過compute.SNN)咒劲，通過計(jì)算每個(gè)細(xì)胞最近鄰之間的鄰域重疊(Jaccard索引)和其鄰近的k.param來構(gòu)造SNN。
combined.data <- FindNeighbors(combined.data, dims = 1:30)
具體在計(jì)算細(xì)胞之間的距離的時(shí)候呢诫隅，用得到的KNN算法腐魂，即鄰近算法。
但是逐纬，這個(gè)算法我也不太懂蛔屹，但是其中有個(gè)事情還蠻有趣。
就是判定鄰近的模塊有兩個(gè)：annoy豁生、rann
其中這個(gè)annoy全稱又叫做Approximate Nearest Neighbors Oh Yeah兔毒，名字還蠻可愛的。
下面這位博主對(duì)于這個(gè)算法有非常詳細(xì)的解釋甸箱，有興趣的家人們自行前往觀看育叁。
FindNeighbors {Seurat} - 簡(jiǎn)書 (jianshu.com)
七、FindClusters()
就是在已經(jīng)計(jì)算完細(xì)胞之間的距離之后芍殖，對(duì)這些細(xì)胞進(jìn)行分類豪嗽。
可以指定分為幾類細(xì)胞。
但是很多參考資料里面最重要強(qiáng)調(diào)的都只是一個(gè)參數(shù)：resolution豌骏。
resolution這個(gè)參數(shù)設(shè)置的大小決定了細(xì)胞類型的多少龟梦，值越大細(xì)胞類型越多。
具體分析的時(shí)候很多時(shí)候就會(huì)問：到底多少合理肯适？到底應(yīng)該分為幾群变秦？
其實(shí)對(duì)于測(cè)序的人來講，很多時(shí)候確實(shí)也不太清楚到底有多少種細(xì)胞框舔。
那就有兩種解決辦法蹦玫。
1）直接看tSNE的圖，物理距離就是判斷的一種方法刘绣。當(dāng)物理距離很近的一群細(xì)胞被拆開了樱溉，那就說明可能沒拆開之前是合理的。但是纬凤，這種方法呢就簡(jiǎn)單粗暴一些福贞。
2）有另外一個(gè)包c(diǎn)lustree，可以對(duì)你的分群數(shù)據(jù)進(jìn)行判斷停士。
如下圖中挖帘，當(dāng)分為4完丽、5和6群的時(shí)候，細(xì)胞之間沒有太多的交叉拇舀，都是在進(jìn)一步細(xì)分逻族。
但是再往下，那些半透明的箭頭就表示骄崩，從一個(gè)細(xì)胞群跑到了另外一個(gè)細(xì)胞群聘鳞，說明就不太合適。
所以要拂，參照上述兩種判斷方法抠璃，可以得到結(jié)果。
八脱惰、RunUMAP()搏嗡、RunTSNE()
上述兩種方法其實(shí)都是對(duì)數(shù)據(jù)進(jìn)行降維。
為什么需要降維枪芒？
如果將單個(gè)細(xì)胞看作一個(gè)數(shù)據(jù)點(diǎn)彻况，那么檢測(cè)的基因數(shù)就是其對(duì)應(yīng)的變量數(shù)，也就是我們所說的維數(shù)舅踪。一般一個(gè)人類scrna-seq 能檢測(cè)到～25k基因纽甘，于是基因表達(dá)數(shù)據(jù)也被稱為高維數(shù)據(jù)。我們知道抽碌， 并不是所有的基因都會(huì)表達(dá)悍赢，并且由于單細(xì)胞測(cè)序技術(shù)的限制，不是所有的轉(zhuǎn)錄分子都能被成功捕獲货徙，再加上測(cè)序深度的差異左权， 每個(gè)細(xì)胞中約能檢測(cè)到10%～50%的轉(zhuǎn)錄分子，這導(dǎo)致了許多基因計(jì)數(shù)為0. 雖然通過數(shù)據(jù)的預(yù)處理我們能過濾掉零計(jì)數(shù)基因痴颊，剩下的基因維數(shù)仍然可以高達(dá)～15k赏迟。并且基因之間的高度相關(guān)性也使得表達(dá)數(shù)據(jù)包含了許多冗余信息， 從而掩蓋了真正的生物學(xué)差異蠢棱。因此锌杀，我們需要通過降維抽取數(shù)據(jù)的概要（component），減少數(shù)據(jù)噪音泻仙，并提高下游分析速度糕再。降維的另外一個(gè)好處就是可以在低維中實(shí)現(xiàn)數(shù)據(jù)的可視化。
常見的降維方法有PCA (principle component analysis)玉转、MDS (multi-dimensional scaling)突想、Sammon mapping、Isomap 、t-SNE （t-distributed stochastic neighbor embedding）和UMAP（Uniform Approximation and Projection method）猾担。
其中PCA袭灯， t-SNE和UMAP在scRNA-seq中使用非常普遍。
之前的分析中已經(jīng)用PCA降維了绑嘹，為什么這里還要降維妓蛮？
PCA輸入的是@scale.data，是屬于線性降維
RunUMAP()和RunTSNE()輸入的數(shù)據(jù)是PCA之后的數(shù)據(jù)（PCA@cell.embedding）圾叼，屬于非線性降維。輸出的結(jié)果保存在各自的slots里面捺癞。
這兩種方法有什么區(qū)別夷蚊？
摘抄自：（跟著小魚頭學(xué)單細(xì)胞測(cè)序-scRNA-seq數(shù)據(jù)的降維和可視化 - 云+社區(qū) - 騰訊云 (tencent.com)
）
t-SNE和UMAP是另外兩種非線性的降維方法，由于其漂亮的可視化效果髓介，這兩種方法在單細(xì)胞數(shù)據(jù)教程中非常受歡迎惕鼓。與PCA不同，這兩種方法不僅限于線性變換唐础，也不受制于準(zhǔn)確表示遠(yuǎn)距離種群之間的距離箱歧，因此使得它們?cè)诘臀豢臻g如何對(duì)細(xì)胞進(jìn)行排列具有更大的自由度，從而在其二維的可視化效果圖中能夠?qū)?fù)雜的細(xì)胞群分成許多不同的簇一膨，使得效果圖比PCA更加直觀和容易解釋呀邢。
t-SNE是以犧牲整體結(jié)構(gòu)（global structure）為代價(jià)，著重于捕獲局部結(jié)構(gòu)（local similarly）豹绪，因此它可能會(huì)夸大細(xì)胞群體之間的差異而忽略群體之間的潛在聯(lián)系价淌。如果簡(jiǎn)單的通過t-SNE結(jié)果圖來解讀細(xì)胞簇之間的關(guān)系，產(chǎn)生的結(jié)論可能會(huì)被圖中的集群大小和位置所誤導(dǎo)瞒津。在效果圖中蝉衣，t-SNE傾向于將密集的簇膨脹，并且壓縮稀疏的簇巷蚪，因此我們不能簡(jiǎn)單的通過圖上的大小來衡量細(xì)胞群的差異病毡。并且由于t-sne無法保證能保留距離較遠(yuǎn)的簇的相對(duì)位置，我們也不能簡(jiǎn)單通過圖中的位置來確定遠(yuǎn)距離細(xì)胞簇之間的關(guān)系屁柏。相比之下啦膜，UMAP與t-SNE類似，同時(shí)在低維空間保留了高維空間細(xì)胞間的關(guān)系前联，因此UMAP更好的保留并反映了細(xì)胞群潛在的拓?fù)浣Y(jié)構(gòu)功戚，對(duì)于細(xì)胞軌跡推斷（trajectory inference）分析來說更實(shí)用 [4]。
在下圖中似嗤，我們可以看到UMAP的可視化更趨向于一種緊湊的視覺效果啸臀，群簇之間的空間更大，也保留了更多的global structure，因此在選擇可視化圖中大家可以根據(jù)具體的需要來選擇乘粒。根據(jù)小編的實(shí)戰(zhàn)經(jīng)驗(yàn)豌注，如果對(duì)軌跡推斷沒有要求，只是看細(xì)胞簇群灯萍，UMAP的效果會(huì)更干凈轧铁，但是在細(xì)胞數(shù)目非常大情況下，由于UMAP最大限度的保留了全局結(jié)構(gòu)旦棉，這也使得每個(gè)簇群的分辨率降低齿风，可能會(huì)使得簇群重疊，從而遮蓋一些小的簇群绑洛，而t-SNE通常能將所有簇群盡可能的“鋪開”救斑，所以這種情況建議大家兩種都畫，然后比較一下真屯。從速度方面來說脸候，同一個(gè)數(shù)據(jù)UMAP的速度要比t-SNE快，這也是UMAP變得更受歡迎的重要原因绑蔫。
九运沦、降維之后的可視化DimPlot、featureplot等等配深。
處理完上述數(shù)據(jù)之后携添，可視化就很簡(jiǎn)單了。
具體的參數(shù)自行查閱凉馆。
轉(zhuǎn)自：http://www.reibang.com/p/9e0fb7c2c43f
