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

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

第一步，數據集導入

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

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

library(dplyr)library(Seurat)# Load the PBMC datasetpbmc.data # Initialize the Seurat object with the raw (non-normalized data).pbmc 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"]] # Visualize QC metrics as a violin plotVlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

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

這里，我們過濾掉上圖所示的離群點，并對基因的表達量進行log轉換。

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

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

pbmc # Identify the 10 most highly variable genestop10 # plot variable features with and without labelsplot1 plot2 CombinePlots(plots = list(plot1, plot2))

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

第四步，細胞分類

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

#Scaling the dataall.genes pbmc all.genes)#Perform linear dimensional reductionpbmc # 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) 定義數據集的“維度”

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

# 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 pbmc JackStrawPlot(pbmc, dims = 1:15)ElbowPlot(pbmc)

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

3)?細胞分類

pbmc pbmc ## 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 # 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 head(cluster1.markers, n = 5)

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

top10 % 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與已發表的細胞類型特意的基因表達marker，可以定義我們劃分出來的細胞類群。

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

new.cluster.ids "NK", "DC", "Platelet")names(new.cluster.ids) pbmc DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

以上就是小編學習到的基本的單細胞轉錄組數據分析流程，整理出來供大家參考，本文主要參考了Seurat官網給出的單細胞轉錄組數據分析教程，若文中有小編解讀錯誤之處，還請廣大讀者批評指正。

Seurat官網地址

https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

?

歡迎關注生信人轉錄組?|?甲基化?|?重測序?|?單細胞?|?m6A|多組學?cytoscape|?limma|?WGCNA|水熊蟲傳奇|linux電泳?|?PCR?|?測序簡史?|?核型?|?NIPT?|?基礎實驗?基因|?2019-nCoV|?富集分析|?聯合分析?|微環境瘟疫追兇|?思路匯總|?學者|?科研?|?撤稿?|?讀博|基因

總結

以上是生活随笔為你收集整理的seurat提取表达矩阵_单细胞数据分析神器——Seurat的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。