前面我們有介紹了利用10x Space Ranger軟件分析空間轉(zhuǎn)錄組原始數(shù)據(jù)得到可用于下游分析的矩陣和鏡像文件甘有。今天來介紹一下怎么利用Space Ranger的結(jié)果文件進(jìn)行后續(xù)分析诉儒，這里主要使用Seurat在進(jìn)行下游分析。
先來回顧一下跑完Space Ranger得到哪些結(jié)果文件：
Outputs:
- Run summary HTML: /opt/sample345/outs/web_summary.html
- Outputs of spatial pipeline: /opt/sample345/outs/spatial
- Run summary CSV: /opt/sample345/outs/metrics_summary.csv
- BAM: /opt/sample345/outs/possorted_genome_bam.bam
- BAM index: /opt/sample345/outs/possorted_genome_bam.bam.bai
- Filtered feature-barcode matrices MEX: /opt/sample345/outs/filtered_feature_bc_matrix
- Filtered feature-barcode matrices HDF5: /opt/sample345/outs/filtered_feature_bc_matrix.h5
- Unfiltered feature-barcode matrices MEX: /opt/sample345/outs/raw_feature_bc_matrix
- Unfiltered feature-barcode matrices HDF5: /opt/sample345/outs/raw_feature_bc_matrix.h5
- Secondary analysis output CSV: /opt/sample345/outs/analysis
- Per-molecule read information: /opt/sample345/outs/molecule_info.h5
- Loupe Browser file: /opt/sample345/outs/cloupe.cloupe
用seurat進(jìn)行下游分析主要用到兩個(gè)結(jié)果文件亏掀。一個(gè)是filtered_feature_bc_matrix.h5文件忱反，一個(gè)是spatial鏡像結(jié)果目錄。
安裝R包
由于seurat分析空間轉(zhuǎn)錄組的R包satijalab-seurat是在GitHub 上的滤愕，如果我們需要直接安裝的話需要先安裝R包devtools温算，然后利用devtools工具中的install_github來安裝GitHub上的R包。
安裝devtools
install.packages('devtools')
安裝satijalab-seurat
devtools::install_github("satijalab/seurat",ref = "spatial")
考慮到直接安裝github上的R包速度是很慢的间影，非趁渍撸考驗(yàn)網(wǎng)速，可能需要多次才能安裝成功宇智，我們也可以直接下載安裝包蔓搞，本地安裝。
#下載
https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial
#安裝
install.packages("satijalab-seurat-v3.1.5-351-g85610bc.tar.gz", repos = NULL, type = "source")
注意該R包還在開發(fā)中随橘，不要和之前安裝的seurat包沖突喂分。
數(shù)據(jù)準(zhǔn)備
這里使用從10x官網(wǎng)下載的小鼠腦組織樣本MouseBrain Serial Section 1 (Sagittal-Posterior)。
下載網(wǎng)址：https://support.10xgenomics.com/spatial-gene-expression/datasets
點(diǎn)擊選擇的樣本机蔗，下載兩個(gè)數(shù)據(jù)就行
cell matrix HDF5 (filtered)和Spatial imaging data
導(dǎo)入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的輸出信息我們可以知道蒲祈，這個(gè)樣本包含3355個(gè)spot點(diǎn)、32285個(gè)基因萝嘁。
基礎(chǔ)統(tǒng)計(jì)作圖
##UMI統(tǒng)計(jì)畫圖
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)
UMI數(shù)大部分集中到10000-20000區(qū)間梆掸，最高不超過80000，并且組織中高UMI數(shù)的區(qū)域主要集中在左下角牙言。后面可以特意性關(guān)注一下左下角區(qū)域的基因的表達(dá)和主要的細(xì)胞類型酸钦。
##gene數(shù)目統(tǒng)計(jì)畫圖
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)
基因數(shù)目大部分處于2500-7500之間，結(jié)合UMI數(shù)據(jù)的分布可以發(fā)現(xiàn)UMI數(shù)目高的區(qū)域基因數(shù)也高咱枉，說明基因數(shù)和UMI數(shù)基本上是呈正相關(guān)的卑硫。
#線粒體統(tǒng)計(jì)
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[-]
總體來說，這個(gè)樣本的線粒體比例不高蚕断，左邊中上區(qū)域有一處線粒體比例稍微高一點(diǎn)欢伏，后面也可以仔細(xì)研究一下這一塊區(qū)域到底是特定的細(xì)胞類型引起的還是組織活性的差異引起的。不過從這張圖我們還可以發(fā)現(xiàn)一個(gè)有意思的現(xiàn)象亿乳，基因和UMI高表達(dá)的區(qū)域往往線粒體比例更低硝拧。
數(shù)據(jù)過濾
做單細(xì)胞RNAseq我們都會(huì)根據(jù)UMI、基因數(shù)、線粒體比例等進(jìn)行過濾障陶，那么做空間轉(zhuǎn)錄組數(shù)據(jù)分析其實(shí)我們也可以按這樣的方式來過濾滋恬。具體的過濾條件需要根據(jù)具體樣本數(shù)據(jù)來定，沒有固定的標(biāo)準(zhǔn)咸这。
比如這個(gè)樣本我們可以設(shè)置過濾條件：
① 基因數(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個(gè)spot點(diǎn)媳维。過濾后我們?cè)诶L制一下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)
那么現(xiàn)在問題來了遏暴，過濾之后組織圖像里面缺了幾塊侄刽，顯得特別丑。那么我們到底應(yīng)不應(yīng)該過濾呢朋凉？過濾數(shù)據(jù)可以減少利群的點(diǎn)州丹，減少對(duì)后面聚類結(jié)果的影響，不過濾數(shù)據(jù)可以讓組織圖像保持完整性杂彭，繪圖更好看一點(diǎn)墓毒，所以這個(gè)還真不好決斷。
數(shù)據(jù)歸一化
Seurat官方推薦使用sctransform 歸一化方法亲怠，它構(gòu)建了基因表達(dá)的正則化負(fù)二項(xiàng)模型所计，以便在保留生物差異的同時(shí)考慮技術(shù)因素。
Sctransform函數(shù)同時(shí)實(shí)現(xiàn)了NormalizeData团秽、ScaleData主胧、FindVariableFeatures三個(gè)函數(shù)的功能。
mydata <- SCTransform(mydata, assay = "Spatial", verbose = FALSE)
基因表達(dá)可視化
Seurat的SpatialFeaturePlot功能擴(kuò)展了FeaturePlot习勤，可以將表達(dá)數(shù)據(jù)覆蓋在組織組織上踪栋。這里展示的Hpca基因是一個(gè)強(qiáng)的海馬marker ，Ttr是一個(gè)脈絡(luò)叢marker 图毕∫亩迹可以通過基因的表達(dá)分布來初步判斷一下海馬區(qū)和脈絡(luò)叢區(qū)處于組織切片的哪個(gè)位置。
從結(jié)果的展示來看予颤，這兩個(gè)marker基因的分布還是挺集中的损肛，這也說明理由空間轉(zhuǎn)錄組數(shù)據(jù)來分析小鼠腦的不同區(qū)域的表達(dá)差異應(yīng)該還是比較準(zhǔn)確的。另外荣瑟，海馬區(qū)的分布可以大概分成3大塊治拿，從上之下第一塊弧形區(qū)域似乎處于線粒體高表達(dá)區(qū)域，而最下面一塊弧形區(qū)處于基因高表達(dá)區(qū)笆焰。后面可以把這三個(gè)不同區(qū)域的數(shù)據(jù)進(jìn)行差異基因和功能的比較也許會(huì)發(fā)現(xiàn)一些有意思的東西劫谅。
降維、聚類和可視化
先進(jìn)行PCA降維，再選擇前30個(gè)維度進(jìn)行聚類和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)
mydata <- RunTSNE(mydata, reduction = "pca",dims = 1:30)
tsne展示結(jié)果：
Umap展示結(jié)果：
tsne和umap兩種展示方式在這次分析里差別不是特別大，tsne相對(duì)來說群于群之間分的更開贯城，而umap則單個(gè)亞群位置更集中熊楼。這個(gè)時(shí)候我們也可以結(jié)合前面marker基因的表達(dá)分布圖來大概判斷一下每個(gè)亞群大概處于小鼠腦的那個(gè)區(qū)。
由于亞群的顏色比較接近能犯，有時(shí)候不太好判斷鲫骗，我們可以是cells.highlight來標(biāo)記特定的亞群。
SpatialDimPlot(mydata, cells.highlight = CellsByIdentities(object = mydata, idents = c(1, 2, 3, 4,
5, 6)), facet.highlight = TRUE, ncol = 3)
當(dāng)然也可以利用LinkedDimPlot和LinkedFeaturePlot函數(shù)進(jìn)行交互式可視化踩晶，這里不再展示执泰。
掃碼關(guān)注公眾號(hào)，內(nèi)容更精彩哦渡蜻！
