空間轉(zhuǎn)錄組與單細(xì)胞轉(zhuǎn)錄組關(guān)聯(lián)分析方法

今天給大家介紹的是BMKMANU S1000空間轉(zhuǎn)錄組數(shù)據(jù)與單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的關(guān)聯(lián)分析。在我們的關(guān)聯(lián)分析中可以使用不同分辨率的空間數(shù)據(jù)分別與單細(xì)胞數(shù)據(jù)進(jìn)行關(guān)聯(lián)分析，本次我們使用spatialDWLS進(jìn)行100μm關(guān)聯(lián)分析，下周我們則會(huì)推出多級(jí)分辨率下的關(guān)聯(lián)分析。詳細(xì)分析流程如下所示：

spatialDWLS

spatialDWLS首先使用細(xì)胞類型富集分析方法來確定哪些類型的細(xì)胞在每個(gè)位置具有較高的概率，然后使用阻尼加權(quán)最小二乘法（DWLS）的擴(kuò)展來確定指定位置的細(xì)胞類型的組成。與現(xiàn)有的去卷積方法相比，關(guān)鍵區(qū)別在于SpatialDWLS包含額外的過濾步驟，以去除不相關(guān)的細(xì)胞類型，從而增強(qiáng)特異性。

1、單細(xì)胞數(shù)據(jù)格式要求

（1）data 1 聚類和注釋結(jié)果（sc_meta.txt）

（2）data 2 基因表達(dá)矩陣（all_cells_count_matrix_filtered.tsv）

2、空間數(shù)據(jù)格式要求

（1）data 1 空間spots位置文件（loc.txt）

（2）data 2 基因表達(dá)矩陣（exp.txt）

#安裝Giotto

library(devtools) ?# if not installed: install.packages(‘devtools’)

library(remotes) ?# if not installed: install.packages(‘remotes’)

remotes::install_github(“RubD/Giotto”)

#加載R包

library(Seurat)

library(SingleR)

library(celldex)

library(Giotto)

library(ggplot2)

library(scatterpie)

library(data.table)

library(tibble)

#單細(xì)胞數(shù)據(jù)分析

#讀入單細(xì)胞數(shù)據(jù)

sc_meta<-read.table(“../sc_meta.txt”,header?=?F,row.names?=?1,sep=”\t”)

sc_data<-read.table(“../all_cells_count_matrix_filtered.tsv”,header?=?T,row.names?=?1)

#設(shè)置python路徑

my_python_path?=?“/xxx/software/miniconda3/bin/python”

instrs?=?createGiottoInstructions(python_path?=?my_python_path)

#創(chuàng)建Giotto 對(duì)象

sc?<-?createGiottoObject(raw_exprs?=?sc_data,instructions?=?instrs)

#表達(dá)數(shù)據(jù)標(biāo)準(zhǔn)分析

sc?<-?normalizeGiotto(gobject?=?sc, scalefactor?=?6000, verbose?=?T)

sc?<-?calculateHVG(gobject?=?sc, save_plot?=?TRUE,save_param=list(save_name=”sc_HVGplot”))

gene_metadata?=?fDataDT(sc)

featgenes?=?gene_metadata[hvg?==?‘yes’]$gene_ID

sc?<-?runPCA(gobject?=?sc, genes_to_use?=?featgenes, scale_unit?=?F)

signPCA(sc, genes_to_use?=?featgenes, scale_unit?=?F, save_plot?=?TRUE, ?save_param=list(save_name=”sc_signPCAplot”))

#使用Giotto中的DEG函數(shù)計(jì)算用于去卷積的Sig

sc@cell_metadata$leiden_clus?<-?as.character(sc_meta$V3)

scran_markers_subclusters?=?findMarkers_one_vs_all(gobject?=?sc,

method?=?‘scran’,

expression_values?=?‘normalized’,

cluster_column?=?‘leiden_clus’)

Sig_scran?<-?unique(scran_markers_subclusters$genes[which(scran_markers_subclusters$ranking?<=?100)])

#計(jì)算每種細(xì)胞類型中特征基因的表達(dá)中值

norm_exp<-2^(sc@norm_expr)-1

id<-sc@cell_metadata$leiden_clus

ExprSubset<-norm_exp[Sig_scran,]

Sig_exp<-NULL

for?(i?in?unique(id)){

Sig_exp<-cbind(Sig_exp,(apply(ExprSubset,1,function(y) mean(y[which(id==i)]))))

}

colnames(Sig_exp)<-unique(id)

#空間數(shù)據(jù)分析

#讀入數(shù)據(jù)

spatial_loc<-read.table(file=”loc.txt”,header?=?F)

spatial_exp<-read.table(file=”exp.txt”,header?=?T,row.names?=?1)

#創(chuàng)建Giotto對(duì)象

st?<-?createGiottoObject(raw_exprs?=?spatial_exp,spatial_locs?=?spatial_loc,

instructions?=?instrs)

#對(duì)空間數(shù)據(jù)進(jìn)行過濾，均一化和聚類

st?<-?filterGiotto(gobject?=?st,

expression_threshold?=?1,

gene_det_in_min_cells?=?5,

min_det_genes_per_cell?=?100,

expression_values?=?c(‘raw’),

verbose?=?T)

st?<-?normalizeGiotto(gobject?=?st)

st?<-?calculateHVG(gobject?=?st, save_plot?=?TRUE, save_param=list(save_name=”st_HVGplot”))

gene_metadata?=?fDataDT(st)

featgenes?=?gene_metadata[hvg?==?‘yes’]$gene_ID

st?<-?runPCA(gobject?=?st, genes_to_use?=?featgenes, scale_unit?=?F)

signPCA(st, genes_to_use?=?featgenes, scale_unit?=?F, save_plot?=?TRUE,save_param=list(save_name=”st_signPCAplot”))

st?<-?createNearestNetwork(gobject?=?st, dimensions_to_use?=?1:10, k?=?10)

st?<-?doLeidenCluster(gobject?=?st, resolution?=?0.5, n_iterations?=?1000)

#基于單細(xì)胞數(shù)據(jù)的基因表達(dá)特征和空間數(shù)據(jù)的Giotto對(duì)象進(jìn)行去卷積分析

st?<-?runDWLSDeconv(gobject?=?st, sign_matrix?=?Sig_exp)

#可視化展示

plot_data?<-?as.data.frame(st@spatial_enrichment$DWLS)[-1]

plot_col?<-?colnames(plot_data)

plot_data$x?<-?as.numeric(as.character(st@spatial_locs$sdimx))

plot_data$y?<-?as.numeric(as.character(st@spatial_locs$sdimy))

min_x?<-?min(plot_data$x)

plot_data$radius?<-?7

df?<-?data.frame()

p?<-?ggplot(df) +?geom_point() +?xlim(min(plot_data$x)-7, max(plot_data$x)+7) +?ylim(min(plot_data$y)-7, max(plot_data$y)+7)

p?<-?p?+?geom_scatterpie(aes(x=x, y=y, r=radius), data=plot_data, cols=plot_col, color=NA, alpha=.8) +scale_y_reverse()+

geom_scatterpie_legend(plot_data$radius, x=1, y=1) +

scale_fill_manual(values?=?c(“#F56867″,”#FEB915″,”#C798EE”,”#59BE86″,”#7495D3″,”#D1D1D1″,”#6D1A9C”,”#15821E”,”#3A84E6″,”#997273″,”#787878″,”#DB4C6C”,”#9E7A7A”,”#554236″))+

theme_classic()

ggsave(p, filename?=?“plot.pdf”,height=8,width=5)