來源
nature原文：Cluster-level quality control was performed after
the standard Seurat clustering pipeline using the following functions
in order: FindVariableFeatures with 2,000 genes, ScaleData, RunPCA,
FindNeighbors with the first 20 PCs and FindClusters with resolution 1,
otherwise default settings
nature原文：Highly variable genes were identified by fitting the mean variance relationship for each sample, to avoid selecting for genes with highly variable between-sample effects and to prioritize those with high within-sample variation
目的：避免選擇的基因有批次效應(yīng)
0. 降維
- scRNA分析就是在于數(shù)據(jù)打交道。本來看不見摸不著的基因斗遏，現(xiàn)在體現(xiàn)在數(shù)據(jù)上，雖然還是摸不著梅猿，但能看得到氓辣。每個(gè)基因現(xiàn)在都作為數(shù)據(jù)的一個(gè)維度存在。假如我們有一個(gè)scRNA數(shù)據(jù)集袱蚓，其中只有兩個(gè)基因钞啸，其中兩個(gè)坐標(biāo)軸表示兩個(gè)基因的表達(dá)量，圖中的點(diǎn)就表示細(xì)胞喇潘，每個(gè)細(xì)胞在圖上都由xy軸確定体斩。也就是說，基因的表達(dá)決定了細(xì)胞的分布颖低。如果有兩個(gè)基因絮吵，那么細(xì)胞的位置就由兩個(gè)維度決定;三個(gè)基因，就由三維空間決定;而我們有2萬多個(gè)基因忱屑，那么就存在2萬多個(gè)維度（高維空間)蹬敲。
- 單細(xì)胞分析的一個(gè)重點(diǎn)就是: 用盡可能少的維度來展示數(shù)據(jù)真實(shí)的結(jié)構(gòu)。降維的過程其實(shí)就是去繁（基因的變化休航洹）存簡伴嗡，每個(gè)基因的變化都對(duì)整體有影響，對(duì)細(xì)胞來講都是一個(gè)變化維度从铲。但這么多維度我們看不了也處理不了闹究，需要盡可能保證真實(shí)差異的前提下減少維度的數(shù)量，因此需要挑出那些更能代表整體差異的基因食店。
-
如何挑選有價(jià)值的基因(feature)
一般來說渣淤，如果一個(gè)基因在細(xì)胞群體中變化幅度很大，就是受關(guān)注對(duì)象吉嫩，我們也會(huì)認(rèn)為是生物因素導(dǎo)致了這么大的差異价认。這樣的基因叫做高度變化基因(highly variable genes ,HyGs)，y降維和聚類分群也是常用的操作自娩，它們都是依賴于基因表達(dá)量（例如綜合細(xì)胞中各個(gè)基因的表達(dá)量用踩，然后把表達(dá)模式相近的細(xì)胞聚在一起)。因此挑選哪些基因進(jìn)行聚類分群和降維分析忙迁，這是非常關(guān)鍵的脐彩，挑選的基因一定要有代表性，盡可能多的包含生物信息而剔除其他技術(shù)噪音姊扔。 - 官方推薦使用2000個(gè)高變基因做后續(xù)的降維
1. 挑選高變基因
FindVariableFeatures(function)
scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst", nfeatures = 3000)
#nfeature的小提琴圖中惠奸，大部分細(xì)胞的基因表達(dá)在2000左右，高變基因閾值越高恰梢，丟失的信息會(huì)越少（維度越多佛南，信息量越多）
#但是選的太多梗掰，降維幅度不大，挑選高變基因意義也不會(huì)太大
top10 <- head(VariableFeatures(scRNA1), 10)
#把top10高變基因挑選出來
plot1 <- VariableFeaturePlot(scRNA1)
#畫出不帶標(biāo)簽的高變基因圖
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5)
#把top10的基因加到圖中
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom")
plot
-
結(jié)果
橫坐標(biāo)：基因的平均表達(dá)量嗅回；縱坐標(biāo)：每個(gè)基因的變異程度
2.scaledata
#如果內(nèi)存足夠最好對(duì)所有基因進(jìn)行中心化
scale.genes <- rownames(scRNA1)
scRNA1 <- ScaleData(scRNA1, features = scale.genes)
##如果內(nèi)存不夠,可以只對(duì)高變基因進(jìn)行標(biāo)準(zhǔn)化（淺推薦）及穗，但是后續(xù)熱圖只能查看高變基因的變化
scale.genes <- VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
-
什么是中心化處理？
Scaling the data(數(shù)據(jù)比例伸縮)
ScaleData函數(shù):
改變每個(gè)基因的表達(dá)绵载，使細(xì)胞間的平均表達(dá)為0
測量每個(gè)基因的表達(dá)埂陆，使細(xì)胞間的差異為1
我們一般將平均值為0，方差值為1的數(shù)據(jù)認(rèn)為是標(biāo)準(zhǔn)數(shù)據(jù)娃豹。改變每個(gè)基因的表達(dá)焚虱，使得跨細(xì)胞的平均表達(dá)為0
縮放每個(gè)基因的表達(dá)，以便跨細(xì)胞的方差為1
該步驟在下游分析中給予相同的權(quán)重培愁，因此高表達(dá)的基因不占優(yōu)勢 scRNA對(duì)象中原始表達(dá)矩陣經(jīng)過標(biāo)準(zhǔn)化和中心化之后，已經(jīng)產(chǎn)生了三套基因表達(dá)數(shù)據(jù)缓窜，可以通過以下命令獲得
GetAssayData(scRNA,slot="counts",assay="RNA")
#原始表達(dá)矩陣
GetAssayData(scRNA,slot= "data",assay="RNA")
#標(biāo)準(zhǔn)化之后的表達(dá)矩陣
GetAssayData(scRNA,slot="scale.data",assay="RNA")
#中心化之后的表達(dá)矩陣
Normalize函數(shù)：會(huì)使用count數(shù)據(jù)（原始數(shù)據(jù)）定续，將歸一化后的結(jié)果放入data數(shù)據(jù)
-
什么是稀疏矩陣（dgcMatrix）
0是用點(diǎn)進(jìn)行表現(xiàn)的（節(jié)省內(nèi)存），查看此類數(shù)據(jù)結(jié)構(gòu)時(shí)應(yīng)該用as.dataframe函數(shù)
3.檢查細(xì)胞周期對(duì)降維聚類的影響
細(xì)胞周期回歸: 上一步找到的高變基因禾锤，常常會(huì)包含一些細(xì)胞周期相關(guān)基因私股。它們會(huì)導(dǎo)致細(xì)胞聚類發(fā)生一定的偏移，即相同類型的細(xì)胞在聚類時(shí)會(huì)因?yàn)榧?xì)胞周期的不同而分開恩掷。
- ?CaseMatch——Match the case of character vectors
#cc.genes包內(nèi)自帶數(shù)據(jù)倡鲸；s期的基因有哪些（s.genes），g2m期的基因有哪些（g2m.genes）
cc.genes
#將cc.genes的數(shù)據(jù)組成一個(gè)新的character黄娘，與高變基因match（PCA降維只用了高變基因）
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA1))
#g2m基因與高變基因匹配
g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match =rownames(scRNA1))
#s基因與高變基因匹配
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA1))
#cellcyclescoring對(duì)每個(gè)細(xì)胞的細(xì)胞周期進(jìn)行打分
scRNA1 <- CellCycleScoring(object=scRNA1, g2m.features=g2m_genes, s.features=s_genes)
#得到周期基因的PCA降維結(jié)果
scRNAa <- RunPCA(scRNA1, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
p
#scRNAb <- ScaleData(scRNA1, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA1))
#vars.to.regress 可消除周期對(duì)降維的全部影響
結(jié)果：
4.PCA降維
-
PCA
PCA (Principal Component Analysis优床，主成分分析)是一種降維方法，致力于解決三類問題:
降維可以緩解維度災(zāi)難問題;
降維可以在壓縮數(shù)據(jù)的同時(shí)讓信息損失最小化;
理解幾百個(gè)維度的數(shù)據(jù)結(jié)構(gòu)很困難誓焦，兩三個(gè)維度的數(shù)據(jù)通過可視化更容易理解胆敞。 - 隨著維度的增加，數(shù)據(jù)的稀疏性會(huì)越來越高杂伟。在高維向量空間中探索同樣的數(shù)據(jù)集比在同樣稀疏的數(shù)據(jù)集中探索更加困難移层。主成分分析也稱為卡爾胡寧-勒夫變換(Karhunen-Loeve Transform)，是一種用于探索高維數(shù)據(jù)結(jié)構(gòu)的技術(shù)赫粥。PCA通常用于高維數(shù)據(jù)集的探索與可視化观话，還可以用于數(shù)據(jù)壓縮，數(shù)據(jù)預(yù)處理等越平。PCA可以把可能具有相關(guān)性的高維變量合成線性無關(guān)的低維變量匪燕，稱為主成分蕾羊。新的低維數(shù)據(jù)集會(huì)盡可能的保留原始數(shù)據(jù)的變量。
- PCA把相關(guān)的基因合并到“metagene”或主成分(PC)中帽驯。PC1解釋最大的數(shù)據(jù)差異龟再，具有最大的標(biāo)準(zhǔn)差（例如對(duì)于一個(gè)實(shí)驗(yàn)，細(xì)胞之間30%的差異由定義了PC1的基因解釋)尼变，PC2解釋了數(shù)據(jù)的第二大部分差異(例如利凑，細(xì)胞之間20%的差異可歸因于PC2中的基因，而8%則歸因于PC3中的基因)嫌术，然后依此類推哀澈，簡單來說PC的排名就是解釋數(shù)據(jù)差異貢獻(xiàn)的順序，其中PC1是排名最高的PC度气，同時(shí)也說明PC排名越低割按，對(duì)解釋數(shù)據(jù)差異的貢獻(xiàn)就越小。
5.聚類
nature原文：Cluster-level quality control was performed after the standard Seurat clustering pipeline using the following functions in order ： FindVariableFeatures with 2,000 genes , ScaleData , RunPCA , FindNleighbors with the first 20 PCs and FindClusters with resolution 1/otherwise defaultsettings.
-
第一步：選擇合適的PCs數(shù)
為了克服scRNA序列數(shù)據(jù)單一特征中的廣泛技術(shù)噪音磷籍，Seurat根據(jù)其PCA分?jǐn)?shù)對(duì)細(xì)胞進(jìn)行聚類适荣，每個(gè)PC基本上表示一個(gè)“元特征”，該特征結(jié)合了相關(guān)特征集上的信息院领。因此弛矛，最主要的主成分代表數(shù)據(jù)集的強(qiáng)大壓縮。但是, 我們應(yīng)該選擇多少個(gè)PC? 10個(gè)?20?還是100? - 利用elbow point選擇
ElbowPlot(scRNA1, ndims=20, reduction="pca")
elbow point就是在它之前變化幅度很大比然，之后變化幅度很小丈氓，屬于一個(gè)轉(zhuǎn)折點(diǎn)。啟發(fā)式評(píng)估方法生成一個(gè)Elbow plot圖强法。在圖中展示了每個(gè)主成分對(duì)數(shù)據(jù)方差的解釋情況(百分比表示)万俗，并進(jìn)行排序。根據(jù)自己需要選擇主成分饮怯。
- 每個(gè)PCs都能捕獲一些生物差異，而且前面的PC比后面的PC包含的差異信息更多硕淑，更有價(jià)值课竣。就像"二八準(zhǔn)則"，僅有20%的變因操縱著80%的局面置媳，PCA中也一樣于樟。于是當(dāng)我們從頭到尾觀察每個(gè)PC的差異貢獻(xiàn)時(shí)，會(huì)有一個(gè)明顯的落差拇囊，落差之前的那些PC就是"二八準(zhǔn)則"中20%的變因迂曲。
- 如果選的PC多，可以避免丟棄一些生物信號(hào)寥袭，但同時(shí)又增加了噪音的風(fēng)險(xiǎn)路捧。很多有經(jīng)驗(yàn)的人會(huì)將PC數(shù)設(shè)置在一個(gè)合理的區(qū)間关霸，例如10-50之間
- 曲線下面積代表總的貢獻(xiàn)度
- 注意: 鑒別數(shù)據(jù)的真實(shí)維度不是件容易的事情，除了上面兩種方法杰扫，Seurat官方文檔還建議將主成分(數(shù)據(jù)異質(zhì)性的相關(guān)來源有關(guān))與GSEA分析相結(jié)合队寇。如Dendritic cell和NKaficionados可能識(shí)別的基因與主成分12和13相關(guān)，定義了罕見的免疫亞群(i.e.MZB1 is a marker for plasmacytoid DCs)章姓。如果不是事先知道的情況下佳遣，很難發(fā)現(xiàn)這些問題。
Seurat官方文檔因此鼓勵(lì)用戶使用不同數(shù)量的PC(10凡伊、15甚至50)重復(fù)下游分析零渐。其實(shí)觀察到的結(jié)果通常沒有顯著差異。因此系忙，在選擇此參數(shù)時(shí)诵盼，可以盡量選大一點(diǎn)的維度，維度太小的話會(huì)對(duì)結(jié)果產(chǎn)生不好的影響银还。
第二步风宁、聚類分析
- 常見的聚類如層次聚類、k-means等见剩，以及為單細(xì)胞專門開發(fā)的SC3聚類杀糯。但是Seurat采用的是譜聚類(Spectral Clustering ) 扫俺。Seurat中的譜聚類基于共享最近鄰圖(SNN）和模塊化優(yōu)化的聚類算法來識(shí)別細(xì)胞簇苍苞。它首先計(jì)算k最近鄰(k-nearestneighbors)并構(gòu)造SNN圖，然后優(yōu)化模塊化功能以確定具體類群狼纬。這個(gè)算法最早在2013年發(fā)表在物理學(xué)的期刊The European Physical Journal B上羹呵。
第一個(gè)難點(diǎn)是譜聚類是一種基于圖論方法
圖論方法中的"圖“和圖像的"圖""不一樣，它屬于離散數(shù)學(xué)疗琉。譜聚類是從圖論中演化出來的算法冈欢，它的主要思想是把所有的細(xì)胞看為高維空間中的點(diǎn)，這些點(diǎn)之間可以用邊連接起來盈简。距離較遠(yuǎn)的兩個(gè)點(diǎn)之間的邊權(quán)重值較低; 距離較近的兩個(gè)點(diǎn)之間的邊權(quán)重值較高凑耻。然后通過對(duì)所有數(shù)據(jù)點(diǎn)組成的圖進(jìn)行"切圖"，讓切圖后不同的子圖間邊權(quán)重和盡可能的低柠贤，而子圖內(nèi)的邊權(quán)重和盡可能的高香浩，從而達(dá)到聚類的目的。第二個(gè)難點(diǎn)是如何把我們的單細(xì)胞的數(shù)據(jù)構(gòu)建成為無向有權(quán)圖
譜聚類中的圖采用的是無向有權(quán)圖臼勉，也就是說邻吭，每個(gè)細(xì)胞（圖中的點(diǎn)）之間的連接是沒有方向性的。但是每條邊上是有權(quán)重信息的宴霸。類比一下地圖囱晴，兩地之間的路是沒有方向(雙向的)膏蚓，但是地點(diǎn)與地點(diǎn)之間路的遠(yuǎn)近距離是不同的。我們手中最初只有細(xì)胞的表達(dá)量矩陣畸写，這個(gè)每個(gè)細(xì)胞中各個(gè)基因的表達(dá)情況驮瞧，屬于“節(jié)點(diǎn)信息"。
那么哪些細(xì)胞之間需要有邊相連艺糜，邊的權(quán)重又該如何計(jì)算呢?
- 這里引入了鄰接矩陣(W)的概念剧董，我們首先定義一種計(jì)算細(xì)胞之間距離度量的方法，來計(jì)算每條邊之間的權(quán)重wij破停，又稱相似性sij翅楼，然后獲得整個(gè)的鄰接矩陣W。
- 終于到了“切圖"的環(huán)節(jié)真慢，實(shí)際就是將圖切成相互沒有連接的k個(gè)子圖毅臊。目標(biāo)就是讓子圖內(nèi)的點(diǎn)權(quán)重和高，子圖間的點(diǎn)權(quán)重和低黑界。實(shí)際就是一個(gè)最優(yōu)化問題管嬉。可以理解為通過D朗鸠，L這兩個(gè)矩陣所給的信息蚯撩，在圖的部分邊進(jìn)行切分震贵，讓完整的圖變成許多子圖蛋哭，每一個(gè)子圖就構(gòu)成了一個(gè)類群。
-
resolution參數(shù)決定類群的“粒度”型型，官網(wǎng)建議將此參數(shù)設(shè)置在0.4-1.2之間忆家，通秤坦剑可為包含大約3000個(gè)cell的數(shù)據(jù)集會(huì)返回良好的聚類結(jié)果。對(duì)于較大的數(shù)據(jù)集芽卿，最佳resolution值通常會(huì)增加揭芍。
（resolution越高，切的刀數(shù)越多）
nature原文：FindVariableFeatures with 2,000 genes, ScaleData, RunPCA, FindNeighbors with the first 20 PCs and FindClusters with resolution 1, otherwise default settings
pc.num=1:20#確定維度
#根據(jù)20個(gè)維度的信息構(gòu)建無向有權(quán)圖
scRNA1 <- FindNeighbors(scRNA1, dims = pc.num)
#resolution=1.0進(jìn)行切圖卸例，目的是尋找cluster称杨，cluster的結(jié)果會(huì)存放到graph中
scRNA1 <- FindClusters(scRNA1, resolution = 1.0)
系統(tǒng)發(fā)育分析(Phylogenetic Analysis of ldentity Classes)
【但是似乎很少有paper放這個(gè)】
nature原文：
Constructs a phylogenetic tree relating the ’average‘ cell from each identity class.
Tree is estimated based on a distance matrix constructed in either gene expression space or PCA spac
scRNA1<-BuildClusterTree(scRNA1)
PlotClusterTree(scRNA1)
6.可視化
- 降維可以定義為將高維數(shù)據(jù)降低為低維同時(shí)盡可能保證數(shù)據(jù)的全局性筷转，降維技術(shù)是機(jī)器學(xué)習(xí)中大數(shù)據(jù)預(yù)處理的一種重要技術(shù)姑原。
- 降維的算法主要分為兩類: (1)在數(shù)據(jù)中保持距離結(jié)構(gòu)的算法——線性降維(PCA、MDS和Sammon映射)旦装，(2)在全局距離上保持局部距離的算法（流形學(xué)習(xí))——非線性降維（t-SNE页衙、Jsomap、LargeVis、Laplacian特征映射和difusion映射)店乐。
- 非線性降維——這個(gè)目的是為了可視化艰躺，而不是特征提取（PCA)眨八，雖然它也可以用來做特征提取腺兴。tSNE和UMAP目前主要用于可視化。用于可視化的降維必然涉及信息丟失并改變細(xì)胞之間的距離廉侧。因此tSNE/UMAP圖應(yīng)僅只用于解釋或傳達(dá)基于更精確的页响、更多維度的定量分析結(jié)果。這樣可以保證分析充分利用了壓縮到二維空間時(shí)丟失的信息段誊。假如二維圖上呈現(xiàn)的細(xì)胞分布與使用更多數(shù)目的PC進(jìn)行聚類獲得的結(jié)果之間存在差異闰蚕，應(yīng)傾向于相信前者（聚類）的結(jié)果。
tSNE和UMAP的區(qū)別
t-分布式隨機(jī)鄰居嵌入(t-SNE: t-Distributed stochastic neighbor embedding)是一種常見的可視化方法连舍。它使用機(jī)器學(xué)習(xí)的算法來降低維數(shù)没陡，非常適合將高維數(shù)據(jù)放到二維或三維空間中可視化展示，并且不會(huì)丟失細(xì)胞之間的相對(duì)距離的信息索赏。
例如盼玄，如果發(fā)現(xiàn)用七個(gè)PC可以很好地表示細(xì)胞的多樣性，就得需要七個(gè)軸或維度來展示細(xì)胞的空間分布潜腻。t-SNE能維持細(xì)胞在七維空間的關(guān)系并在二維圖上展示細(xì)胞埃儿，所以在七維圖上相鄰的細(xì)胞在二維圖上仍然相鄰。同時(shí)PCA分析是線性的融涣。t-SNE是非線性降維方法童番。
統(tǒng)一流形逼近與投影(UMAP：Uniform Manifold Approximation and Projection)是一種新的降維技術(shù)，其理論基礎(chǔ)是黎曼幾何和代數(shù)拓?fù)浔┬摹Ｏ鄬?duì)于T-SNE降維妓盲，UMAP的優(yōu)點(diǎn)在于:(1)能夠盡可能多的保留全局結(jié)構(gòu)（2）耗時(shí)更短 (3）對(duì)嵌入維數(shù)沒有限制杂拨，可以擴(kuò)展到更大的維度的數(shù)據(jù)集专普。而T-SNE對(duì)數(shù)據(jù)的維度有所限制（Hinton的T-SNE實(shí)驗(yàn)中先用PCA降到50,再進(jìn)一步的使用T-SNE，密集的數(shù)據(jù)才能達(dá)到較好的可視化效果)弹沽。T-SNE是建立在二維的基礎(chǔ)上檀夹，而UMAP則是建立在三維的黎曼流行上，其相關(guān)定義策橘、計(jì)算方法都會(huì)發(fā)生一定的變化炸渡。
- t-SNE圖
#先選擇多少個(gè)維度
scRNA1 = RunTSNE(scRNA1, dims = pc.num)
embed_tsne <- Embeddings(scRNA1, 'tsne')
write.csv(embed_tsne,'embed_tsne.csv')
plot1 = DimPlot(scRNA1, reduction = "tsne")
plot1
#label = TRUE把注釋展現(xiàn)在圖中
DimPlot(scRNA1, reduction = "tsne",label = TRUE)
#把cluster都標(biāo)記在圖中
ggsave("tSNE.pdf", plot = plot1, width = 8, height = 7)
#保存
結(jié)果：
- UMAP圖
scRNA1 <- RunUMAP(scRNA1, dims = pc.num)
embed_umap <- Embeddings(scRNA1, 'umap')
write.csv(embed_umap,'embed_umap.csv')
plot2 = DimPlot(scRNA1, reduction = "umap")
結(jié)果：
