植物空間轉(zhuǎn)錄組分析1：Seurat基本流程 - 簡(jiǎn)書 (jianshu.com)
植物空間轉(zhuǎn)錄組分析2：STEEL+Seurat - 簡(jiǎn)書 (jianshu.com)
植物空間轉(zhuǎn)錄組分析 3：hdWGCNA - 簡(jiǎn)書 (jianshu.com)
空間轉(zhuǎn)錄組是2022生命科學(xué)十大創(chuàng)新產(chǎn)品名單，因此將來會(huì)在生物學(xué)領(lǐng)域有非常大的應(yīng)用空間，目前植物類的相關(guān)文章較少荞雏，在本次的系列教程中，將使用復(fù)旦大學(xué)戚繼團(tuán)隊(duì)蘭花空間轉(zhuǎn)錄組的數(shù)據(jù)晓殊，希望大家一起學(xué)習(xí)，掌握植物空間轉(zhuǎn)錄組基本的分析方法(*^▽^*)
數(shù)據(jù)連接OSF | Spatiotemporal atlas of organogenesis in development of orchid flowers,與單細(xì)胞的數(shù)據(jù)結(jié)構(gòu)基本一致伤提，多了spatial這個(gè)文件夾巫俺，主要包含的就是切片信息和spot定位信息
1：前言
在本篇文章中，作者一共制作了蘭花空間轉(zhuǎn)錄組切片肿男，并使用了STEEL方法進(jìn)行聚類介汹，出于學(xué)習(xí)考慮，本次分析先使用Seurat的常見流程分析其中一個(gè)切片舶沛，后續(xù)推文中將使用文章中的STELL方法進(jìn)行聚類并使用Seurat的其他函數(shù)進(jìn)行后續(xù)分析
本次只分析Slide1嘹承，感興趣的可以試試其他切片
2：數(shù)據(jù)載入
和單細(xì)胞一樣，只是多了spatial這個(gè)文件夾需要輸入
##=============================1.安裝依賴包=====================================
##BiocManager::install(c('Seurat','ggplot2','cowplot','dplyr'))
##BiocManager::install("monocle",force = TRUE)
library(Seurat)
library(ggplot2)
library(cowplot)
library(dplyr)
library(magrittr)
library(gtools)
library(stringr)
library(Matrix)
library(tidyverse)
library(patchwork)
##=============================2.載入數(shù)據(jù)=======================================
orc1<- Load10X_Spatial("Slide_1",filename = "filtered_feature_bc_matrix.h5",assay = "spatial")
3: 質(zhì)量控制
不只是能體現(xiàn)小提琴圖如庭，還能將每個(gè)spot的測(cè)序質(zhì)量體現(xiàn)在切片上赶撰，幫我們確定哪些組織可能會(huì)出現(xiàn)問題
##=============================3.質(zhì)量控制=======================================
dir.create("QC")
##UMI統(tǒng)計(jì)畫圖
plot1 <- VlnPlot(orc1, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(orc1, features = "nCount_Spatial") + theme(legend.position = "right")
pearplot <- plot_grid(plot1,plot2)
ggsave("QC/Slide1_UMI.pdf", plot = pearplot, width = 12, height = 5)
##Feature/Gene統(tǒng)計(jì)畫圖
plot1 <- VlnPlot(orc1, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(orc1, features = "nFeature_Spatial") + theme(legend.position = "right")
pearplot <- plot_grid(plot1,plot2)
ggsave("QC/Slide1_Feature.pdf", plot = pearplot, width = 12, height = 5)
可以明顯的看到，不管是UMI還是Feature柱彻，在spot中差異主要是有與切片組織的細(xì)胞數(shù)目以及細(xì)胞類型導(dǎo)致的。比如餐胀，在花瓣中可以看到非常少的UMI以及Feature哟楷，在中間的分生組織中，UMI和Feature的數(shù)目非常多否灾。因此卖擅，單細(xì)胞的標(biāo)準(zhǔn)方法(如LogNormalize函數(shù))可能會(huì)有問題，在空間轉(zhuǎn)錄組分析中墨技，一般使用其他方法進(jìn)行標(biāo)準(zhǔn)化惩阶。
4. 數(shù)據(jù)標(biāo)準(zhǔn)化
目前空間轉(zhuǎn)錄組數(shù)據(jù)推薦使用sctransform，具體的原理我還沒有看懂扣汪，所以先根據(jù)流程走一遍断楷，看看效果
為了探究標(biāo)準(zhǔn)化方法的不同,sctransform和log規(guī)范化結(jié)果如何與UMIs的數(shù)量相關(guān)。為了進(jìn)行比較崭别，首先運(yùn)行sctransform來存儲(chǔ)所有基因的值冬筒，并通過NormalizeData運(yùn)行一個(gè)log規(guī)范化過程恐锣。
##============================4.數(shù)據(jù)標(biāo)準(zhǔn)化======================================
## 目前空間轉(zhuǎn)錄組推薦使用SCTransform路操，集合了單細(xì)胞標(biāo)準(zhǔn)化的NormalizeData()用踩，ScaleData()博敬，F(xiàn)indVariableFeatures()谒主，
dir.create("Normalize")
orc1 <- SCTransform(orc1, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
## 比較SCTransform與log Normalization的區(qū)別
orc1 <- NormalizeData(orc1, verbose = FALSE, assay = "Spatial")
## 計(jì)算UMI與Feature的相關(guān)性店量，分為5組
orc1 <- GroupCorrelation(orc1, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
orc1 <- GroupCorrelation(orc1, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)
p1 <- GroupCorrelationPlot(orc1, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") +
theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(orc1, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") +
theme(plot.title = element_text(hjust = 0.5))
p3 <- plot_grid(p1, p2)
ggsave("Normalize/Normalization_cor_Slide1.pdf", plot = p3, width = 12, height = 5)
對(duì)于上面的箱形圖滥朱，計(jì)算每個(gè)特征(基因)與UMIs數(shù)量(這里的nCount_Spatial變量)的相關(guān)性鞍历。然后乌询，根據(jù)基因的平均表達(dá)將它們分組呀打，并生成這些相關(guān)性的箱形圖矢赁。
同時(shí)這里需要注意，SCT標(biāo)準(zhǔn)化后有個(gè)單獨(dú)的assay(SCT)聚磺，里面包含三個(gè)矩陣坯台，分別為count,data以及scale.data，log標(biāo)準(zhǔn)化也有一個(gè)assay(Spatial)瘫寝，里面也包含三個(gè)矩陣蜒蕾，但是scale.data的數(shù)據(jù)為0
標(biāo)準(zhǔn)化之后進(jìn)行高變基因的篩選，這部分沒什么可說的焕阿，與單細(xì)胞是一致的
top10 <- head(VariableFeatures(orc1), 10)
plot1 <- VariableFeaturePlot(orc1)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5)
ggsave("Normalize/VariableFeature_Slide1.pdf", plot = plot2, width = 9, height = 6)
5.數(shù)據(jù)降維與聚類
由于文章中使用了STEEL方法咪啡，所以這次的分析只是走走流程，之后會(huì)按照文章中的STEEL方法進(jìn)行測(cè)試
##============================5.數(shù)據(jù)降維與聚類==================================
dir.create("cluster")
## PCA降維并提取主成分
orc1 <- RunPCA(orc1, features = VariableFeatures(orc1_new))
plot1 <- DimPlot(orc1, reduction = "pca", group.by="orig.ident")
plot2 <- ElbowPlot(orc1, ndims=40, reduction="pca")
pearplot <- plot_grid(plot1,plot2)
ggsave("cluster/PCA_Slide1.pdf", plot = pearplot, width = 13, height = 6)
pdf("cluster/PC_heatmap_Slide1.pdf",height = 12,width = 12)
DimHeatmap(orc1, dims = 1:9, nfeatures=10, cells = 200, balanced = TRUE)
dev.off()
## 細(xì)胞聚類
## 此步利用 細(xì)胞-PC值 矩陣計(jì)算細(xì)胞之間的距離暮屡，
## 然后利用距離矩陣來聚類撤摸。其中有兩個(gè)參數(shù)需要人工選擇，
## 第一個(gè)是FindNeighbors()函數(shù)中的dims參數(shù)褒纲，需要指定哪些pc軸用于分析准夷，選擇依據(jù)是之前介紹的cluster/pca.png文件中的右圖。
## 第二個(gè)是FindClusters()函數(shù)中的resolution參數(shù)莺掠，需要指定一個(gè)數(shù)值衫嵌，用于決定clusters的相對(duì)數(shù)量，數(shù)值越大cluters越多彻秆。
orc1 <- FindNeighbors(orc1, dims = 1:20)
orc1 <- FindClusters(orc1, resolution = 0.9)
metadata <- orc1@meta.data
table(orc1@meta.data$seurat_clusters)
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster_Slide1.csv',row.names = F)
使用resolution = 0.9只生成了12個(gè)聚類楔绞，原文中有40個(gè)，大家在設(shè)置的時(shí)候可以把resolution調(diào)到2甚至3左右看看效果
## 非線性降維
## tsne
orc1 <- RunTSNE(orc1, dims = 1:20)
embed_tsne <- Embeddings(orc1, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne_Slide1.csv')
p1 <- DimPlot(orc1, reduction = "tsne", label = TRUE)
p2 <- SpatialDimPlot(orc1, label.size = 3, pt.size.factor = 1.3)
pearplot <- plot_grid(p1,p2)
ggsave("cluster/tsne_Slide1.pdf", plot = p2, width = 6, height = 6)
## UMAP
orc1 <- RunUMAP(orc1,dims = 1:20)
embed_umap <- Embeddings(orc1, 'umap')
write.csv(embed_umap,'cluster/embed_umap_Slide1.csv')
p1 <- DimPlot(orc1, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(orc1, label = TRUE, label.size = 3, pt.size.factor = 1.3)
pearplot <- plot_grid(p1,p2)
ggsave("cluster/umap_Slide1.pdf", plot = pearplot, width = 13, height = 6)
6. 部分基因展示
此處是我比較糾結(jié)的地方唇兑，因?yàn)橐婚_始得到的數(shù)據(jù)和原文中有點(diǎn)不符酒朵，后來發(fā)現(xiàn)可能的原因是文章中使用的展示數(shù)據(jù)是counts
## featureplot
dir.create("featureplot")
p1 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="counts",pt.size.factor = 1.6,min.cutoff = 0.1)
p2 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="data",pt.size.factor = 1.6,min.cutoff = 0.1)
p3 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="scale.data",pt.size.factor = 1.6,min.cutoff = 0.1)
pearplot <- plot_grid(p1,p2,p3,align = "v",axis = "b",ncol = 3)
ggsave("featureplot/SCTSlide1.pdf", plot = pearplot, width = 18, height = 7)
總結(jié)
目前先把前期的工作的進(jìn)行了一下，原文中進(jìn)行了兩個(gè)擬時(shí)分析扎附，但因?yàn)榫垲惒灰粯铀韵炔贿M(jìn)行分析蔫耽，研究完STEEL之后再進(jìn)行擬時(shí)分析。
總體來講與單細(xì)胞變化不大留夜，有優(yōu)勢(shì)也有劣勢(shì)吧针肥，對(duì)于細(xì)胞分型來說饼记，空間轉(zhuǎn)錄組有先天優(yōu)勢(shì)，可以直接根據(jù)切片信息判斷細(xì)胞類型慰枕，但是空間轉(zhuǎn)錄組最小的單位spot具则，不能代表單個(gè)細(xì)胞，分辨率可能還存在問題
轉(zhuǎn)載請(qǐng)注明>>>周小釗的博客, 打賞推薦博客
