【空间转录组】MIA分析

之前讲过一篇空间转录组的文献，里面首次提出了Multimodal intersection analysis(MIA)的空间转录组分析思路。

讲解视频在B站

MIA分析可以用来评估空间上某个region或者cluster中富集的细胞类型。需要单细胞和空间转录组两种组学数据，数据最好配对。

MIA原理

上图是示例图，一个region是否富含某一种细胞类型，看的是一个region高表达的基因和一个celltype高表达的基因是不是有足够多的重叠。

原理不算难，跟常规的「差异基因」做「富集分析」原理是一样的，详细原理之前写过一篇帖子【富集分析的原理与实现】，感兴趣可以看看。

下文的分析会复现这两张小图：原文的fig.2e和fig.2h

1.单细胞数据分析流程

library(Seurat)
library(tidyverse)

pdacA.list=readRDS("PDAC-A.list.rds")
# length(pdacA.list)
# 4
# names(pdacA.list)
# "scRNA_count" "scRNA_anno"  "ST_count"    "ST_coord"
# 该列表包含单细胞count矩阵，单细胞数据集的注释，空间转录组count矩阵，空间转录组spot的坐标
# 由公粽号【TOP生物信息】整理

### 0.根据原文献fig.2h，将小注释合并 ########################################
pdacA.list$scRNA_anno$celltype[str_detect(pdacA.list$scRNA_anno$celltype,"^Ductal")] = "Ductal"
pdacA.list$scRNA_anno$celltype[str_detect(pdacA.list$scRNA_anno$celltype,"^Macrophages")] = "Macrophages"
pdacA.list$scRNA_anno$celltype[str_detect(pdacA.list$scRNA_anno$celltype,"^mDCs")] = "mDCs"

### 1.单细胞转录组流程 #####################################################
pdacA.seu = CreateSeuratObject(counts = pdacA.list$scRNA_count)
pdacA.seu <- NormalizeData(pdacA.seu, normalization.method = "LogNormalize", scale.factor = 10000)
pdacA.seu <- FindVariableFeatures(pdacA.seu, selection.method = "vst", nfeatures = 2000)
pdacA.seu <- ScaleData(pdacA.seu)
pdacA.seu <- RunPCA(pdacA.seu, npcs = 50, verbose = FALSE)

### 添加注释信息
pdacA.seu@meta.data$CB=rownames(pdacA.seu@meta.data)
pdacA.seu@meta.data=pdacA.seu@meta.data%>%inner_join(pdacA.list$scRNA_anno,by="CB")
rownames(pdacA.seu@meta.data)=pdacA.seu@meta.data$CB

### 降维聚类
pdacA.seu <- FindNeighbors(pdacA.seu, dims = 1:20)
pdacA.seu <- FindClusters(pdacA.seu, resolution = 0.5)
pdacA.seu <- RunUMAP(pdacA.seu, dims = 1:20)
pdacA.seu <- RunTSNE(pdacA.seu, dims = 1:20)
DimPlot(pdacA.seu,group.by = "celltype",reduction = "umap",pt.size = 1,label = T,repel = T,label.size = 4)

### 找差异基因
# 控制三个阈值：logfc.threshold p_val_adj d
Idents(pdacA.seu)="celltype"
maintype_marker=FindAllMarkers(pdacA.seu,logfc.threshold = 0.5,only.pos = T)
maintype_marker=maintype_marker%>%filter(p_val_adj < 1e-05)
maintype_marker$d=maintype_marker$pct.1 - maintype_marker$pct.2
maintype_marker=maintype_marker%>%filter(d > 0.2)
maintype_marker=maintype_marker%>%arrange(cluster,desc(avg_log2FC))
maintype_marker=as.data.frame(maintype_marker)

2.空间转录组流程

注意：

• 这个文献只能下载到空转矩阵和spot坐标，没有图像。因此下面的流程和单细胞转录组一模一样；
• 若是完整的空转数据集，则有一些步骤是不一样的，具体流程可以参考seurat官网，https://satijalab.org/seurat/articles/spatial_vignette.html

pdacA.V.seu = CreateSeuratObject(counts = pdacA.list$ST_count)
pdacA.V.seu <- SCTransform(pdacA.V.seu, verbose = FALSE)
pdacA.V.seu <- RunPCA(pdacA.V.seu, assay = "SCT", verbose = FALSE)
pdacA.V.seu <- FindNeighbors(pdacA.V.seu, reduction = "pca", dims = 1:30)
pdacA.V.seu <- FindClusters(pdacA.V.seu, verbose = FALSE)
pdacA.V.seu <- RunUMAP(pdacA.V.seu, reduction = "pca", dims = 1:30)
pdacA.V.seu <- RunTSNE(pdacA.V.seu, reduction = "pca", dims = 1:30)

### 整合坐标、region
#（此处把坐标类比单细胞转录组分析中的注释信息）
# region信息由公粽号【TOP生物信息】整理
region.info=read.table("PDAC_A.region.txt",header = T,sep = "\t",stringsAsFactors = F)
pdacA.V.seu@meta.data$SB=rownames(pdacA.V.seu@meta.data)
pdacA.V.seu@meta.data=pdacA.V.seu@meta.data%>%inner_join(region.info,by = "SB")
rownames(pdacA.V.seu@meta.data)=pdacA.V.seu@meta.data$SB

### 画图看看
library(RColorBrewer)
library(scales)
color_region=c("#666666","#766fb1","#e42a88","#189d77")
names(color_region)=c("Cancer region","Pancreatic tissue","Duct epithelium","Stroma")

pdacA.V.seu@meta.data%>%ggplot(aes(x=x_coord,y=y_coord,fill=region))+
  geom_tile(color="white")+
  scale_fill_manual("Cluster assignments",values = color_region)+
  theme_void()+
  theme(
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 14),
    legend.position = "bottom",
    legend.direction = "vertical"
  )+
  guides(fill = guide_legend(override.aes = list(size=10)))
ggsave("fig.2e.pdf",width = 10,height = 14,units = "cm")

### 找region特异基因
Idents(pdacA.V.seu)="region"
region_marker=FindAllMarkers(pdacA.V.seu,logfc.threshold = 0,only.pos = T)
region_marker=region_marker%>%filter(p_val_adj < 0.1)
region_marker$d=region_marker$pct.1 - region_marker$pct.2
region_marker=region_marker%>%filter(d > 0.05)
region_marker=region_marker%>%arrange(cluster,desc(avg_log2FC))
region_marker=as.data.frame(region_marker)

说明：

• 1.上述找两个DEG数据框的方法不唯一，阈值也不唯一
• 2.第二个DEG数据框也可以是空间cluster的marker
• 3.MIA分析模式在单细胞和空间转录组场景都可以应用，空转场景是看细胞亚群的富集程度，单细胞场景是做细胞亚群注释

3.MIA分析

region_specific=region_marker[,c("cluster","gene")]
colnames(region_specific)[1]="region"

celltype_specific=maintype_marker[,c("cluster","gene")]
colnames(celltype_specific)[1]="celltype"

N=length(union(rownames(pdacA.seu),rownames(pdacA.V.seu)))

library(RColorBrewer)
library(scales)
color_region=c("#666666","#766fb1","#e42a88","#189d77")
names(color_region)=c("Cancer region","Pancreatic tissue","Duct epithelium","Stroma")

source("syMIA.R")
miares=syMIA(region_specific,celltype_specific,N,color_region)

之后会返回一个数据框，以及自动生成一张图

可以看到cancer region同时富集了cancer cell和fibroblast，这和文章结果是「一致」的。

当然上图和原文有些差别，原因是我没有采用原文找差异基因的方法，而是选用了seurat中的常规方法，这个灵活选择。

本文数据整理和代码编写花费大量时间，故不无偿提供，有需要的朋友可以关注公粽号，然后回复2022B