版本信息：
Seurat v2.0不是3.0嘱函！現(xiàn)在Seurat更新了3.0版本，下載也是默認(rèn)的3.0埂蕊，這篇記錄只適用于用2.0的往弓。
梗概
- 將Cellranger中的基因表達(dá)矩陣filtered_gene_bc_matrices用于分析。
- 進(jìn)行質(zhì)量控制（QC）蓄氧，以刪除異常細(xì)胞函似；
- 標(biāo)準(zhǔn)化與歸一化，消除技術(shù)噪音與批次效應(yīng)匀们；
- 主成分分析（PCA）與挑選
- t-SNE聚類
參考網(wǎng)站：https://satijalab.org/seurat/pbmc3k_tutorial.html
（注意＝闪堋Ｗ几泄朴！現(xiàn)在這個(gè)網(wǎng)站會(huì)自動(dòng)跳轉(zhuǎn)到3.0版本）
Seurat的安裝：R中運(yùn)行install.packages("Seurat")
上次結(jié)果：
經(jīng)過(guò)Cellranger的數(shù)據(jù)整理之后重抖，得到：
- Filtered gene-barcode matrices MEX: /outs/filtered_gene_bc_matrices
此輸出結(jié)果應(yīng)為基因-細(xì)胞的表達(dá)矩陣，用Seurat包進(jìn)行后續(xù)分析祖灰。
Seurat是一種R包钟沛，設(shè)計(jì)用于QC，分析和探索單細(xì)胞RNA-seq數(shù)據(jù)局扶。 Seurat旨在使用戶能夠從單細(xì)胞轉(zhuǎn)錄組測(cè)量中識(shí)別和解釋異質(zhì)性來(lái)源恨统，并整合不同類型的單細(xì)胞數(shù)據(jù)。
運(yùn)行R三妈，并且加載這兩個(gè)包
library(Seurat)
library(dplyr)
讀取數(shù)據(jù)
spleen.data <- Read10X(data.dir = '/GRCh38/')
dim(spleen.data)
[1] 33694 1960
原始數(shù)據(jù)的基因數(shù)為33694畜埋，細(xì)胞數(shù)為1960.
比較普通與疏松矩陣的內(nèi)存使用：
> dense.size <- object.size(x = as.matrix(x = spleen.data))
> dense.size
530488272 bytes
#轉(zhuǎn)化為疏松矩陣，查看大小
> sparse.size <- object.size(x = spleen.data)
> sparse.size
45955656 bytes
> dense.size/sparse.size
11.5 bytes
初始化Seurat對(duì)象：
命令CreateSeuratObject
輸入數(shù)據(jù)spleen.data
留下所有在>=3個(gè)細(xì)胞中表達(dá)的基因min.cells = 3畴蒲；
留下所有檢測(cè)到>=200個(gè)基因的細(xì)胞min.genes = 200悠鞍。
(為了除去一些質(zhì)量差的細(xì)胞)
spleen <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")
spleen
An object of class seurat in project 10X_spleen
15655 genes across 1959 samples.
剩下15655 基因和 1959 個(gè)細(xì)胞
質(zhì)量控制
以下步驟包括Seurat中scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理工作流程。這些代表了Seurat對(duì)象的創(chuàng)建模燥，基于QC指標(biāo)的細(xì)胞選擇和過(guò)濾咖祭，數(shù)據(jù)標(biāo)準(zhǔn)化和縮放，以及高度可變基因的檢測(cè)蔫骂。
mito.genes <- grep(pattern = "^MT-", x = rownames(x = spleen@data), value = TRUE)
percent.mito <- Matrix::colSums(spleen@raw.data[mito.genes, ])/Matrix::colSums(spleen@raw.data)
spleen <- AddMetaData(object = spleen, metadata = percent.mito, col.name = "percent.mito")
VlnPlot(object = spleen, features.plot = c("nGene", "nUMI", "percent.mito"), nCol = 3)
> par(mfrow = c(1, 2))
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "percent.mito")
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "nGene")
過(guò)濾細(xì)胞么翰，根據(jù)上面的兩幅圖，去除異常值辽旋，這里選擇基因數(shù)從300-5000浩嫌，線粒體基因占比大于0.1的細(xì)胞。（主要看小提琴圖1和圖3）
spleen <- FilterCells(spleen, subset.names = c("nGene", "percent.mito"), low.thresholds = c(300, -Inf), high.thresholds = c(5000,0.10))
查看過(guò)濾掉剩下多少細(xì)胞：
spleen
An object of class seurat in project 10X_spleen
15655 genes across 1940 samples.
剩下15655個(gè)基因补胚，1940個(gè)細(xì)胞固该。
數(shù)據(jù)標(biāo)準(zhǔn)化
加個(gè)log：
spleen <- NormalizeData(object=spleen, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
spleen <- FindVariableGenes(object = spleen, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
Calculating gene means
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
TEXT_SHOW_BACKTRACE environmental variable.
> length(x=spleen@var.genes)
[1] 1829
縮放數(shù)據(jù)并刪除不需要的變體來(lái)源
您的單細(xì)胞數(shù)據(jù)集可能包含“不感興趣”的變異來(lái)源。這不僅包括技術(shù)噪音糖儡，還包括批次效應(yīng)伐坏，甚至包括生物變異來(lái)源（細(xì)胞周期階段）。正如(Buettner, et al NBT握联，2015)中所建議的那樣桦沉，從分析中回歸這些信號(hào)可以改善下游維數(shù)減少和聚類。為了減輕這些信號(hào)的影響金闽，Seurat構(gòu)建線性模型以基于用戶定義的變量預(yù)測(cè)基因表達(dá)纯露。這些模型的縮放得分殘差存儲(chǔ)在Scale.data槽中，用于降維和聚類代芜。
我們可以消除由批次（如果適用）驅(qū)動(dòng)的基因表達(dá)中的細(xì)胞 - 細(xì)胞變異埠褪，細(xì)胞比對(duì)率（由Drop-seq數(shù)據(jù)的Drop-seq工具提供），檢測(cè)到的分子數(shù)量和線粒體基因表達(dá)。對(duì)于循環(huán)細(xì)胞钞速，我們還可以學(xué)習(xí)“細(xì)胞周期”評(píng)分（參見(jiàn)此處的示例）并對(duì)其進(jìn)行回歸贷掖。在這個(gè)有絲分裂后血細(xì)胞的簡(jiǎn)單例子中，我們回歸了每個(gè)細(xì)胞檢測(cè)到的分子數(shù)量以及線粒體基因含量百分比渴语。
spleen <-ScaleData(spleen, vars.to.regress = c("nUMI","percent.mito"))
Regressing out: nUMI, percent.mito
|=========================================================================================| 100%
Time Elapsed: 18.0711550712585 secs
Scaling data matrix
|=========================================================================================| 100%
PCA分析
主成分分析是什么苹威？
主成分分析，是考察多個(gè)變量間相關(guān)性一種多元統(tǒng)計(jì)方法驾凶，研究如何通過(guò)少數(shù)幾個(gè)主成分來(lái)揭示多個(gè)變量間的內(nèi)部結(jié)構(gòu)牙甫，即從原始變量中導(dǎo)出少數(shù)幾個(gè)主成分，使它們盡可能多地保留原始變量的信息调违，且彼此間互不相關(guān).通常數(shù)學(xué)上的處理就是將原來(lái)P個(gè)指標(biāo)作線性組合窟哺，作為新的綜合指標(biāo)。
將數(shù)據(jù)集降維技肩，利用低階的變量去反應(yīng)整體的結(jié)果脏答。
spleen <- RunPCA(spleen, pc.genes = spleen@var.genes, do.print = TRUE, pcs.print = 1:5, genes.print = 5)
[1] "PC1"
[1] "CD69" "CD79A" "TRAC" "CD3D" "MS4A1"
[1] ""
[1] "FCN1" "LYZ" "SERPINA1" "CSTA" "RP11-1143G9.4"
[1] ""
[1] ""
[1] "PC2"
[1] "CD79A" "MS4A1" "VPREB3" "CD79B" "HLA-DQB1"
[1] ""
[1] "NKG7" "CST7" "GZMA" "CD7" "CCL5"
[1] ""
[1] ""
[1] "PC3"
[1] "TRDC" "KLRF1" "MS4A1" "CD79B" "IRF8"
[1] ""
[1] "IL7R" "TRAC" "CD3D" "CD2" "CD3G"
[1] ""
[1] ""
[1] "PC4"
[1] "GIMAP7" "GZMB" "FGFBP2" "SPON2" "PRF1"
[1] ""
[1] "BAG3" "HSPD1" "FKBP4" "DNAJA1" "ZFAND2A"
[1] ""
[1] ""
[1] "PC5"
[1] "UBE2C" "TYMS" "MKI67" "TOP2A" "AURKB"
[1] ""
[1] "FCGR3A" "FGFBP2" "SPON2" "GNLY" "GZMB"
[1] ""
[1] ""
PCElbowPlot(spleen)
選擇了前10個(gè)PC成分
spleen <- FindClusters(spleen, reduction.type = "pca", dims.use = 1:10, resolution = 0.6, print.output = 0, save.SNN = TRUE)
PrintFindClustersParams(spleen)
Parameters used in latest FindClusters calculation run on: 2018-10-01 21:59:55
=============================================================================
Resolution: 0.6
-----------------------------------------------------------------------------
Modularity Function Algorithm n.start n.iter
1 1 100 10
-----------------------------------------------------------------------------
Reduction used k.param prune.SNN
pca 30 0.0667
-----------------------------------------------------------------------------
Dims used in calculation
=============================================================================
1 2 3 4 5 6 7 8 9 10
細(xì)胞聚類
spleen <- RunTSNE(spleen, dims.use = 1:10, do.fast= TRUE)
TSNEPlot(spleen)
> saveRDS(spleen, file = "/spleen_1.rds")
將R變量保存，利于后續(xù)的分析亩鬼。
一些補(bǔ)充：
過(guò)濾低質(zhì)量細(xì)胞：
在 scRNA-seq 分析中殖告，有些細(xì)胞質(zhì)量比較低，比如細(xì)胞處于凋亡狀態(tài),細(xì)胞中 RNA 發(fā)生降解等,這些細(xì)胞的存在會(huì)影響分析雳锋，因此我們第一步需要對(duì)細(xì)胞進(jìn)行過(guò)濾黄绩。主要可分為三類:
①利用細(xì)胞檢測(cè)到的基因數(shù)或者是 reads 比對(duì)率來(lái)判斷技術(shù)噪音。
但不管是基因檢測(cè)數(shù)目還是比對(duì)率都跟實(shí)驗(yàn)方法有很大相關(guān)性玷过。 如果比對(duì)率太低,表明 RNA 可能發(fā)生了降解,或者文庫(kù)有污染或者細(xì)胞裂解不完全爽丹。
②如果實(shí)驗(yàn)中加入了 spike-ins（本實(shí)驗(yàn)沒(méi)有），可以通過(guò)計(jì)算比對(duì)到內(nèi)源性 RNA 和外源性 RNA(spike-ins)的 reads 比例來(lái)過(guò)濾低質(zhì)量細(xì)胞辛蚊。
比值偏低表明細(xì)胞中的 RNA 數(shù)量較低粤蝎，細(xì)胞可丟棄。但是也需要注意其實(shí)當(dāng)細(xì)胞狀態(tài)不一樣袋马，比如處于不同細(xì)胞周期時(shí)初澎，細(xì)胞的 RNA 數(shù)量是具有很大差異的。不過(guò)我們依然認(rèn)為在一大群細(xì)胞中虑凛，spike-ins比例特別高的細(xì)胞在很大概率上應(yīng)該被排除在外碑宴。軟件 SinQC (Single-cell RNA-seq Quality Control)可以根據(jù)比對(duì)率和檢測(cè)到的基因數(shù)來(lái)過(guò)濾細(xì)胞。
③根據(jù)整體的基因表達(dá)譜來(lái)定義技術(shù)噪音桑谍。
比如對(duì)細(xì)胞進(jìn)行聚類分析延柠，PCA 分析等，將 outlier 細(xì)胞刪除锣披，或者細(xì)胞表達(dá)中位值低于某一設(shè)定閾值時(shí)將該細(xì)胞過(guò)濾掉贞间。當(dāng)然這種方法也存在誤刪具有真正生物學(xué)差異的細(xì)胞,因此在刪除細(xì)胞時(shí)需要小心贿条，可與上述另外兩種方法連用。
如果你的數(shù)據(jù)量過(guò)大增热，使用Seurat時(shí)內(nèi)存不足整以，請(qǐng)看
實(shí)驗(yàn)記錄11：海量scRNA-seq數(shù)據(jù)的質(zhì)量控制、PCA钓葫、聚類
