單細胞數據分析神器——Seurat

yjt2004us 2020-02-28

展開全文

近年來，單細胞技術日益火熱，并且有著愈演愈烈的趨勢。在2015年至2017年，甚至對某細胞群體或組織進行單細胞測序，解析其細胞成分就能發(fā)一篇CNS級別的文章。近兩三年，單細胞技術從最開始的基因組，轉錄組測序，發(fā)展成現在的單細胞DNA甲基化，單細胞ATAC-seq等等。測序手段也從早期的10X Genomics、 Drop-seq等，發(fā)展為現在的多種多樣個性化的方法。研究內容更不僅僅局限于解析細胞群體的成分，而是向研究細胞功能和生物學特性發(fā)展。今天小編向大家簡單一個實用并且易上手的單細胞數據分析軟件——Seurat,大家躺在床上為國家做貢獻的同時也能get新技能。

介紹一下今天的主角，Seurat是由New York Genome Center, Satija Lab開發(fā)的單細胞數據分析集成軟件包。其功能不僅包含基本的數據分析流程，如質控，細胞篩選，細胞類型鑒定，特征基因選擇，差異表達分析，數據可視化等等。同時也包括一些高級功能，如時序單細胞數據分析，不同組學單細胞數據整合分析等。今天，小編以官網中提供的單細胞基因表達數據為例，為大家簡單介紹一下Seurat軟件包中的基礎分析流程，希望能夠拋磚引玉，祝大家在科研的道路上越走越遠。

第一步，數據集導入

在本教程中，我們將分析從10X基因組學免費獲得的外周血單個核細胞（PBMC）數據集，來源于Illumina NextSeq 500測得的2700個單細胞轉錄組數據。首先，我們需要把數據集存儲成Seurat可識別的數據格式，

讀入的數據可以是一個矩陣，行代表基因，列代表細胞。

library(dplyr)library(Seurat)# Load the PBMC datasetpbmc.data <- Read10X(data.dir = '../data/pbmc3k/filtered_gene_bc_matrices/hg19/')# Initialize the Seurat object with the raw (non-normalized data).pbmc <- CreateSeuratObject(counts = pbmc.data, project = 'pbmc3k', min.cells = 3, min.features = 200)pbmc## An object of class Seurat ## 13714 features across 2700 samples within 1 assay## Active assay: RNA (13714 features)

數據導入成功以后，我們可以看到pbmc對象中包含了一個13714（基因數）X 2700（細胞數）的矩陣，其實在數據導入的時候，數據集中測到的少于200個基因的細胞（min.features = 200），和少于3個細胞覆蓋的基因（min.cells = 3），就已經被過濾掉了。

第二步，數據質控

質控的參數主要有兩個：1.每個細胞測到的unique feature數目（unique feature代表一個細胞檢測到的基因的數目，可以根據數據的質量進行調整）2.每個細胞檢測到的線粒體基因的比例，理論上線粒體基因組與核基因組相比，只占很小一部分。所以線粒體基因表達比例過高的細胞會被過濾。

在質控前，我們需要目測一下數據集的基本情況。

# The [[ operator can add columns to object metadata. This is a great place to stash QC statspbmc[['percent.mt']] <- PercentageFeatureSet(pbmc, pattern = '^MT-')# Visualize QC metrics as a violin plotVlnPlot(pbmc, features = c('nFeature_RNA', 'nCount_RNA', 'percent.mt'), ncol = 3)

代碼運行后會得到下圖所示結果，nFeature_RNA代表每個細胞測到的基因數目，nCount代表每個細胞測到所有基因的表達量之和，percent.mt代表測到的線粒體基因的比例。

這里，我們過濾掉上圖所示的離群點，并對基因的表達量進行l(wèi)og轉換。

第三步，鑒定細胞間表達量高變的基因（feature selection）

這一步的目的是鑒定出細胞與細胞之間表達量相差很大的基因，用于后續(xù)鑒定細胞類型，我們使用默認參數，即“vst”方法選取2000個高變基因。

pbmc <- FindVariableFeatures(pbmc, selection.method = 'vst', nfeatures = 2000)# Identify the 10 most highly variable genestop10 <- head(VariableFeatures(pbmc), 10)# plot variable features with and without labelsplot1 <- VariableFeaturePlot(pbmc)plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)CombinePlots(plots = list(plot1, plot2))

代碼運行后會得到下圖所示的結果。

第四步，細胞分類

1) 分類前首先要對數據集進行降維

#Scaling the dataall.genes <- rownames(pbmc)pbmc <- ScaleData(pbmc, features = all.genes)#Perform linear dimensional reductionpbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))# Examine and visualize PCA results a few different waysprint(pbmc[['pca']], dims = 1:5, 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 ## PC_ 4 ## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 ## Negative: VIM, IL7R, S100A6, IL32, S100A8 ## PC_ 5 ## Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY ## Negative: LTB, IL7R, CKB, VIM, MS4A7

官網中提供了不同的方法可視化降維的結果，此處不再贅述。感興趣的讀者可以自行探索降維的結果，代碼如下。

VizDimLoadings(pbmc, dims = 1:2, reduction = 'pca')DimPlot(pbmc, reduction = 'pca')DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)

2) 定義數據集的“維度”

這個小步驟就比較關鍵了，我們需要選擇出主成分的數目，用于后續(xù)細胞分類。值得注意的是，這里定義的“維度”并不代表細胞類型的數目，而是對細胞分類時需要用到的一個參數。

# NOTE: This process can take a long time for big datasets, comment out for expediency. More# approximate techniques such as those implemented in ElbowPlot() can be used to reduce# computation timepbmc <- JackStraw(pbmc, num.replicate = 100)pbmc <- ScoreJackStraw(pbmc, dims = 1:20)JackStrawPlot(pbmc, dims = 1:15)ElbowPlot(pbmc)

JackStraw（左圖）和Elbow（右圖）都可以決定數據的“維度”。但是Elbow比較直觀，我們選擇Elbow結果進行解讀?？梢钥吹剑鞒煞郑≒C）7到10之間，數據的標準差基本不在下降。所以我們需要在7到10之間進行選擇，為了尊重官網的建議，我們選取10，即前10個主成分用于細胞的分類。

3) 細胞分類

pbmc <- FindNeighbors(pbmc, dims = 1:10)pbmc <- FindClusters(pbmc, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck## ## Number of nodes: 2638## Number of edges: 96033## ## Running Louvain algorithm...## Maximum modularity in 10 random starts: 0.8720## Number of communities: 9## Elapsed time: 0 seconds

這里我們設置了dims = 1:10 即選取前10個主成分來分類細胞。分類的結果如下,可以看到，細胞被分為9個類別。

# Look at cluster IDs of the first 5 cellshead(Idents(pbmc), 5)## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG## 1 3 1 2 6## Levels: 0 1 2 3 4 5 6 7 8

4) 可視化分類結果

TSNE和UMAP兩種方法經常被用于可視化細胞類別。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')pbmc <- RunUMAP(pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clustersDimPlot(pbmc, reduction = 'umap')
saveRDS(pbmc, file = '../output/pbmc_tutorial.rds')

第五步，提取各個細胞類型的marker gene

利用 FindMarkers 命令，可以找到找到各個細胞類型中于其他類別的差異表達基因，作為該細胞類型的生物學標記基因。其中ident.1參數設置待分析的細胞類別，min.pct表示該基因表達數目占該類細胞總數的比例

# find all markers of cluster 1cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)head(cluster1.markers, n = 5)

利用 DoHeatmap 命令可以可視化marker基因的表達

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)DoHeatmap(pbmc, features = top10$gene) + NoLegend()

第六步，探索感興趣的基因

Seurat提供了小提琴圖和散點圖兩種方法，使我們能夠方便的探索感興趣的基因在各個細胞類型中的表達情況

VlnPlot(pbmc, features = c('MS4A1', 'CD79A'))

我們能夠看到，MS4A1和CD79A兩個基因在細胞群體3中特異性表達。

FeaturePlot(pbmc, features = c('MS4A1', 'GNLY', 'CD3E', 'CD14', 'FCER1A', 'FCGR3A', 'LYZ', 'PPBP', 'CD8A'))

這種展示方法把基因表達量映射到UMAP結果中，同樣可以直觀的看到基因表達的特異性。

第七步，利用先驗知識定義細胞類型

通過對比我們鑒定的marker gene與已發(fā)表的細胞類型特意的基因表達marker，可以定義我們劃分出來的細胞類群。

最后，給我們定義好的細胞類群加上名稱

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()