06Mapping.r

#!/usr/bin/env Rscript

#### 脚本使用说明
#1) mamba activate r452
#2) 应用: 多个子object映射到一个主object

### 加载包
source("/data/med-wangcq/01CondaEnv/02Git_repo/scRNA/00Functions/00_loadPackages_v4.r")
source("/data/med-wangcq/01CondaEnv/02Git_repo/scRNA/00Functions/00_ColorPanel.r")

cols <- c("#8494FF", "#FDC086", "#7FC97F", "#386CB0", "#F0027F", "#BF5B17", "#1F78B4",
                  "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00", "#CAB2D6",
                  "#6A3D9A", "#FFFF99", "#B15928", "#A6CEE3", "#1F9CEF", "#B2DFDB", "#A676A6",
                  "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999", "#E69F00", "#56B4E9",
                  "#009E73", "#F0E442", "#0072B2", "#D55E00")

### 基本参数
option_list <- list(
    make_option(c("-r", "--rdsList"), type = "character",action = "store", help = "三列, tab分隔tsv/csv/xls/xlsx文件,第一列指定rds文件类型[main/sub],第二列提供rds文件路径[.qs/.rds],第三列指定映射使用列"),
    make_option(c("-m", "--rmCluster"), type = "character",action = "store", help = "默认'no', 对mapping ident后的主脚本处理, subset需要的细胞类型出图", default = "no"),
    make_option(c("-o", "--outdir"), type = "character",action = "store", help = "默认'./', 输出文件夹路径", default = "./"),
    make_option(c("-p", "--prefix"), type = "character",action = "store", help = "默认'no', 输出文件前缀", default = "no"),
    make_option(c("-c", "--ct_ord"), type = "character",action = "store", help = "默认'default', 输出细胞类型顺序, 逗号分隔", default = 'default'),
    make_option(c("-a", "--assay"), type = "character",action = "store", help = "默认'RNA', 指定使用assay", default = "RNA"), 
    make_option(c("-d", "--reduction"), type = "character",action = "store", help = "默认'umap', 降维坐标矩阵", default = "umap"),
    make_option(c("-t", "--topn"), type = "integer",action = "store", help = "默认'10', 输出top N差异基因", default = "10")
    )

# 解析参数
opt = parse_args(OptionParser(option_list = option_list, 
                              usage = "Usage: %prog [options] \nDescription: This Script is for Cluster analysis of scRNAseq!"))

# rdsList <- "/data/med-wangcq/01CondaEnv/02Git_repo/scRNA/01DataDemo/Mapping_DemoData.xlsx"
# rmCluster <- "rmCluster"
# outdir <- "/scratch/2026-07-27/med-wangcq/Others/Dong/02Result/"
# prefix <- "main_303132_PC"
# ct_ord <- "default"
# assay <- "RNA"
# reduction <- "umap"
# topn <- 10

rdsList <- opt$rdsList
rmCluster <- opt$rmCluster
outdir <- opt$outdir
prefix <- opt$prefix
ct_ord <- opt$ct_ord
assay <- opt$assay
reduction <- opt$reduction
topn <- opt$topn

### 参数处理
#1) 建立新文件夹
outdir1 <- ifelse(prefix == "no", file.path(outdir,"05Mapping"), file.path(outdir,"05Mapping",prefix))
outrds <- file.path(outdir,"05Mapping","00rds")
if(!file.exists(outdir1)) dir.create(outdir1, recursive=TRUE)
if(!file.exists(outrds)) dir.create(outrds, recursive=TRUE)

cat("输出文件路径",outdir1,"\n")
cat("rds存储路径",outrds,"\n")

#2) 提取rds文件路径
ext <- tools::file_ext(rdsList)
cat("文件后缀: ",ext,"\n")
if(ext == "xlsx"){
    data <- read_xlsx(path=rdsList,trim_ws=TRUE,sheet=1,col_names=FALSE) %>% as.data.frame()
}else if(ext %in% c("tsv","xls")){
    data <- read.table(rdsList,header=FALSE,sep="\t")
}else if(ext == "csv"){
    data <- read.csv(rdsList,header=FALSE,sep=",")
}
colnames(data) <- c("Type","Path","ident")

#3) 合并映射
# 建立list存储rds结果
rds_list <- list()
rds_list[["main"]]$rds <- c()
rds_list[["main"]]$ident <- c()
rds_list[["sub"]]$rds <- c()
rds_list[["sub"]]$ident <- c()

for (i in 1:nrow(data)){
    cat("Processing: ",data$Type[i],"\n")
    cat("File Path: ",data$Path[i],"\n")
    cat("Using ident: ",data$ident[i],"\n")
    ext <- tools::file_ext(data$Path[i])
    if (ext == "qs"){
        rds_list[[data$Type[i]]]$rds <- c(rds_list[[data$Type[i]]]$rds, qread(data$Path[i]))
        rds_list[[data$Type[i]]]$ident <- c(rds_list[[data$Type[i]]]$ident, data$ident[i])
    }else if(ext == "rds"){
        rds_list[[data$Type[i]]]$rds <- c(rds_list[[data$Type[i]]]$rds, readRDS(data$Path[i]))
        rds_list[[data$Type[i]]]$ident <- c(rds_list[[data$Type[i]]]$ident, data$ident[i])
    }
}

# 映射细胞类型
cat("Processing Mapping ...\n")
main <- rds_list[["main"]]$rds[[1]]
meta <- main@meta.data
meta$Celltype_mapping <- meta[,rds_list[["main"]]$ident] %>% as.character()
meta$Mapping_source <- "main"
for(j in 1:length(rds_list[["sub"]]$ident)){
    cat("Processing sub Data ...\n")
    id <- rds_list[["sub"]]$ident[j]
    cat("Using ident: ",id,"\n")
    meta_sub <- rds_list[["sub"]]$rds[[j]]@meta.data
    print(unique(meta_sub[,id]))
    meta[rownames(meta_sub),"Celltype_mapping"] <- as.character(meta_sub[,id])
    meta[rownames(meta_sub),"Mapping_source"] <- paste0("sub_",j)
}
print(unique(meta$Celltype_mapping))
main@meta.data <- meta

#4) sub细胞类型
if(rmCluster != "no"){
    rmCluster <- unlist(strsplit(rmCluster,","))
    sub <- subset(main, subset=Celltype_mapping %in% rmCluster, invert=TRUE)
}else{
    sub <- main
}
print(levels(sub))

if(ct_ord != "default"){
    ct_ord <- unlist(strsplit(ct_ord,","))
}else{
    ct_ord <- unique(sub@meta.data$Celltype_mapping)
}
sub@meta.data$Celltype_mapping <- factor(sub@meta.data$Celltype_mapping, levels=ct_ord)
Idents(sub) <- sub@meta.data$Celltype_mapping
print(levels(sub))

rds_path <- file.path(outrds, ifelse(prefix == "no", "Mapping.qs", paste0("Mapping_",prefix,".qs")))
qsave(sub, rds_path)            # 输出存储文件
cat("rds file stores in path: ",rds_path,"\n")

#5) 输出画图
#### 作图：
print("Plotting")
sp_ord <- levels(sub$Sample)
gp_ord <- levels(sub$Group)
sp_cols <- cols[1:length(levels(sub$Sample))]
names(sp_cols) <- levels(sub$Sample)
gp_cols <- cols[1:length(levels(sub$Group))]
names(gp_cols) <- levels(sub$Group)
ct_cols <- cols[1:length(unique(sub$Celltype_mapping))]
names(ct_cols) <- levels(sub$Celltype_mapping)
print(sp_cols)
print(gp_cols)
print(ct_cols)

# 1)UMAP图
p <- DimPlot(sub, reduction = reduction, pt.size=0.1,cols=ct_cols,raster=FALSE)
pdf(file.path(outdir1,"UMAPplot.pdf"),width=6.5,height=5)
print(p)
dev.off()

p <- DimPlot(sub, reduction = reduction, pt.size=0.1,cols=ct_cols,label=T,raster=FALSE) + NoLegend()
pdf(file.path(outdir1,"UMAPplot_label.pdf"),width=5,height=5)
print(p)
dev.off()

p <- DimPlot(sub, reduction = reduction, split.by="Sample",pt.size=0.1,cols=ct_cols,raster=FALSE)
pdf(file.path(outdir1,"UMAPplot_splitbySample.pdf"),width=length(sp_ord)*5,height=5)
print(p)
dev.off()

p <- DimPlot(sub, reduction = reduction, split.by="Group",pt.size=0.1,cols=ct_cols,raster=FALSE)
pdf(file.path(outdir1,"UMAPplot_splitbyGroup.pdf"),width=length(gp_ord)*5,height=5)
print(p)
dev.off()

# 3)占比图bar图:
#分样本：
cellCount <- table(as.character(sub$Celltype_mapping),as.character(sub$Sample))
cellCount <- cellCount[ct_ord,sp_ord]
cellProp <- prop.table(cellCount,margin=1)
write.csv(cellCount,file.path(outdir1,"cellCount_Sample.csv"))
write.csv(cellProp,file.path(outdir1,"cellProp_SampleinCelltype.csv"))

#利用geom_bar()绘制填充式条形图——'fill'
plotdata <- as.data.frame(cellProp)
colnames(plotdata) <- c("Celltype","Sample","Prop")
plotdata$Sample <- factor(plotdata$Sample,levels=sp_ord)
plotdata$Celltype <- factor(plotdata$Celltype,levels=ct_ord)
p <- ggplot(plotdata, aes(Celltype, weight = Prop, fill = Sample)) +
        geom_bar(color = "black", width = 0.7, position = 'fill') +
        labs( y = 'Relative abundance (%)') +
        scale_fill_manual(values = sp_cols)+
        scale_y_continuous(expand = c(0,0)) +
        theme_classic() +
        theme(panel.background = element_blank(),
        panel.grid = element_blank(), 
        #axis.line = element_line(colour = "#000000",size = 2),
        #axis.text = element_text(colour = "#000000" ,size = 27),
        axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75,size = 10) ##就是这里
        #axis.ticks = element_line(colour = "#000000" ,size = 2) ,
        #axis.ticks.length = unit(2,'mm'),
        #plot.margin = unit(c(0.5,0,0,0),"cm"),
        #axis.title.y = element_text(size = 27),
        #axis.title.x = element_blank(),
        #plot.title = element_text(hjust = 0.5),
        #legend.position = "none"
        )

pdf(file.path(outdir1,"Barplot_SampleinCelltype.pdf"),width=6,height=5)
print(p)
dev.off()

# 细胞类型在各样本中占比：
cellProp <- prop.table(cellCount,margin=2)
write.csv(cellProp,file.path(outdir1,"cellProp_CelltypeinSample.csv"))
plotdata <- as.data.frame(cellProp)
colnames(plotdata) <- c("Celltype","Sample","Prop")
plotdata$Sample <- factor(plotdata$Sample,levels=sp_ord)
plotdata$Celltype <- factor(plotdata$Celltype,levels=ct_ord)
p <- ggplot(plotdata, aes(Sample, weight = Prop, fill = Celltype)) +
        geom_bar(color = "black", width = 0.7, position = 'fill') +
        labs( y = 'Relative abundance (%)') +
        scale_fill_manual(values = ct_cols)+
        scale_y_continuous(expand = c(0,0)) +
        theme_classic() +
        theme(panel.background = element_blank(),
        panel.grid = element_blank(), 
        #axis.line = element_line(colour = "#000000",size = 2),
        #axis.text = element_text(colour = "#000000" ,size = 27),
        axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75,size = 10) ##就是这里
        #axis.ticks = element_line(colour = "#000000" ,size = 2) ,
        #axis.ticks.length = unit(2,'mm'),
        #plot.margin = unit(c(0.5,0,0,0),"cm"),
        #axis.title.y = element_text(size = 27),
        #axis.title.x = element_blank(),
        #plot.title = element_text(hjust = 0.5),
        #legend.position = "none"
        )

pdf(file.path(outdir1,"Barplot_CelltypeinSample.pdf"),width=6,height=5)
print(p)
dev.off()

#分组：
if (length(gp_ord) >= 2) {
    cellCount <- table(as.character(sub$Celltype_mapping),as.character(sub$Group))
    cellCount <- cellCount[ct_ord,gp_ord]
    cellProp <- prop.table(cellCount,margin=1)
    write.csv(cellCount,file.path(outdir1,"cellCount_Group.csv"))
    write.csv(cellProp,file.path(outdir1,"cellProp_GroupinCelltype.csv"))
    #1) 画热图
    dd <- as.data.frame.matrix(cellProp)
    pdf(file.path(outdir1,"Heatmap_GroupinCelltype.pdf"),width=6,height=6)
    print(pheatmap(dd,
                    cellwidth = 20,
                    cellheight = 20,
                    cluster_rows = TRUE,
                    cluster_cols = TRUE,
                    color = colorRampPalette(c('#ffffff','#ff0000'))(100),
                    main = "Group in Celltype"))
    dev.off()
    
    ## 细胞类型在各组中占比：
    cellProp <- prop.table(cellCount,margin=2)
    write.csv(cellProp,file.path(outdir1,"cellProp_CelltypeinGroup.csv"))
    #1) 画热图
    dd <- as.data.frame.matrix(cellProp)
    pdf(file.path(outdir1,"Heatmap_CelltypeinGroup.pdf"),width=6,height=6)
    print(pheatmap(dd,
                    cellwidth = 20,
                    cellheight = 20,
                    cluster_rows = TRUE,
                    cluster_cols = TRUE,
                    color = colorRampPalette(c('#ffffff','#ff0000'))(100),
                    main = "Celltype in Group"))
    dev.off()
}
#2) 画堆积图
#利用geom_bar()绘制填充式条形图——'fill'
plotdata <- as.data.frame(cellProp)
colnames(plotdata) <- c("Celltype","Group","Prop")
plotdata$Group <- factor(plotdata$Group,levels=gp_ord)
plotdata$Celltype <- factor(plotdata$Celltype,levels=ct_ord)

if (length(gp_ord) == 2) {
    perc <- table(as.character(sub$Group)) %>% as.data.frame()
    linecept <- perc$Freq[perc$Var1 == gp_ord[2]]/sum(perc$Freq)
    p <- ggplot(plotdata, aes(Celltype, weight = Prop, fill = Group)) +
            geom_bar(color = "black", width = .7, position = 'fill') +
            geom_hline(yintercept=linecept, linetype='dashed', col = 'black',size = 1) +
            labs( y = 'Relative abundance (%)') +
            #scale_fill_brewer(palette = "Set3")+
            scale_fill_manual(values = gp_cols)+
            scale_y_continuous(expand = c(0,0)) +
            theme_classic() +
            theme(panel.background = element_blank(),
            panel.grid = element_blank(), 
            axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.8,size = 10))

    pdf(file.path(outdir1,"Barplot_GroupinCelltype.pdf"),width=6,height=5)
    print(p)
    dev.off()
}else if (length(gp_ord) > 2) {
    colnames(plotdata) <- c("Celltype","Group","Prop")
    plotdata$Group <- factor(plotdata$Group,levels=gp_ord)
    plotdata$Celltype <- factor(plotdata$Celltype,levels=ct_ord)
    p <- ggplot(plotdata, aes(Celltype, weight = Prop, fill = Group)) +
            geom_bar(color = "black", width = .7, position = 'fill') +
            labs( y = 'Relative abundance (%)') +
            #scale_fill_brewer(palette = "Set3")+
            scale_fill_manual(values = gp_cols)+
            scale_y_continuous(expand = c(0,0)) +
            theme_classic() +
            theme(panel.background = element_blank(),
            panel.grid = element_blank(), 
            axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.8,size = 10))

    pdf(file.path(outdir1,"Barplot_GroupinCelltype.pdf"),width=6,height=5)
    print(p)
    dev.off()
}

#2) 画堆积图
plotdata <- as.data.frame(cellProp)
colnames(plotdata) <- c("Celltype","Group","Prop")
plotdata$Group <- factor(plotdata$Group,levels=gp_ord)
plotdata$Celltype <- factor(plotdata$Celltype,levels=ct_ord)

p <- ggplot(plotdata, aes(Group, weight = Prop, fill = Celltype)) +
            geom_bar(color = "black", width = .7, position = 'fill') +
            labs( y = 'Relative abundance (%)') +
            #scale_fill_brewer(palette = "Set3")+
            scale_fill_manual(values = ct_cols)+
            scale_y_continuous(expand = c(0,0)) +
            theme_classic() +
            theme(panel.background = element_blank(),
            panel.grid = element_blank(), 
            axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.8,size = 10))

pdf(file.path(outdir1,"Barplot_CelltypeinGroup.pdf"),width=6,height=5)
print(p)
dev.off()

#### 输出每个注释后细胞类型的topn差异基因：
count <- table(as.character(sub$Celltype_mapping)) %>% as.data.frame()
clusters <- levels(sub)
print(clusters)
topGenePC <- c()
topGeneFC <- c()
### 这里前面注释掉，前面已经进行了SCT标准化：
for (ct in clusters) {
    print(ct)
    DEGs <- FindMarkers(sub, ident.1 = ct, verbose = FALSE,logfc.threshold = 0.25,min.cells.group = 5)
    #data <- DEGs[order(DEGs$avg_log2FC,decreasing=T),]
    DEGs$pct1_pct2 <- DEGs$pct.1-DEGs$pct.2
    data <- DEGs[order(DEGs$pct1_pct2,decreasing=T),]
    write.csv(data,file.path(outdir1,paste0("Cluster",ct,"_DEGs_all.csv")),row.names=T)
    data2 <- data[which(data$p_val_adj < 0.05),]
    data2 <- data2[which(!grepl("^mt-",rownames(data2),ignore.case = TRUE) & 
                        !grepl("^Gm\\d",rownames(data2),ignore.case = TRUE) &
                        !grepl("\\dRik$",rownames(data2),ignore.case = TRUE)),] # &
                        #!grepl("^AC\\d",rownames(data2),ignore.case = TRUE) &
                        #!grepl("^AL\\d",rownames(data2),ignore.case = TRUE) &
                        #!grepl("^RP\\d",rownames(data2),ignore.case = TRUE)),]
    write.csv(data2,file.path(outdir1,paste0("Cluster",ct,"_DEGs_Clean.csv")),row.names=T)
    tpgPC <- rownames(data2)[1:topn]
    print(tpgPC)
    topGenePC <- c(topGenePC,tpgPC)
    data3 <- data2[order(data2$avg_log2FC,decreasing=T),]
    tpgFC <- rownames(data3)[1:topn]
    print(tpgFC)
    topGeneFC <- c(topGeneFC,tpgFC)
}
topGenePC <- unique(topGenePC)
avg_dataPC <- AverageExpression(sub,features=topGenePC,assays = assay)
plotdatPC <- avg_dataPC[[assay]]
write.csv(plotdatPC,file.path(outdir1,paste0("Top",topn,"_DEGs_plot_pct.csv")),row.names=T)

pdf(file.path(outdir1,paste0("Top",topn,"_DotPlot_pct_flip.pdf")),height=16,width=length(levels(sub))*0.8+2)
print(DotPlot(sub,features=topGenePC, assay = assay) + 
            coord_flip()+ #翻转
            theme(panel.grid = element_blank(), 
                    axis.text.x=element_text(angle = 90, hjust = 0.5,vjust=0.5)) + #轴标签
            labs(x=NULL,y=NULL) + 
            guides(size = guide_legend("Percent Expression") ) + #legend
            scale_color_gradientn(colours = c('#330066','#336699','#66CC66','#FFCC33'))
            )
dev.off()

pdf(file.path(outdir1,paste0("Top",topn,"_DotPlot_pct.pdf")),height=length(levels(sub))*0.6+2,width=16)
print(DotPlot(sub,features=topGenePC, assay = assay) + 
            theme(panel.grid = element_blank(), 
                    axis.text.x=element_text(angle = 90, hjust = 0.5,vjust=0.5)) + #轴标签
            labs(x=NULL,y=NULL) + 
            guides(size = guide_legend("Percent Expression") ) + #legend
            scale_color_gradientn(colours = c('#330066','#336699','#66CC66','#FFCC33'))
            )
dev.off()

sub <- ScaleData(sub,features=topGenePC,assay = assay, model.use="negbinom")
#pdf(file.path(outdir1,paste0("Top",topn,"_Heatmap_avrg_pct.pdf")),height=20,width=10)
#print(pheatmap(plotdatPC,cluster_cols=TRUE,cluster_rows=TRUE,scale="row"))
#dev.off()
# https://zhuanlan.zhihu.com/p/652376957
pdf(file.path(outdir1,paste0("Top",topn,"_Heatmap_pct.pdf")),height=16,width=10)
print(DoHeatmap(sub,
        features=topGenePC,
        size = 4,
        assay = assay,
        angle = 90,
        group.colors = ct_cols,
        raster = FALSE,
        draw.lines = FALSE)+ #设置组别颜色
    scale_fill_gradientn(colors = c("navy","white","firebrick3")))
dev.off()

topGeneFC <- unique(topGeneFC)
avg_dataFC <- AverageExpression(sub,features=topGeneFC,assays = assay)
plotdatFC <- avg_dataFC[[assay]]
write.csv(plotdatFC,file.path(outdir1,paste0("Top",topn,"_DEGs_plot_lgFC.csv")),row.names=T)

pdf(file.path(outdir1,paste0("Top",topn,"_DotPlot_lgFC_flip.pdf")),height=16,width=length(levels(sub))*0.5+2)
print(DotPlot(sub,features=topGeneFC, assay = assay) +
            coord_flip() + #翻转
            theme(panel.grid = element_blank(), 
                    axis.text.x=element_text(angle = 90, hjust = 0.5,vjust=0.5)) + #轴标签
            labs(x=NULL,y=NULL) + 
            guides(size = guide_legend("Percent Expression") ) + #legend
            scale_color_gradientn(colours = c('#330066','#336699','#66CC66','#FFCC33')))
dev.off()

pdf(file.path(outdir1,paste0("Top",topn,"_DotPlot_lgFC.pdf")),height=length(levels(sub))*0.6+1,width=16)
print(DotPlot(sub,features=topGeneFC, assay = assay) +
            theme(panel.grid = element_blank(), 
                    axis.text.x=element_text(angle = 90, hjust = 0.5,vjust=0.5)) + #轴标签
            labs(x=NULL,y=NULL) + 
            guides(size = guide_legend("Percent Expression") ) + #legend
            scale_color_gradientn(colours = c('#330066','#336699','#66CC66','#FFCC33')))
dev.off()

sub <- ScaleData(sub,features=topGeneFC,assay = assay, model.use="negbinom")
#pdf(file.path(outdir1,paste0("Top",topn,"_Heatmap_avrg_lgFC.pdf")),height=20,width=10)
#print(pheatmap(plotdatFC,cluster_cols=TRUE,cluster_rows=TRUE,scale="row"))
#dev.off()
pdf(file.path(outdir1,paste0("Top",topn,"_Heatmap_lgFC.pdf")),height=16,width=10)
print(DoHeatmap(sub,
        features=topGeneFC,
        size = 4,
        assay = assay,
        angle = 90,
        group.colors = ct_cols,
        raster = FALSE,
        draw.lines = FALSE)+ #设置组别颜色
    scale_fill_gradientn(colors = c("navy","white","firebrick3")))
dev.off()