Seurat 分析Visium空間轉(zhuǎn)錄組

一、Seurat v3.2 對空間轉(zhuǎn)錄組Visium的結(jié)果分析中實現(xiàn)的功能：

  • 歸一化
  • 降維和聚類
  • 檢測空間可變特征(spatially-variable features)
  • 互動式的可視化
  • 與單細胞RNA-seq數(shù)據(jù)整合
  • 處理多個切片

二、安裝Seurat v3.2 ：

#For learning Seurat v2.3
#Author:Robin 2019.12.22
# Enter commands in R (or R studio, if installed)

# Install the devtools package from Hadley Wickham
install.packages("devtools")

#devtools::install_github("satijalab/seurat", ref = "spatial")
devtools::install_github('satijalab/seurat-data')

成功安裝！

> library(Seurat)
> packageVersion("Seurat")
[1] ‘3.2.2’

加載包：

library(Seurat)
library(SeuratData)
library(ggplot2)
library(cowplot)
library(dplyr)

三、下載數(shù)據(jù)集

  • 從10X Visium流程中生成結(jié)果文件的方法：
    先來回顧一下跑完Space Ranger得到哪些結(jié)果文件：
├── analysis
│   ├── clustering
│   ├── diffexp
│   ├── pca
│   ├── tsne
│   └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5  #實際上需要的文件①
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
├── spatial  # 實際上需要的文件②；記錄空間信息:用戶提供的原始全分辨率brightfield圖像的下采樣版本。下采樣是通過box濾波實現(xiàn)的，它對全分辨率圖像中像素塊的RGB值進行平均，得到下采樣圖像中一個像素點的RGB值。
│   ├── aligned_fiducials.jpg 
│   ├── detected_tissue_image.jpg
│   ├── scalefactors_json.json
│   ├── tissue_hires_image.png
│   ├── tissue_lowres_image.png
│   └── tissue_positions_list.csv
└── web_summary.html

9 directories, 20 files

用seurat進行下游分析主要用到兩個結(jié)果文件。一個是filtered_feature_bc_matrix.h5文件，一個是spatial鏡像結(jié)果目錄。

  • 從10x官網(wǎng)下載測試數(shù)據(jù)

這里使用從10x官網(wǎng)下載的小鼠腦組織樣本MouseBrain Serial Section 1 (Sagittal-Posterior)。

下載網(wǎng)址：https://support.10xgenomics.com/spatial-gene-expression/datasets

點擊選擇的樣本，下載兩個數(shù)據(jù)就行
cell matrix HDF5 (filtered)和Spatial imaging data

image.png

導入R包讀取文件

library("Seurat")
library("ggplot2")
library("cowplot")
library("dplyr")
library("hdf5r")
##讀取矩陣文件
name='Posterior1'
expr <- "/data/Sagittal-Posterior1/V1_Mouse_Brain_Sagittal_Posterior_filtered_feature_bc_matrix.h5"
expr.mydata <- Seurat::Read10X_h5(filename =  expr )
mydata <- Seurat::CreateSeuratObject(counts = expr.mydata, project = 'Posterior1', assay = 'Spatial')
mydata$slice <- 1
mydata$region <- 'Posterior1' #命名
#讀取鏡像文件
imgpath <- "/data/Sagittal-Posterior1/spatial"
img <- Seurat::Read10X_Image(image.dir = imgpath)
Seurat::DefaultAssay(object = img) <- 'Spatial'
img <- img[colnames(x = mydata)]
mydata[['image']] <- img
mydata  #查看數(shù)據(jù)
An object of class Seurat
32285 features across 3355 samples within 1 assay
Active assay: Spatial (32285 features)

從mydata的輸出信息我們可以知道，這個樣本包含3355個spot點、32285個基因。

也可以直接加載自己的數(shù)據(jù)：

brain <-  Load10X_Spatial(data.dir = a,
                  filename = 'V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.h5',
                  assay = "Spatial",
                  slice = 'spatial',
                  filter.matrix = T)

四、基礎統(tǒng)計作圖 

##UMI統(tǒng)計畫圖
plot1 <- VlnPlot(mydata, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "nCount_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)
image.png

UMI數(shù)大部分集中到10000-20000區(qū)間，最高不超過80000，并且組織中高UMI數(shù)的區(qū)域主要集中在左下角。后面可以特意性關注一下左下角區(qū)域的基因的表達和主要的細胞類型。

##gene數(shù)目統(tǒng)計畫圖
plot1 <- VlnPlot(mydata, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "nFeature_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)
image.png

基因數(shù)目大部分處于2500-7500之間，結(jié)合UMI數(shù)據(jù)的分布可以發(fā)現(xiàn)UMI數(shù)目高的區(qū)域基因數(shù)也高，說明基因數(shù)和UMI數(shù)基本上是呈正相關的。

#線粒體統(tǒng)計
mydata[["percent.mt"]] <- PercentageFeatureSet(mydata, pattern = "^mt[-]")
plot1 <- VlnPlot(mydata, features = "percent.mt", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "percent.mt") + theme(legend.position = "right")
plot_grid(plot1, plot2)

注意如果是人的數(shù)據(jù)pattern = "^mt[-]改成pattern = "^MT[-]

image.png

總體來說，這個樣本的線粒體比例不高，左邊中上區(qū)域有一處線粒體比例稍微高一點，后面也可以仔細研究一下這一塊區(qū)域到底是特定的細胞類型引起的還是組織活性的差異引起的。不過從這張圖我們還可以發(fā)現(xiàn)一個有意思的現(xiàn)象，基因和UMI高表達的區(qū)域往往線粒體比例更低。

五、數(shù)據(jù)過濾

做單細胞RNAseq我們都會根據(jù)UMI、基因數(shù)、線粒體比例等進行過濾，那么做空間轉(zhuǎn)錄組數(shù)據(jù)分析其實我們也可以按這樣的方式來過濾。具體的過濾條件需要根據(jù)具體樣本數(shù)據(jù)來定，沒有固定的標準。
比如這個樣本我們可以設置過濾條件：
① 基因數(shù)大于200，小于7500
② UMI數(shù)大于1000，小于60000
③ 線粒體比例小于25%

mydata2 <- subset(mydata, subset = nFeature_Spatial > 200 & nFeature_Spatial <7500 & nCount_Spatial > 1000 & nCount_Spatial < 60000 & percent.mt < 25)
mydata2
An object of class Seurat
32285 features across 2977 samples within 1 assay
Active assay: Spatial (32285 features)

過濾后還剩2977個spot點。過濾后我們在繪制一下UMI分布圖。

plot1 <- VlnPlot(mydata2, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata2, features = "nCount_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)
image.png

那么現(xiàn)在問題來了，過濾之后組織圖像里面缺了幾塊，顯得特別丑。那么我們到底應不應該過濾呢？過濾數(shù)據(jù)可以減少利群的點，減少對后面聚類結(jié)果的影響，不過濾數(shù)據(jù)可以讓組織圖像保持完整性，繪圖更好看一點，所以這個還真不好決斷。

六、數(shù)據(jù)歸一化

Seurat官方推薦使用sctransform 歸一化方法，它構(gòu)建了基因表達的正則化負二項模型，以便在保留生物差異的同時考慮技術因素。
Sctransform函數(shù)同時實現(xiàn)了NormalizeData、ScaleData、FindVariableFeatures三個函數(shù)的功能。

有關sctransform的更多信息，請參見 here的預印和here的Seurat教程。sctransform將數(shù)據(jù)歸一化，檢測高方差特征high-variance features，并將數(shù)據(jù)存儲在SCTassay中。

mydata <- SCTransform(mydata, assay = "Spatial", verbose = FALSE)

七、基因表達可視化

Seurat的SpatialFeaturePlot功能擴展了FeaturePlot，可以將表達數(shù)據(jù)覆蓋在組織組織上。這里展示的Hpca基因是一個強的海馬marker ，Ttr是一個脈絡叢marker ?？梢酝ㄟ^基因的表達分布來初步判斷一下海馬區(qū)和脈絡叢區(qū)處于組織切片的哪個位置。

SpatialFeaturePlot(mydata, features = c("Hpca", "Ttr"))
image.png

從結(jié)果的展示來看，這兩個marker基因的分布還是挺集中的，這也說明理由空間轉(zhuǎn)錄組數(shù)據(jù)來分析小鼠腦的不同區(qū)域的表達差異應該還是比較準確的。另外，海馬區(qū)的分布可以大概分成3大塊，從上之下第一塊弧形區(qū)域似乎處于線粒體高表達區(qū)域，而最下面一塊弧形區(qū)處于基因高表達區(qū)。后面可以把這三個不同區(qū)域的數(shù)據(jù)進行差異基因和功能的比較也許會發(fā)現(xiàn)一些有意思的東西。

7.1 Seurat中的默認參數(shù)強調(diào)分子數(shù)據(jù)的可視化。

但是，您還可以通過更改以下參數(shù)來調(diào)整斑點的大?。捌渫该鞫龋?，以改善組織學圖像的可視化：

  • pt.size.factor-這將縮放斑點的大小。默認值為1.6
  • alpha-最小和最大透明度。默認值為c（1，1）。
  • 嘗試設置為alphac（0.1，1），以降低具有較低表達的點的透明度
p1 <- SpatialFeaturePlot(mydata, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(mydata, features = "Ttr", alpha = c(0.1, 1))
plot_grid(p1, p2)
image.png

八、降維、聚類和可視化

先進行PCA降維，再選擇前30個維度進行聚類和umap降維。

降維、聚類和可視化

mydata <- RunPCA(mydata, assay = "SCT", verbose = FALSE)
mydata <- FindNeighbors(mydata, reduction = "pca", dims = 1:30)
mydata <- FindClusters(mydata, verbose = FALSE)
mydata <- RunUMAP(mydata, reduction = "pca", dims = 1:30)

再然后，我們可以在UMAP空間(使用DimPlot)或使用SpatialDimPlot將聚類clustering的結(jié)果顯示在圖像上

p1 <- DimPlot(mydata, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(mydata, label = TRUE, label.size = 3)
plot_grid(p1, p2)
image.png

由于顏色太多，因此可視化哪個立體像素voxel屬于哪個cluster可能是一個挑戰(zhàn)。
解決方案：

  • 設置label參數(shù)會在每個cluster的中間處放置一個彩色框（請參見上圖），
  • 而SpatialDimPlot的do.hover參數(shù)允許交互式地查看當前點的標識。
# move your mouse
>SpatialDimPlot(mydata, do.hover = TRUE)
Error in SpatialDimPlot(mydata, do.hover = TRUE) : 
  參數(shù)沒有用(do.hover = TRUE)

操作失敗


image

由于亞群的顏色比較接近，有時候不太好判斷，我們可以是cells.highlight來標記特定的亞群。

SpatialDimPlot(mydata, cells.highlight = CellsByIdentities(object = mydata, idents = c(1, 2, 3, 4, 
   5, 6)), facet.highlight = TRUE, ncol = 3)
image.png

此外，LinkedDimPlot和LinkedFeaturePlot函數(shù)支持交互式可視化。這些圖將UMAP表示與組織圖像表示聯(lián)系起來，并允許交互選擇。例如，您可以在UMAP圖中選擇一個區(qū)域，圖像表示中相應的點將突出顯示。

LinkedDimPlot(brain)
image

tsne和umap兩種展示方式在這次分析里差別不是特別大，tsne相對來說群于群之間分的更開，而umap則單個亞群位置更集中。這個時候我們也可以結(jié)合前面marker基因的表達分布圖來大概判斷一下每個亞群大概處于小鼠腦的那個區(qū)。

九、亞群marker基因分析

同時分析所有亞群的marker基因

data.markers <- FindAllMarkers(mydata, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25,test.use = "wilcox")

參數(shù)解析：
logfc.threshold: logfc閾值，閾值以上的基因才輸出，默認0.25
only.pos： 是否只輸出正的marker，也就是只輸出在該亞群中高表達的marker，默認是FALSE
test.use： 算法選擇，indAllMarkers提供8種差異分析算法，默認是wilcox
min.pct: 最低表達比例，指基因在兩組細胞中至少有一組中有表達的細胞的比例等于或超過該閾值

查看每個亞群top3差異基因

topgene<-data.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
topgene
# A tibble: 48 x 7
# Groups:   cluster [16]
       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene 
       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
 1 2.04e-245     2.03  1     0.911 3.73e-241 0       Plp1 
 2 3.04e-238     1.69  1     0.829 5.55e-234 0       Trf  
 3 1.13e-232     1.66  1     0.998 2.07e-228 0       Mbp  
 4 3.47e- 83     0.670 1     0.998 6.34e- 79 1       Mbp  
 5 3.81e- 77     0.643 1     0.818 6.96e- 73 1       Mobp 
 6 5.76e- 60     0.653 0.98  0.686 1.05e- 55 1       Agt  
 7 1.50e-161     2.14  1     0.655 2.74e-157 2       Pcp2 
 8 3.85e-155     1.81  1     0.808 7.04e-151 2       Pvalb
 9 6.11e-155     1.70  1     0.856 1.12e-150 2       Itpr1
10 3.04e-163     1.85  0.944 0.261 5.56e-159 3       Igf2 
# ... with 38 more rows

也可以單獨兩個或多個亞群比較分析marker基因分析

markers <- FindMarkers(mydata, ident.1 = 2,, ident.2 = 6)

繪制每個亞群top3差異基因熱圖

DoHeatmap(mydata, features = topgene$gene,size = 2) + NoLegend()
image.png

繪制marker基因小提琴圖

VlnPlot(mydata,features = c('Plp1','Trf'), group.by = "seurat_clusters",pt.size = 0)
image.png

繪制marker基因tsne圖

FeaturePlot(mydata, features = c('Plp1','Trf'), cols = c("grey", "red"),reduction = "tsne")
image.png

參考文章

作者：鄧老師呦
鏈接：http://www.itdecent.cn/p/137376261b13
作者：Geekero
鏈接：http://www.itdecent.cn/p/07593e4d99a9
來源：簡書
著作權歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權，非商業(yè)轉(zhuǎn)載請注明出處。

最后編輯于 ：
?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容