Liger_DimPlot.r

#!/usr/bin/env Rscript

#### 脚本说明
#1) 运行环境: mamba activate r452
#2) 应用: 对Seurat对象重新绘制单细胞降维图 (基于Liger[https://github.com/welch-lab/liger.git])

### 调用包
library(Seurat)
library(ggplot2)
library(dplyr)
library(patchwork)
library(qs)

### 调用Liger原生函数
liger_path = "/data/med-wangcq/01CondaEnv/01Demo_Script/00Git_Clone/liger/R/"
# invisible(lapply(list.files(liger_path, pattern = "\\.r$", full.names = TRUE), source))
source(file.path(liger_path,"ggplotting.R"))
source(file.path(liger_path,"util.R"))
source(file.path(liger_path,"aaa.R"))

# ### 基本参数
# inRDS <- "/scratch/2026-08-19/med-wangcq/Others/AnJQ_source/AnalysisData/Sam68OECM_scFASTseq/02Result/01scRNA/03LIGER/02_Annotation/00rds/Annotation_rmCluster_main_Seurat.qs"
# useDimRed <- "liger_UMAP"
# slot <- "cellMeta"
# colorBy <- "Phase"
# shapeBy <- NULL
# titles <- NULL

# ### 基本处理
# object <- qread(inRDS)

### 整理成Liger 绘图包需要的格式
to_ligerDF <- function(obj, addCol = NULL, reduc = "umap"){
    ## 参数说明：
    #1) obj: Seurat对象
    #2) addCol: 作图分组信息列, 向量, NULL的时候使用默认idents 
    #3) reduc: 降维坐标使用
    # 获取降维坐标矩阵
    reDuc_emb <- Embeddings(obj, reduction = reduc)
    # 添加注释信息
    if(!is.null(addCol)){
        anno <- obj@meta.data[, addCol]
    }else{
        # 默认 Idents
        anno <- data.frame(Ident_use=Idents(obj))
    }
    plotDF <- cbind(reDuc_emb, anno)
    return(plotDF)
}

### Liger 降维图 函数
Liger_ggScatter <- function(plotDF,
                            x,
                            y,
                            colorBy = NULL,
                            shapeBy = NULL,
                            dotOrder = c("shuffle", "ascending", "descending"),
                            dotSize = getOption("ligerDotSize"),
                            dotAlpha = 0.9,
                            trimHigh = NULL,
                            trimLow = NULL,
                            zeroAsNA = TRUE,
                            raster = NULL,
                            labelBy = colorBy,
                            labelText = TRUE,
                            labelTextSize = 4,
                            ggrepelLabelTick = FALSE,
                            seed = 1) {
        dotOrder <- match.arg(dotOrder)
        set.seed(seed)
        raster <- .checkRaster(nrow(plotDF), raster)
        if (!is.null(colorBy)) {
            if (is.character(plotDF[[colorBy]]) ||
                is.logical(plotDF[[colorBy]])) {
                plotDF[[colorBy]] <- factor(plotDF[[colorBy]])
            }
            if (dotOrder == "shuffle") {
                # Always put NA at bottom layer, i.e. plot them first
                isNA <- which(is.na(plotDF[[colorBy]]))
                nonNA <- which(!is.na(plotDF[[colorBy]]))
                idx <- sample(nonNA)
                plotDF <- plotDF[c(isNA, idx), ]
            } else if (dotOrder == "ascending") {
                plotDF <- plotDF[order(plotDF[[colorBy]], decreasing = FALSE,
                                    na.last = FALSE),]
            } else {
                plotDF <- plotDF[order(plotDF[[colorBy]], decreasing = TRUE,
                                    na.last = FALSE),]
            }
            if (!is.factor(plotDF[[colorBy]]) &&
                !is.character(plotDF[[colorBy]]) &&
                !is.logical(plotDF[[colorBy]])) {
                if (!is.null(trimHigh))
                    plotDF[[colorBy]][plotDF[[colorBy]] > trimHigh] <- trimHigh
                if (!is.null(trimLow))
                    plotDF[[colorBy]][plotDF[[colorBy]] < trimLow] <- trimLow
                if (isTRUE(zeroAsNA))
                    plotDF[[colorBy]][plotDF[[colorBy]] == 0] <- NA
            }
        }

        # TODO: Any way to avoid using these many conditions?
        if (!is.null(colorBy)) {
            if (!is.null(shapeBy))
                p <- ggplot2::ggplot(plotDF,
                                    ggplot2::aes(x = .data[[x]],
                                                y = .data[[y]],
                                                color = .data[[colorBy]],
                                                shape = .data[[shapeBy]]))
            else
                p <- ggplot2::ggplot(plotDF,
                                    ggplot2::aes(x = .data[[x]],
                                                y = .data[[y]],
                                                color = .data[[colorBy]]))
        } else {
            if (!is.null(shapeBy))
                p <- ggplot2::ggplot(plotDF,
                                    ggplot2::aes(x = .data[[x]],
                                                y = .data[[y]],
                                                shape = .data[[shapeBy]]))
            else
                p <- ggplot2::ggplot(plotDF,
                                    ggplot2::aes(x = .data[[x]],
                                                y = .data[[y]]))
        }

        if (isTRUE(raster)) {
            p <- p + scattermore::geom_scattermore(pointsize = dotSize,
                                                alpha = dotAlpha)
        } else {
            p <- p + ggplot2::geom_point(size = dotSize, stroke = 0,
                                        alpha = dotAlpha)
        }

        # For categorical grouping
        if (!is.null(labelBy) &&
            (is.factor(plotDF[[labelBy]]) || is.character(plotDF[[labelBy]]))) {
            if (isTRUE(labelText)) {
                textData <- dplyr::group_by(plotDF, .data[[labelBy]]) %>%
                    dplyr::summarise(
                        x = stats::median(.data[[x]]),
                        y = stats::median(.data[[y]])
                    ) %>%
                    dplyr::filter(!is.na(.data[[labelBy]]), !is.na(x), !is.na(y))
                if (!requireNamespace("ggrepel", quietly = TRUE)) {
                    p <- p + ggplot2::annotate(
                        "text", x = textData$x, y = textData$y,
                        label = textData[[labelBy]], color = "black",
                        size = labelTextSize
                    )
                } else {
                    if (isTRUE(ggrepelLabelTick)) {
                        p <- p + ggrepel::geom_text_repel(
                            data = textData,
                            mapping = ggplot2::aes(x = .data[["x"]],
                                                y = .data[["y"]],
                                                label = .data[[labelBy]]),
                            color = "black", size = labelTextSize, inherit.aes = FALSE,
                            bg.colour = "white", bg.r = .2,
                            force = 1,
                            min.segment.length = 0,
                            nudge_y = 1
                        )
                    } else {
                        p <- p + ggrepel::geom_text_repel(
                            data = textData,
                            mapping = ggplot2::aes(x = .data[["x"]],
                                                y = .data[["y"]],
                                                label = .data[[labelBy]]),
                            color = "black", size = labelTextSize, inherit.aes = FALSE,
                            bg.colour = "white", bg.r = .2, hjust = 0.5, vjust = 0.5
                        )
                    }
                }

                # Important to have `inherit.aes = F` above, otherwise
                # `geom_text_repel` looks for "shapeBy" setting which this newly
                # generated label coordinate table just doesn't have.
            }
        }
        p <- .ggplotLigerTheme(p)

        return(p)
}


### 主调用脚本
plot_main <- function(obj,
                      colorBy=NULL,
                      split.by=NULL,
                      shapeBy = NULL, 
                      useDimRed="umap", 
                      ncol = 1, 
                      pt.size = 0.2,
                      pt.alpha = 0.5,
                      sWidth = 6, 
                      sHeight = 5, 
                      gap = 0.5,
                      outdir = "./",
                      PlotName = "DimPlot"){
        # sWidth, sHeight, gap单位英寸, 指定小图尺寸
        plotDF <- to_ligerDF(obj, addCol = c(colorBy,shapeBy,split.by), reduc = useDimRed)
        if(!is.null(split.by)){
            df_list <- split(plotDF, plotDF[,split.by])
            plot_list <- lapply(names(df_list), function(gp){
                            pDF <- df_list[[gp]]
                            p <- Liger_ggScatter(plotDF = pDF, 
                                        x = colnames(pDF)[1], 
                                        y = colnames(pDF)[2], 
                                        colorBy = colorBy,
                                        shapeBy = shapeBy,
                                        dotOrder = c("shuffle", "ascending", "descending"),
                                        dotSize = pt.size,
                                        # getOption("ligerDotSize"),
                                        dotAlpha = pt.alpha,
                                        trimHigh = NULL,
                                        trimLow = NULL,
                                        zeroAsNA = TRUE,
                                        raster = NULL,
                                        labelBy = colorBy,
                                        labelText = TRUE,
                                        labelTextSize = 4,
                                        ggrepelLabelTick = FALSE,
                                        seed = 1)
                            return(p)}
                            )
            # 拼图并控制布局
            nrow <- ceiling(length(plot_list)/ncol)
            final_plot <- wrap_plots(plot_list, ncol = ncol, nrow = nrow, byrow = TRUE, guides = "collect") +
                                        plot_layout(widths = rep(1,ncol), heights = rep(1, nrow))
            # 自动计算整体尺寸
            total_width <- ncol * sWidth + (ncol - 1) * gap
            total_height <- nrow * sHeight + (nrow - 1) * gap
        }else{
            final_plot <- Liger_ggScatter(plotDF = plotDF, 
                                    x = colnames(plotDF)[1], 
                                    y = colnames(plotDF)[2], 
                                    colorBy = colorBy,
                                    shapeBy = shapeBy,
                                    dotOrder = c("shuffle", "ascending", "descending"),
                                    dotSize = pt.size,
                                    # getOption("ligerDotSize"),
                                    dotAlpha = pt.alpha,
                                    trimHigh = NULL,
                                    trimLow = NULL,
                                    zeroAsNA = TRUE,
                                    raster = NULL,
                                    labelBy = colorBy,
                                    labelText = TRUE,
                                    labelTextSize = 4,
                                    ggrepelLabelTick = FALSE,
                                    seed = 1)
            ncol <- 1
            nrow <- 1
            # 自动计算整体尺寸
            total_width <- sWidth + (ncol - 1) * gap
            total_height <- sHeight + (nrow - 1) * gap
        }
        ggsave(file.path(outdir,paste0(PlotName,".pdf")), final_plot, width = total_width, height = total_height, units = "in", dpi = 300)
        # 或者
        # pdf(file.path(outdir,paste0(PlotName,".pdf")), width = total_width, height = total_height)
        # print(final_plot)
        # dev.off()
}