课前准备---空间基因梯度2

作者，Evil Genius

目前我们的课程已经上了20节课了，预计还有3周，我们本年度的课程就要全部结束了，大家都是硕士博士，并且身处国家重点单位、高校，多学点空间分析技能，自然对数据的分析很有帮助，并且输出一些高质量的分析文章，对这个社会、科研等发展，都有很正向的作用。

不知不觉都已经走了很多年了，其实我自己也没什么拿得出手的成果，唯一让我觉得自己走过这段过程的印记，大概就是经常写关于单细胞空间的“日记”吧。

感觉课程上完，应该可以突破200万字了。虽然比起大佬们很寒酸啊，没有什么好的文章发表，也“水”过低分的文章，基本上不了台面。

至于空间基因梯度，主要是有如下的运用。

对伤口损伤的小鼠大脑皮层（损伤后3天）进行空间转录组梯度分析中，空间梯度分析从损伤核心（深红色点）向周边（浅粉色）的区域内进行分析，其中还这涉及到基因调控的一些内容。

我们需要实现的目标

完整的脚本如下

#! usr/R
###zhaoyunfei
###20240801

library(Seurat)
library(Matrix)
library(RcppML)  
library(ggplot2)
library(dplyr)
library(LSGI)
library(msigdbr)
library(hypeR)
library(magick)

sample = 'liver2'

img <- Read10X_Image(image.dir = "/home/samples/DB/Spatial/data/ST/ST-liver2/",image.name = "tissue_lowres_image.png", filter.matrix = TRUE)

data <- Load10X_Spatial(data.dir = "/home/samples/DB/Spatial/data/ST/ST-liver2/",filename = "filtered_feature_bc_matrix.h5",assay = "Spatial", slice = sample, filter.matrix = TRUE, to.upper = FALSE,image = img)

data <- NormalizeData(data)

###NMF
scan.nmf.mse <- function(obj, ranks = seq(1, 30, 2), tol = 1e-04) {
    # users can customize the scan by changing 'ranks'
    dat <- obj@assays$Spatial@data
    errors <- c()
    ranks <- seq(1, 30, 2)
    for (i in ranks) {
        # cat('rank: ', i, '\n')
        mod <- RcppML::nmf(dat, i, tol = 1e-04, verbose = F)
        mse_i <- mse(dat, mod$w, mod$d, mod$h)
        errors <- c(errors, mse_i)
    }
    results <- data.frame(rank = ranks, MSE = errors)
    return(results)
}

sr.nmf <- function(obj, k = 10, tol = 1e-06, assay = "Spatial") {
    dat <- obj@assays$Spatial@data
    nmf_model <- RcppML::nmf(dat, k = k, tol = tol, verbose = F)
    embeddings <- t(nmf_model$h)
    rownames(embeddings) <- colnames(obj)
    colnames(embeddings) <- paste0("nmf_", 1:k)
    loadings <- nmf_model$w
    rownames(loadings) <- rownames(obj)
    obj@reductions$nmf <- CreateDimReducObject(embeddings = embeddings,
        loadings = loadings, key = "nmf_", assay = assay)
    return(obj)
}

scan.nmf.res <- scan.nmf.mse(obj = data)

res = 300

png('NMF.k.png',width = 7 * res ,height = 8 * res, res = 300 ,type = c("cairo"))

ggplot(scan.nmf.res, aes(x = rank, y = MSE)) + geom_point(size = 0.7) + geom_smooth(method = "loess", span = 0.2, color = "black",linewidth = 1, se = F) + labs(x = "NMF rank", y = "MSE") + theme_classic() + scale_y_continuous(expand = c(0.01, 0)) + theme(aspect.ratio = 1)

dev.off()

data <- sr.nmf(obj = data, k = 10, tol = 1e-05)

###Prepare input data for LSGI

spatial_coords <- data@images[[sample]]@coordinates[, c(4, 5)]

colnames(spatial_coords) <- c("X", "Y")

embeddings <- data@reductions$nmf@cell.embeddings

# neighborhood
lsgi.res <- local.traj.preprocessing(spatial_coords = spatial_coords, n.grids.scale = 2, embeddings = embeddings, n.cells.per.meta = 25)

####注释
get.nmf.info <- function(obj, top.n = 50) {
    feature.loadings <- as.data.frame(obj@reductions$nmf@feature.loadings)

    top.gene.list <- list()
    for (i in 1:ncol(feature.loadings)) {
        o <- order(feature.loadings[, i], decreasing = T)[1:top.n]
        features <- rownames(feature.loadings)[o]
        top.gene.list[[colnames(feature.loadings)[i]]] <- features
    }
    nmf.info <- list(feature.loadings = feature.loadings, top.genes = top.gene.list)
    return(nmf.info)
}

nmf_info <- get.nmf.info(data)

mdb_h <- msigdbr(species = "Homo sapiens", category = "H")

gene.set.list <- list()
for (gene.set.name in unique(mdb_h$gs_name)) {
    gene.set.list[[gene.set.name]] <- mdb_h[mdb_h$gs_name %in%
        gene.set.name, ]$gene_symbol
}

# run hypeR test
mhyp <- hypeR(signature = nmf_info$top.genes, genesets = gene.set.list,
    test = "hypergeometric", background = rownames(nmf_info[["feature.loadings"]]))
hyper.data <- mhyp$data
hyper.res.list <- list()
for (nmf.name in names(hyper.data)) {
    res <- as.data.frame(hyper.data[[nmf.name]]$data)
    hyper.res.list[[nmf.name]] <- res
}

png('NMF.enrichment.png',width = 7 * res ,height = 8 * res, res = 300 ,type = c("cairo"))
ggplot(hyper.res.list[[1]][1:5, ], aes(x = reorder(label, -log10(fdr)),
    y = overlap/signature, fill = -log10(fdr))) + geom_bar(stat = "identity",
    show.legend = T) + xlab("Gene Set") + ylab("Gene Ratio") +
    viridis::scale_fill_viridis() + theme_classic() + coord_flip() +
    theme(axis.text.x = element_text(color = "black", size = 12,
        angle = 45, hjust = 1), axis.text.y = element_text(color = "black",
        size = 8, angle = 0), panel.border = element_rect(colour = "black",
        fill = NA, size = 1))

dev.off()

plot.overlay.factor <- function(object, info, sel.factors = NULL,
    r_squared_thresh = 0.6, minimum.fctr = 20) {
    scf <- object@images[[sample]]@scale.factors[["lowres"]]
    object <- subset(object, cells = rownames(info$spatial_coords))
    print(identical(rownames(object@meta.data), rownames(info$spatial_coords)))
    object <- rotateSeuratImage(object, rotation = "L90")
    object@meta.data <- cbind(object@meta.data, info$embeddings)

    lin.res.df <- get.ind.rsqrs(info)
    lin.res.df <- na.omit(lin.res.df)
    lin.res.df <- lin.res.df[lin.res.df$rsquared > r_squared_thresh,
        ]
    if (!is.null(sel.factors)) {
        lin.res.df <- lin.res.df[lin.res.df$fctr %in% sel.factors,
            ]
    }
    lin.res.df <- lin.res.df %>%
        group_by(fctr) %>%
        filter(n() >= minimum.fctr) %>%
        ungroup()

    spatial_coords <- info$spatial_coords
    spatial_coords$X <- spatial_coords$X * scf
    spatial_coords$Y <- spatial_coords$Y * scf

    lin.res.df$Xend <- lin.res.df$X + lin.res.df$vx.u
    lin.res.df$Yend <- lin.res.df$Y + lin.res.df$vy.u

    lin.res.df$X <- lin.res.df$X * scf
    lin.res.df$Xend <- lin.res.df$Xend * scf
    lin.res.df$Y <- lin.res.df$Y * scf
    lin.res.df$Yend <- lin.res.df$Yend * scf

    p <- SpatialFeaturePlot(object, features = NULL, alpha = c(0)) +
        NoLegend() + geom_segment(data = as.data.frame(lin.res.df),
        aes(x = X, y = Y, xend = Xend, yend = Yend, color = fctr,
            fill = NULL), linewidth = 0.4, arrow = arrow(length = unit(0.1,
            "cm"))) + scale_color_brewer(palette = "Paired") +
        # scale_fill_gradient(low='lightgrey', high='navy')
        # +
    theme_classic() + theme(axis.text.x = element_text(face = "bold",
        color = "black", size = 12, angle = 0, hjust = 1), axis.text.y = element_text(face = "bold",
        color = "black", size = 12, angle = 0))

    return(p)
}

# Adapted from this link:
# https://github.com/satijalab/seurat/issues/2702#issuecomment-1626010475
rotateSeuratImage <- function(seuratVisumObject, slide = sample,
    rotation = "Vf") {
    if (!(rotation %in% c("180", "Hf", "Vf", "L90", "R90"))) {
        cat("Rotation should be either 180, L90, R90, Hf or Vf\n")
        return(NULL)
    } else {
        seurat.visium <- seuratVisumObject
        ori.array <- (seurat.visium@images)[[slide]]@image
        img.dim <- dim(ori.array)[1:2]/(seurat.visium@images)[[slide]]@scale.factors$lowres
        new.mx <- c()
        # transform the image array
        for (rgb_idx in 1:3) {
            each.mx <- ori.array[, , rgb_idx]
            each.mx.trans <- rotimat(each.mx, rotation)
            new.mx <- c(new.mx, list(each.mx.trans))
        }

        # construct new rgb image array
        new.X.dim <- dim(each.mx.trans)[1]
        new.Y.dim <- dim(each.mx.trans)[2]
        new.array <- array(c(new.mx[[1]], new.mx[[2]], new.mx[[3]]),
            dim = c(new.X.dim, new.Y.dim, 3))

        # swap old image with new image
        seurat.visium@images[[slide]]@image <- new.array

        ## step4: change the tissue pixel-spot index
        img.index <- (seurat.visium@images)[[slide]]@coordinates

        # swap index
        if (rotation == "Hf") {
            seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[2] -
                img.index$imagecol
        }

        if (rotation == "Vf") {
            seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[1] -
                img.index$imagerow
        }

        if (rotation == "180") {
            seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[1] -
                img.index$imagerow
            seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[2] -
                img.index$imagecol
        }

        if (rotation == "L90") {
            seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[2] -
                img.index$imagecol
            seurat.visium@images[[slide]]@coordinates$imagecol <- img.index$imagerow
        }

        if (rotation == "R90") {
            seurat.visium@images[[slide]]@coordinates$imagerow <- img.index$imagecol
            seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[1] -
                img.index$imagerow
        }

        return(seurat.visium)
    }
}

rotimat <- function(foo, rotation) {
    if (!is.matrix(foo)) {
        cat("Input is not a matrix")
        return(foo)
    }
    if (!(rotation %in% c("180", "Hf", "Vf", "R90", "L90"))) {
        cat("Rotation should be either L90, R90, 180, Hf or Vf\n")
        return(foo)
    }
    if (rotation == "180") {
        foo <- foo %>%
            .[, dim(.)[2]:1] %>%
            .[dim(.)[1]:1, ]
    }
    if (rotation == "Hf") {
        foo <- foo %>%
            .[, dim(.)[2]:1]
    }

    if (rotation == "Vf") {
        foo <- foo %>%
            .[dim(.)[1]:1, ]
    }
    if (rotation == "L90") {
        foo = t(foo)
        foo <- foo %>%
            .[dim(.)[1]:1, ]
    }
    if (rotation == "R90") {
        foo = t(foo)
        foo <- foo %>%
            .[, dim(.)[2]:1]
    }
    return(foo)
}

png('Spatial.gredient.png',width = 7 * res ,height = 8 * res, res = 300 ,type = c("cairo"))

plot.overlay.factor(object = data, info = lsgi.res, sel.factors = NULL)

dev.off()

生活很好，有你更好