|
1. 將Cellranger中的基因表達(dá)矩陣filtered_gene_bc_matrices用于分析。 2. 進(jìn)行質(zhì)量控制（QC），以刪除異常細(xì)胞； 3. 標(biāo)準(zhǔn)化與歸一化，消除技術(shù)噪音與批次效應(yīng)； 4. 主成分分析（PCA）與挑選 5. t-SNE聚類
參考網(wǎng)站：https:///seurat/pbmc3k_tutorial.html
Seurat的安裝：R中運(yùn)行install.packages('Seurat') 經(jīng)過Cellranger的數(shù)據(jù)整理之后，得到：
Seurat是一種R包，設(shè)計用于QC，分析和探索單細(xì)胞RNA-seq數(shù)據(jù)。 Seurat旨在使用戶能夠從單細(xì)胞轉(zhuǎn)錄組測量中識別和解釋異質(zhì)性來源，并整合不同類型的單細(xì)胞數(shù)據(jù)。 運(yùn)行R，并且加載這兩個包 library(Seurat)
library(dplyr)
讀取數(shù)據(jù) spleen.data <- Read10X(data.dir = '/data1/zll/project/deepBase3/HCA/filtered_gene_bc_matrices/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對象： 命令CreateSeuratObject
輸入數(shù)據(jù)spleen.data
留下所有在>=3個細(xì)胞中表達(dá)的基因min.cells = 3；
留下所有檢測到>=200個基因的細(xì)胞min.genes = 200。
(為了除去一些) 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 個細(xì)胞 以下步驟包括Seurat中scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理工作流程。這些代表了Seurat對象的創(chuàng)建，基于QC指標(biāo)的細(xì)胞選擇和過濾，數(shù)據(jù)標(biāo)準(zhǔn)化和縮放，以及高度可變基因的檢測。
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)
VlnPlot_of_spleen.png > par(mfrow = c(1, 2))
> GenePlot(object = spleen, gene1 = 'nUMI', gene2 = 'percent.mito')
> GenePlot(object = spleen, gene1 = 'nUMI', gene2 = 'nGene')
過濾細(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))
查看過濾掉剩下多少細(xì)胞： spleen
An object of class seurat in project 10X_spleen
15655 genes across 1940 samples.
剩下15655個基因，1940個細(xì)胞。 數(shù)據(jù)標(biāo)準(zhǔn)化加個log： spleen <- NormalizeData(object=spleen, normalization.method = 'LogNormalize', scale.factor = 10000) Performing log-normalization0% 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 means0% 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
您的單細(xì)胞數(shù)據(jù)集可能包含“不感興趣”的變異來源。這不僅包括技術(shù)噪音，還包括批次效應(yīng)，甚至包括生物變異來源（細(xì)胞周期階段）。正如(Buettner, et al NBT，2015)中所建議的那樣，從分析中回歸這些信號可以改善下游維數(shù)減少和聚類。為了減輕這些信號的影響，Seurat構(gòu)建線性模型以基于用戶定義的變量預(yù)測基因表達(dá)。這些模型的縮放得分殘差存儲在Scale.data槽中，用于降維和聚類。
我們可以消除由批次（如果適用）驅(qū)動的基因表達(dá)中的細(xì)胞 - 細(xì)胞變異，細(xì)胞比對率（由Drop-seq數(shù)據(jù)的Drop-seq工具提供），檢測到的分子數(shù)量和線粒體基因表達(dá)。對于循環(huán)細(xì)胞，我們還可以學(xué)習(xí)“細(xì)胞周期”評分（參見此處的示例）并對其進(jìn)行回歸。在這個有絲分裂后血細(xì)胞的簡單例子中，我們回歸了每個細(xì)胞檢測到的分子數(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%
主成分分析是什么？
主成分分析，是考察多個變量間相關(guān)性一種多元統(tǒng)計方法，研究如何通過少數(shù)幾個主成分來揭示多個變量間的內(nèi)部結(jié)構(gòu)，即從原始變量中導(dǎo)出少數(shù)幾個主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關(guān).通常數(shù)學(xué)上的處理就是將原來P個指標(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] 'FCN1''LYZ''SERPINA1''CSTA''RP11-1143G9.4' [1] 'PC2' [1] 'CD79A''MS4A1''VPREB3''CD79B''HLA-DQB1' [1] 'NKG7' 'CST7' 'GZMA' 'CD7' 'CCL5' [1] 'PC3' [1] 'TRDC' 'KLRF1' 'MS4A1' 'CD79B' 'IRF8' [1] 'IL7R' 'TRAC' 'CD3D' 'CD2' 'CD3G' [1] 'PC4' [1] 'GIMAP7' 'GZMB''FGFBP2' 'SPON2''PRF1' [1] 'BAG3''HSPD1''FKBP4''DNAJA1''ZFAND2A' [1] 'PC5'[1] 'UBE2C' 'TYMS''MKI67' 'TOP2A' 'AURKB' [1] 'FCGR3A' 'FGFBP2' 'SPON2' 'GNLY''GZMB' 選擇了前10個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 spleen <- RunTSNE(spleen, dims.use = 1:10, do.fast= TRUE)
TSNEPlot(spleen) TSNE.jpeg > saveRDS(spleen, file = '/Users/shinianyike/Desktop/zll/Seurat/spleen_results/spleen_1.rds')
將R變量保存，利于后續(xù)的分析。 過濾低質(zhì)量細(xì)胞：
在 scRNA-seq 分析中，有些細(xì)胞質(zhì)量比較低，比如細(xì)胞處于凋亡狀態(tài),細(xì)胞中 RNA 發(fā)生降解等,這些細(xì)胞的存在會影響分析，因此我們第一步需要對細(xì)胞進(jìn)行過濾。主要可分為三類:
①利用細(xì)胞檢測到的基因數(shù)或者是 reads 比對率來判斷技術(shù)噪音。
但不管是基因檢測數(shù)目還是比對率都跟實驗方法有很大相關(guān)性。 如果比對率太低,表明 RNA 可能發(fā)生了降解,或者文庫有污染或者細(xì)胞裂解不完全。 ②如果實驗中加入了 spike-ins（本實驗沒有），可以通過計算比對到內(nèi)源性 RNA 和外源性 RNA(spike-ins)的 reads 比例來過濾低質(zhì)量細(xì)胞。 比值偏低表明細(xì)胞中的 RNA 數(shù)量較低，細(xì)胞可丟棄。但是也需要注意其實當(dāng)細(xì)胞狀態(tài)不一樣，比如處于不同細(xì)胞周期時，細(xì)胞的 RNA 數(shù)量是具有很大差異的。不過我們依然認(rèn)為在一大群細(xì)胞中，spike-ins比例特別高的細(xì)胞在很大概率上應(yīng)該被排除在外。軟件 SinQC (Single-cell RNA-seq Quality Control)可以根據(jù)比對率和檢測到的基因數(shù)來過濾細(xì)胞。 ③根據(jù)整體的基因表達(dá)譜來定義技術(shù)噪音。 比如對細(xì)胞進(jìn)行聚類分析，PCA 分析等，將 outlier 細(xì)胞刪除，或者細(xì)胞表達(dá)中位值低于某一設(shè)定閾值時將該細(xì)胞過濾掉。當(dāng)然這種方法也存在誤刪具有真正生物學(xué)差異的細(xì)胞,因此在刪除細(xì)胞時需要小心，可與上述另外兩種方法連用。
作者：MC學(xué)公衛(wèi)
鏈接：https://www.jianshu.com/p/866a2f0097fe
來源：簡書
|
