From a1b3d98f16686648a84e7d3b093e8e0d55512571 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Mon, 22 Dec 2025 11:18:37 +0800 Subject: [PATCH 01/19] merge and rename fq.gz file --- script/merge_and_rename_fq.sh | 142 ++++++++++++++++++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 script/merge_and_rename_fq.sh diff --git a/script/merge_and_rename_fq.sh b/script/merge_and_rename_fq.sh new file mode 100644 index 0000000..3fcacde --- /dev/null +++ b/script/merge_and_rename_fq.sh @@ -0,0 +1,142 @@ +#!/bin/bash +set -euo pipefail + +# ====================== +# Help document +# ====================== +usage() { +cat <_1.fq.gz + _2.fq.gz + + Mode "dir": + raw_dir/OldSampleID/*.fq.gz + fastq names can be arbitrary but must contain R1 / R2 or 1 / 2 + +Usage: + bash merge_and_rename_fq.sh -m metadata.txt -i raw_fq_dir -o output_dir [-M name|dir] + +Options: + -m Metadata file (tab- or space-delimited, first line is header) + -i Directory containing raw fastq files (recursive search) + -o Output directory for renamed fastq symbolic links + -M Match mode: + name (default) : match fastq file names + dir : match by directory name + -h Show this help message and exit + +Output: + output_dir/ + ├── SampleID_R1.fq.gz + └── SampleID_R2.fq.gz + +EOF +exit 0 +} + +# ====================== +# Parse arguments +# ====================== +mode="name" + +while getopts "m:i:o:M:h" opt; do + case $opt in + m) metadata=$OPTARG ;; + i) raw_dir=$OPTARG ;; + o) out_dir=$OPTARG ;; + M) mode=$OPTARG ;; + h) usage ;; + *) usage ;; + esac +done + +# ====================== +# Check arguments +# ====================== +[[ -z "${metadata:-}" || -z "${raw_dir:-}" || -z "${out_dir:-}" ]] && usage + +mkdir -p "$out_dir" + +echo "[INFO] Metadata : $metadata" +echo "[INFO] Raw dir : $raw_dir" +echo "[INFO] Out dir : $out_dir" +echo "[INFO] Mode : $mode" +echo "[INFO] Start processing..." +echo + +# ====================== +# Main loop +# ====================== +tail -n +2 "$metadata" | while read -r new old; do + # 去掉 Windows 回车 + new=$(echo "$new" | tr -d '\r') + old=$(echo "$old" | tr -d '\r') + + [[ -z "$new" || -z "$old" ]] && continue + + echo "[INFO] Processing: $old → $new" + + R1="" + R2="" + + if [[ "$mode" == "name" ]]; then + # ========= 按 fastq 文件名匹配 ========= + R1=$(find "$raw_dir" -type f -name "${old}_1.fq.gz" | head -n 1) + R2=$(find "$raw_dir" -type f -name "${old}_2.fq.gz" | head -n 1) + + elif [[ "$mode" == "dir" ]]; then + # ========= 按目录名匹配 ========= + sample_dir=$(find "$raw_dir" -type d -name "$old" | head -n 1) + + if [[ -z "$sample_dir" ]]; then + echo "[WARNING] Sample directory not found: $old" + continue + fi + + R1=$(find "$sample_dir" -type f \( \ + -name "*_R1_*.fq.gz" -o \ + -name "*_R1.fq.gz" -o \ + -name "*_1.fq.gz" \ + \) | head -n 1) + + R2=$(find "$sample_dir" -type f \( \ + -name "*_R2_*.fq.gz" -o \ + -name "*_R2.fq.gz" -o \ + -name "*_2.fq.gz" \ + \) | head -n 1) + else + echo "[ERROR] Unknown mode: $mode" + exit 1 + fi + + if [[ -z "$R1" || -z "$R2" ]]; then + echo "[WARNING] Missing fastq files for $old, skip" + continue + fi + + ln -sf "$(realpath "$R1")" "${out_dir}/${new}_R1.fq.gz" + ln -sf "$(realpath "$R2")" "${out_dir}/${new}_R2.fq.gz" + + echo "[OK] ${new}_R1.fq.gz ${new}_R2.fq.gz" +done + +echo +echo "[DONE] All samples processed successfully." From ca8c20c173c8c7ff791d59ad01b3afddb3267f5c Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Thu, 25 Dec 2025 15:03:31 +0800 Subject: [PATCH 02/19] venn infor extract and bracken batch Boxplots --- script/BrackenBatchBoxplots.R | 181 +++++++++++++++++++++++++++++ script/venn_info_extract.R | 213 ++++++++++++++++++++++++++++++++++ 2 files changed, 394 insertions(+) create mode 100644 script/BrackenBatchBoxplots.R create mode 100644 script/venn_info_extract.R diff --git a/script/BrackenBatchBoxplots.R b/script/BrackenBatchBoxplots.R new file mode 100644 index 0000000..e8dce13 --- /dev/null +++ b/script/BrackenBatchBoxplots.R @@ -0,0 +1,181 @@ + +#!/usr/bin/env Rscript + +# ============================================================ +# plot_species_boxplot.R +# ============================================================ + +# =========================== +# 1. 定义需要的包 +# =========================== +packages <- c( + "optparse", + "data.table", + "dplyr", + "ggplot2", + "ggpubr", + "cowplot" +) + +# =========================== +# 2. 检查是否安装，没有就安装 +# =========================== +for (pkg in packages) { + if (!requireNamespace(pkg, quietly = TRUE)) { + message(paste0("⚠️ Installing package: ", pkg)) + install.packages(pkg, repos = "https://cran.r-project.org") + } +} + +# =========================== +# 3. 加载包 +# =========================== +suppressPackageStartupMessages({ + lapply(packages, library, character.only = TRUE) +}) + + +# ====================== +# 1. 命令行参数 +# ====================== +option_list <- list( + make_option(c("-e", "--exp"), type = "character", + help = "Abundance table (e.g. bracken.S.txt)", metavar = "file"), + make_option(c("-m", "--meta"), type = "character", + help = "Metadata table", metavar = "file"), + make_option(c("-s", "--species"), type = "character", + help = "Species list, comma-separated (e.g. \"A,B,C\")"), + make_option(c("--groups"), type = "character", + default = "VfNnVa5Rs,ViNnVa5Rs,VfNnVb5Rs,ViNnVb5Rs", + help = "Group order, comma-separated"), + make_option(c("-o", "--out"), type = "character", + default = "Species_boxplot.pdf", + help = "Output PDF [default: %default]"), + make_option(c("--sample_col"), type = "character", + default = "SampleID", + help = "Sample column name in metadata [default: %default]"), + make_option(c("--group_col"), type = "character", + default = "GroupID", + help = "Group column name in metadata [default: %default]"), + make_option(c("--ncol"), type = "integer", + default = 4, + help = "Number of columns in plot layout [default: %default]"), + make_option(c("--pdf_width"), type = "double", default = 14, + help = "PDF width [default: %default], unit = cm"), + make_option(c("--pdf_height"), type = "double", default = 10, + help = "PDF height [default: %default], unit = cm") +) + +opt <- parse_args(OptionParser(option_list = option_list)) + +if (is.null(opt$exp) || is.null(opt$meta) || is.null(opt$species)) { + stop("❌ Required arguments: -e, -m, --species\n", call. = FALSE) +} + +# ====================== +# 2. 解析 species（命令行输入） +# ====================== +species <- unlist(strsplit(opt$species, ",")) +species <- trimws(species) +if (length(species) == 0) stop("❌ No species provided.", call. = FALSE) + +# ====================== +# 3. 解析 groups（命令行输入） +# ====================== +group_levels <- unlist(strsplit(opt$groups, ",")) +group_levels <- trimws(group_levels) + +# ====================== +# 4. 读入丰度表 +# ====================== +Exp <- fread(opt$exp, sep = "\t", fill = TRUE, data.table = FALSE) +if (!"Taxonomy" %in% colnames(Exp)) { + stop("❌ Abundance table must contain 'Taxonomy' column.", call. = FALSE) +} +rownames(Exp) <- Exp$Taxonomy +Exp <- Exp[, -1] + +# ====================== +# 5. 每列相对丰度 +# ====================== +Exp <- Exp %>% + mutate(across(where(is.numeric), ~ .x / sum(.x))) + +Exp2 <- as.data.frame(t(Exp)) +Exp2$SampleID <- rownames(Exp2) +colnames(Exp2) <- gsub("\\.", " ", colnames(Exp2)) +Exp2 <- Exp2[, c(ncol(Exp2), 1:(ncol(Exp2) - 1))] + +# ====================== +# 6. 物种存在性检查 +# ====================== +missing_sp <- setdiff(species, colnames(Exp2)) +if (length(missing_sp) > 0) { + warning("⚠️ Species not found:\n", paste(missing_sp, collapse = ", ")) +} +species <- intersect(species, colnames(Exp2)) +if (length(species) == 0) { + stop("❌ None of the species found in abundance table.", call. = FALSE) +} +Exp_plot <- Exp2[, c("SampleID", species)] + +# ====================== +# 7. 读入 metadata 并合并 +# ====================== +meta <- read.table(opt$meta, header = TRUE, sep = "\t", stringsAsFactors = FALSE) +if (!opt$sample_col %in% colnames(meta) || !opt$group_col %in% colnames(meta)) { + stop("❌ sample_col or group_col not found in metadata.", call. = FALSE) +} + +Exp_plot <- Exp_plot[Exp_plot$SampleID %in% meta[[opt$sample_col]], ] +Exp_plot <- Exp_plot[match(meta[[opt$sample_col]], Exp_plot$SampleID), ] + +Exp_plot$group <- meta[[opt$group_col]] +Exp_plot$group <- factor(Exp_plot$group, levels = group_levels) +Exp_plot <- Exp_plot[!is.na(Exp_plot$group), ] + +# ====================== +# 8. 作图 +# ====================== +colors <- c("#EB746A", "#7AA82C", "#1EB5B8", "#A07DB7") +comparisons <- combn(group_levels, 2, simplify = FALSE) + +plist <- list() +for (i in seq_along(species)) { + bar_tmp <- Exp_plot[, c(species[i], "group")] + colnames(bar_tmp) <- c("Relative abundance", "group") + + p <- ggboxplot( + bar_tmp, + x = "group", + y = "Relative abundance", + color = "group", + add = "jitter", + palette = colors + ) + + ggtitle(species[i]) + + theme_bw() + + theme( + axis.title = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1), + plot.title = element_text(hjust = 0.5, face = "bold"), + legend.position = "none" + ) + + stat_compare_means( + method = "wilcox.test", + comparisons = comparisons, + label = "p.signif" + ) + + plist[[i]] <- p +} + +# ====================== +# 9. 输出 PDF +# ====================== +#pdf(opt$out, width = 14, height = 10) +pdf(opt$out, width = opt$pdf_width, height = opt$pdf_height) +plot_grid(plotlist = plist, ncol = opt$ncol) +dev.off() + +cat("✅ Done! Output:", opt$out, "\n") diff --git a/script/venn_info_extract.R b/script/venn_info_extract.R new file mode 100644 index 0000000..1879d04 --- /dev/null +++ b/script/venn_info_extract.R @@ -0,0 +1,213 @@ +#!/usr/bin/env Rscript + +# ------------------------------- +# 自动检查并安装缺失包 +# ------------------------------- +packages <- c("optparse", "VennDiagram", "gtools", "grid") + +for (pkg in packages) { + if (!requireNamespace(pkg, quietly = TRUE)) { + install.packages(pkg, repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN") + } + suppressPackageStartupMessages(library(pkg, character.only = TRUE)) +} + + +############################ +## 参数定义 +############################ + +option_list <- list( + make_option(c("-t", "--taxonomy"), + type = "character", + help = "taxonomy.spf file"), + + make_option(c("-m", "--metadata"), + type = "character", + help = "metadata.txt file"), + + make_option(c("-g", "--groups"), + type = "character", + help = "GroupID list, comma separated"), + + make_option(c("-k", "--k_common"), + type = "character", + default = "2", + help = "k groups for common species, e.g. 2 or 2,3"), + + make_option(c("-o", "--outdir"), + type = "character", + default = "results", + help = "Output directory"), + + make_option(c("--max_venn"), + type = "integer", + default = 4, + help = "Max group number for Venn plot") +) + +opt <- parse_args(OptionParser(option_list = option_list)) + +############################ +## 参数解析 +############################ + +groups <- unlist(strsplit(opt$groups, ",")) +k_common <- as.integer(unlist(strsplit(opt$k_common, ","))) + +dir.create(opt$outdir, showWarnings = FALSE, recursive = TRUE) + +############################ +## 读入数据 +############################ + +data_s <- read.table(opt$taxonomy, + header = TRUE, + sep = "\t", + stringsAsFactors = FALSE) + +metadata <- read.table(opt$metadata, + header = TRUE, + sep = "\t", + stringsAsFactors = FALSE) + +############################ +## 筛选 metadata +############################ + +sub_metadata <- metadata[metadata$GroupID %in% groups, ] +rownames(sub_metadata) <- sub_metadata$SampleID + +############################ +## 物种矩阵处理 +############################ + +# data_s2 <- data_s[, c(7, 9:ncol(data_s))] + +if ("Species" %in% colnames(data_s)) { + # Species 数据 + data_s2 <- data_s[, c(7, 9:ncol(data_s))] +} else { + # KO 数据 + data_s2 <- data_s +} + +# 统一第一列列名为 Feature +colnames(data_s2)[1] <- "Feature" + +# data_s3 <- aggregate(. ~ Species, +# data = data_s2, +# FUN = sum) + +data_s3 <- aggregate(. ~ Feature, + data = data_s2, + FUN = sum) + +#rownames(data_s3) <- data_s3$Species +rownames(data_s3) <- data_s3$Feature +data_s3 <- data_s3[, -1] + +data_s4 <- data_s3[, colnames(data_s3) %in% rownames(sub_metadata)] + +############################ +## 分组求平均 +############################ + +data_s4_t <- as.data.frame(t(data_s4)) +data_s4_t$GroupID <- sub_metadata[rownames(data_s4_t), "GroupID"] + +group_mean <- aggregate(. ~ GroupID, + data = data_s4_t, + FUN = mean) + +rownames(group_mean) <- group_mean$GroupID +group_mean <- group_mean[, -1] + +############################ +## 二值化 +############################ + +bin_mat <- as.data.frame(t(group_mean)) +bin_mat[] <- as.integer(bin_mat > 0) + +############################ +## Venn 图（<= max_venn） +############################ + +if (ncol(bin_mat) <= opt$max_venn) { + + venn_list <- lapply(colnames(bin_mat), function(g) { + rownames(bin_mat)[bin_mat[, g] == 1] + }) + names(venn_list) <- colnames(bin_mat) + + venn.plot <- venn.diagram( + x = venn_list, + filename = NULL, + fill = c("#4DBBD5", "#E64B35", "#00A087", "#3C5488")[1:ncol(bin_mat)], + alpha = 0.5, + cex = 1.3, + cat.cex = 1.2, + cat.pos = seq(-30, 30, length.out = ncol(bin_mat)), + cat.dist = rep(0.05, ncol(bin_mat)), + main = "Feature Presence Venn" + ) + + pdf(file.path(opt$outdir, "Venn.pdf"), width = 6, height = 6) + grid.draw(venn.plot) + dev.off() +} + +############################ +## 特有 & 共有物种 +############################ + +group_species <- lapply(colnames(bin_mat), function(g) { + rownames(bin_mat)[bin_mat[, g] == 1] +}) +names(group_species) <- colnames(bin_mat) + +## 特有物种 +for (g in names(group_species)) { + + others <- setdiff(names(group_species), g) + + unique_sp <- setdiff( + group_species[[g]], + unique(unlist(group_species[others])) + ) + + write.csv( + data.frame(Feature = unique_sp), + file = file.path(opt$outdir, paste0("Unique_", g, ".csv")), + row.names = FALSE + ) +} + +## 指定 k 组共有 +for (k in k_common) { + + if (k > length(group_species)) next + + combs <- combinations( + n = length(group_species), + r = k, + v = names(group_species) + ) + + for (i in 1:nrow(combs)) { + + gs <- combs[i, ] + + common_sp <- Reduce(intersect, group_species[gs]) + + write.csv( + data.frame(Feature = common_sp), + file = file.path( + opt$outdir, + paste0("Common_", paste(gs, collapse = "_"), ".csv") + ), + row.names = FALSE + ) + } +} From 81e0ee235d2cd1c7a48b142513038dae8add6fef Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Fri, 23 Jan 2026 20:05:45 +0800 Subject: [PATCH 03/19] Add files via upload --- script/stack_data_prepare.sh | 112 ++++++++++++++++++++++++----------- 1 file changed, 77 insertions(+), 35 deletions(-) diff --git a/script/stack_data_prepare.sh b/script/stack_data_prepare.sh index bf4b856..1943a0a 100644 --- a/script/stack_data_prepare.sh +++ b/script/stack_data_prepare.sh @@ -1,35 +1,77 @@ -#!/bin/bash - -levels=("Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species" "Strain") -input_file="metaphlan4/taxonomy.spf" -header=$(head -n 1 "$input_file") -sample_start=9 - -for i in {1..8}; do - echo "Processing: ${levels[$i-1]}" - awk -v col=$i -v sample_start="$sample_start" -v header="$header" 'BEGIN {OFS="\t"} - NR == 1 { - printf "Taxonomy\t" - for (j = sample_start; j <= NF; j++) { - printf "%s%s", (j == sample_start ? "" : OFS), $j - } - print ""; - next; - } - { - key = $col - for (j = sample_start; j <= NF; j++) { - sums[key, j] += $j - } - keys[key] = 1; - } - END { - for (key in keys) { - printf "%s", key - for (j = sample_start; j <= NF; j++) { - printf "\t%.5f", sums[key, j] - } - print "" - } - }' "$input_file" > "metaphlan4/${levels[$i-1]}.txt" -done + +#!/bin/bash +# set -euo pipefail + +# ====================== +# 参数解析 +# ====================== +usage() { + echo "Usage: $0 -i -o " + echo "" + echo " -i MetaPhlAn taxonomy.spf file" + echo " -o Output directory" + exit 1 +} + +while getopts "i:o:" opt; do + case ${opt} in + i ) input_file=$OPTARG ;; + o ) outdir=$OPTARG ;; + * ) usage ;; + esac +done + +# 参数检查 +if [[ -z "${input_file:-}" || -z "${outdir:-}" ]]; then + usage +fi + +if [[ ! -f "$input_file" ]]; then + echo "ERROR: input file not found: $input_file" + exit 1 +fi + +mkdir -p "$outdir" + +# ====================== +# 主逻辑 +# ====================== +levels=("Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species" "Strain") +sample_start=9 + +for i in {1..8}; do + level=${levels[$i-1]} + echo "Processing: $level" + + awk -v col="$i" -v sample_start="$sample_start" 'BEGIN {OFS="\t"} + NR == 1 { + printf "Taxonomy" + for (j = sample_start; j <= NF; j++) { + printf "\t%s", $j + } + print "" + next + } + { + key = $col + for (j = sample_start; j <= NF; j++) { + sums[key, j] += $j + } + keys[key] = 1 + } + END { + for (key in keys) { + printf "%s", key + for (j = sample_start; j <= NF; j++) { + printf "\t%.5f", sums[key, j] + } + print "" + } + } + ' "$input_file" > "${outdir}/${level}.txt" +done + +echo "All levels processed successfully." + + + From b4ee0108c16a46e9f04dd9f7ba346166f1652be3 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Fri, 22 May 2026 17:02:44 +0800 Subject: [PATCH 04/19] Delete script/compare_valcano.R --- script/compare_valcano.R | 190 --------------------------------------- 1 file changed, 190 deletions(-) delete mode 100644 script/compare_valcano.R diff --git a/script/compare_valcano.R b/script/compare_valcano.R deleted file mode 100644 index c0b0ef5..0000000 --- a/script/compare_valcano.R +++ /dev/null @@ -1,190 +0,0 @@ -#!/usr/bin/env Rscript - -# Copyright 2024-2026 Defeng Bai - -# If used this script, please cited: -# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 - -# 手动运行脚本请，需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 - -# 更新 -# 2024/11/12：增加利用MaAsLin2差异分析结果绘制火山图代码 -# 2025/11/27：规范代码 - -# 1.1 程序功能描述和主要步骤 - -# 程序功能：火山图展示组间差异 -# Functions: Species difference illustrated using volcano plot - -options(warn = -1) # Turn off warning - - -## 设置输入输出文件和参数 - -# 修改下面`default=`后面的文件和参数。 -# -# 输入文件为MaAsLin2组间差异分析得到的两个结果文件(metaphlan4/MaAsLin2_overall_difference.csv)+趋势信息(metaphlan4/MaAsLin2_enriched_depleted.csv) -# -# 输入文件"-i", "--input"，metaphlan4/MaAsLin2_overall_difference.csv; MaAsLin2差异分析结果汇总信息； -# -# 实验设计"-d", "--trend"，默认`metaphlan4/MaAsLin2_enriched_depleted.csv`，差异分析结果趋势信息； -# -# 图片宽"-w", "--width"，默认89 mm，根据图像布局可适当增大或缩小 -# -# 图片高"-e", "--height"，默认59 mm，根据图像布局可适当增大或缩小 -# -# 分组列名"-o", "--output"，默认为输出目录； - - -# 1.2 依赖包安装 - -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -a = rownames(installed.packages()) - -# install CRAN -install_CRAN <- c("ggplot2", "BiocManager", "optparse", "dplyr", "magrittr") -for (i in install_CRAN) { - if (!i %in% a) - install.packages(i, repos = site) - require(i,character.only=T) - a = rownames(installed.packages()) -} - -# install bioconductor -install_bioc <- c("ggplot2", "multcompView") -for (i in install_bioc) { - if (!i %in% a) - BiocManager::install(i, update = F) - a = rownames(installed.packages()) -} - -# install github -if (!"amplicon" %in% a){ - devtools::install_github("microbiota/amplicon") -} - -if (!"ggvolcano" %in% a){ - devtools::install_github("BioSenior/ggvolcano", force = TRUE) -} - - -# 1.3 解析命令行 -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -# site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="metaphlan4/MaAsLin2_overall_difference.csv", - help="Unfiltered OTU table [default %default]"), - make_option(c("-d", "--trend"), type="character", default="metaphlan4/MaAsLin2_enriched_depleted.csv", - help="metadata file or metadata [default %default]"), - make_option(c("-o", "--output"), type="character", default="metaphlan4/", - help="Output quantile value for filter feature table [default %default]"), - make_option(c("-w", "--width"), type="numeric", default=89, - help="Width of figure [default %default]"), - make_option(c("-e", "--height"), type="numeric", default=59, - help="Height of figure [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# # Install related packages -# if (FALSE){ -# source("https://bioconductor.org/biocLite.R") -# biocLite(c("ggplot2","magrittr","dplyr")) -# } -# # load related packages -# suppressWarnings(suppressMessages(library("ggplot2"))) -# suppressWarnings(suppressMessages(library("dplyr"))) -# suppressWarnings(suppressMessages(library("magrittr"))) - -# 2. 依赖关系检查、安装和加载 - -# 依赖包列表 -package_list <- c( - "ggplot2", "BiocManager", "optparse", "dplyr", "magrittr", "ggvolcano" -) - -# 批量安装和加载 -for (p in package_list) { - # 如果未安装，则安装 - if (!requireNamespace(p, quietly = TRUE)) { - install.packages(p, repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/") - } - - # 批量加载，抑制警告和消息 - suppressWarnings( - suppressMessages( - library(p, character.only = TRUE) - ) - ) -} - - -# 读取输入文件 - -# 读取趋势表格 -data_species01 <- read.csv(opts$trend, row.names = 1) -data_species01$species <- rownames(data_species01) -data_species01 <- data_species01[, -c(1:4)] - -# 读取差异结果汇总表 -data_MWAS <- read.csv(opts$input, row.names = 1) -data_MWAS$species <- rownames(data_MWAS) - -data_species02 <- merge(data_species01, data_MWAS, by = "species") - -# devtools::install_github("BioSenior/ggvolcano", force = TRUE) -# library(ggVolcano) -data_vol <- data_species02 -data_vol = as.data.frame(data_vol) - -data_vol2 <- data_vol -data_vol2$padj2 <- -log10(data_vol2$FDR) - -# logFC = 0.5 -# P.Value = 0.05 -# library(ggplot2) -p_volcano1 <- ggplot(data = data_vol2, aes(x = -Beta, y = padj2)) + - geom_point(alpha = 0.4, size = 2.0, aes(color = NPC.association)) + - ylab("-log10(Pvalue)") + - scale_color_manual(values = c("#74add1","#a60026", "grey")) + - geom_vline(xintercept = 0, lty = 4, col = "black", lwd = 0.4) + - geom_hline(yintercept = -log10(0.75), lty = 4, col = "black", lwd = 0.4) + - labs(x = "Coef. (by MaAsLin2)", y= bquote(atop(-Log[10]~italic(FDR))))+ - theme_bw() - -# add labels -# library(dplyr) -# select top 5 enriched species -up_data1 <- filter(data_vol2, data_vol2$NPC.association == "enriched") -up_data2 <- arrange(up_data1, desc(up_data1$padj2)) -up_data_5 <- up_data2[1:5, ] - -# select top 25 depleted species -down_data1 <- filter(data_vol2, data_vol2$NPC.association == "depleted") -down_data2 <- arrange(down_data1, desc(down_data1$padj2)) -down_data_5 <- down_data2[1:5, ] - -# using geom_text_repel() to add labels -library(ggrepel) -p_volcano2 <- p_volcano1 + - geom_text_repel(data = up_data_5, aes(x = -Beta, - y = padj2, - label = up_data_5$Feature), size = 2, max.overlaps = 20) + - geom_text_repel(data = down_data_5, aes(x = -Beta, - y = padj2, - label = down_data_5$Feature), size = 2, max.overlaps = 20)+ - theme(legend.position = c(0.84, 0.85),panel.grid = element_blank()) -#ggsave(paste0(opts$output,"Difference_valcano.pdf"), p_volcano2, width=59 * 1.5, height=50 * 1.5, unit='mm') -ggsave(paste0(opts$output,"Difference_valcano.pdf"), p_volcano2, width=opts$width, height=opts$height, unit='mm') From a601d04939b20fe8818065ab1a111082965b98ad Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Fri, 22 May 2026 17:03:47 +0800 Subject: [PATCH 05/19] Update scripts Update compare_volcano.R and compare_MaAsLin2.R scripts --- script/compare_MaAsLin2.R | 182 ++++++-------------- script/compare_volcano.R | 342 +++++++++++++++++++++----------------- 2 files changed, 243 insertions(+), 281 deletions(-) diff --git a/script/compare_MaAsLin2.R b/script/compare_MaAsLin2.R index d2fad4c..b00f025 100644 --- a/script/compare_MaAsLin2.R +++ b/script/compare_MaAsLin2.R @@ -10,6 +10,7 @@ # 更新 # 2024/11/12：增加MaAsLin2方法用于组间差异比较 # 2025/11/27：规范代码 +# 2026/5/22：更新软件包，规范代码 # Clean enviroment object @@ -54,45 +55,16 @@ options(warn = -1) # Turn off warning # 分组列名"-o", "--output"，默认为输出目录，输出差异分析结果文件; -# 1.2 依赖包安装 - -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -a = rownames(installed.packages()) - -# install CRAN -install_CRAN <- c("optparse", "dplyr", "reshape2", "readxl", "tibble", "openxlsx", - "foreach", "data.table", "gridExtra", "scales", "ggplot2", "ggh4x", - "ggfortify", "ggvenn", "ggrepel", "vegan", "pairwiseCI", "vcd", - "igraph", "sampling", "CVXR", "DescTools", "phyloseq") -for (i in install_CRAN) { - if (!i %in% a) - install.packages(i, repos = site) -} - -# install bioconductor -install_bioc <- c( "phyloseq", "ANCOMBC", "Maaslin2") -for (i in install_bioc) { - if (!i %in% a) - BiocManager::install(i, update = F) # , site_repository=site - a = rownames(installed.packages()) -} - -# install github -if (!"amplicon" %in% a){ - devtools::install_github("microbiota/amplicon") -} - - -# 1.3 解析命令行 +# 1. 解析命令行 # 解析参数-h显示帮助信息 suppressMessages( library("optparse") ) if (TRUE){ option_list = list( - make_option(c("-i", "--input"), type="character", default="metaphlan4/Species.txt", + make_option(c("-i", "--input"), type="character", default="result/metaphlan4/Species.txt", help="Unfiltered OTU table [default %default]"), - make_option(c("-g", "--metadata"), type="character", default="metadata.txt", + make_option(c("-g", "--metadata"), type="character", default="result/metadata.txt", help="metadata file or metadata [default %default]"), make_option(c("-a", "--min_abundance"), type="numeric", default="0", help="Minimum abundance used for analysis [default %default]"), @@ -112,7 +84,7 @@ if (TRUE){ help="Fixed effects [default %default]"), make_option(c("-c", "--correction"), type="character", default="BH", help="Fixed effects [default %default]"), - make_option(c("-o", "--output"), type="character", default="metaphlan4/", + make_option(c("-o", "--output"), type="character", default="result/metaphlan4/", help="Output quantile value for filter feature table [default %default]") ) opts = parse_args(OptionParser(option_list=option_list)) @@ -123,93 +95,58 @@ if (TRUE){ # 2. 依赖关系检查、安装和加载 -suppress <- function(x){invisible(capture.output(suppressMessages(suppressWarnings(x))))} +cran_repo <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 依赖包列表 -package_list <- c( - "optparse","dplyr", "reshape2", "readxl", "tibble", "openxlsx", - "foreach", "data.table", "gridExtra", "scales", "ggplot2", "ggh4x", - "ggfortify", "ggvenn", "ggrepel", "vegan", "pairwiseCI", "vcd", - "igraph", "sampling", "CVXR", "DescTools", "phyloseq" +# CRAN 包 +cran_pkgs <- c( + "optparse", "dplyr", "reshape2", "readxl", "tibble", "openxlsx", + "foreach", "data.table", "gridExtra", "scales", "ggplot2", "ggh4x", + "ggfortify", "ggvenn", "ggrepel", "vegan", "pairwiseCI", "vcd", + "igraph", "sampling", "CVXR", "DescTools", "devtools" ) -# 批量安装和加载 -for (p in package_list) { - # 如果未安装，则安装 - if (!requireNamespace(p, quietly = TRUE)) { - install.packages(p, repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/") +# Bioconductor 包 +bioc_pkgs <- c("phyloseq", "Maaslin2") + +# GitHub 包 +github_pkgs <- c("amplicon" = "microbiota/amplicon") + + +# ---------- 安装 CRAN ---------- +for (pkg in cran_pkgs) { + if (!requireNamespace(pkg, quietly = TRUE)) { + install.packages(pkg, repos = cran_repo) + } +} + +# ---------- 安装 Bioconductor ---------- +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager", repos = cran_repo) +} + +for (pkg in bioc_pkgs) { + if (!requireNamespace(pkg, quietly = TRUE)) { + BiocManager::install(pkg, update = FALSE, ask = FALSE) + } +} + +# ---------- 安装 GitHub ---------- +for (pkg in names(github_pkgs)) { + if (!requireNamespace(pkg, quietly = TRUE)) { + devtools::install_github(github_pkgs[pkg]) } - - # 批量加载，抑制警告和消息 - suppressWarnings( - suppressMessages( - library(p, character.only = TRUE) - ) - ) } +# ---------- 批量加载 ---------- +all_pkgs <- c(cran_pkgs, bioc_pkgs, names(github_pkgs)) -# # Install related packages -# # 基于CRAN安装R包，检测没有则安装 Installing R packages based on CRAN and installing them if they are not detected -# p_list = c("dplyr", "reshape2", "readxl", "tibble", "openxlsx", -# "foreach", "data.table", "gridExtra", "scales", "ggplot2", "ggh4x", -# "ggfortify", "ggvenn", "ggrepel", "vegan", "pairwiseCI", "vcd", -# "igraph", "sampling", "CVXR", "DescTools") -# for(p in p_list){if (!requireNamespace(p)){install.packages(p)} -# library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)} -# -# #### LOAD REQUIRED R PACKAGES #### -# -# # options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor") -# if (!requireNamespace("BiocManager", quietly = TRUE)) -# install.packages("BiocManager", repos = site) -# a = rownames(installed.packages()) -# -# # install CRAN -# install_CRAN <- c("ggplot2", "reshape2", "readxl", "tibble","openxlsx", "foreach", -# "data.table", "gridExtra","scales", "ggh4x", "ggfortify", "ggvenn", -# "ggrepel", "vegan", "pairwiseCI", "vcd", "igraph") -# for (i in install_CRAN) { -# if (!i %in% a) -# install.packages(i, repos = site) -# } -# -# # install bioconductor -# install_bioc <- c( "phyloseq", "ANCOMBC", "Maaslin2") -# for (i in install_bioc) { -# if (!i %in% a) -# BiocManager::install(i, update = F) # , site_repository=site -# a = rownames(installed.packages()) -# } -# -# # install github -# if (!"amplicon" %in% a){ -# devtools::install_github("microbiota/amplicon") -# } - -# suppress <- function(x){invisible(capture.output(suppressMessages(suppressWarnings(x))))} -# suppress(library(dplyr)) -# suppress(library(reshape2)) -# suppress(library(readxl)) -# suppress(library(phyloseq)) -# suppress(library(tibble)) -# suppress(library(openxlsx)) -# suppress(library(foreach)) -# suppress(library(data.table)) -# suppress(library(gridExtra)) -# suppress(library(scales)) -# suppress(library(ggplot2)) -# suppress(library(ggh4x)) -# suppress(library(ggfortify)) -# suppress(library(ggvenn)) -# suppress(library(ggrepel)) -# suppress(library(vegan)) -# suppress(library(pairwiseCI)) -# suppress(library(vcd)) -# # library(CVXR) -# suppress(library(ANCOMBC)) -# suppress(library(Maaslin2)) -# suppress(library(igraph)) +invisible( + lapply(all_pkgs, function(pkg) { + suppressPackageStartupMessages( + library(pkg, character.only = TRUE) + ) + }) +) # 3. 读取输入文件 @@ -263,25 +200,6 @@ library(sampling) sam_data = as.data.frame(sample_data(ps)) common_columns = colnames(sam_data)[colnames(sam_data) %in% colnames(metadata)] input_metadata <- data.frame(sam_data[, common_columns]) -# capt<- capture.output(fits <- suppressWarnings(Maaslin2(input_data, input_metadata, -# output='temp_directory', -# min_prevalence=0.05, -# #min_variance, -# #normalization='NONE', -# normalization='CLR', -# #normalization= 'TSS', -# #transform='LOG', -# transform = "NONE", -# correction = "BH", -# analysis_method='LM', -# max_significance=0.05, -# fixed_effects = c('Group'), -# #standardize=FALSE, -# standardize=TRUE, -# #standardize=FALSE, -# plot_heatmap=TRUE, -# plot_scatter=TRUE -# ))) capt<- capture.output(fits <- suppressWarnings(Maaslin2(input_data, input_metadata, output='temp_directory', diff --git a/script/compare_volcano.R b/script/compare_volcano.R index cbd38f3..dd219fd 100644 --- a/script/compare_volcano.R +++ b/script/compare_volcano.R @@ -1,149 +1,193 @@ -#!/usr/bin/env Rscript - -# Copyright 2016-2023 Yong-Xin Liu - -# If used this script, please cited: -# Jingying Zhang, Yong-Xin Liu, et. al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology 37, 676-684, doi:10.1038/s41587-019-0104-4 (2019). - -# 手动运行脚本请，需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 - -#----1. 参数 Parameters#---- - -#----1.1 功能描述 Function description#---- - -# 程序功能：差异比较火山图绘制 -# Functions: Volcano plot - -options(warn = -1) # Turn off warning - - -## 设置输入输出文件和参数 - -# 修改下面`default=`后面的文件和参数。 -# -# 输入文件为特征表(A-B.txt)+分组信息(metadata.tsv) -# -# 输入文件"-i", "--input"，otutab.txt; 特征表 -# -# 实验设计"-d", "--design"，默认`metadata.tsv`，可手动修改文件位置； -# -# 分组列名"-n", "--group"，默认将metadata.tsv中的Group列作为分组信息，可修改为任意列名； -# -# 分组列名"-c", "--compare_pair"，默认将比较metadata.tsv中的Group列的前两个值，建议手动设定； -# -# 分组列名"-t", "--threhold"，丰度筛选阈值，默认千分之1 -# -# 图片宽"-w", "--width"，默认89 mm，根据图像布局可适当增大或缩小 -# -# 图片高"-e", "--height"，默认59 mm，根据图像布局可适当增大或缩小 - - -#----1.2 参数缺少值 Default values#---- -# 设置清华源加速下载 -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析是否安装，安装并加载 -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="result/compare/KO-WT.txt", - help="Feature table [default %default]"), - make_option(c("-t", "--threshold"), type="numeric", default=0.1, - help="Threshold of relative abundance 0.1% [default %default]"), - make_option(c("-d", "--design"), type="character", default="result/metadata.tsv", - help="Design file or metadata [default %default]"), - make_option(c("-n", "--group"), type="character", default="Group", - help="Group name [default %default]"), - make_option(c("-c", "--compare"), type="character", default="KO-WT", - help="Groups comparison [default %default]"), - make_option(c("-m", "--method"), type="character", default="wilcox", - help="Compare method, alternative edgeR or t.test [default %default]"), - make_option(c("-p", "--pvalue"), type="numeric", default=0.05, - help="Threshold of P-value [default %default]"), - make_option(c("-f", "--fdr"), type="numeric", default=0.1, - help="Threshold of FDR [default %default]"), - make_option(c("-o", "--output"), type="character", default="", - help="Output directory; name according to input [default %default]"), - make_option(c("-w", "--width"), type="numeric", default=89, - help="Figure width [default %default]"), - make_option(c("-e", "--height"), type="numeric", default=59, - help="Figure heidth [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) - # suppressWarnings(dir.create(opts$output)) -} -# 设置输出文件缺省值，如果为空，则为输入+pcoa.pdf -if(opts$output == ""){ - opts$output = paste0(opts$input, ".volcano.pdf")} - - -#----1.3. 加载包 Load packages#---- - -suppressWarnings(suppressMessages(library(amplicon))) -# 5.23修改 加载包 -package_list <- c('cowplot') -# 判断R包加载是否成功来决定是否安装后再加载 -for(pac in package_list){ - if(!suppressWarnings(suppressMessages(require(pac, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){ - install.packages(pac, repos=site) - suppressWarnings(suppressMessages(library(pac, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))) - } -} - -#----2. 读取文件 Read files#---- - -#----2.1 读取差异比较结果#---- -diff = read.table(opts$input, header=T, row.names=1, sep="\t", comment.char="") -diffN = as.data.frame(table(diff$level)) -# 5.23修改 对区间赋值 -diff$log2FC[diff$log2FC> 4]=4 -diff$log2FC[diff$log2FC< -4]=-4 -breaks= c(-4,-3,-2,-1,0,1,2,3,4) - - -#----3. 统计绘图#---- - -main_theme = theme(panel.background = element_blank(), - panel.grid = element_blank(), - axis.line.x = element_line(size = .5, colour = "black"), - axis.line.y = element_line(size = .5, colour = "black"), - axis.ticks = element_line(size = .5, color = "black"), - axis.text = element_text(size = 7, color = "black"), - legend.position = "right", - legend.background = element_blank(), - legend.key = element_blank(), - legend.text = element_text(size = 7), - text = element_text(family = "sans", size = 7)) - -# p =ggplot(diff[diff$log2FC < 4 & diff$log2FC > -4, ], aes(x=log2FC, y=log2CPM, color=level)) + -p =ggplot(diff, aes(x=log2FC, y=log2CPM, color=level)) + - geom_point()+scale_colour_manual(values=c("green","red","grey"))+ - theme_classic()+ - #加抖动 - # geom_jitter(data = subset(diff, log2FC == 4), width = 0.2, height = 0.2, color = "red")+ - # geom_jitter(data = subset(diff, log2FC == -4), width = 0.2, height = 0.2, color = "green")+ - # 5.23修改 重设刻度间隔 - scale_x_continuous(breaks=seq(-4,4,1))+ - theme_classic()+ - labs(x="log2(fold change)", y="log2(count per million)", - title=gsub(".txt","",basename(opts$input)))+ - annotate("text",x=-3,y=15,label=table(diff$level)[1])+ - annotate("text",x=3,y=15,label=table(diff$level)[2]) + - main_theme -# 5.23修改 加图例 -# legend= ggplot(diff, aes(x=log2FC, y=log2CPM, color=level)) + -# geom_point() + theme_classic()+ -# scale_colour_manual(values=c("green","red","grey")) + -# main_theme -# -# p1=plot_grid(p, -# plot_grid(get_legend(legend), -# ncol = 1, -# align = "hv"), -# ncol = 2,align = "hv",rel_widths=c(9,1)) -p -ggsave(opts$output, p, width = opts$width, height = opts$height, units = "mm") - +#!/usr/bin/env Rscript + +# Copyright 2024-2026 Defeng Bai + +# If used this script, please cited: +# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 + +# 手动运行脚本请，需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 + +# 更新 +# 2024/11/12：增加利用MaAsLin2差异分析结果绘制火山图代码 +# 2025/11/27：规范代码 +# 2026/5/22：更新依赖包安装 + +# 1.1 程序功能描述和主要步骤 + +# 程序功能：火山图展示组间差异 +# Functions: Species difference illustrated using volcano plot + +options(warn = -1) # Turn off warning + + +## 设置输入输出文件和参数 + +# 修改下面`default=`后面的文件和参数。 +# +# 输入文件为MaAsLin2组间差异分析得到的两个结果文件(metaphlan4/MaAsLin2_overall_difference.csv)+趋势信息(metaphlan4/MaAsLin2_enriched_depleted.csv) +# +# 输入文件"-i", "--input"，metaphlan4/MaAsLin2_overall_difference.csv; MaAsLin2差异分析结果汇总信息； +# +# 实验设计"-d", "--trend"，默认`metaphlan4/MaAsLin2_enriched_depleted.csv`，差异分析结果趋势信息； +# +# 图片宽"-w", "--width"，默认89 mm，根据图像布局可适当增大或缩小 +# +# 图片高"-e", "--height"，默认59 mm，根据图像布局可适当增大或缩小 +# +# 分组列名"-o", "--output"，默认为输出目录； + + +# 1.2 依赖包安装 + +# CRAN 镜像 +site <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN" + +# 安装函数 +install_if_missing <- function(pkg, type = "CRAN", repo = NULL) { + if (!requireNamespace(pkg, quietly = TRUE)) { + + if (type == "CRAN") { + install.packages(pkg, repos = site) + + } else if (type == "BIOC") { + if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager", repos = site) + } + BiocManager::install(pkg, update = FALSE, ask = FALSE) + + } else if (type == "GITHUB") { + if (!requireNamespace("devtools", quietly = TRUE)) { + install.packages("devtools", repos = site) + } + devtools::install_github(repo, force = TRUE) + } + } +} + +# ===== CRAN 包 ===== +cran_pkgs <- c("ggplot2", "BiocManager", "optparse", "dplyr", "magrittr") +for (pkg in cran_pkgs) install_if_missing(pkg, "CRAN") + +# ===== Bioconductor 包 ===== +bioc_pkgs <- c("multcompView") +for (pkg in bioc_pkgs) install_if_missing(pkg, "BIOC") + +# ===== GitHub 包 ===== +install_if_missing("amplicon", "GITHUB", repo = "microbiota/amplicon") +install_if_missing("ggvolcano", "GITHUB", repo = "BioSenior/ggvolcano") + + +# 1.3 解析命令行 +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +# site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-i", "--input"), type="character", default="metaphlan4/MaAsLin2_overall_difference.csv", + help="Unfiltered OTU table [default %default]"), + make_option(c("-d", "--trend"), type="character", default="metaphlan4/MaAsLin2_enriched_depleted.csv", + help="metadata file or metadata [default %default]"), + make_option(c("-o", "--output"), type="character", default="metaphlan4/", + help="Output quantile value for filter feature table [default %default]"), + make_option(c("-w", "--width"), type="numeric", default=89, + help="Width of figure [default %default]"), + make_option(c("-e", "--height"), type="numeric", default=59, + help="Height of figure [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + + +# 2. 依赖关系检查、安装和加载 + +# 依赖包列表 +package_list <- c( + "ggplot2", "BiocManager", "optparse", "dplyr", "magrittr", "ggvolcano" +) + +# 批量安装和加载 +for (p in package_list) { + # 如果未安装，则安装 + if (!requireNamespace(p, quietly = TRUE)) { + install.packages(p, repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/") + } + + # 批量加载，抑制警告和消息 + suppressWarnings( + suppressMessages( + library(p, character.only = TRUE) + ) + ) +} + + +# 读取输入文件 + +# 读取趋势表格 +data_species01 <- read.csv(opts$trend, row.names = 1) +data_species01$species <- rownames(data_species01) +data_species01 <- data_species01[, -c(1:4)] + +# 读取差异结果汇总表 +data_MWAS <- read.csv(opts$input, row.names = 1) +data_MWAS$species <- rownames(data_MWAS) + +data_species02 <- merge(data_species01, data_MWAS, by = "species") + +data_vol <- data_species02 +data_vol = as.data.frame(data_vol) + +data_vol2 <- data_vol +data_vol2$padj2 <- -log10(data_vol2$FDR) + +# logFC = 0.5 +# P.Value = 0.05 +# library(ggplot2) +p_volcano1 <- ggplot(data = data_vol2, aes(x = -Beta, y = padj2)) + + geom_point(alpha = 0.4, size = 2.0, aes(color = NPC.association)) + + ylab("-log10(Pvalue)") + + scale_color_manual(values = c("#74add1","#a60026", "grey")) + + geom_vline(xintercept = 0, lty = 4, col = "black", lwd = 0.4) + + geom_hline(yintercept = -log10(0.75), lty = 4, col = "black", lwd = 0.4) + + labs(x = "Coef. (by MaAsLin2)", y= bquote(atop(-Log[10]~italic(FDR))))+ + theme_bw() + +# add labels +# library(dplyr) +# select top 5 enriched species +up_data1 <- filter(data_vol2, data_vol2$NPC.association == "enriched") +up_data2 <- arrange(up_data1, desc(up_data1$padj2)) +up_data_5 <- up_data2[1:5, ] + +# select top 25 depleted species +down_data1 <- filter(data_vol2, data_vol2$NPC.association == "depleted") +down_data2 <- arrange(down_data1, desc(down_data1$padj2)) +down_data_5 <- down_data2[1:5, ] + +# using geom_text_repel() to add labels +library(ggrepel) +p_volcano2 <- p_volcano1 + + geom_text_repel(data = up_data_5, aes(x = -Beta, + y = padj2, + label = up_data_5$Feature), size = 2, max.overlaps = 20) + + geom_text_repel(data = down_data_5, aes(x = -Beta, + y = padj2, + label = down_data_5$Feature), size = 2, max.overlaps = 20)+ + theme( + legend.position = c(0.20, 0.75), + legend.text = element_text(size = 6), + legend.title = element_text(size = 7), + legend.background = element_rect(fill = "transparent", colour = NA), # 图例背景透明 + legend.key = element_rect(fill = "transparent", colour = NA), # 图例色块背景透明 + panel.grid = element_blank() + ) +ggsave(paste0(opts$output,"Difference_valcano.pdf"), p_volcano2, width=opts$width, height=opts$height, unit='mm') From 157b23bcbf61cda4eccc627b5457612ffbe67d8f Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:32:42 +0800 Subject: [PATCH 06/19] Update scripts for 2StatPlot.sh --- script/FastSpar_input01.R | 313 ++++++++++++++++++++++++++++++++++++ script/GTDBTK_tax.sh | 61 +++++++ script/remove_otu_header.sh | 38 +++++ script/run_FastSpar.sh | 140 ++++++++++++++++ script/run_sparcc3.sh | 119 ++++++++++++++ script/txt2tsv.sh | 55 +++++++ 6 files changed, 726 insertions(+) create mode 100644 script/FastSpar_input01.R create mode 100644 script/GTDBTK_tax.sh create mode 100644 script/remove_otu_header.sh create mode 100644 script/run_FastSpar.sh create mode 100644 script/run_sparcc3.sh create mode 100644 script/txt2tsv.sh diff --git a/script/FastSpar_input01.R b/script/FastSpar_input01.R new file mode 100644 index 0000000..a2c9e20 --- /dev/null +++ b/script/FastSpar_input01.R @@ -0,0 +1,313 @@ +#!/usr/bin/env Rscript + +############################# +## Install & Load packages +############################# + +# 检查并安装软件包 +packages <- c( + "optparse", + "Hmisc", + "igraph" +) + +for(pkg in packages){ + + if(!require(pkg, + character.only = TRUE, + quietly = TRUE)){ + + install.packages( + pkg, + repos="https://cloud.r-project.org" + ) + + library(pkg, + character.only = TRUE) + + } + +} + +# 定义输入参数 +option_list <- list( + + make_option("--abundance", + type="character", + help="MAG abundance table"), + + make_option("--taxonomy", + type="character", + help="taxonomy table"), + + make_option("--metadata", + type="character", + help="metadata table"), + + make_option("--group1", + default="Centenarians", + type="character"), + + make_option("--group2", + default="Young", + type="character"), + + make_option("--abundance_cutoff", + default=0.0007, + type="double"), + + make_option("--prevalence_cutoff", + default=0, + type="double"), + + make_option("--outdir", + default="result", + type="character") + +) + +opt <- parse_args( + OptionParser(option_list=option_list) +) + +dir.create( + opt$outdir, + showWarnings = FALSE, + recursive = TRUE +) + +# 读取数据 + +MAG_table_raw <- + read.table( + opt$abundance, + header=TRUE, + row.names=1, + sep="\t", + check.names=FALSE + ) + +tax_table_raw <- + read.table( + opt$taxonomy, + header=TRUE, + row.names=1, + sep="\t", + check.names=FALSE + ) + +meta_table_raw <- + read.table( + opt$metadata, + header=TRUE, + row.names=1, + sep="\t", + check.names=FALSE + ) + + +# 按组提取样本 + +MAG_table_group1 <- + MAG_table_raw[, + rownames( + meta_table_raw[ + meta_table_raw$Group==opt$group1,]) + ] + +MAG_table_group2 <- + MAG_table_raw[, + rownames( + meta_table_raw[ + meta_table_raw$Group==opt$group2,]) + ] + +MAG_tab_group1 <- t(MAG_table_group1) + +MAG_tab_group2 <- t(MAG_table_group2) + + +# 计算Prevalence + +prev_group1 <- colMeans(MAG_tab_group1>0) + +prev_group2 <- colMeans(MAG_tab_group2>0) + +write.table( + data.frame(prevalence=prev_group1), + file=file.path(opt$outdir, + "Prevalence_group1.tsv"), + sep="\t", + quote=FALSE +) + +write.table( + data.frame(prevalence=prev_group2), + file=file.path(opt$outdir, + "Prevalence_group2.tsv"), + sep="\t", + quote=FALSE +) + + +# 绘制 Prevalence 曲线数据 + +calculate_prevalence <- function(data,group){ + + threshold <- seq(0.01,1,0.01) + + prevalence <- + sapply( + threshold, + function(x) + sum(data>x) + ) + + data.frame( + threshold, + prevalence, + group + ) + +} + +curve1 <- + calculate_prevalence( + prev_group1, + opt$group1) + +curve2 <- + calculate_prevalence( + prev_group2, + opt$group2) + +curve <- + rbind( + curve1, + curve2) + +write.table( + curve, + file=file.path( + opt$outdir, + "Prevalence_curve.tsv"), + sep="\t", + quote=FALSE, + row.names=FALSE +) + +# 计算 Relative abundance +total_relative_abundance <- + colSums(t(MAG_table_raw))/ + sum(colSums(t(MAG_table_raw))) + +write.table( + data.frame( + Relative_abundance= + total_relative_abundance), + file=file.path( + opt$outdir, + "Relative_abundance.tsv"), + sep="\t", + quote=FALSE +) + +# Relative abundance 曲线 +calculate_rela <- function(data){ + + threshold <- seq( + 0, + max(data), + by=0.00025 + ) + + species <- + sapply( + threshold, + function(x) + sum(data>x) + ) + + data.frame( + threshold, + species + ) + +} + +rela <- + calculate_rela( + total_relative_abundance) + +write.table( + rela, + file=file.path( + opt$outdir, + "Relative_abundance_curve.tsv"), + sep="\t", + quote=FALSE, + row.names=FALSE +) + +# MAG筛选 +is.above.prev <- + prev_group1>=opt$prevalence_cutoff & + prev_group2>=opt$prevalence_cutoff + +is.above.abundance <- + total_relative_abundance>= + opt$abundance_cutoff + +target_MAG <- + rownames(MAG_table_raw)[ + is.above.prev & + is.above.abundance + ] + +cat( + "Selected MAG:", + length(target_MAG), + "\n" +) + +# 输出筛选后的MAG +target_MAG_table <- + MAG_table_raw[target_MAG,] + +write.table( + target_MAG_table, + file=file.path( + opt$outdir, + "Selected_MAG.tsv"), + sep="\t", + quote=FALSE +) + +# 输出FastSpar输入文件 + +write.table( + MAG_table_group1[target_MAG,], + file=file.path( + opt$outdir, + paste0( + "MAG_for_FastSpar_", + opt$group1, + ".txt")), + sep="\t", + quote=FALSE +) + +write.table( + MAG_table_group2[target_MAG,], + file=file.path( + opt$outdir, + paste0( + "MAG_for_FastSpar_", + opt$group2, + ".txt")), + sep="\t", + quote=FALSE +) + + + diff --git a/script/GTDBTK_tax.sh b/script/GTDBTK_tax.sh new file mode 100644 index 0000000..2a9148f --- /dev/null +++ b/script/GTDBTK_tax.sh @@ -0,0 +1,61 @@ +# #!/bin/bash +# +# awk -F'\t' ' +# BEGIN{ +# OFS="\t" +# print "user_genome","classification","ID" +# } +# NR>1{ +# species=$8 +# sub(/^s__/,"",species) +# +# print $1,species,species" "$1 +# } +# ' itol/tax01.txt > itol/GTDBTK_npc.txt + + +#!/bin/bash +# ============================================ +# Convert GTDB-Tk taxonomy to iTOL annotation +# Usage: +# bash gtdbtk2itol.sh -i input.txt -o output.txt +# ============================================ + +usage() { + echo "Usage: $0 -i -o " + exit 1 +} + +while getopts "i:o:h" opt; do + case $opt in + i) INPUT="$OPTARG" ;; + o) OUTPUT="$OPTARG" ;; + h) usage ;; + *) usage ;; + esac +done + +[ -z "$INPUT" ] && usage +[ -z "$OUTPUT" ] && usage + +if [ ! -f "$INPUT" ]; then + echo "Error: Input file '$INPUT' not found." + exit 1 +fi + +awk -F'\t' ' +BEGIN{ + OFS="\t" + print "user_genome","classification","ID" +} +NR>1{ + species=$8 + sub(/^s__/,"",species) + + print $1,species,species" "$1 +} +' "$INPUT" > "$OUTPUT" + +echo "Done!" +echo "Input : $INPUT" +echo "Output: $OUTPUT" \ No newline at end of file diff --git a/script/remove_otu_header.sh b/script/remove_otu_header.sh new file mode 100644 index 0000000..d141b27 --- /dev/null +++ b/script/remove_otu_header.sh @@ -0,0 +1,38 @@ +#!/bin/bash +############################################################################### +# Script Name : remove_otu_header.sh +# Description : Remove "#OTU ID" from the first line of a FastSpar output file. +# +# Usage: +# bash remove_otu_header.sh input.tsv output.txt +# +# Example: +# bash remove_otu_header.sh \ +# Centenarians_correlation.tsv \ +# R_Centenarians2.txt +############################################################################### + +# 检查参数 +if [ $# -ne 2 ]; then + echo "Usage: bash $0 " + exit 1 +fi + +INPUT=$1 +OUTPUT=$2 + +# 判断输入文件是否存在 +if [ ! -f "$INPUT" ]; then + echo "Error: Input file '$INPUT' not found!" + exit 1 +fi + +# 创建输出目录 +mkdir -p "$(dirname "$OUTPUT")" + +# 删除第一行的 #OTU ID +sed '1s/^#OTU ID//' "$INPUT" > "$OUTPUT" + +echo "Done!" +echo "Input : $INPUT" +echo "Output: $OUTPUT" \ No newline at end of file diff --git a/script/run_FastSpar.sh b/script/run_FastSpar.sh new file mode 100644 index 0000000..896f9e7 --- /dev/null +++ b/script/run_FastSpar.sh @@ -0,0 +1,140 @@ +#!/bin/bash +set -euo pipefail +############################################## +## Default parameters +############################################## + +THREADS=20 +BOOTSTRAP=1000 +ITERATIONS=1000 + +############################################## +## Parse parameters +############################################## + +usage(){ +echo " +Usage: +bash run_FastSpar.sh \ + -c MAG_for_FastSpar_Centenarians_fixed.tsv \ + -y MAG_for_FastSpar_Young_fixed.tsv \ + -o FastSpar_result \ + -t 20 \ + -b 1000 \ + -i 1000 +" +exit 1 +} + +while getopts "c:y:o:t:b:i:h" opt +do + case $opt in + c) CENT=$OPTARG;; + y) YOUNG=$OPTARG;; + o) OUTDIR=$OPTARG;; + t) THREADS=$OPTARG;; + b) BOOTSTRAP=$OPTARG;; + i) ITERATIONS=$OPTARG;; + h) usage;; + *) usage;; + esac +done + +############################################## +## Check parameters +############################################## + +[ -z "${CENT:-}" ] && usage +[ -z "${YOUNG:-}" ] && usage +[ -z "${OUTDIR:-}" ] && usage + +############################################## +## Check software +############################################## + +for exe in fastspar fastspar_bootstrap fastspar_pvalues parallel +do + command -v ${exe} >/dev/null 2>&1 || { + echo "ERROR: ${exe} not found!" + exit 1 + } +done + +mkdir -p ${OUTDIR} + +############################################## +## Function +############################################## + +run_fastspar(){ + +GROUP=$1 +INPUT=$2 + +echo "========================================" +echo "Processing ${GROUP}" +echo "========================================" + +mkdir -p ${OUTDIR}/${GROUP} +mkdir -p ${OUTDIR}/${GROUP}/bootstrap +mkdir -p ${OUTDIR}/${GROUP}/bootstrap_cor +mkdir -p ${OUTDIR}/${GROUP}/bootstrap_cov + +############################################## +## Step1 correlation +############################################## + +fastspar \ + --otu_table ${INPUT} \ + --correlation ${OUTDIR}/${GROUP}/${GROUP}_correlation.tsv \ + --covariance ${OUTDIR}/${GROUP}/${GROUP}_covariance.tsv \ + --iterations ${ITERATIONS} \ + --threads ${THREADS} + +############################################## +## Step2 bootstrap +############################################## + +fastspar_bootstrap \ + --otu_table ${INPUT} \ + --number ${BOOTSTRAP} \ + --prefix ${OUTDIR}/${GROUP}/bootstrap/${GROUP}_otu_ + +############################################## +## Step3 bootstrap correlation +############################################## + +parallel -j ${THREADS} \ +fastspar \ + --otu_table {} \ + --correlation ${OUTDIR}/${GROUP}/bootstrap_cor/{/.}_cor.tsv \ + --covariance ${OUTDIR}/${GROUP}/bootstrap_cov/{/.}_cov.tsv \ + --iterations 5 \ +::: ${OUTDIR}/${GROUP}/bootstrap/*.tsv + +############################################## +## Step4 pvalue +############################################## + +fastspar_pvalues \ + --otu_table ${INPUT} \ + --correlation ${OUTDIR}/${GROUP}/${GROUP}_correlation.tsv \ + --prefix ${OUTDIR}/${GROUP}/bootstrap_cor/${GROUP}_otu_ \ + --permutations ${BOOTSTRAP} \ + --outfile ${OUTDIR}/${GROUP}/${GROUP}_pvalues.tsv + +echo "${GROUP} finished." + +} + +############################################## +## Run +############################################## + +run_fastspar Centenarians ${CENT} + +run_fastspar Young ${YOUNG} + +echo "=====================================" +echo "All FastSpar analyses completed!" +echo "=====================================" diff --git a/script/run_sparcc3.sh b/script/run_sparcc3.sh new file mode 100644 index 0000000..c9269ab --- /dev/null +++ b/script/run_sparcc3.sh @@ -0,0 +1,119 @@ +#!/bin/bash +# Copyright 2026 Defeng Bai + +set -euo pipefail + +############################################################### +# Usage +# +# bash run_sparcc3.sh input.txt R.txt P.txt +# +############################################################### + +if [ $# -ne 3 ]; then + echo "Usage:" + echo "bash $0 input_abundance.txt correlation.txt pvalue.txt" + exit 1 +fi + +############################# +# Input +############################# + +INPUT=$(realpath "$1") +R_OUT=$(realpath "$2") +P_OUT=$(realpath "$3") + +############################# +# Parameters +############################# + +ITER=20 +BOOTSTRAP=100 + +############################# +# Working directory +############################# + +WD=$(pwd) + +############################# +# Download SparCC3 +############################# + +if [ ! -d "${WD}/SparCC3" ]; then + echo "Downloading SparCC3..." + git clone https://github.com/JCSzamosi/SparCC3.git +else + echo "SparCC3 already exists." +fi + +cd "${WD}/SparCC3" + +mkdir -p data +mkdir -p example/basis_corr +mkdir -p example/pvals + +############################# +# Copy input +############################# + +cp "${INPUT}" data/input.txt + +############################# +# Step1 Correlation +############################# + +python SparCC.py \ + data/input.txt \ + -i ${ITER} \ + --cor_file=example/basis_corr/correlation.tsv \ + > example/basis_corr/sparcc.log + +############################# +# Step2 Bootstrap +############################# + +python MakeBootstraps.py \ + data/input.txt \ + -n ${BOOTSTRAP} \ + -t bootstrap_#.txt \ + -p example/pvals/ \ + >> example/basis_corr/sparcc.log + +############################# +# Step3 Bootstrap correlations +############################# + +for ((n=0;n> example/basis_corr/sparcc.log +done + +############################# +# Step4 P values +############################# + +python PseudoPvals.py \ + example/basis_corr/correlation.tsv \ + example/pvals/bootstrap_cor_#.txt \ + ${BOOTSTRAP} \ + -o example/pvals/pvalue.tsv \ + -t two_sided \ + >> example/basis_corr/sparcc.log + +############################# +# Output +############################# + +cp example/basis_corr/correlation.tsv "${R_OUT}" +cp example/pvals/pvalue.tsv "${P_OUT}" + +echo +echo "Finished!" +echo "Correlation : ${R_OUT}" +echo "P-value : ${P_OUT}" diff --git a/script/txt2tsv.sh b/script/txt2tsv.sh new file mode 100644 index 0000000..d9b94b7 --- /dev/null +++ b/script/txt2tsv.sh @@ -0,0 +1,55 @@ +#!/bin/bash + +############################################################################### +# Script Name : txt2tsv.sh +# Description : Convert whitespace-delimited txt to tab-delimited tsv +# and rename the first header to "#OTU ID". +# +# Usage: +# bash txt2tsv.sh input.txt output.tsv +# +# Example: +# bash txt2tsv.sh \ +# ./coverm/MAG_for_FastSpar_Centenarians.txt \ +# ./coverm/MAG_for_FastSpar_Centenarians_fixed.tsv +############################################################################### + +# 检查参数 +if [ $# -ne 2 ]; then + echo "Usage: bash $0 " + exit 1 +fi + +INPUT=$1 +OUTPUT=$2 + +# 判断输入文件是否存在 +if [ ! -f "$INPUT" ]; then + echo "Error: Input file '$INPUT' does not exist!" + exit 1 +fi + +# 创建输出目录 +mkdir -p "$(dirname "$OUTPUT")" + +# 转换 +awk ' +BEGIN{ + OFS="\t" +} +{ + # 将连续空格压缩为TAB + $1=$1 + + # 第一行为表头 + if(NR==1){ + $1="#OTU ID" + } + + print +} +' "$INPUT" > "$OUTPUT" + +echo "Done!" +echo "Input : $INPUT" +echo "Output: $OUTPUT" \ No newline at end of file From a51c2f7239c81a29f7a7695b3b66ec8efd55c5fe Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:03:51 +0800 Subject: [PATCH 07/19] Update scripts for 2StatPlot.sh --- script/SparCC_data_processing.R | 420 ++++++++++++++++++++++---------- script/TCG_FromNetwork.R | 340 ++++++++++++++++++++++++++ script/extract_eggnog_anno01.sh | 46 ++++ script/extract_pathways.sh | 66 +++++ script/get_graphml.R | 163 +++++++++++++ 5 files changed, 913 insertions(+), 122 deletions(-) create mode 100644 script/TCG_FromNetwork.R create mode 100644 script/extract_eggnog_anno01.sh create mode 100644 script/extract_pathways.sh create mode 100644 script/get_graphml.R diff --git a/script/SparCC_data_processing.R b/script/SparCC_data_processing.R index abaf5c0..cfb6505 100644 --- a/script/SparCC_data_processing.R +++ b/script/SparCC_data_processing.R @@ -1,122 +1,298 @@ -#!/usr/bin/env Rscript - -# Copyright 2024 De-feng Bai - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 - -# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 -# Set Work Directory - Choose Directory / To Source File Location - -# Clean enviroment object -rm(list=ls()) - -# 1.1 简介 Introduction #---- - -# 程序功能：处理Metaphlan4得到的物种组成数据作为SparCC计算的输入 -# Functions: Processing of species composition data obtained from Metaphlan4 as input for SparCC calculations - - -options(warn = -1) # Turn off warning - -# 1.2 参数 Parameters #---- -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="result12/metaphlan4/Species.txt", - help="Metaphlan4 species table"), - make_option(c("-g", "--group"), type="character", default="result12/metadata.txt", - help="Unfiltered OTU table [default %default]"), - make_option(c("-o", "--output"), type="character", default="result12/metaphlan4/", - help="Output file for SparCC analysis in different groups [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# Install related packages -if (FALSE){ - source("https://bioconductor.org/biocLite.R") - biocLite(c("reshape2","ggplot2","ggprism","dplyr","plyr", - "igraph")) # ,"vegan" -} -# load related packages -suppressWarnings(suppressMessages(library("reshape2"))) -suppressWarnings(suppressMessages(library("ggplot2"))) -suppressWarnings(suppressMessages(library("ggprism"))) -# suppressWarnings(suppressMessages(library("vegan"))) -suppressWarnings(suppressMessages(library("dplyr"))) -suppressWarnings(suppressMessages(library("plyr"))) -suppressWarnings(suppressMessages(library("igraph"))) - - -# metadata -#design <- read.table(file = "result12/metaphlan4/metadata.txt", sep = "\t", header = T, row.names=1) -design <- read.table(opts$group, sep = "\t", header = T, row.names=1) -#df3 <- read.table(file = "result12/metaphlan4/Species.txt", sep = "\t", header = T, check.names = FALSE) -df3 <- read.table(opts$input, sep = "\t", header = T, check.names = FALSE) -df_species <- df3 - -# sum of Species -data<-aggregate(.~ Taxonomy,data=df_species,sum) -rownames(data) = data$Taxonomy -data = data[, -1] - -#If species are too much, please select part of them according to relative abundace -#species_selected <- read.table(file = "data/Species_selected.txt", sep = "\t", header = T, check.names = FALSE) -#rownames(species_selected) <- species_selected$ID -#data <- data[rownames(data)%in%rownames(species_selected), ] - -data3 = apply(data, 2, function(x)x/100) -#sum(data3$H001F) - -# Check data -dim(data3) -data3 = data3 * 100000 -OTU.table.filtered.colnames <- colnames(data3) -OTU.table.filtered.sparcc <- cbind(rownames(data3), data3) -colnames(OTU.table.filtered.sparcc) <- c("OTU_id", OTU.table.filtered.colnames) -OTU.table.filtered.sparcc2 <- t(OTU.table.filtered.sparcc) -OTU.table.filtered.sparcc2 <- OTU.table.filtered.sparcc2[-1,] -OTU.table.filtered.sparcc2 <- as.data.frame(OTU.table.filtered.sparcc2) -#OTU.table.filtered.sparcc2$group <- rownames(OTU.table.filtered.sparcc2) -#OTU.table.filtered.sparcc2$group = as.character(OTU.table.filtered.sparcc2$group) -#OTU.table.filtered.sparcc2$group = gsub("[0-9]","", OTU.table.filtered.sparcc2$group) -otutab <- as.data.frame(t(OTU.table.filtered.sparcc2)) - -# Select by manual set group -# NPC group -if (TRUE){ - sub_design = subset(design, Group %in% c("Cancer")) - sub_design$Group = factor(sub_design$Group, levels=c("Cancer")) -} -idx = rownames(sub_design) %in% colnames(otutab) -sub_design_Cancer = sub_design[idx,] -sub_otutab_Cancer = otutab[,rownames(sub_design_Cancer)] -sub_otutab_Cancer <- as.data.frame(sub_otutab_Cancer) -#write.table(sub_otutab_npc, file = "results/Species_sparcc_npc_p01_11R_adjusted_1031.txt", row.names = T, sep = "\t", quote = T, col.names = T) -write.table(sub_otutab_Cancer, file = paste(opts$output, "Cancer_sparcc.txt", sep=""), row.names = T, sep = "\t", quote = T, col.names = T) - -# Healthy group -if (TRUE){ - sub_design = subset(design, Group %in% c("Normal")) - sub_design$Group = factor(sub_design$Group, levels=c("Normal")) -} -idx = rownames(sub_design) %in% colnames(otutab) -sub_design_Normal = sub_design[idx,] -sub_otutab_Normal = otutab[,rownames(sub_design_Normal)] -sub_otutab_Normal = as.data.frame(sub_otutab_Normal) -sub_otutab_Normal = as.data.frame(sub_otutab_Normal) -#write.table(sub_otutab_healthy, file = "results/Species_sparcc_npc_h01_11R_adjusted_1031.txt", row.names = T, sep = "\t", quote = T, col.names = T) -write.table(sub_otutab_Normal, file = paste(opts$output, "Normal_sparcc.txt", sep=""), row.names = T, sep = "\t", quote = T, col.names = T) +#!/usr/bin/env Rscript + +# Copyright 2024 De-feng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# 1.1 简介 Introduction #---- + +# 程序功能：处理得到的物种组成数据作为SparCC计算的输入 +# Functions: Processing Metaphlan4 species abundance table +# as SparCC input files for each group + +# Clean environment +rm(list = ls()) + +options(warn = -1) + +# ========================================================= +# Parameters +# ========================================================= + +site = "https://mirrors.tuna.tsinghua.edu.cn/CRAN" + +# Install optparse if needed +if (!suppressWarnings( + suppressMessages( + require("optparse", + character.only = TRUE, + quietly = TRUE, + warn.conflicts = FALSE) + ) +)) { + + install.packages("optparse", repos = site) + + require(optparse) +} + +# Parse arguments +option_list = list( + + make_option( + c("-i", "--input"), + type = "character", + default = "result12/metaphlan4/Species.txt", + help = "Metaphlan4 species table" + ), + + make_option( + c("-g", "--group"), + type = "character", + default = "result12/metadata.txt", + help = "Metadata file" + ), + + make_option( + c("-o", "--output"), + type = "character", + default = "result12/metaphlan4/", + help = "Output directory" + ), + + make_option( + c("--groups"), + type = "character", + default = NULL, + help = "Groups to extract, separated by comma. Example: Cancer,Normal" + ) +) + +opts = parse_args(OptionParser(option_list = option_list)) + +cat("=====================================\n") +cat("Parameters:\n") +print(opts) +cat("=====================================\n") + +# ========================================================= +# Load packages +# ========================================================= + +packages <- c( + "reshape2", + "ggplot2", + "ggprism", + "dplyr", + "plyr", + "igraph" +) + +for (p in packages) { + + suppressWarnings( + suppressMessages( + library(p, + character.only = TRUE) + ) + ) +} + +# ========================================================= +# Create output directory +# ========================================================= + +if (!dir.exists(opts$output)) { + + dir.create(opts$output, + recursive = TRUE) +} + +# ========================================================= +# Read metadata +# ========================================================= + +design <- read.table( + opts$group, + sep = "\t", + header = TRUE, + row.names = 1, + check.names = FALSE +) + +# Check Group column +if (!"Group" %in% colnames(design)) { + + stop("ERROR: metadata file must contain a column named 'Group'") +} + +# ========================================================= +# Read species abundance table +# ========================================================= + +df3 <- read.table( + opts$input, + sep = "\t", + header = TRUE, + check.names = FALSE +) + +df_species <- df3 + +# ========================================================= +# Aggregate taxonomy +# ========================================================= + +# 获取第一列列名 +feature_col <- colnames(df_species)[1] + +cat("Feature column:", feature_col, "\n") + +data <- aggregate( + df_species[, -1], + by = list(df_species[[feature_col]]), + FUN = sum +) + +colnames(data)[1] <- feature_col + +rownames(data) <- data[[feature_col]] + +data <- data[, -1] + +# ========================================================= +# Convert abundance +# ========================================================= + +data3 <- apply(data, + 2, + function(x) x / 100) + +cat("Data dimension:\n") +print(dim(data3)) + +# SparCC recommended scaling +# data3 <- data3 * 100000 + +# ========================================================= +# Format OTU table +# ========================================================= + +# keep numeric matrix +data3 <- as.data.frame(data3) + +# species as rownames +rownames(data3) <- rownames(data) + +# samples as columns +colnames(data3) <- colnames(data) + +# transpose: +# rows = species +# cols = samples +otutab <- as.data.frame(data3) + +# ensure numeric +otutab[] <- lapply(otutab, as.numeric) + +# check +cat("OTU table dimension:\n") +print(dim(otutab)) + +cat("OTU table preview:\n") +print(head(otutab[,1:min(5,ncol(otutab))])) + + +# ========================================================= +# Determine groups +# ========================================================= + +if (is.null(opts$groups)) { + + group_list <- unique(design$Group) + +} else { + + group_list <- unlist( + strsplit(opts$groups, ",") + ) +} + +cat("Groups to process:\n") +print(group_list) + +# ========================================================= +# Function for extracting group +# ========================================================= + +extract_group <- function(group_name, + design, + otutab, + output_dir) { + + cat("\nProcessing group:", group_name, "\n") + + # subset metadata + sub_design <- subset( + design, + Group %in% c(group_name) + ) + + sub_design$Group <- factor( + sub_design$Group, + levels = c(group_name) + ) + + # sample intersection + idx <- rownames(sub_design) %in% colnames(otutab) + + sub_design <- sub_design[idx, ] + + # subset otu table + sub_otutab <- otutab[, rownames(sub_design)] + + sub_otutab <- as.data.frame(sub_otutab) + + # output file + outfile <- paste0( + output_dir, + "/", + group_name, + "_sparcc.txt" + ) + + write.table( + sub_otutab, + file = outfile, + row.names = TRUE, + sep = "\t", + quote = FALSE, + col.names = TRUE + ) + + cat("Output:", outfile, "\n") + cat("Samples:", ncol(sub_otutab), "\n") + cat("Species:", nrow(sub_otutab), "\n") +} + +# ========================================================= +# Run all groups +# ========================================================= + +for (g in group_list) { + + extract_group( + group_name = g, + design = design, + otutab = otutab, + output_dir = opts$output + ) +} + +cat("\nAll done!\n") + + diff --git a/script/TCG_FromNetwork.R b/script/TCG_FromNetwork.R new file mode 100644 index 0000000..ae325bb --- /dev/null +++ b/script/TCG_FromNetwork.R @@ -0,0 +1,340 @@ +#!/usr/bin/env Rscript + +# Copyright 2026 De-feng Bai + +if (!requireNamespace("optparse", quietly = TRUE)) { + install.packages("optparse", + repos = "https://cloud.r-project.org") +} + +library(optparse) + +option_list <- list( + + make_option( + c("-g", "--gtdb"), + type = "character", + help = "GTDBTK annotation file" + ), + + make_option( + "--group1_r", + type = "character", + help = "Group1 FastSpar correlation matrix" + ), + + make_option( + "--group1_p", + type = "character", + help = "Group1 FastSpar pvalue matrix" + ), + + make_option( + "--group2_r", + type = "character", + help = "Group2 FastSpar correlation matrix" + ), + + make_option( + "--group2_p", + type = "character", + help = "Group2 FastSpar pvalue matrix" + ), + + make_option(c("-o", "--output"), type="character", default="coverm/", + help="Output directory [default %default]") +) + +opts <- parse_args( + OptionParser( + option_list = option_list + ) +) + + +# 检查并安装软件包 + +options( + repos = c( + CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/" + ) +) + +options( + BioC_mirror = + "https://mirrors.tuna.tsinghua.edu.cn/bioconductor" +) + +# ===================================================== +# Install and load packages +# ===================================================== + +cran_packages <- c( + "reshape2", + "ggplot2", + "dplyr", + "purrr", + "igraph" +) + +bioc_packages <- c( + "impute", + "preprocessCore", + "WGCNA" +) + +# ---------- CRAN ---------- +for(pkg in cran_packages){ + + if(!requireNamespace(pkg, quietly = TRUE)){ + + cat("Installing CRAN package:", pkg, "\n") + + install.packages( + pkg, + repos = "https://cloud.r-project.org" + ) + } + + suppressPackageStartupMessages( + library(pkg, + character.only = TRUE) + ) +} + +# ---------- Bioconductor ---------- +if(length(bioc_packages) > 0){ + + if(!requireNamespace("BiocManager", + quietly = TRUE)){ + + install.packages( + "BiocManager", + repos = "https://cloud.r-project.org" + ) + } + + for(pkg in bioc_packages){ + + if(!requireNamespace(pkg, + quietly = TRUE)){ + + cat("Installing Bioconductor package:", pkg, "\n") + + BiocManager::install( + pkg, + ask = FALSE, + update = FALSE + ) + } + + suppressPackageStartupMessages( + library(pkg, + character.only = TRUE) + ) + } +} + +# library(reshape2) +# library(ggplot2) +# library(dplyr) +# library(purrr) +# library(ggplot2) +# library(igraph) +# library(WGCNA) + +# if (!requireNamespace("BiocManager", quietly = TRUE)) +# install.packages("BiocManager") + +#BiocManager::install("impute") + +#BiocManager::install("preprocessCore") + +#---- +#0. local function +FastsparToABC=function(Fastspar_R,Fastspar_P) +{ + tmpdata=Fastspar_R + tmpdata[lower.tri(tmpdata,diag = T)]=NA #set lower.tri and diag = NA + tmpdata=melt(tmpdata) # melt the matrix + tmpdata$node2=rep(colnames(Fastspar_R,length(colnames(Fastspar_R)))) # add node2 id + tmpdata=tmpdata[which(is.na(tmpdata$value)==FALSE),] + tmpR_abc=tmpdata + tmpdata=Fastspar_P + tmpdata[lower.tri(tmpdata,diag = T)]=NA #set lower.tri and diag = NA + tmpdata=melt(tmpdata) # melt the matrix + tmpdata$node2=rep(colnames(Fastspar_P,length(colnames(Fastspar_P)))) # add node2 id + tmpdata=tmpdata[which(is.na(tmpdata$value)==FALSE),] + tmpP_abc=tmpdata + myresult=data.frame(Source=tmpP_abc$variable,Target=tmpP_abc$node2,FastsparR=tmpR_abc$value,FastsparP=tmpP_abc$value,stringsAsFactors = F) + myresult$Source=as.character(myresult$Source) + return(myresult) +} + +FastsparABC_furprocess=function(myFastsparABC) +{ + myFastsparABC$Adjp=p.adjust(myFastsparABC$FastsparP,method = "BH") + myFastsparABC$NorP="P" + myFastsparABC$NorP[myFastsparABC$FastsparR<0]="N" + myFastsparABC$R2=myFastsparABC$FastsparR^2 + for(i in 1:nrow(myFastsparABC)) + { + if(myFastsparABC$Source[i]=0.05)]="U" + ind=which(mytestedge[[i]]$Adjp<0.05) + print(length(ind)) + tmpedge=mytestedge[[i]][ind,] + print(length(unique(c(tmpedge$Source,tmpedge$Target)))) + mysigedge=c(mysigedge,tmpedge$link_R) + mysigedge.list[[i]]=tmpedge + alledge.list[[i]]=mytestedge[[i]] +} +mysigedge_count=table(mysigedge) +table(mysigedge_count) + +c=alledge.list%>%reduce(full_join,by="link") +d=c[,c("link","NorP.x","NorP.y")] +d[is.na(d)]="U" +e=table(paste(d$NorP.x,d$NorP.y,d$NorP,sep = "")) +e=data.frame(e) +e$Freq_ratio=e$Freq/sum(e$Freq)*100 +e=e[order(e$Freq),] +e$Var1=factor(e$Var1,levels = e$Var1) +ggplot(data = e,aes(x=Var1,y=Freq))+ + geom_bar(stat = "identity",aes(fill=Var1))+ + coord_flip()+ + labs(x="",y="Number")+ + theme_bw()+ + theme(text=element_text(size=40))+theme(legend.position = "none") #--> figure S2 + +#---- +#4 get stable correlations +a=names(mysigedge_count)[which(mysigedge_count==2)] # =2 -> in all 2 groups +alink=gsub("=.*","",a) +mysigedge_stable=data.frame( + Source=gsub("/.*","",alink), + Target=gsub(".*/","",alink), + NorP=gsub(".*=","",a), + link=alink +) +mysigedge_stable$NorPv2=ifelse(mysigedge_stable$NorP=="P",1,-1) +mysignode_stable=data.frame(node=unique(c(mysigedge_stable$Source,mysigedge_stable$Target))) +mysignode_stable$GTDBTK=myallGTDBTK$classification[match(mysignode_stable$node,myallGTDBTK$ID)] + +#write.table(mysigedge_stable,"372stable.edge2.txt",sep = "\t",quote = F,row.names = F) +#write.table(mysignode_stable,"372stable.node2.txt",sep = "\t",quote = F,row.names = F) + +# write.table(mysigedge_stable, +# file = "stable.edge2.txt", +# output_dir = opts$output) +# +# write.table(mysignode_stable, +# file = "stable.node2.txt", +# output_dir = opts$output) + +write.table( + mysigedge_stable, + file = file.path(opts$output, "stable.edge2.txt"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) + +write.table( + mysignode_stable, + file = file.path(opts$output, "stable.node2.txt"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) + + diff --git a/script/extract_eggnog_anno01.sh b/script/extract_eggnog_anno01.sh new file mode 100644 index 0000000..67a993d --- /dev/null +++ b/script/extract_eggnog_anno01.sh @@ -0,0 +1,46 @@ +#!/bin/bash + +# 固定输入输出文件 +input="../temp/eggnog/output" +output="eggnog/data_venn1.txt" + +awk -F'\t' ' +BEGIN{ + OFS="\t" +} + +# 跳过注释行（eggNOG文件常见） +/^#/ {next} + +# 读取表头 +header==0{ + for(i=1;i<=NF;i++){ + if($i=="query") q=i + if($i=="COG_category") cog=i + if($i=="KEGG_ko") kegg=i + if($i=="GOs") go=i + if($i=="EC") ec=i + } + + if(!q || !cog || !kegg || !go || !ec){ + print "Error: required columns not found!" > "/dev/stderr" + exit 1 + } + + print "query","COGs","KEGG","GOs","ECs" + header=1 + next +} + +# 数据行处理 +{ + c = ($cog != "-" && $cog != "") ? 1 : 0 + k = ($kegg != "-" && $kegg != "") ? 1 : 0 + g = ($go != "-" && $go != "") ? 1 : 0 + e = ($ec != "-" && $ec != "") ? 1 : 0 + + print $q,c,k,g,e +} +' "$input" > "$output" + +echo "Done! Output written to $output" \ No newline at end of file diff --git a/script/extract_pathways.sh b/script/extract_pathways.sh new file mode 100644 index 0000000..c4dde67 --- /dev/null +++ b/script/extract_pathways.sh @@ -0,0 +1,66 @@ +#!/bin/bash + + + +# 检查参数 +if [ $# -eq 0 ]; then + echo "Usage:" + echo "bash $0 \"Pathway1\" \"Pathway2\" ..." + exit 1 +fi + +# 输入文件 +ANNOTATION="humann4/KEGG_KO_Annotation.tsv" +KO_FILE="humann4/ko_stratified.tsv" + +# 输出文件 +SELECTED="humann4/KEGG_Selected_Pathways.tsv" +FILTERED="humann4/KEGG_Selected_Pathways_filtered.tsv" +SANKEY="humann4/data_sankey2.txt" + +# 构建正则表达式 +PATTERN=$(printf "%s|" "$@") +PATTERN=${PATTERN%|} + +echo "Selected pathways:" +printf " %s\n" "$@" + +############################ +# Step 1. 提取目标通路 +############################ +awk -F'\t' -v pattern="$PATTERN" ' +NR==1 || $9 ~ ("^(" pattern ")$") +' "$ANNOTATION" > "$SELECTED" + +############################ +# Step 2. 保留样本中存在的KO +############################ +awk -F'\t' ' +NR==FNR{ + if(FNR>1){ + split($1,a,"|") + ko[a[1]]=1 + } + next +} +FNR==1 || ($1 in ko) +' "$KO_FILE" "$SELECTED" > "$FILTERED" + +############################ +# Step 3. 生成 Sankey 输入文件 +############################ +awk -F'\t' ' +BEGIN{ + OFS="\t" + print "Level2_Name","Level3_Name","KO" +} +NR>1{ + print $7,$9,$1 +} +' "$FILTERED" > "$SANKEY" + +echo "Done!" +echo "Output files:" +echo " $SELECTED" +echo " $FILTERED" +echo " $SANKEY" diff --git a/script/get_graphml.R b/script/get_graphml.R new file mode 100644 index 0000000..a183178 --- /dev/null +++ b/script/get_graphml.R @@ -0,0 +1,163 @@ +#!/usr/bin/env Rscript + +# Copyright 2026 Defeng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# Clean enviroment object +rm(list=ls()) + +# 1.1 简介 Introduction #---- + +# 程序功能：SparCC可视化 +# Functions: Sparcc Visualization + + +options(warn = -1) # Turn off warning + +# 1.2 参数 Parameters #---- +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-R", "--Correlation"), type="character", default="metaphlan4/R_Centenarians.txt", + help="Metaphlan4 species table"), + make_option(c("-P", "--Pvalue"), type="character", default="metaphlan4/P_Centenarians.txt", + help="Unfiltered OTU table [default %default]"), + make_option("--Group", type="character"), + make_option(c("-r", "--output"), type="character", default="metaphlan4/", + help="Output file for SparCC analysis in different groups [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + +options(repos = c(CRAN = "https://cloud.r-project.org")) + +packages <- c( + "reshape2", + "ggplot2", + "ggprism", + "dplyr", + "plyr", + "igraph", + "tidyverse", + "ggraph", + "magrittr" +) + +installed <- rownames(installed.packages()) +missing_pkgs <- packages[!(packages %in% installed)] + +if(length(missing_pkgs) > 0){ + install.packages(missing_pkgs, dependencies = TRUE) +} + +if (!requireNamespace("BiocManager", quietly = TRUE)) + install.packages("BiocManager") + +pkgs <- c( + "clusterProfiler", + "DOSE", + "enrichplot", + "ComplexHeatmap", + "pathview", + "KEGGREST" +) + +for(pkg in pkgs){ + if(!requireNamespace(pkg, quietly = TRUE)){ + BiocManager::install(pkg, update = FALSE, ask = FALSE) + } +} + +# if (!requireNamespace("clusterProfiler", quietly = TRUE)) +# BiocManager::install("clusterProfiler", update = FALSE, ask = FALSE) + +#library(clusterProfiler) + +# load related packages +suppressWarnings(suppressMessages(library("reshape2"))) +suppressWarnings(suppressMessages(library("ggplot2"))) +suppressWarnings(suppressMessages(library("ggprism"))) +# suppressWarnings(suppressMessages(library("vegan"))) +suppressWarnings(suppressMessages(library("dplyr"))) +suppressWarnings(suppressMessages(library("plyr"))) +suppressWarnings(suppressMessages(library("igraph"))) +suppressWarnings(suppressMessages(library("tidyverse"))) +suppressWarnings(suppressMessages(library("ggraph"))) +suppressWarnings(suppressMessages(library("magrittr"))) +suppressWarnings(suppressMessages(library("clusterProfiler"))) + + +#r.cor <- read.table("data3/r.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) +r.cor <- read.table(opts$Correlation, sep="\t", header=T, check.names=F,row.names = 1) +#p.cor <- read.table("data3/p.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) +p.cor <- read.table(opts$Pvalue, sep="\t", header=T, check.names=F,row.names = 1) + +r.cor[p.cor>0.05] <- 0 + +# Build network connection attributes and node attributes +# Convert data to long format for merging and add connection properties +r.cor$from = rownames(r.cor) +p.cor$from = rownames(p.cor) +p_value <- p.cor %>% + gather(key = "to", value = "p", -from) %>% + data.frame() +p_value$FDR <- p.adjust(p_value$p,"BH") +p_value <- p_value[, -3] + +cor.data<- r.cor %>% + gather(key = "to", value = "r", -from) %>% + data.frame() %>% + left_join(p_value, by=c("from","to")) +cor.data <- as.data.frame(cor.data) +cor.data <- cor.data[cor.data$FDR <= 0.05 & cor.data$from != cor.data$to, ] +cor.data <- cor.data[abs(cor.data$r) >= 0.6 & cor.data$from != cor.data$to, ] +cor.data <- cor.data %>% + plyr::mutate( + linecolor = ifelse(r > 0,"positive","negative"), + linesize = abs(r) + ) + +# Set node properties +vertices <- c(as.character(cor.data$from),as.character(cor.data$to)) %>% + as_tibble() %>% + group_by(value) %>% + clusterProfiler::summarise() +colnames(vertices) <- "name" + +# Build graph data structure and add network basic attributes, save data +# Building a graph data structure +graph <- graph_from_data_frame(cor.data, vertices = vertices, directed = FALSE) +E(graph)$weight <- abs(E(graph)$r) +V(graph)$label <- V(graph)$name + +# 利用“louvain”算法进行进行聚类群划分 +# calculate community membership and modularity of networks +patients.clusters <- cluster_louvain(graph) +V(graph)$Cluster <- patients.clusters$membership + +# save data +#write_graph(graph, file = paste(opts$output, "Centenarians_01.graphml", sep=""), format="graphml") +write_graph( + graph, + file = file.path( + opts$output, + paste0(opts$Group, "_01.graphml") + ), + format = "graphml" +) From 0c382dec2a504c9dba9a365b461af61032cfa60b Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:07:25 +0800 Subject: [PATCH 08/19] Add files via upload --- script/SparCC_visualization.R | 446 ++++++++++++++------------ script/cluster_from_Cytoscape.R | 542 ++++++++++++++++++++++++++++++++ 2 files changed, 791 insertions(+), 197 deletions(-) create mode 100644 script/cluster_from_Cytoscape.R diff --git a/script/SparCC_visualization.R b/script/SparCC_visualization.R index 2f3c4ef..e7bf210 100644 --- a/script/SparCC_visualization.R +++ b/script/SparCC_visualization.R @@ -1,197 +1,249 @@ -#!/usr/bin/env Rscript - -# Copyright 2024 De-feng Bai - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 - -# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 -# Set Work Directory - Choose Directory / To Source File Location - -# Clean enviroment object -rm(list=ls()) - -# 1.1 简介 Introduction #---- - -# 程序功能：SparCC可视化 -# Functions: Sparcc Visualization - - -options(warn = -1) # Turn off warning - -# 1.2 参数 Parameters #---- -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--Correlation"), type="character", default="metaphlan4/sxtr_cov_mat_Centenarians.tsv", - help="Metaphlan4 species table"), - make_option(c("-P", "--Pvalue"), type="character", default="metaphlan4/sxtr_pvals_Centenarians.two_sided.tsv", - help="Unfiltered OTU table [default %default]"), - make_option(c("-r", "--output"), type="character", default="metaphlan4/", - help="Output file for SparCC analysis in different groups [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# Install related packages -if (FALSE){ - source("https://bioconductor.org/biocLite.R") - biocLite(c("reshape2","ggplot2","ggprism","dplyr","plyr", - "igraph","tidyverse","ggraph")) # ,"vegan" -} -# load related packages -suppressWarnings(suppressMessages(library("reshape2"))) -suppressWarnings(suppressMessages(library("ggplot2"))) -suppressWarnings(suppressMessages(library("ggprism"))) -# suppressWarnings(suppressMessages(library("vegan"))) -suppressWarnings(suppressMessages(library("dplyr"))) -suppressWarnings(suppressMessages(library("plyr"))) -suppressWarnings(suppressMessages(library("igraph"))) -suppressWarnings(suppressMessages(library("tidyverse"))) -suppressWarnings(suppressMessages(library("ggraph"))) - - -#r.cor <- read.table("data3/r.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) -r.cor <- read.table(opts$Correlation, sep="\t", header=T, check.names=F,row.names = 1) -#p.cor <- read.table("data3/p.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) -p.cor <- read.table(opts$Pvalue, sep="\t", header=T, check.names=F,row.names = 1) - -r.cor[p.cor>0.05] <- 0 - -# Build network connection attributes and node attributes -# Convert data to long format for merging and add connection properties -r.cor$from = rownames(r.cor) -p.cor$from = rownames(p.cor) -p_value <- p.cor %>% - gather(key = "to", value = "p", -from) %>% - data.frame() -p_value$FDR <- p.adjust(p_value$p,"BH") -p_value <- p_value[, -3] - -cor.data<- r.cor %>% - gather(key = "to", value = "r", -from) %>% - data.frame() %>% - left_join(p_value, by=c("from","to")) -cor.data <- as.data.frame(cor.data) -cor.data <- cor.data[cor.data$FDR <= 0.05 & cor.data$from != cor.data$to, ] -cor.data <- cor.data[abs(cor.data$r) >= 0.6 & cor.data$from != cor.data$to, ] -cor.data <- cor.data %>% - plyr::mutate( - linecolor = ifelse(r > 0,"positive","negative"), - linesize = abs(r) - ) - -# Set node properties -vertices <- c(as.character(cor.data$from),as.character(cor.data$to)) %>% - as_tibble() %>% - group_by(value) %>% - clusterProfiler::summarise() -colnames(vertices) <- "name" - -# Build graph data structure and add network basic attributes, save data -# Building a graph data structure -graph <- graph_from_data_frame(cor.data, vertices = vertices, directed = FALSE) -E(graph)$weight <- abs(E(graph)$r) -V(graph)$label <- V(graph)$name - -# 利用“louvain”算法进行进行聚类群划分 -# calculate community membership and modularity of networks -patients.clusters <- cluster_louvain(graph) -V(graph)$Cluster <- patients.clusters$membership - -# save data -write_graph(graph, file = paste(opts$output, "Centenarians_01.graphml", sep=""), format="graphml") - -# 可视化方式1：基于Gephi软件进行可视化 https://gephi.org/ -# Visualized in Gephi software -# The same procedure for healthy group - -# healthy组和patients组相同 -# healthy.clusters <- cluster_louvain(healthy.igraph.s) -# V(healthy.igraph.s)$Cluster <- healthy.clusters$membership - -# 可视化方式2：利用igraph进行可视化 -g <- graph -# 准备网络图布局数据 -# Preparing network diagram layout data。 -layout1 <- layout_in_circle(g) -layout5 <- layout_with_graphopt(g) - -# 设置绘图颜色 -# Setting the drawing color -#color <- c("#879b56","#ce77ad","#5ea6c2") - -color = c("#d2da93","#5196d5","#00ceff","#ff630d","#35978b", - "#e5acd7","#77aecd","#ec8181","#dfc6a5","#e50719", - "#d27e43","#8a4984","#fe5094","#8d342e","#f94e54", - "#ffad00","#36999d","#00fc8d","#b64aa0","#9b82e1", - "#fec44f","#e0f3db","#fa9fb5","#c994c7","#807dba", - "#ec7014","#a6bddb","#ef3b2c","#fe9929","#a1d99b") - -names(color) <- unique(V(g)$Cluster) -V(g)$point.col <- color[match(V(g)$Cluster,names(color))] - -# 边颜色按照相关性正负设置 -# The edge color is set according to the positive or negative correlation -E(g)$color <- ifelse(E(g)$linecolor == "positive","#ff878c","#5ea6c2") - -pdf(file=paste(opts$output, "network_group_Centenarians_cluster.pdf", sep=""), width=10, height=12) -par(mar=c(5,2,1,2)) -plot.igraph(g, layout=layout5, - vertex.color=V(g)$point.col, - vertex.border=V(g)$point.col, - vertex.size=6, - vertex.frame.color="white", - vertex.label=g$name, - vertex.label.cex=0.8, - vertex.label.dist=0, - vertex.label.degree = pi/2, - vertex.label.col="black", - edge.arrow.size=0.5, - edge.width=abs(E(g)$r)*15, - edge.curved = FALSE -) - -# 设置图例 -legend( - title = "Cluster", - list(x = min(layout1[,1])-0.05, - y = min(layout1[,2])-0.05), - legend = c(unique(V(g)$Cluster)), - fill = color, - #pch=1 -) - -legend( - title = "|r-value|", - list(x = min(layout1[,1])+0.6, - y = min(layout1[,2])-0.05), - legend = c(0.2,0.4,0.6,0.8,1.0), - col = "black", - lty=1, - lwd=c(0.2,0.4,0.6,0.8,1.0)*4, -) - -legend( - title = "Correlation (±)", - list(x = min(layout1[,1])+1.0, - y = min(layout1[,2])-0.05), - legend = c("positive","negative"), - col = c("#ff878c",rgb(0,147,0,maxColorValue = 255)), - lty=1, - lwd=1 -) -dev.off() - +#!/usr/bin/env Rscript + +# Copyright 2024 De-feng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# Clean enviroment object +rm(list=ls()) + +# 1.1 简介 Introduction #---- + +# 程序功能：SparCC可视化 +# Functions: Sparcc Visualization + + +options(warn = -1) # Turn off warning + +# 1.2 参数 Parameters #---- +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-R", "--Correlation"), type="character", default="metaphlan4/sxtr_cov_mat_Centenarians.tsv", + help="Metaphlan4 species table"), + make_option(c("-P", "--Pvalue"), type="character", default="metaphlan4/sxtr_pvals_Centenarians.two_sided.tsv", + help="Unfiltered OTU table [default %default]"), + make_option("--Group", type="character"), + make_option(c("-r", "--output"), type="character", default="metaphlan4/", + help="Output file for SparCC analysis in different groups [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + + +# Install related packages +options(repos = c(CRAN = "https://cloud.r-project.org")) + +packages <- c( + "reshape2", + "ggplot2", + "ggprism", + "dplyr", + "plyr", + "igraph", + "tidyverse", + "ggraph", + "magrittr" +) + +installed <- rownames(installed.packages()) +missing_pkgs <- packages[!(packages %in% installed)] + +if(length(missing_pkgs) > 0){ + install.packages(missing_pkgs, dependencies = TRUE) +} + +if (!requireNamespace("BiocManager", quietly = TRUE)) + install.packages("BiocManager") + +pkgs <- c( + "clusterProfiler", + "DOSE", + "enrichplot", + "ComplexHeatmap", + "pathview", + "KEGGREST" +) + +for(pkg in pkgs){ + if(!requireNamespace(pkg, quietly = TRUE)){ + BiocManager::install(pkg, update = FALSE, ask = FALSE) + } +} + +# if (!requireNamespace("clusterProfiler", quietly = TRUE)) +# BiocManager::install("clusterProfiler", update = FALSE, ask = FALSE) + +#library(clusterProfiler) + +# load related packages +suppressWarnings(suppressMessages(library("reshape2"))) +suppressWarnings(suppressMessages(library("ggplot2"))) +suppressWarnings(suppressMessages(library("ggprism"))) +# suppressWarnings(suppressMessages(library("vegan"))) +suppressWarnings(suppressMessages(library("dplyr"))) +suppressWarnings(suppressMessages(library("plyr"))) +suppressWarnings(suppressMessages(library("igraph"))) +suppressWarnings(suppressMessages(library("tidyverse"))) +suppressWarnings(suppressMessages(library("ggraph"))) +suppressWarnings(suppressMessages(library("magrittr"))) +suppressWarnings(suppressMessages(library("clusterProfiler"))) + + +#r.cor <- read.table("data3/r.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) +r.cor <- read.table(opts$Correlation, sep="\t", header=T, check.names=F,row.names = 1) +#p.cor <- read.table("data3/p.cor_count_Y.txt", sep="\t", header=T, check.names=F,row.names = 1) +p.cor <- read.table(opts$Pvalue, sep="\t", header=T, check.names=F,row.names = 1) + +r.cor[p.cor>0.05] <- 0 + +# Build network connection attributes and node attributes +# Convert data to long format for merging and add connection properties +r.cor$from = rownames(r.cor) +p.cor$from = rownames(p.cor) +p_value <- p.cor %>% + gather(key = "to", value = "p", -from) %>% + data.frame() +p_value$FDR <- p.adjust(p_value$p,"BH") +p_value <- p_value[, -3] + +cor.data<- r.cor %>% + gather(key = "to", value = "r", -from) %>% + data.frame() %>% + left_join(p_value, by=c("from","to")) +cor.data <- as.data.frame(cor.data) +cor.data <- cor.data[cor.data$FDR <= 0.05 & cor.data$from != cor.data$to, ] +cor.data <- cor.data[abs(cor.data$r) >= 0.6 & cor.data$from != cor.data$to, ] +cor.data <- cor.data %>% + plyr::mutate( + linecolor = ifelse(r > 0,"positive","negative"), + linesize = abs(r) + ) + +# Set node properties +vertices <- c(as.character(cor.data$from),as.character(cor.data$to)) %>% + as_tibble() %>% + group_by(value) %>% + clusterProfiler::summarise() +colnames(vertices) <- "name" + +# Build graph data structure and add network basic attributes, save data +# Building a graph data structure +graph <- graph_from_data_frame(cor.data, vertices = vertices, directed = FALSE) +E(graph)$weight <- abs(E(graph)$r) +V(graph)$label <- V(graph)$name + +# 利用“louvain”算法进行进行聚类群划分 +# calculate community membership and modularity of networks +patients.clusters <- cluster_louvain(graph) +V(graph)$Cluster <- patients.clusters$membership + +# save data +# write_graph(graph, file = paste(opts$output, "Centenarians_01.graphml", sep=""), format="graphml") + +# 可视化方式1：基于Gephi软件进行可视化 https://gephi.org/ +# Visualized in Gephi software +# The same procedure for healthy group + +# healthy组和patients组相同 +# healthy.clusters <- cluster_louvain(healthy.igraph.s) +# V(healthy.igraph.s)$Cluster <- healthy.clusters$membership + +# 可视化方式2：利用igraph进行可视化 +g <- graph +# 准备网络图布局数据 +# Preparing network diagram layout data。 +layout1 <- layout_in_circle(g) +layout5 <- layout_with_graphopt(g) + +# 设置绘图颜色 +# Setting the drawing color +#color <- c("#879b56","#ce77ad","#5ea6c2") + +color = c("#d2da93","#5196d5","#00ceff","#ff630d","#35978b", + "#e5acd7","#77aecd","#ec8181","#dfc6a5","#e50719", + "#d27e43","#8a4984","#fe5094","#8d342e","#f94e54", + "#ffad00","#36999d","#00fc8d","#b64aa0","#9b82e1", + "#fec44f","#e0f3db","#fa9fb5","#c994c7","#807dba", + "#ec7014","#a6bddb","#ef3b2c","#fe9929","#a1d99b") + +names(color) <- unique(V(g)$Cluster) +V(g)$point.col <- color[match(V(g)$Cluster,names(color))] + +# 边颜色按照相关性正负设置 +# The edge color is set according to the positive or negative correlation +E(g)$color <- ifelse(E(g)$linecolor == "positive","#ff878c","#5ea6c2") + +#pdf(file=paste(opts$output, "network_group_Centenarians_cluster.pdf", sep=""), width=10, height=12) +pdf( + file = paste0( + opts$output, + "/network_group_", + opts$Group, + "_cluster.pdf" + ), + width = 10, + height = 12 +) +par(mar=c(5,2,1,2)) +plot.igraph(g, layout=layout5, + vertex.color=V(g)$point.col, + vertex.border=V(g)$point.col, + vertex.size=6, + vertex.frame.color="white", + vertex.label=g$name, + vertex.label.cex=0.8, + vertex.label.dist=0, + vertex.label.degree = pi/2, + vertex.label.col="black", + edge.arrow.size=0.5, + edge.width=abs(E(g)$r)*15, + edge.curved = FALSE +) + +# 设置图例 +legend( + title = "Cluster", + list(x = min(layout1[,1])-0.05, + y = min(layout1[,2])-0.05), + legend = c(unique(V(g)$Cluster)), + fill = color, + #pch=1 +) + +legend( + title = "|r-value|", + list(x = min(layout1[,1])+0.6, + y = min(layout1[,2])-0.05), + legend = c(0.2,0.4,0.6,0.8,1.0), + col = "black", + lty=1, + lwd=c(0.2,0.4,0.6,0.8,1.0)*4, +) + +legend( + title = "Correlation (±)", + list(x = min(layout1[,1])+1.0, + y = min(layout1[,2])-0.05), + legend = c("positive","negative"), + col = c("#ff878c",rgb(0,147,0,maxColorValue = 255)), + lty=1, + lwd=1 +) +dev.off() + diff --git a/script/cluster_from_Cytoscape.R b/script/cluster_from_Cytoscape.R new file mode 100644 index 0000000..5a977f5 --- /dev/null +++ b/script/cluster_from_Cytoscape.R @@ -0,0 +1,542 @@ +# #!/usr/bin/env Rscript +# +# if (!requireNamespace("optparse", quietly = TRUE)) { +# install.packages("optparse", +# repos = "https://cloud.r-project.org") +# } +# +# library(optparse) +# +# +# # 检查并安装软件包 +# +# options( +# repos = c( +# CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/" +# ) +# ) +# +# options( +# BioC_mirror = +# "https://mirrors.tuna.tsinghua.edu.cn/bioconductor" +# ) +# +# # ===================================================== +# # Install and load packages +# # ===================================================== +# +# cran_packages <- c( +# "reshape2", +# "ggplot2", +# "dplyr", +# "purrr", +# "igraph" +# ) +# +# bioc_packages <- c( +# "impute", +# "preprocessCore", +# "WGCNA" +# ) +# +# # ---------- CRAN ---------- +# for(pkg in cran_packages){ +# +# if(!requireNamespace(pkg, quietly = TRUE)){ +# +# cat("Installing CRAN package:", pkg, "\n") +# +# install.packages( +# pkg, +# repos = "https://cloud.r-project.org" +# ) +# } +# +# suppressPackageStartupMessages( +# library(pkg, +# character.only = TRUE) +# ) +# } +# +# # ---------- Bioconductor ---------- +# if(length(bioc_packages) > 0){ +# +# if(!requireNamespace("BiocManager", +# quietly = TRUE)){ +# +# install.packages( +# "BiocManager", +# repos = "https://cloud.r-project.org" +# ) +# } +# +# for(pkg in bioc_packages){ +# +# if(!requireNamespace(pkg, +# quietly = TRUE)){ +# +# cat("Installing Bioconductor package:", pkg, "\n") +# +# BiocManager::install( +# pkg, +# ask = FALSE, +# update = FALSE +# ) +# } +# +# suppressPackageStartupMessages( +# library(pkg, +# character.only = TRUE) +# ) +# } +# } +# +# +# #---- +# #5.load ccluster from Cytoscape +# mysignode_stable_ccluster=read.table("./coverm/ccCluster.txt",header = T,sep = "\t",check.names = F,row.names = 1) +# mysignode_stable$CCcluster=mysignode_stable_ccluster$ccCluster[match(mysignode_stable$node,rownames(mysignode_stable_ccluster))] +# mysignode_stable$CCcluster=paste("C",mysignode_stable$CCcluster,sep = "") +# #5.1 find clusters in C1 +# C1genome=mysignode_stable$node[which(mysignode_stable$CCcluster=="C1")] +# ind=which(mysigedge_stable$Source%in%C1genome&mysigedge_stable$Target%in%C1genome) +# tmpdata=graph.data.frame(mysigedge_stable[ind,],directed = F) +# tmpdata=get.adjacency(tmpdata,attr="NorPv2",sparse = F) +# tmpdist=1-tmpdata +# tmpcluster=hclust(as.dist(tmpdist),method = "average") +# hc=hclust(as.dist(tmpdist),method = "average") +# b=plot(hc,hang=-1,cex=0.8) +# mynamicmods=cutreeDynamic(dendro = hc,distM=tmpdist) +# table(mynamicmods) +# +# mynamicmods=data.frame(ID=colnames(tmpdata),g=mynamicmods) +# mynamicmods$g=plyr::mapvalues(mynamicmods$g,from=c(0:25),to=LETTERS) +# mynamicmods$g=paste("C1",mynamicmods$g,sep = "") +# mynamicmods$ID_raw=rownames(myallGTDBTK)[match(mynamicmods$ID,myallGTDBTK$ID)] + + + +#!/usr/bin/env Rscript + +# ===================================================== +# Install and load packages +# ===================================================== + +if (!requireNamespace("optparse", quietly = TRUE)) { + install.packages( + "optparse", + repos = "https://cloud.r-project.org" + ) +} + +library(optparse) + +options( + repos = c( + CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/" + ) +) + +options( + BioC_mirror = + "https://mirrors.tuna.tsinghua.edu.cn/bioconductor" +) + +cran_packages <- c( + "reshape2", + "ggplot2", + "dplyr", + "purrr", + "igraph", + "plyr" +) + +bioc_packages <- c( + "impute", + "preprocessCore", + "WGCNA" +) + +# ---------- CRAN ---------- +for(pkg in cran_packages){ + + if(!requireNamespace(pkg, quietly = TRUE)){ + + cat("Installing CRAN package:", pkg, "\n") + + install.packages( + pkg, + repos = "https://cloud.r-project.org" + ) + } + + suppressPackageStartupMessages( + library( + pkg, + character.only = TRUE + ) + ) +} + +# ---------- Bioconductor ---------- +if(!requireNamespace( + "BiocManager", + quietly = TRUE +)){ + install.packages( + "BiocManager", + repos = "https://cloud.r-project.org" + ) +} + +for(pkg in bioc_packages){ + + if(!requireNamespace( + pkg, + quietly = TRUE + )){ + + cat( + "Installing Bioconductor package:", + pkg, + "\n" + ) + + BiocManager::install( + pkg, + ask = FALSE, + update = FALSE + ) + } + + suppressPackageStartupMessages( + library( + pkg, + character.only = TRUE + ) + ) +} + +# ===================================================== +# Parameters +# ===================================================== + +option_list <- list( + + make_option( + c("--ccCluster"), + type = "character", + default = "./coverm/ccCluster.txt", + help = "ccCluster.txt" + ), + + make_option( + c("--node"), + type = "character", + default = "./coverm/372stable.node.txt", + help = "stable.node2.txt" + ), + + make_option( + c("--edge"), + type = "character", + default = "./coverm/372stable.edge.txt", + help = "stable.edge2.txt" + ), + + make_option( + c("--gtdb"), + type = "character", + default = "./coverm/GTDBTK.txt", + help = "myallGTDBTK.txt" + ), + + make_option( + c("--cluster"), + type = "character", + default = "C1", + help = "Cluster ID (default=C1)" + ), + + make_option( + c("-o","--output"), + type = "character", + default = "./coverm", + help = "Output directory" + ) +) + +opts <- parse_args( + OptionParser( + option_list = option_list + ) +) + +dir.create( + opts$output, + recursive = TRUE, + showWarnings = FALSE +) + +# ===================================================== +# Read files +# ===================================================== + +cat("Reading files ...\n") + +mysignode_stable_ccluster <- read.table( + opts$ccCluster, + header = TRUE, + sep = "\t", + check.names = FALSE, + row.names = 1 +) + +mysignode_stable <- read.table( + opts$node, + header = TRUE, + sep = "\t", + check.names = FALSE +) + +mysigedge_stable <- read.table( + opts$edge, + header = TRUE, + sep = "\t", + check.names = FALSE +) + +myallGTDBTK <- read.table( + opts$gtdb, + header = TRUE, + sep = "\t", + check.names = FALSE, + row.names = 1 +) + +# ===================================================== +# Merge cluster information +# ===================================================== + +mysignode_stable$CCcluster <- + mysignode_stable_ccluster$ccCluster[ + match( + mysignode_stable$node, + rownames(mysignode_stable_ccluster) + ) + ] + +mysignode_stable$CCcluster <- + paste0( + "C", + mysignode_stable$CCcluster + ) + +# ===================================================== +# Extract target cluster +# ===================================================== + +cat( + "Processing cluster:", + opts$cluster, + "\n" +) + +target_genome <- + mysignode_stable$node[ + mysignode_stable$CCcluster == + opts$cluster + ] + +if(length(target_genome) < 2){ + + stop( + "Less than 2 genomes found in ", + opts$cluster + ) +} + +ind <- which( + mysigedge_stable$Source %in% + target_genome & + mysigedge_stable$Target %in% + target_genome +) + +if(length(ind) == 0){ + + stop( + "No edges found in ", + opts$cluster + ) +} + +# ===================================================== +# Build adjacency matrix +# ===================================================== + +cat("Building network ...\n") + +tmpgraph <- graph.data.frame( + mysigedge_stable[ind,], + directed = FALSE +) + +tmpdata <- get.adjacency( + tmpgraph, + attr = "NorPv2", + sparse = FALSE +) + +tmpdist <- 1 - tmpdata + +# ===================================================== +# Hierarchical clustering +# ===================================================== + +cat("Running hclust ...\n") + +hc <- hclust( + as.dist(tmpdist), + method = "average" +) + +# ===================================================== +# Save clustering tree +# ===================================================== + +pdf( + file.path( + opts$output, + paste0( + opts$cluster, + "_hclust_tree.pdf" + ) + ), + width = 10, + height = 8 +) + +plot( + hc, + hang = -1, + cex = 0.8, + main = paste( + opts$cluster, + "Hierarchical Clustering" + ) +) + +dev.off() + +# ===================================================== +# Dynamic Tree Cut +# ===================================================== + +cat( + "Running dynamic tree cut ...\n" +) + +mynamicmods <- cutreeDynamic( + dendro = hc, + distM = tmpdist +) + +print( + table(mynamicmods) +) + +# ===================================================== +# Generate module table +# ===================================================== + +module_letters <- c( + LETTERS, + paste0( + rep(LETTERS, each = 26), + LETTERS + ) +) + +module_ids <- + sort(unique(mynamicmods)) + +module_map <- setNames( + module_letters[ + seq_along(module_ids) + ], + module_ids +) + +mynamicmods <- data.frame( + ID = colnames(tmpdata), + Module = module_map[ + as.character(mynamicmods) + ], + stringsAsFactors = FALSE +) + +mynamicmods$Module <- + paste0( + opts$cluster, + "_", + mynamicmods$Module + ) + +if("ID" %in% colnames(myallGTDBTK)){ + + mynamicmods$ID_raw <- + rownames(myallGTDBTK)[ + match( + mynamicmods$ID, + myallGTDBTK$ID + ) + ] +} + +# ===================================================== +# Save results +# ===================================================== + +write.table( + mynamicmods, + file = file.path( + opts$output, + paste0( + opts$cluster, + "_dynamic_modules.txt" + ) + ), + sep = "\t", + quote = FALSE, + row.names = FALSE +) + +cat( + "\nFinished!\n" +) + +cat( + "Tree : ", + file.path( + opts$output, + paste0( + opts$cluster, + "_hclust_tree.pdf" + ) + ), + "\n" +) + +cat( + "Module table : ", + file.path( + opts$output, + paste0( + opts$cluster, + "_dynamic_modules.txt" + ) + ), + "\n" +) + From cdfebf6414aa7c536f2a574d2194feaca261baeb Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:16:07 +0800 Subject: [PATCH 09/19] Update scripts for 2StatPlot.sh --- script/KEGG_Orthology_data_extract.sh | 98 ++++++ script/function_sankey.R | 218 ++++++++------ script/merged_tax.sh | 82 ++--- script/phylogenetic_anno_data.sh | 123 ++++++++ script/phylogenetic_tree.R | 379 ++++++++++++------------ script/phylogenetic_tree_env.R | 411 +++++++++++++++----------- 6 files changed, 824 insertions(+), 487 deletions(-) create mode 100644 script/KEGG_Orthology_data_extract.sh create mode 100644 script/phylogenetic_anno_data.sh diff --git a/script/KEGG_Orthology_data_extract.sh b/script/KEGG_Orthology_data_extract.sh new file mode 100644 index 0000000..7194a9b --- /dev/null +++ b/script/KEGG_Orthology_data_extract.sh @@ -0,0 +1,98 @@ +#!/bin/bash + +OUTDIR="humann4" + +mkdir -p ${OUTDIR} + +echo "Downloading KEGG hierarchy..." + +wget -q -O ${OUTDIR}/ko00001.keg \ +"https://www.kegg.jp/kegg-bin/download_htext?htext=ko00001.keg&format=htext" + +echo "Parsing hierarchy..." + +awk ' +BEGIN{ + OFS="\t" + + print "KO", + "Gene", + "Description", + "Level1_ID", + "Level1_Name", + "Level2_ID", + "Level2_Name", + "Level3_ID", + "Level3_Name" +} + +/^A/{ + line=$0 + sub(/^A/,"",line) + + match(line,/^[0-9]+/) + L1ID=substr(line,RSTART,RLENGTH) + + sub(/^[0-9]+ /,"",line) + L1NAME=line +} + +/^B/{ + line=$0 + sub(/^B[ ]+/,"",line) + + match(line,/^[0-9]+/) + L2ID=substr(line,RSTART,RLENGTH) + + sub(/^[0-9]+ /,"",line) + L2NAME=line +} + +/^C/{ + line=$0 + sub(/^C[ ]+/,"",line) + + match(line,/^[0-9]+/) + L3ID=substr(line,RSTART,RLENGTH) + + sub(/^[0-9]+ /,"",line) + sub(/ \[PATH:ko[0-9]+\]/,"",line) + + L3NAME=line +} + +/^D/{ + line=$0 + sub(/^D[ ]+/,"",line) + + split(line,a," ") + + KO=a[1] + + sub(/^K[0-9]+[ ]+/,"",line) + + split(line,b,";") + + Gene=b[1] + + Desc=b[2] + gsub(/^[ ]+/,"",Desc) + + print KO, + Gene, + Desc, + L1ID, + L1NAME, + L2ID, + L2NAME, + L3ID, + L3NAME +} +' ${OUTDIR}/ko00001.keg > ${OUTDIR}/KEGG_KO_Annotation.tsv + +echo "" +echo "Finished." +echo "" +echo "Generated files:" +echo "${OUTDIR}/ko00001.keg" +echo "${OUTDIR}/KEGG_KO_Annotation.tsv" diff --git a/script/function_sankey.R b/script/function_sankey.R index 99655f4..3dee4b6 100644 --- a/script/function_sankey.R +++ b/script/function_sankey.R @@ -1,92 +1,126 @@ -#!/usr/bin/env Rscript - -# Copyright 2024 De-feng Bai - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 - -# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 -# Set Work Directory - Choose Directory / To Source File Location - -# Clean enviroment object -rm(list=ls()) - -# 1.1 简介 Introduction #---- - -# 程序功能：不同层级功能通路和KO基因桑基图展示 -# Functions: Sankey diagram display of functional pathways and KO genes at different levels - - - -options(warn = -1) # Turn off warning - -# 1.2 参数 Parameters #---- -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} - -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="result/humann3/data_sankey7.txt", - help="MAGs taxonomy tables [default %default]"), - make_option(c("-o", "--output"), type="character", default="result/humann3/", - help="Output Sankey diagram display of functional pathways and KO genes at different levels [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# Install related packages -if (FALSE){ - source("https://bioconductor.org/biocLite.R") - biocLite(c("ggsankey","ggplot2","ggalluvial")) -} -# load related packages -suppressWarnings(suppressMessages(library("ggsankey"))) -suppressWarnings(suppressMessages(library("ggplot2"))) -suppressWarnings(suppressMessages(library("ggalluvial"))) - - -#df01 <- read.table(file = "data/data_sankey7.txt", sep = "\t", header = T, check.names = FALSE) -df01 <- read.table(file = opts$input, sep = "\t", header = T, check.names = FALSE) -data <- df01 -df <- to_lodes_form(data[,1:ncol(data)], - axes = 1:ncol(data), - id = "value") - -# Set color -col<- rep(c('#3690c0', '#f16913', '#238b45', '#ff6f81', '#fc9272', '#ffc2c0','#8c96c6', - '#bfd3e6', '#fae6f0', '#eb6fa6', '#ff88b5', '#00b1a5',"#ffa68f","#ffca75","#ccebc5","#7bccc4", - "#6baed6","#2171b5","#c6dbef","#448c99","#67a9cf","#b8d8c9","#88419d","#d4b9da","#fee6ce", - "#8f9898","#bfcfcb"), 6) - - -# Sankey diagram -p3 <- ggplot(df, aes(x = x, fill=stratum, label=stratum, - stratum = stratum, alluvium = value), width = 0.1)+ - geom_flow(width = 0.1, - curve_type = "sine", - alpha = 0.6, - color = 'white', - size = 0.05)+ - geom_stratum(width = 0.1, color = "white")+ - geom_text(stat = 'stratum', size = 3.5, color = 'black')+ - scale_fill_manual(values = col)+ - theme_void()+ - theme(legend.position = 'none', - text = element_text(size = 18)) -#p3 - -pdf(file=paste(opts$output, "MAGs_KEGG_Sankey.pdf", sep=""), height = 5.2, width = 11) -p3 -dev.off() - +#!/usr/bin/env Rscript + +# Copyright 2024 De-feng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# Clean enviroment object +rm(list=ls()) + +# 1.1 简介 Introduction #---- + +# 程序功能：不同层级功能通路和KO基因桑基图展示 +# Functions: Sankey diagram display of functional pathways and KO genes at different levels + + + +options(warn = -1) # Turn off warning + +# 1.2 参数 Parameters #---- +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} + +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-i", "--input"), type="character", default="result/humann3/data_sankey7.txt", + help="MAGs taxonomy tables [default %default]"), + make_option(c("-o", "--output"), type="character", default="result/humann3/", + help="Output Sankey diagram display of functional pathways and KO genes at different levels [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + + +# CRAN镜像 +site <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN" + +# CRAN包 +cran_pkgs <- c( + "ggplot2", + "ggalluvial", + "remotes" +) + +# 检查并安装CRAN包 +for(pkg in cran_pkgs){ + + if(!requireNamespace(pkg, quietly = TRUE)){ + + install.packages(pkg, repos = site) + + } + +} + +# ggsankey（GitHub安装） +if(!requireNamespace("ggsankey", quietly = TRUE)){ + + remotes::install_github("davidsjoberg/ggsankey") + +} + +# 加载包 +pkg_list <- c( + "ggsankey", + "ggplot2", + "ggalluvial" +) + +for(pkg in pkg_list){ + + suppressWarnings( + suppressMessages( + library(pkg, character.only = TRUE) + ) + ) + +} + + +df01 <- read.table(file = opts$input, sep = "\t", header = T, check.names = FALSE) +data <- df01 +df <- to_lodes_form(data[,1:ncol(data)], + axes = 1:ncol(data), + id = "value") + +# Set color +col<- rep(c('#3690c0', '#f16913', '#238b45', '#ff6f81', '#fc9272', '#ffc2c0','#8c96c6', + '#bfd3e6', '#fae6f0', '#eb6fa6', '#ff88b5', '#00b1a5',"#ffa68f","#ffca75","#ccebc5","#7bccc4", + "#6baed6","#2171b5","#c6dbef","#448c99","#67a9cf","#b8d8c9","#88419d","#d4b9da","#fee6ce", + "#8f9898","#bfcfcb"), 6) + + +# Sankey diagram +p3 <- ggplot(df, aes(x = x, fill=stratum, label=stratum, + stratum = stratum, alluvium = value), width = 0.1)+ + geom_flow(width = 0.1, + curve_type = "sine", + alpha = 0.6, + color = 'white', + size = 0.05)+ + geom_stratum(width = 0.1, color = "white")+ + geom_text(stat = 'stratum', size = 3.5, color = 'black')+ + scale_fill_manual(values = col)+ + theme_void()+ + theme(legend.position = 'none', + text = element_text(size = 18)) +#p3 + +pdf(file=paste(opts$output, "MAGs_KEGG_Sankey.pdf", sep=""), height = 5.2, width = 11) +p3 +dev.off() + diff --git a/script/merged_tax.sh b/script/merged_tax.sh index 51da07f..8850852 100644 --- a/script/merged_tax.sh +++ b/script/merged_tax.sh @@ -1,41 +1,41 @@ -#!/bin/bash - -# 创建输出文件并添加表头 -echo -e "Name\tCompleteness\tContamination\tGenome_Size\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies" > result/checkm2/parsed_taxonomy_filled.tsv - -# 处理数据并生成最终结果 -join -t $'\t' -1 1 -2 1 \ - <(tail -n +2 result/checkm2/quality_report.tsv | cut -f1-3,9 | sort -k1,1) \ - <(tail -n +2 temp/gtdb_classify/tax.bac120.summary.tsv | cut -f1-2 | sort -k1,1) | \ -awk -F '\t' -v OFS="\t" ' -{ - # 默认分类为 unclassified - domain="unclassified" - phylum="unclassified" - class="unclassified" - order="unclassified" - family="unclassified" - genus="unclassified" - species="unclassified" - - # 将分类信息分割到数组中 - split($5, tax, ";") - for (i in tax) { - if (tax[i] ~ /^d__/) domain=substr(tax[i], 4) - if (tax[i] ~ /^p__/) phylum=substr(tax[i], 4) - if (tax[i] ~ /^c__/) class=substr(tax[i], 4) - if (tax[i] ~ /^o__/) order=substr(tax[i], 4) - if (tax[i] ~ /^f__/) family=substr(tax[i], 4) - if (tax[i] ~ /^g__/) genus=substr(tax[i], 4) - if (tax[i] ~ /^s__/) species=substr(tax[i], 4) - } - # 打印最终结果 - print $1, $2, $3, $4, domain, phylum, class, order, family, genus, species -}' >> result/checkm2/parsed_taxonomy_filled.tsv - -# 将空白字段替换为 "unclassified" -awk -v OFS="\t" '{for(i=5;i<=11;i++) if($i=="") $i="unclassified"; print}' result/checkm2/parsed_taxonomy_filled.tsv > result/checkm2/taxonomy_merge.txt - -rm result/checkm2/parsed_taxonomy_filled.tsv -# result/checkm2/taxonomy_merge.csv - +#!/bin/bash + +# 创建输出文件并添加表头 +echo -e "Name\tCompleteness\tContamination\tGenome_Size\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies" > checkm2/parsed_taxonomy_filled.tsv + +# 处理数据并生成最终结果 +join -t $'\t' -1 1 -2 1 \ + <(tail -n +2 checkm2/quality_report.tsv | cut -f1-3,9 | sort -k1,1) \ + <(tail -n +2 ../temp/gtdb_classify/tax.bac120.summary.tsv | cut -f1-2 | sort -k1,1) | \ +awk -F '\t' -v OFS="\t" ' +{ + # 默认分类为 unclassified + domain="unclassified" + phylum="unclassified" + class="unclassified" + order="unclassified" + family="unclassified" + genus="unclassified" + species="unclassified" + + # 将分类信息分割到数组中 + split($5, tax, ";") + for (i in tax) { + if (tax[i] ~ /^d__/) domain=substr(tax[i], 4) + if (tax[i] ~ /^p__/) phylum=substr(tax[i], 4) + if (tax[i] ~ /^c__/) class=substr(tax[i], 4) + if (tax[i] ~ /^o__/) order=substr(tax[i], 4) + if (tax[i] ~ /^f__/) family=substr(tax[i], 4) + if (tax[i] ~ /^g__/) genus=substr(tax[i], 4) + if (tax[i] ~ /^s__/) species=substr(tax[i], 4) + } + # 打印最终结果 + print $1, $2, $3, $4, domain, phylum, class, order, family, genus, species +}' >> checkm2/parsed_taxonomy_filled.tsv + +# 将空白字段替换为 "unclassified" +awk -v OFS="\t" '{for(i=5;i<=11;i++) if($i=="") $i="unclassified"; print}' checkm2/parsed_taxonomy_filled.tsv > checkm2/taxonomy_merge.txt + +rm checkm2/parsed_taxonomy_filled.tsv +# result/checkm2/taxonomy_merge.csv + diff --git a/script/phylogenetic_anno_data.sh b/script/phylogenetic_anno_data.sh new file mode 100644 index 0000000..dde86ff --- /dev/null +++ b/script/phylogenetic_anno_data.sh @@ -0,0 +1,123 @@ +#!/bin/bash + +set -e + +# ========================================== +# 输入文件 +# ========================================== +GTDB="../temp/gtdb_classify/tax.bac120.summary.tsv" +CHECKM="checkm2/quality_report.tsv" +ABUND="coverm/abundance.tsv" + +# ========================================== +# 输出目录 +# ========================================== +OUTDIR="../temp/gtdb_infer" +mkdir -p ${OUTDIR} + +OUTFILE="${OUTDIR}/annotation.txt" + +# ========================================== +# 去除Windows换行符 +# ========================================== +sed -i 's/\r$//' ${GTDB} +sed -i 's/\r$//' ${CHECKM} +sed -i 's/\r$//' ${ABUND} + +# ========================================== +# 写表头 +# ========================================== +ABUND_HEADER=$(head -n 1 ${ABUND} | cut -f2-) + +echo -e "ID\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\tscore\tcompleteness\tcontamination\tstrain_heterogeneity\tsize\tN50\tcluster_members\t${ABUND_HEADER}" \ +> ${OUTFILE} + +# ========================================== +# 主程序 +# ========================================== +awk -F '\t' -v OFS='\t' ' + +######################################################## +# 读取 checkm +######################################################## +FNR==NR { + + if(NR==1) next + + id=$1 + + completeness=$2 + contamination=$3 + + score=completeness-(5*contamination) + + size=$9 + N50=$7 + + checkm[id]=score"\t"completeness"\t"contamination"\t0\t"size"\t"N50"\t1" + + next +} + +######################################################## +# 读取 abundance +######################################################## +FILENAME==ARGV[2] { + + if(FNR==1) next + + id=$1 + + abund="" + + for(i=2;i<=NF;i++){ + abund=abund"\t"$i + } + + abundance[id]=substr(abund,2) + + next +} + +######################################################## +# 读取 GTDB +######################################################## +FILENAME==ARGV[3] { + + if(FNR==1) next + + id=$1 + + split($2,a,";") + + domain=a[1] + phylum=a[2] + class=a[3] + order=a[4] + family=a[5] + genus=a[6] + species=a[7] + + print \ + id,\ + domain,\ + phylum,\ + class,\ + order,\ + family,\ + genus,\ + species,\ + checkm[id],\ + abundance[id] + +} +' \ +${CHECKM} \ +${ABUND} \ +${GTDB} \ +>> ${OUTFILE} + +echo "======================================" +echo "annotation file generated:" +echo "${OUTFILE}" +echo "======================================" diff --git a/script/phylogenetic_tree.R b/script/phylogenetic_tree.R index 31e9b42..16e6046 100644 --- a/script/phylogenetic_tree.R +++ b/script/phylogenetic_tree.R @@ -1,183 +1,196 @@ -#!/usr/bin/env Rscript - -# Copyright 2024 De-feng Bai - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 - -# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 -# Set Work Directory - Choose Directory / To Source File Location - -# Clean enviroment object -rm(list=ls()) - -# 1.1 简介 Introduction #---- - -# 程序功能：微生物物种系统发育树 -# Functions: Phylogenetic analysis for microbiota species - - -options(warn = -1) # Turn off warning - -# 1.2 参数 Parameters #---- -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} - -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-p", "--input"), type="character", default="result/gtdb_95/tax.unrooted.tree", - help="MAGs unrooted tree [default %default]"), - make_option(c("-a", "--annotation"), type="character", default="result/gtdb_95/annotation2.txt", - help="MAGs annotations [default %default]"), - make_option(c("-o", "--output"), type="character", default="result/gtdb_95/", - help="Output Phylogenetic tree for microbiota species [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# Install related packages -if (FALSE){ - source("https://bioconductor.org/biocLite.R") - biocLite(c("ggtreeExtra","ggtree","treeio","tidytree","ggstar","ggplot2", - "ggnewscale","TDbook")) -} -# load related packages -suppressWarnings(suppressMessages(library("ggtreeExtra"))) -suppressWarnings(suppressMessages(library("ggtree"))) -suppressWarnings(suppressMessages(library("treeio"))) -suppressWarnings(suppressMessages(library("tidytree"))) -suppressWarnings(suppressMessages(library("ggstar"))) -suppressWarnings(suppressMessages(library("ggplot2"))) -suppressWarnings(suppressMessages(library("ggnewscale"))) -suppressWarnings(suppressMessages(library("TDbook"))) - - - -# Load tree data -#tree <- read.tree("tax.unrooted.tree") -tree <- read.tree(opts$input) - -#data_anno <- read.table("annotation2.txt", header = TRUE, sep = "\t") -data_anno <- read.table(opts$annotation, header = TRUE, sep = "\t") - -data1 <- subset(data_anno, select = c(ID, Phylum)) -data1 <- data1[order(data1$Phylum), ] -color_palette <-c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0") -unique_phyla <- unique(data1$Phylum) -num_unique_phyla <- length(unique_phyla) -if (num_unique_phyla > length(color_palette)) { - color_mapping <- setNames(rep(color_palette, length.out = num_unique_phyla), unique_phyla) -} else { - color_mapping <- setNames(color_palette[1:num_unique_phyla], unique_phyla) -} -data1$color <- color_mapping[data1$Phylum] -data1 <- data1[order(data1$ID), ] -data2 <- subset(data_anno, select = c(ID, MY1, MY2, MY3, FY1, FY2, FY3)) -data2$Prevalence<- rowMeans(data2[, c("MY1", "MY2", "MY3","FY1", "FY2", "FY3")]) -data2$Group <- "Young" -data3 <- subset(data_anno, select = c(ID, ME1, ME2, ME3, FE1, FE2, FE3 )) -data3$Prevalence<- rowMeans(data3[, c("ME1", "ME2", "ME3", "FE1","FE2","FE3")]) -data3$Group <- "Elderly" -data4 <- subset(data_anno, select = c(ID, MC1, MC2, MC3, FC1, FC2, FC3)) -data4$Prevalence<- rowMeans(data4[, c("MC1", "MC2", "MC3", "FC1","FC2","FC3")]) -data4$Group <- "Centenarian" -#如果数据相差太大，可以做Min-Max Normalization - -min_max_normalize <- function(x) { - return((x - min(x)) / (max(x) - min(x))) -} -data2$Prevalence <- min_max_normalize(data2$Prevalence) -data3$Prevalence <- min_max_normalize(data3$Prevalence) -data4$Prevalence <- min_max_normalize(data4$Prevalence) -# Prepare the tree -data1$label <- data1$ID -tree2 <- full_join(tree, data1, by = "label") -# Create the circular layout tree -#p <- ggtree(tree2, layout = "fan", size = 0.15, open.angle = 0) -#p - -p2<- ggtree(tree2, - aes(color=Phylum),#支长颜色按照分组进行着色 - layout="fan",#进化树类型 - open.angle=0,#开口角度 - linewidth=0.75,#分支线条粗细 - show.legend = F)+ - #geom_tiplab(aes(color = label %in% df_label),#设定标签颜色根据筛选条件突出显示特定标签 - # size=3.5,#字体大小 - # align = T,#使用虚线连接标签与分支 - # linetype = 3,linewidth = 0.4,offset = 12.5,show.legend = F) - new_scale_fill() + - geom_fruit( - data=data2,#数据 - geom = geom_col,#绘图类型 - mapping = aes(y=ID, x= Prevalence, fill = Group), - offset = 0.1, - pwidth = 0.1, - width=0.5, - #show.legend = FALSE - )+ - scale_fill_manual( - #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), - values=c("#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), - #values=c("red","black"), - guide=guide_legend(keywidth=1, keyheight=1, order=2), - name="Cooling")+ - new_scale_fill() + - geom_fruit( - data=data3,#数据 - geom = geom_col,#绘图类型 - mapping = aes(y=ID, x= Prevalence, fill = Group), - offset = 0.01, - pwidth = 0.1, - width=0.5, - #show.legend = FALSE - )+ - scale_fill_manual( - #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), - values=c("#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), - #values=c("red","black"), - guide=guide_legend(keywidth=1, keyheight=1, order=2), - name="Mature")+ - new_scale_fill() + - geom_fruit( - data=data4,#数据 - geom = geom_col,#绘图类型 - mapping = aes(y=ID, x= Prevalence, fill = Group), - offset = 0.01, - pwidth = 0.1, - width=0.5, - #show.legend = FALSE - )+ - #scale_color_manual(values=c("black","#1aafd0","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238"))+ - scale_color_manual(values=c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0"))+ - new_scale_fill() + - geom_fruit( - data=data1, - geom=geom_tile, - mapping=aes(y=ID, fill=data1$Phylum), - #color="grey10", - width=0.05, - offset=0.1#, - #show.legend = FALSE - )+ - scale_fill_manual( - #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), - values=c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0"), - #values=c("red","black"), - guide=guide_legend(keywidth=1, keyheight=1, order=2), - name="Taxonomy") -pdf(file=paste(opts$output, "MAGs_phylogenetic_tree01.pdf", sep=""), height = 5.2, width = 11) -p2 -dev.off() - +#!/usr/bin/env Rscript + +# Copyright 2024 De-feng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# Clean enviroment object +rm(list=ls()) + +# 1.1 简介 Introduction #---- + +# 程序功能：微生物物种系统发育树 +# Functions: Phylogenetic analysis for microbiota species + + +options(warn = -1) # Turn off warning + +# 1.2 参数 Parameters #---- +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} + +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-p", "--input"), type="character", default="result/gtdb_95/tax.unrooted.tree", + help="MAGs unrooted tree [default %default]"), + make_option(c("-a", "--annotation"), type="character", default="result/gtdb_95/annotation2.txt", + help="MAGs annotations [default %default]"), + make_option(c("-o", "--output"), type="character", default="result/gtdb_95/", + help="Output Phylogenetic tree for microbiota species [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + + +# Install related packages +if (FALSE){ + source("https://bioconductor.org/biocLite.R") + biocLite(c("ggtreeExtra","ggtree","treeio","tidytree","ggstar","ggplot2", + "ggnewscale","TDbook")) +} +# load related packages +suppressWarnings(suppressMessages(library("ggtreeExtra"))) +suppressWarnings(suppressMessages(library("ggtree"))) +suppressWarnings(suppressMessages(library("treeio"))) +suppressWarnings(suppressMessages(library("tidytree"))) +suppressWarnings(suppressMessages(library("ggstar"))) +suppressWarnings(suppressMessages(library("ggplot2"))) +suppressWarnings(suppressMessages(library("ggnewscale"))) +suppressWarnings(suppressMessages(library("TDbook"))) + + + +# Load tree data +#tree <- read.tree("tax.unrooted.tree") +tree <- read.tree(opts$input) + +#data_anno <- read.table("annotation2.txt", header = TRUE, sep = "\t") +data_anno <- read.table(opts$annotation, header = TRUE, sep = "\t") + +data1 <- subset(data_anno, select = c(ID, Phylum)) +data1 <- data1[order(data1$Phylum), ] +color_palette <-c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0") +unique_phyla <- unique(data1$Phylum) +num_unique_phyla <- length(unique_phyla) +if (num_unique_phyla > length(color_palette)) { + color_mapping <- setNames(rep(color_palette, length.out = num_unique_phyla), unique_phyla) +} else { + color_mapping <- setNames(color_palette[1:num_unique_phyla], unique_phyla) +} +data1$color <- color_mapping[data1$Phylum] +data1 <- data1[order(data1$ID), ] +# data2 <- subset(data_anno, select = c(ID, MY1, MY2, MY3, FY1, FY2, FY3)) +# data2$Prevalence<- rowMeans(data2[, c("MY1", "MY2", "MY3","FY1", "FY2", "FY3")]) +# data2$Group <- "Young" + +data2 <- subset(data_anno, select = c(ID, Y1, Y2, Y3)) +data2$Prevalence<- rowMeans(data2[, c("Y1", "Y2", "Y3")]) +data2$Group <- "Young" + +# data3 <- subset(data_anno, select = c(ID, ME1, ME2, ME3, FE1, FE2, FE3 )) +# data3$Prevalence<- rowMeans(data3[, c("ME1", "ME2", "ME3", "FE1","FE2","FE3")]) +# data3$Group <- "Elderly" + +# data4 <- subset(data_anno, select = c(ID, MC1, MC2, MC3, FC1, FC2, FC3)) +# data4$Prevalence<- rowMeans(data4[, c("MC1", "MC2", "MC3", "FC1","FC2","FC3")]) +# data4$Group <- "Centenarian" + +data4 <- subset(data_anno, select = c(ID, C1, C2, C3)) +data4$Prevalence<- rowMeans(data4[, c("C1", "C2", "C3")]) +data4$Group <- "Centenarian" + +#如果数据相差太大，可以做Min-Max Normalization + +min_max_normalize <- function(x) { + return((x - min(x)) / (max(x) - min(x))) +} +data2$Prevalence <- min_max_normalize(data2$Prevalence) +#data3$Prevalence <- min_max_normalize(data3$Prevalence) +data4$Prevalence <- min_max_normalize(data4$Prevalence) +# Prepare the tree +data1$label <- data1$ID +tree2 <- full_join(tree, data1, by = "label") +# Create the circular layout tree +#p <- ggtree(tree2, layout = "fan", size = 0.15, open.angle = 0) +#p + +p2<- ggtree(tree2, + aes(color=Phylum),#支长颜色按照分组进行着色 + layout="fan",#进化树类型 + open.angle=0,#开口角度 + linewidth=0.75,#分支线条粗细 + show.legend = F)+ + #geom_tiplab(aes(color = label %in% df_label),#设定标签颜色根据筛选条件突出显示特定标签 + # size=3.5,#字体大小 + # align = T,#使用虚线连接标签与分支 + # linetype = 3,linewidth = 0.4,offset = 12.5,show.legend = F) + new_scale_fill() + + geom_fruit( + data=data2,#数据 + geom = geom_col,#绘图类型 + mapping = aes(y=ID, x= Prevalence, fill = Group), + offset = 0.1, + pwidth = 0.1, + width=0.5, + #show.legend = FALSE + )+ + scale_fill_manual( + #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), + values=c("#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), + #values=c("red","black"), + guide=guide_legend(keywidth=1, keyheight=1, order=2), + name="Cooling")+ + # new_scale_fill() + + # geom_fruit( + # data=data3,#数据 + # geom = geom_col,#绘图类型 + # mapping = aes(y=ID, x= Prevalence, fill = Group), + # offset = 0.01, + # pwidth = 0.1, + # width=0.5, + # #show.legend = FALSE + # )+ + # scale_fill_manual( + # #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), + # values=c("#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), + # #values=c("red","black"), + # guide=guide_legend(keywidth=1, keyheight=1, order=2), + # name="Mature")+ + new_scale_fill() + + geom_fruit( + data=data4,#数据 + geom = geom_col,#绘图类型 + mapping = aes(y=ID, x= Prevalence, fill = Group), + offset = 0.01, + pwidth = 0.1, + width=0.5, + #show.legend = FALSE + )+ + #scale_color_manual(values=c("black","#1aafd0","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238"))+ + scale_color_manual(values=c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0"))+ + new_scale_fill() + + geom_fruit( + data=data1, + geom=geom_tile, + mapping=aes(y=ID, fill=data1$Phylum), + #color="grey10", + width=0.05, + offset=0.1#, + #show.legend = FALSE + )+ + scale_fill_manual( + #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), + values=c("#bf812d", "#c7eae5", "#80cdc1", "#01665e", "#de77ae", "#d73027", "#4393c3","#e0e0e0"), + #values=c("red","black"), + guide=guide_legend(keywidth=1, keyheight=1, order=2), + name="Taxonomy") +pdf(file=paste(opts$output, "MAGs_phylogenetic_tree01.pdf", sep=""), height = 5.2, width = 11) +p2 +dev.off() + + + diff --git a/script/phylogenetic_tree_env.R b/script/phylogenetic_tree_env.R index 8e7e638..38981fd 100644 --- a/script/phylogenetic_tree_env.R +++ b/script/phylogenetic_tree_env.R @@ -1,171 +1,240 @@ -#!/usr/bin/env Rscript - -# Copyright 2024 De-feng Bai - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 - -# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 -# Set Work Directory - Choose Directory / To Source File Location - -# Clean enviroment object -rm(list=ls()) - -# 1.1 简介 Introduction #---- - -# 程序功能：微生物物种系统发育树 -# Functions: Phylogenetic analysis for microbiota species - - -options(warn = -1) # Turn off warning - -# 1.2 参数 Parameters #---- -# 设置清华源加速下载 -# (Optional) Set up Tsinghua Mirror to speed up download -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析包是否安装，没安装则安装，然后加载 -# Determine whether the command line parsing package is installed, install it if it is not installed, then load -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} - -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-p", "--input"), type="character", default="result/gtdb_95/tax.unrooted.tree", - help="MAGs unrooted tree [default %default]"), - make_option(c("-a", "--annotation"), type="character", default="result/gtdb_95/annotation2.txt", - help="MAGs annotations [default %default]"), - make_option(c("-o", "--output"), type="character", default="result/gtdb_95/", - help="Output Phylogenetic tree for microbiota species [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -print("You are using the following parameters:") -print(opts) - - -# Install related packages -if (FALSE){ - source("https://bioconductor.org/biocLite.R") - biocLite(c("ggtreeExtra","ggtree","treeio","tidytree","ggstar","ggplot2", - "ggnewscale","TDbook")) -} -# load related packages -suppressWarnings(suppressMessages(library("ggtreeExtra"))) -suppressWarnings(suppressMessages(library("ggtree"))) -suppressWarnings(suppressMessages(library("treeio"))) -suppressWarnings(suppressMessages(library("tidytree"))) -suppressWarnings(suppressMessages(library("ggstar"))) -suppressWarnings(suppressMessages(library("ggplot2"))) -suppressWarnings(suppressMessages(library("ggnewscale"))) -suppressWarnings(suppressMessages(library("TDbook"))) - - -# Load tree data -#tree <- read.tree("tax.unrooted.tree") -tree <- read.tree(opts$input) - -#data_anno <- read.table("annotation2.txt", header = TRUE, sep = "\t") -data_anno <- read.table(opts$annotation, header = TRUE, sep = "\t") - -#各个物种占比注释 - -data3 <- data_anno[, c("ID", "Phylum")] - -# 创建一个数据框存储门的名称和相应的百分比#后续可手动将比例尺加到图上 -phylum_counts <- table(data3$Phylum) -phylum_proportion <- prop.table(phylum_counts) -phylum_percentage <- phylum_proportion * 100 -phylum_percentage_df <- data.frame( - Phylum = names(phylum_percentage), - Percentage = as.numeric(phylum_percentage) -) - -phylum_percentage_sorted <- phylum_percentage_df[order(-phylum_percentage_df$Percentage), ] -# 提取占比多的Phylum -top_phylum <- phylum_percentage_sorted$Phylum[1:15] -# 创建新的列，将其他替换为 'others' -data4 <- data3 -data4$Phylum <- ifelse(data4$Phylum %in% top_phylum, data4$Phylum, "Others") -# 检查结果 -head(data4) - -data1 <- subset(data4, select = c(ID, Phylum)) -data1 <- data1[order(data1$Phylum), ] -color_palette <- c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", - "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806","#737373") -unique_phyla <- unique(data1$Phylum) -num_unique_phyla <- length(unique_phyla) -if (num_unique_phyla > length(color_palette)) { - color_mapping <- setNames(rep(color_palette, length.out = num_unique_phyla), unique_phyla) -} else { - color_mapping <- setNames(color_palette[1:num_unique_phyla], unique_phyla) -} -data1$color <- color_mapping[data1$Phylum] -data1 <- data1[order(data1$ID), ] -data2 <- subset(data_anno, select = c(ID, CO1, CO2, CO3, CO4,MA1, MA2, MA3, MA4,ME1, ME2, ME3, ME4,TH1, TH2, TH3, TH4)) -data2$Prevalence<- rowMeans(data2[, c("CO1", "CO2", "CO3", "CO4","MA1", "MA2", "MA3", "MA4","ME1", "ME2", "ME3", "ME4","TH1", "TH2", "TH3", "TH4")]) -data2$Group <- "ALL" - -# Prepare the tree -data1$label <- data1$ID -tree2 <- full_join(tree, data1, by = "label") -# Create the circular layout tree -p <- ggtree(tree2, layout = "fan", size = 0.15, open.angle = 0) - -p -top_phylum_sorted <- sort(top_phylum) -data1$Phylum <- factor(data1$Phylum, levels = c(top_phylum_sorted, "Others")) -p2<- ggtree(tree2, - aes(color=Phylum),#支长颜色按照分组进行着色 - layout="fan",#进化树类型 - open.angle=0,#开口角度 - linewidth=0.3,#分支线条粗细 - show.legend = F)+ - #geom_tiplab(aes(color = label %in% df_label),#设定标签颜色根据筛选条件突出显示特定标签 - # size=3.5,#字体大小 - # align = T,#使用虚线连接标签与分支 - # linetype = 3,linewidth = 0.4,offset = 12.5,show.legend = F)+ - #scale_color_manual(values=c("black","#1aafd0","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238"))+ - scale_color_manual(values=c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", - "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806", "#737373"))+ - new_scale_fill() + - geom_fruit( - data=data2,#数据 - geom = geom_col,#绘图类型 - mapping = aes(y=ID, x= Prevalence, fill = Group), - offset = 0.03, - pwidth = 0.5, - width=0.6, - #show.legend = FALSE - )+ - scale_fill_manual( - #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), - values=c("black","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), - #values=c("red","black"), - guide=guide_legend(keywidth=1, keyheight=1, order=2), - name="Cooling")+ - new_scale_fill() + - geom_fruit( - data=data1, - geom=geom_tile, - mapping=aes(y=ID, fill=data1$Phylum), - #color="grey10", - width=0.03, - offset=-0.15#, - #show.legend = FALSE - )+ - scale_fill_manual( - #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), - values=c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", - "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806","#737373"), - #values=c("red","black"), - guide=guide_legend(keywidth=1, keyheight=1, order=2), - name="Host taxonomy") -#p2 -pdf(file=paste(opts$output, "MAGs_phylogenetic_tree02.pdf", sep=""), height = 5.2, width = 11) -p2 -dev.off() - +#!/usr/bin/env Rscript + +# Copyright 2024 De-feng Bai + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, et al. 2021. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein & Cell 12: 315-330. https://doi.org/10.1007/s13238-020-00724-8 + +# 手动运行脚本，使用 Ctrl+Shift+H 或 Session 需要设置工作目录 +# Set Work Directory - Choose Directory / To Source File Location + +# Clean enviroment object +rm(list=ls()) + +# 1.1 简介 Introduction #---- + +# 程序功能：微生物物种系统发育树 +# Functions: Phylogenetic analysis for microbiota species + + +options(warn = -1) # Turn off warning + +# 1.2 参数 Parameters #---- +# 设置清华源加速下载 +# (Optional) Set up Tsinghua Mirror to speed up download +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析包是否安装，没安装则安装，然后加载 +# Determine whether the command line parsing package is installed, install it if it is not installed, then load +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} + +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-p", "--input"), type="character", default="result/gtdb_95/tax.unrooted.tree", + help="MAGs unrooted tree [default %default]"), + make_option(c("-a", "--annotation"), type="character", default="result/gtdb_95/annotation2.txt", + help="MAGs annotations [default %default]"), + make_option(c("-o", "--output"), type="character", default="result/gtdb_95/", + help="Output Phylogenetic tree for microbiota species [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +print("You are using the following parameters:") +print(opts) + + +# Install related packages +# if (FALSE){ +# source("https://bioconductor.org/biocLite.R") +# biocLite(c("ggtreeExtra","ggtree","treeio","tidytree","ggstar","ggplot2", +# "ggnewscale","TDbook")) +# } +# # load related packages +# suppressWarnings(suppressMessages(library("ggtreeExtra"))) +# suppressWarnings(suppressMessages(library("ggtree"))) +# suppressWarnings(suppressMessages(library("treeio"))) +# suppressWarnings(suppressMessages(library("tidytree"))) +# suppressWarnings(suppressMessages(library("ggstar"))) +# suppressWarnings(suppressMessages(library("ggplot2"))) +# suppressWarnings(suppressMessages(library("ggnewscale"))) +# suppressWarnings(suppressMessages(library("TDbook"))) + + +# 设置CRAN镜像 +site <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN" + +# CRAN包 +cran_pkgs <- c( + "ggplot2", + "ggstar", + "ggnewscale", + "BiocManager" +) + +# Bioconductor包 +bioc_pkgs <- c( + "ggtreeExtra", + "ggtree", + "treeio", + "tidytree", + "TDbook" +) + +#-------------------------------------------------- +# 安装CRAN包 +#-------------------------------------------------- +for(pkg in cran_pkgs){ + + if(!requireNamespace(pkg, quietly = TRUE)){ + + install.packages( + pkg, + repos = site + ) + + } + +} + +#-------------------------------------------------- +# 安装Bioconductor包 +#-------------------------------------------------- +for(pkg in bioc_pkgs){ + + if(!requireNamespace(pkg, quietly = TRUE)){ + + BiocManager::install( + pkg, + ask = FALSE, + update = FALSE + ) + + } + +} + +#-------------------------------------------------- +# 加载所有包 +#-------------------------------------------------- +all_pkgs <- c(cran_pkgs, bioc_pkgs) + +for(pkg in all_pkgs){ + + suppressWarnings( + suppressMessages( + library(pkg, character.only = TRUE) + ) + ) + +} + + +# Load tree data +#tree <- read.tree("tax.unrooted.tree") +tree <- read.tree(opts$input) + +#data_anno <- read.table("annotation2.txt", header = TRUE, sep = "\t") +data_anno <- read.table(opts$annotation, header = TRUE, sep = "\t") + +#各个物种占比注释 + +data3 <- data_anno[, c("ID", "Phylum")] + +# 创建一个数据框存储门的名称和相应的百分比#后续可手动将比例尺加到图上 +phylum_counts <- table(data3$Phylum) +phylum_proportion <- prop.table(phylum_counts) +phylum_percentage <- phylum_proportion * 100 +phylum_percentage_df <- data.frame( + Phylum = names(phylum_percentage), + Percentage = as.numeric(phylum_percentage) +) + +phylum_percentage_sorted <- phylum_percentage_df[order(-phylum_percentage_df$Percentage), ] +# 提取占比多的Phylum +top_phylum <- phylum_percentage_sorted$Phylum[1:15] +# 创建新的列，将其他替换为 'others' +data4 <- data3 +data4$Phylum <- ifelse(data4$Phylum %in% top_phylum, data4$Phylum, "Others") +# 检查结果 +head(data4) + +data1 <- subset(data4, select = c(ID, Phylum)) +data1 <- data1[order(data1$Phylum), ] +color_palette <- c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", + "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806","#737373") +unique_phyla <- unique(data1$Phylum) +num_unique_phyla <- length(unique_phyla) +if (num_unique_phyla > length(color_palette)) { + color_mapping <- setNames(rep(color_palette, length.out = num_unique_phyla), unique_phyla) +} else { + color_mapping <- setNames(color_palette[1:num_unique_phyla], unique_phyla) +} +data1$color <- color_mapping[data1$Phylum] +data1 <- data1[order(data1$ID), ] +data2 <- subset(data_anno, select = c(ID, CO1, CO2, CO3, CO4,MA1, MA2, MA3, MA4,ME1, ME2, ME3, ME4,TH1, TH2, TH3, TH4)) +data2$Prevalence<- rowMeans(data2[, c("CO1", "CO2", "CO3", "CO4","MA1", "MA2", "MA3", "MA4","ME1", "ME2", "ME3", "ME4","TH1", "TH2", "TH3", "TH4")]) +data2$Group <- "ALL" + +# Prepare the tree +data1$label <- data1$ID +tree2 <- full_join(tree, data1, by = "label") +# Create the circular layout tree +p <- ggtree(tree2, layout = "fan", size = 0.15, open.angle = 0) + +p +top_phylum_sorted <- sort(top_phylum) +data1$Phylum <- factor(data1$Phylum, levels = c(top_phylum_sorted, "Others")) +p2<- ggtree(tree2, + aes(color=Phylum),#支长颜色按照分组进行着色 + layout="fan",#进化树类型 + open.angle=0,#开口角度 + linewidth=0.3,#分支线条粗细 + show.legend = F)+ + #geom_tiplab(aes(color = label %in% df_label),#设定标签颜色根据筛选条件突出显示特定标签 + # size=3.5,#字体大小 + # align = T,#使用虚线连接标签与分支 + # linetype = 3,linewidth = 0.4,offset = 12.5,show.legend = F)+ + #scale_color_manual(values=c("black","#1aafd0","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238"))+ + scale_color_manual(values=c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", + "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806", "#737373"))+ + new_scale_fill() + + geom_fruit( + data=data2,#数据 + geom = geom_col,#绘图类型 + mapping = aes(y=ID, x= Prevalence, fill = Group), + offset = 0.03, + pwidth = 0.5, + width=0.6, + #show.legend = FALSE + )+ + scale_fill_manual( + #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), + values=c("black","#6a67ce","#ffb900","#fc636b","#aeb6b8","#e53238","lightblue"), + #values=c("red","black"), + guide=guide_legend(keywidth=1, keyheight=1, order=2), + name="Cooling")+ + new_scale_fill() + + geom_fruit( + data=data1, + geom=geom_tile, + mapping=aes(y=ID, fill=data1$Phylum), + #color="grey10", + width=0.03, + offset=-0.15#, + #show.legend = FALSE + )+ + scale_fill_manual( + #values=c("#4285f4", "#34a853", "#fbbc05","#ea4335"), + values=c("#fff7f3", "#deebf7","#9ecae1", "#4292c6", "#2171b5", "#08519c","#bc80bd", "#ffffb3", + "#bebada", "#fb8072", "#fdb462","#b3de69","#fccde5","#8dd3c7","#b35806","#737373"), + #values=c("red","black"), + guide=guide_legend(keywidth=1, keyheight=1, order=2), + name="Host taxonomy") +#p2 +pdf(file=paste(opts$output, "MAGs_phylogenetic_tree02.pdf", sep=""), height = 5.2, width = 11) +p2 +dev.off() + From b975ce2d348573587c1a8de367d719feb8376d72 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:30:28 +0800 Subject: [PATCH 10/19] Add files via upload From 458fd95432109ecbfc70dfdb9294da12490684bc Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:31:07 +0800 Subject: [PATCH 11/19] Add files via upload From c7a3642b3c15c0fc003d3eb14e499309ef57ff20 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:31:45 +0800 Subject: [PATCH 12/19] Add files via upload --- script/tax_stackplot_order.R | 298 +++++++++++++++++++++++------------ 1 file changed, 195 insertions(+), 103 deletions(-) diff --git a/script/tax_stackplot_order.R b/script/tax_stackplot_order.R index 77c8b58..145d628 100644 --- a/script/tax_stackplot_order.R +++ b/script/tax_stackplot_order.R @@ -10,6 +10,7 @@ # 更新 # 2024/11/12：增加按丰度排序堆叠柱状图功能 # 2025/11/27：规范代码 +# 2026/5/22：解决代码报错，规范图例显示 # 1.1 程序功能描述和主要步骤 @@ -101,7 +102,7 @@ if(opts$output==""){opts$output=opts$input} # 依赖包列表 package_list <- c( "ggplot2", "BiocManager", "optparse","patchwork","reshape2","magrittr", - "ggprism", "dplyr", "plyr" + "ggprism", "dplyr", "plyr","scales","grid" ) # 批量安装和加载 @@ -119,142 +120,233 @@ for (p in package_list) { ) } +# ========================================= +# 1. 读取输入文件 +# ========================================= + +# 实验设计 Metadata +metadata <- read.table( + opts$design, + header = TRUE, + row.names = 1, + sep = "\t", + comment.char = "", + stringsAsFactors = FALSE +) -# suppressWarnings(suppressMessages(library(amplicon))) -# suppressWarnings(suppressMessages(library(patchwork))) -# suppressWarnings(suppressMessages(library(reshape2))) -# suppressWarnings(suppressMessages(library(magrittr))) -# suppressWarnings(suppressMessages(library(ggplot2))) -# suppressWarnings(suppressMessages(library(ggprism))) -# suppressWarnings(suppressMessages(library(dplyr))) -# suppressWarnings(suppressMessages(library(plyr))) +# 物种组成矩阵 Taxonomy matrix +taxonomy <- read.table( + opts$input, + header = TRUE, + sep = "\t", + comment.char = "", + quote = "", + stringsAsFactors = FALSE +) +# ========================================= +# 2. 数据处理 +# ========================================= -# 3. 读取输入文件 +# 按 Taxonomy 合并重复物种 +data <- aggregate(. ~ Taxonomy, data = taxonomy, sum) -# 实验设计 Metadata -metadata = read.table(opts$design, header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F) +# 设置行名 +rownames(data) <- data$Taxonomy +data <- data[, -1] -# 物种组成矩阵 Taxonomy matrix -taxonomy = read.table(opts$input, header=T, sep="\t", comment.char="", quote = "") - - -# 4. 分析绘图 -# 计算微生物相对丰度之和，避免有重复统计 -data = taxonomy -data <- aggregate(.~ Taxonomy,data=data,sum) -rownames(data) = data$Taxonomy -data = data[, -1] - -# 计算相对丰度 -# Calculate relative abundance -data = apply(data , 2, function(x) x/sum(x)) -data = as.data.frame(data) - -# Decreased sort by abundance -# 相对丰度按降序排列 -mean_sort = data[(order(-rowSums(data))), ] -mean_sort = as.data.frame(mean_sort) -mean_sort2 = t(mean_sort) -mean_sort2 = mean_sort2[order(-mean_sort2[,1]),] -mean_sort3 = t(mean_sort2) -mean_sort3 = as.data.frame(mean_sort3) - -# Phylum水平展示前5个 -# Top 5 -other = colSums(mean_sort3[opts$legend:dim(mean_sort3)[1], ]) -mean_sort3 = mean_sort3[(opts$legend - 1):1, ] -mean_sort3 = rbind(other,mean_sort3) -rownames(mean_sort3)[1] = c("others") -mean_sort3 = as.data.frame(mean_sort3) - -# Add taxonomy -# 加入微生物分类信息 -mean_sort3$tax = rownames(mean_sort3) -data_all = as.data.frame(reshape2::melt(mean_sort3, id.vars = c("tax"))) -data_all$group = data_all$variable -data_all$group = as.character(data_all$group) -data_all$group = gsub("[0-9]","", data_all$group) - -# 从 metadata 建立映射表：首字母（大写） -> 完整组名（取第一个匹配） +# 相对丰度标准化 +data <- apply(data, 2, function(x) x / sum(x)) +data <- as.data.frame(data) + +# ========================================= +# 3. 丰度排序 +# ========================================= + +# 按总丰度降序排列 +mean_sort <- data[order(-rowSums(data)), ] +mean_sort <- as.data.frame(mean_sort) + +# 样本按第一物种丰度排序 +mean_sort2 <- t(mean_sort) +mean_sort2 <- mean_sort2[order(-mean_sort2[,1]), ] +mean_sort3 <- t(mean_sort2) +mean_sort3 <- as.data.frame(mean_sort3) + +# ========================================= +# 4. Top N物种处理 +# ========================================= + +# Top N + others +other <- colSums(mean_sort3[opts$legend:nrow(mean_sort3), ]) + +mean_sort3 <- mean_sort3[(opts$legend - 1):1, ] +mean_sort3 <- rbind(other, mean_sort3) + +rownames(mean_sort3)[1] <- "others" +mean_sort3 <- as.data.frame(mean_sort3) + +# ========================================= +# 5. 转长表 +# ========================================= + +# 添加 taxonomy 名称 +mean_sort3$tax <- rownames(mean_sort3) + +data_all <- reshape2::melt(mean_sort3, id.vars = "tax") +data_all <- as.data.frame(data_all) + +# group 提取 +data_all$group <- as.character(data_all$variable) +data_all$group <- gsub("[0-9]", "", data_all$group) + +# ========================================= +# 6. Metadata分组映射 +# ========================================= + +# 从 metadata 建立映射表：首字母 -> 完整组名 map <- metadata %>% mutate(prefix = toupper(substr(Group, 1, 1))) %>% - # 若同一首字母有多个完整组，默认保留第一个；如需其它策略可改这里 group_by(prefix) %>% slice(1) %>% ungroup() %>% select(prefix, full_group = Group) -# 将 data_all 的 group 首字母转换并左连接映射表，然后用完整名替换（若找不到则保留原值） +# 替换 group 名称 data_all <- data_all %>% mutate(prefix = toupper(substr(group, 1, 1))) %>% left_join(map, by = "prefix") %>% mutate(group = if_else(!is.na(full_group), full_group, group)) %>% select(-prefix, -full_group) -# 给分组排序 -# Sort for different groups -levels(as.factor(data_all$group)) - +# 保持 group 顺序 data_all2 <- data_all %>% mutate(group = factor(group, levels = unique(group))) +# ========================================= +# 7. 分面作图 +# ========================================= -# 根据样本数量确定每个分面的宽度，图例在顶部 -# Determine the width of each facet based on the number of samples, the legend is at the top plots <- lapply(split(data_all2, data_all2$group), function(df) { + group_name <- unique(df$group) - ggplot(df, aes(x = factor(variable, levels = unique(df$variable)), - y = value, fill = factor(tax, levels = unique(df$tax)))) + + + ggplot( + df, + aes( + x = factor(variable, levels = unique(df$variable)), + y = value, + fill = factor(tax, levels = unique(df$tax)) + ) + ) + geom_bar(stat = "identity", position = "stack", width = 1) + scale_y_continuous(labels = scales::percent, expand = c(0, 0)) + theme_classic() + labs(x = group_name, y = NULL) + - scale_fill_manual(values = c("#e5acd7", "#00ceff", "#ff630d", "#35978b","#d2da93", - "#5196d5", "#77aecd", "#ec8181", "#dfc6a5", "#e50719", - "#d27e43", "#8a4984", "#fe5094", "#8d342e", "#f94e54", - "#ffad00", "#36999d", "#00fc8d", "#b64aa0", "#9b82e1")) + - guides(fill = guide_legend(title = "Phylum")) + scale_fill_manual( + values = c( + "#e5acd7", "#00ceff", "#ff630d", "#35978b", "#d2da93", + "#5196d5", "#77aecd", "#ec8181", "#dfc6a5", "#e50719", + "#d27e43", "#8a4984", "#fe5094", "#8d342e", "#f94e54", + "#ffad00", "#36999d", "#00fc8d", "#b64aa0", "#9b82e1" + ), + ) }) -# 移除后三个分面的所有 y 轴元素和图例 -# Remove all y-axis elements and legends for the last three facets -for (i in 2:2) { - plots[[i]] <- plots[[i]] + theme(axis.text.y = element_blank(), - axis.text.x = element_blank(), - axis.ticks.y = element_blank(), - axis.ticks.x = element_blank(), - axis.title.y = element_blank(), - axis.line.y = element_blank(), - legend.position = "none") +# ========================================= +# 8. 去掉后续分面的坐标轴（不处理legend） +# ========================================= + +if (length(plots) > 1) { + for (i in 2:length(plots)) { + plots[[i]] <- plots[[i]] + + theme( + axis.text.y = element_blank(), + axis.text.x = element_blank(), + axis.ticks.y = element_blank(), + axis.ticks.x = element_blank(), + axis.title.y = element_blank(), + axis.line.y = element_blank(), + legend.position = "none" + ) + } } -# 为第一个分面保留 y 轴标签和图例 -# Keep y-axis label and legend for the first facet -plots[[1]] <- plots[[1]] + - ylab("Percentage (%)") + - theme(axis.text.x = element_blank(), - axis.ticks.x = element_blank(), - legend.position = "left", - legend.justification = c("left", "top")) +# 第一个分面也去掉 legend +plots[[1]] <- plots[[1]] + + ylab("Percentage (%)") + + theme( + axis.text.x = element_blank(), + axis.ticks.x = element_blank(), + legend.position = "none" + ) + +# ========================================= +# 9. patchwork组合（legend自动宽度，不重叠） +# ========================================= -# 每个分面的宽度由样本数量决定 -# The width of each facet is determined by the number of samples sample_counts <- table(data_all2$group) -relative_widths <- sample_counts / sum(sample_counts) - -# 使用 patchwork 组合图形，设置每个分面的宽度并统一图例 -# Use patchwork to combine graphics, set the width of each facet and unify the legend -p <- wrap_plots(plots) + - plot_layout(widths = relative_widths, guides = "collect") & - theme(legend.position = "top", - legend.justification = "center", - legend.direction = "horizontal", - legend.key.size = unit(0.3, "cm"), - legend.text = element_text(size = 8), - legend.spacing.x = unit(0.1, "cm"), - axis.title.y = element_text(size = 10), - plot.margin = unit(c(0.05, 0.05, 0.05, 0.05), "cm")) +relative_widths <- as.numeric(sample_counts) / sum(sample_counts) + +# 主图（无legend） +main_plot <- wrap_plots( + plots, + widths = relative_widths +) + +# 先把 tax 转 factor +data_all2$tax <- factor(data_all2$tax, levels = unique(data_all2$tax)) + +# 单独legend图 +legend_plot <- ggplot( + data_all2, + aes( + x = variable, + y = value, + fill = factor(tax, levels = unique(data_all2$tax)) + ) +) + + geom_bar(stat = "identity") + + scale_fill_manual( + values = c( + "#e5acd7", "#00ceff", "#ff630d", "#35978b", "#d2da93", + "#5196d5", "#77aecd", "#ec8181", "#dfc6a5", "#e50719", + "#d27e43", "#8a4984", "#fe5094", "#8d342e", "#f94e54", + "#ffad00", "#36999d", "#00fc8d", "#b64aa0", "#9b82e1" + ), + name = strsplit(basename(opts$output), "[._]")[[1]][1] + ) + + theme_classic()+ + theme( + legend.position = "right", + legend.box.margin = margin(0, 10, 0, 10), + # 图例标题字体 + legend.title = element_text(size = 8), + # 图例标签字体 + legend.text = element_text(size = 6), + # 图例方块大小 + legend.key.size = unit(0.4, "cm"), + # 图例行间距 + legend.spacing.y = unit(0.1, "cm") + ) + + +# 提取legend +legend <- cowplot::get_legend(legend_plot) + +# 计算legend真实宽度（单位：inch） +legend_width <- sum(convertWidth(legend$widths, "in", valueOnly = TRUE)) + +# 主图给固定宽度，legend按真实宽度占比 +main_width <- 4 # 主图基准宽度（可调整） + +p <- cowplot::plot_grid( + main_plot, + legend, + nrow = 1, + rel_widths = c(main_width, legend_width), + align = "h" +) # 保存 Saving # 大家可以修改图片名称和位置，长宽单位为毫米 From f857dccd2b7b0e40c33bdcd24dc7c81dd0586556 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:32:10 +0800 Subject: [PATCH 13/19] Add files via upload --- script/metaphlan4_boxplot_compare.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/script/metaphlan4_boxplot_compare.R b/script/metaphlan4_boxplot_compare.R index 36e05c2..98e0880 100644 --- a/script/metaphlan4_boxplot_compare.R +++ b/script/metaphlan4_boxplot_compare.R @@ -167,7 +167,10 @@ abu = taxonomy[,9:dim(taxonomy)[2]] merge = cbind(abu, grp) # group_by传变量，前面加".dots=" -mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) +# mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) +mergeTax = merge %>% + group_by(across(all_of(opts$taxonomy))) %>% + summarise(across(everything(), sum), .groups = "drop") # 合并后表格转换为数据框 mergeTax = as.data.frame(mergeTax) From 6dda701c9af737b2179fe745a6168cd0094556da Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:32:40 +0800 Subject: [PATCH 14/19] Add files via upload --- script/metaphlan_boxplot.R | 309 +++++++++++++++++++------------------ 1 file changed, 156 insertions(+), 153 deletions(-) diff --git a/script/metaphlan_boxplot.R b/script/metaphlan_boxplot.R index 07d7a77..e0ce9a8 100644 --- a/script/metaphlan_boxplot.R +++ b/script/metaphlan_boxplot.R @@ -1,153 +1,156 @@ -#!/usr/bin/env Rscript -# -# Copyright 2016-2020 Yong-Xin Liu - -# If used this script, please cited: -# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo & Yang Bai. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein Cell 41, 1-16, doi:10.1007/s13238-020-00724-8 (2020). -# Jingying Zhang, Yong-Xin Liu, et. al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology 37, 676-684, doi:10.1038/s41587-019-0104-4 (2019). - -# 1. 分析前准备：帮助、参数、依赖包和读取文件 - -# 命令行运行为当前目录；Rstudio手动运行脚本前需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 - - -# 1.1 程序功能描述和主要步骤 - -# 程序功能：绘制物种组成热图 -# 主要步骤: -# - 读取Metaphlan2物种组表 result/metaphlan2/taxonomy.spf(己标准化为100) -# - 物种组成按丰度均值排序 -# - 筛选Top N个物种绘制热图 - -# # 程序使用示例USAGE -# # 显示脚本帮助 help -# Rscript db/script/metaphlan_hclust_boxplot.R -h -# # 默认读取result/metaphlan2/taxonomy.spf，按指定列合并、排序并取Top25种绘制热图，输出至输入目录 -# Rscript db/script/metaphlan_hclust_boxplot.R -# # 完整参数：-i输入Metaphlan2文件； -# # -t 分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain，界门纲目科属种株，推荐门，目，属 -# # -n 输出物种数量，默认为25，最大为合并后的数量 -# # -o输出图表前缀，默认根据输入文件、物种级别和数量自动生成； -# Rscript db/script/metaphlan_hclust_boxplot.R \ -# -i result/metaphlan2/taxonomy.spf \ -# -t Species \ -# -n 25 \ -# -o result/metaphlan2/boxplot_Species - - -# 1.2 解析命令行 -# 设置清华源加速下载 -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析是否安装，安装并加载 -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} - - - -# 解析参数-h显示帮助信息 -# 此版本参数在windows中报错显示3，无法正常运行 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="result/kraken2/taxonomy_count.norm.spf", help="Metaphlan2 [default %default]"), - make_option(c("-t", "--taxonomy"), type="character", default="Species", help="Taxonomy level [default %default]"), - make_option(c("-n", "--TopN"), type="numeric", default="25", help="Number of taxonomy showing [default %default]"), - make_option(c("-w", "--width"), type="numeric", default=8, - help="Width of figure [default %default]"), - make_option(c("-e", "--height"), type="numeric", default=5, - help="Height of figure [default %default]"), - make_option(c("-o", "--output"), type="character", default="", help="Output boxplot filename [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) - - # 调置如果无调设置输出，根据其它参数设置默认输出 - prefix = gsub("taxonomy.spf$", "", opts$input, perl = T) - if (opts$output==""){opts$output=paste0(prefix, "boxplot", opts$taxonomy, opts$TopN)} - - # 显示输入输出确认是否正确 - # Metaphlan2物种组成表 - print(paste("The input file: ", opts$input, sep = "")) - # 绘制的分类级别, 默认为种 - print(paste("Taxonomy level: ", opts$taxonomy, ". Default if Species", sep = "")) - # 选择绘制高丰度物种数量，默认30，最大为物种级别非冗余条目数量 - print(paste("Number of taxonomy showing: ", opts$TopN, sep = "")) - # 输出文件名，不填则为输入目录+boxplot+taxonomy - print(paste("Output boxplot filename: ", opts$output, sep = "")) -} - - -# 1.3 安装CRAN来源常用包 -# 依赖包列表：参数解析、数据变换、绘图和开发包安装、安装依赖、ggplot主题 -package_list = c("ggplot2","dplyr","reshape2") -# 判断R包加载是否成功来决定是否安装后再加载 -for(p in package_list){ - if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){ - install.packages(p, repos=site) - suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))) - } -} - - -# 1.4 读取输入文件 - -# 读取metaphlan2文件 -# 默认的quote会跳过2/3的数据，导致行减少产生NA，改默认值为空 -taxonomy = read.table(opts$input, header=T, sep="\t", quote = "", row.names=NULL, comment.char="") -print(paste0("All taxonomy annotations are ", dim(taxonomy)[1], " lines!")) -# 去除NA，否则无法计算 -taxonomy = na.omit(taxonomy) -# 显示样本总数据，有冗余 -# colSums(taxonomy) - - - -# 2. 计算过程 - -## 2.1 按指定组合并 - -grp = taxonomy[, opts$taxonomy, drop=F] -abu = taxonomy[,9:dim(taxonomy)[2]] -merge = cbind(abu, grp) -# group_by传变量，前面加".dots=" -mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) -# 合并后表格转换为数据框 -mergeTax = as.data.frame(mergeTax) -# 按丰度排序 -idx = order(rowMeans(mergeTax[,2:dim(mergeTax)[2]]), decreasing = T) -mergeTax = mergeTax[idx,] -# 添加行名 -rownames(mergeTax)=mergeTax[,1] - -## 2.1 筛选TopN绘图 - -# remove rownames line -Top = mergeTax[,-1] -# normalization to percentage 100 -Top = as.data.frame(t(t(Top)/colSums(Top,na=T) * 100)) -# Select TopN line for plotting -Top = head(Top, n=opts$TopN) -Top = cbind(row.names(Top), Top) -colnames(Top)[1] = opts$taxonomy - -# melt table -data_all = as.data.frame(melt(Top, id.vars = opts$taxonomy)) -# set taxonomy order by abundance, default by alphabet -data_all[[opts$taxonomy]] = factor(as.character(data_all[[opts$taxonomy]]), levels=rownames(Top)) - -# 3. 结果输出 - -# 保存表格 -write.table(mergeTax, file=paste(opts$output, ".txt", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=T) - - -# 保存图片 -p = ggplot(data_all, aes(x=data_all[[opts$taxonomy]], y = value, fill = data_all[[opts$taxonomy]] )) + - geom_boxplot()+ - theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())+ - xlab(opts$taxonomy)+ylab("Percentage (%)")+ theme_classic()+ - theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1),legend.position="none") -p = p + geom_jitter(position=position_jitter(0.17), size=1, alpha=0.7) -p -ggsave(paste(opts$output, ".pdf", sep=""), p, width = opts$width, height = opts$height) - +#!/usr/bin/env Rscript +# +# Copyright 2016-2020 Yong-Xin Liu + +# If used this script, please cited: +# Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo & Yang Bai. A practical guide to amplicon and metagenomic analysis of microbiome data. Protein Cell 41, 1-16, doi:10.1007/s13238-020-00724-8 (2020). +# Jingying Zhang, Yong-Xin Liu, et. al. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology 37, 676-684, doi:10.1038/s41587-019-0104-4 (2019). + +# 1. 分析前准备：帮助、参数、依赖包和读取文件 + +# 命令行运行为当前目录；Rstudio手动运行脚本前需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 + + +# 1.1 程序功能描述和主要步骤 + +# 程序功能：绘制物种组成热图 +# 主要步骤: +# - 读取Metaphlan2物种组表 result/metaphlan2/taxonomy.spf(己标准化为100) +# - 物种组成按丰度均值排序 +# - 筛选Top N个物种绘制热图 + +# # 程序使用示例USAGE +# # 显示脚本帮助 help +# Rscript db/script/metaphlan_hclust_boxplot.R -h +# # 默认读取result/metaphlan2/taxonomy.spf，按指定列合并、排序并取Top25种绘制热图，输出至输入目录 +# Rscript db/script/metaphlan_hclust_boxplot.R +# # 完整参数：-i输入Metaphlan2文件； +# # -t 分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain，界门纲目科属种株，推荐门，目，属 +# # -n 输出物种数量，默认为25，最大为合并后的数量 +# # -o输出图表前缀，默认根据输入文件、物种级别和数量自动生成； +# Rscript db/script/metaphlan_hclust_boxplot.R \ +# -i result/metaphlan2/taxonomy.spf \ +# -t Species \ +# -n 25 \ +# -o result/metaphlan2/boxplot_Species + + +# 1.2 解析命令行 +# 设置清华源加速下载 +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析是否安装，安装并加载 +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} + + + +# 解析参数-h显示帮助信息 +# 此版本参数在windows中报错显示3，无法正常运行 +if (TRUE){ + option_list = list( + make_option(c("-i", "--input"), type="character", default="result/kraken2/taxonomy_count.norm.spf", help="Metaphlan2 [default %default]"), + make_option(c("-t", "--taxonomy"), type="character", default="Species", help="Taxonomy level [default %default]"), + make_option(c("-n", "--TopN"), type="numeric", default="25", help="Number of taxonomy showing [default %default]"), + make_option(c("-w", "--width"), type="numeric", default=8, + help="Width of figure [default %default]"), + make_option(c("-e", "--height"), type="numeric", default=5, + help="Height of figure [default %default]"), + make_option(c("-o", "--output"), type="character", default="", help="Output boxplot filename [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) + + # 调置如果无调设置输出，根据其它参数设置默认输出 + prefix = gsub("taxonomy.spf$", "", opts$input, perl = T) + if (opts$output==""){opts$output=paste0(prefix, "boxplot", opts$taxonomy, opts$TopN)} + + # 显示输入输出确认是否正确 + # Metaphlan2物种组成表 + print(paste("The input file: ", opts$input, sep = "")) + # 绘制的分类级别, 默认为种 + print(paste("Taxonomy level: ", opts$taxonomy, ". Default if Species", sep = "")) + # 选择绘制高丰度物种数量，默认30，最大为物种级别非冗余条目数量 + print(paste("Number of taxonomy showing: ", opts$TopN, sep = "")) + # 输出文件名，不填则为输入目录+boxplot+taxonomy + print(paste("Output boxplot filename: ", opts$output, sep = "")) +} + + +# 1.3 安装CRAN来源常用包 +# 依赖包列表：参数解析、数据变换、绘图和开发包安装、安装依赖、ggplot主题 +package_list = c("ggplot2","dplyr","reshape2") +# 判断R包加载是否成功来决定是否安装后再加载 +for(p in package_list){ + if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){ + install.packages(p, repos=site) + suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))) + } +} + + +# 1.4 读取输入文件 + +# 读取metaphlan2文件 +# 默认的quote会跳过2/3的数据，导致行减少产生NA，改默认值为空 +taxonomy = read.table(opts$input, header=T, sep="\t", quote = "", row.names=NULL, comment.char="") +print(paste0("All taxonomy annotations are ", dim(taxonomy)[1], " lines!")) +# 去除NA，否则无法计算 +taxonomy = na.omit(taxonomy) +# 显示样本总数据，有冗余 +# colSums(taxonomy) + + + +# 2. 计算过程 + +## 2.1 按指定组合并 + +grp = taxonomy[, opts$taxonomy, drop=F] +abu = taxonomy[,9:dim(taxonomy)[2]] +merge = cbind(abu, grp) +# group_by传变量，前面加".dots=" +# mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) +mergeTax = merge %>% + group_by(across(all_of(opts$taxonomy))) %>% + summarise(across(everything(), sum), .groups = "drop") +# 合并后表格转换为数据框 +mergeTax = as.data.frame(mergeTax) +# 按丰度排序 +idx = order(rowMeans(mergeTax[,2:dim(mergeTax)[2]]), decreasing = T) +mergeTax = mergeTax[idx,] +# 添加行名 +rownames(mergeTax)=mergeTax[,1] + +## 2.1 筛选TopN绘图 + +# remove rownames line +Top = mergeTax[,-1] +# normalization to percentage 100 +Top = as.data.frame(t(t(Top)/colSums(Top,na=T) * 100)) +# Select TopN line for plotting +Top = head(Top, n=opts$TopN) +Top = cbind(row.names(Top), Top) +colnames(Top)[1] = opts$taxonomy + +# melt table +data_all = as.data.frame(melt(Top, id.vars = opts$taxonomy)) +# set taxonomy order by abundance, default by alphabet +data_all[[opts$taxonomy]] = factor(as.character(data_all[[opts$taxonomy]]), levels=rownames(Top)) + +# 3. 结果输出 + +# 保存表格 +write.table(mergeTax, file=paste(opts$output, ".txt", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=T) + + +# 保存图片 +p = ggplot(data_all, aes(x=data_all[[opts$taxonomy]], y = value, fill = data_all[[opts$taxonomy]] )) + + geom_boxplot()+ + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())+ + xlab(opts$taxonomy)+ylab("Percentage (%)")+ theme_classic()+ + theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1),legend.position="none") +p = p + geom_jitter(position=position_jitter(0.17), size=1, alpha=0.7) +p +ggsave(paste(opts$output, ".pdf", sep=""), p, width = opts$width, height = opts$height) + From e4f8dcb773d6772adb709e4307e06a6d92b5c248 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:37:30 +0800 Subject: [PATCH 15/19] Add files via upload From 6feacd3da9d42d6dacce2661b26836a1a618724a Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:38:20 +0800 Subject: [PATCH 16/19] Add files via upload --- script/metaphlan_hclust_heatmap.R | 313 +++++++++++++++--------------- 1 file changed, 158 insertions(+), 155 deletions(-) diff --git a/script/metaphlan_hclust_heatmap.R b/script/metaphlan_hclust_heatmap.R index 1956807..2bec90f 100644 --- a/script/metaphlan_hclust_heatmap.R +++ b/script/metaphlan_hclust_heatmap.R @@ -1,156 +1,159 @@ -#!/usr/bin/env Rscript -# -# Copyright 2016-2026 Yong-Xin Liu - -# If used this script, please cited: -# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 - -# 1. 分析前准备：帮助、参数、依赖包和读取文件 - -# 命令行运行为当前目录；Rstudio手动运行脚本前需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 - - -# 1.1 程序功能描述和主要步骤 - -# 程序功能：绘制物种组成热图 -# 主要步骤: -# - 读取Metaphlan2物种组表 result/metaphlan2/taxonomy.spf(己标准化为100) -# - 物种组成按丰度均值排序 -# - 筛选Top N个物种绘制热图 - -# # 程序使用示例USAGE -# # 显示脚本帮助 help -# Rscript db/script/metaphlan_hclust_heatmap.R -h -# # 默认读取result/metaphlan2/taxonomy.spf，按指定列合并、排序并取Top25种绘制热图，输出至输入目录 -# # 完整参数：-i输入Metaphlan2文件； -# # -t 分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain，界门纲目科属种株，推荐门，目，属 -# # -n 输出物种数量，默认为25，最大为合并后的数量 -# # -o输出图表前缀，默认根据输入文件、物种级别和数量自动生成； -# Rscript db/script/metaphlan_hclust_heatmap.R \ -# -i result/metaphlan2/taxonomy.spf \ -# -t Species \ -# -n 25 \ -# -o result/metaphlan2/heatmap_Species - - -# 1.2 解析命令行 -# 设置清华源加速下载 -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -# 判断命令行解析是否安装，安装并加载 -if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { - install.packages("optparse", repos=site) - require("optparse",character.only=T) -} - - -# 解析参数-h显示帮助信息 -# 此版本参数在windows中报错显示3，无法正常运行 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="metaphlan4/taxonomy.spf", help="Metaphlan2 [default %default]"), - make_option(c("-t", "--taxonomy"), type="character", default="Genus", help="Taxonomy level [default %default]"), - make_option(c("-n", "--TopN"), type="numeric", default="25", help="Number of taxonomy showing [default %default]"), - make_option(c("-o", "--output"), type="character", default="", help="Output heatmap filename [default %default]"), - make_option(c("-w", "--width"), type="numeric", default=183, - help="Figure width in mm [default %default]"), - make_option(c("-e", "--height"), type="numeric", default=118, - help="Figure heigth in mm [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) - - # 调置如果无调设置输出，根据其它参数设置默认输出 - prefix = gsub("taxonomy.spf$", "", opts$input, perl = T) - if (opts$output==""){opts$output=paste0(prefix, "Heatmap", opts$taxonomy, opts$TopN)} - - # 显示输入输出确认是否正确 - # Metaphlan2物种组成表 - print(paste("The input file: ", opts$input, sep = "")) - # 绘制的分类级别, 默认为种 - print(paste("Taxonomy level: ", opts$taxonomy, ". Default if Species", sep = "")) - # 选择绘制高丰度物种数量，默认30，最大为物种级别非冗余条目数量 - print(paste("Number of taxonomy showing: ", opts$TopN, sep = "")) - # 输出文件名，不填则为输入目录+heatmap+taxonomy - print(paste("Output heatmap filename: ", opts$output, sep = "")) -} - - -# 1.3 安装CRAN来源常用包 -# 依赖包列表：参数解析、数据变换、绘图和开发包安装、安装依赖、ggplot主题 -package_list = c("pheatmap","ggplot2","dplyr","stringr") -# 判断R包加载是否成功来决定是否安装后再加载 -for(p in package_list){ - if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){ - install.packages(p, repos=site) - suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))) - } -} - - -# 1.4 读取输入文件 - -# 读取metaphlan2文件 -# 默认的quote会跳过2/3的数据，导致行减少产生NA，改默认值为空 -taxonomy = read.table(opts$input, header=T, sep="\t", quote = "", row.names=NULL, comment.char="") -print(paste0("All taxonomy annotations are ", dim(taxonomy)[1], " lines!")) -# 去除NA，否则无法计算 -taxonomy = na.omit(taxonomy) -# 显示样本总数据，有冗余 -# colSums(taxonomy) - - - -# 2. 计算过程 - -## 2.1 按指定组合并 - -grp = taxonomy[, opts$taxonomy, drop=F] - -# 判断metaphlan 8列 或kraken 7列，从第9或8列取 -n=9 -if (str_detect(opts$input,"kraken")) { - n=8 -} -abu = taxonomy[,n:dim(taxonomy)[2]] - -merge = cbind(abu, grp) -# group_by传变量，前面加".dots=" -mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) -# 合并后表格转换为数据框 -mergeTax = as.data.frame(mergeTax) -# 按丰度排序 -idx = order(rowMeans(mergeTax[,2:dim(mergeTax)[2]]), decreasing = T) -mergeTax = mergeTax[idx,] -# 添加行名 -rownames(mergeTax)=mergeTax[,1] - -## 2.1 筛选TopN绘图 - -# remove rownames line -Top = mergeTax[,-1] -# normalization to percentage 100 -Top = as.data.frame(t(t(Top)/colSums(Top,na=T) * 100)) -# Select TopN line for plotting -Top = head(Top, n=opts$TopN) - -pheatmap(Top) - -# 3. 结果输出 - -# 保存表格 -write.table(mergeTax, file=paste(opts$output, ".txt", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=T) - - -# 保存图片 -pheatmap(Top, - filename = paste0(opts$output, ".pdf"),fontsize=7, - width = opts$width/25.4, height = opts$height/25.4, - main = paste("Top", opts$TopN, opts$taxonomy, sep=" ")) - - -# width=dim(Top)[2], height=dim(Top)[1]/4, -# scale = "row", -# cutree_rows=2,cutree_cols = 2, -# annotation_col = anno_col, annotation_row = anno_row, -# annotation_names_row= T,annotation_names_col=T, -# show_rownames=F,show_colnames=T, +#!/usr/bin/env Rscript +# +# Copyright 2016-2026 Yong-Xin Liu + +# If used this script, please cited: +# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 + +# 1. 分析前准备：帮助、参数、依赖包和读取文件 + +# 命令行运行为当前目录；Rstudio手动运行脚本前需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 + + +# 1.1 程序功能描述和主要步骤 + +# 程序功能：绘制物种组成热图 +# 主要步骤: +# - 读取Metaphlan2物种组表 result/metaphlan2/taxonomy.spf(己标准化为100) +# - 物种组成按丰度均值排序 +# - 筛选Top N个物种绘制热图 + +# # 程序使用示例USAGE +# # 显示脚本帮助 help +# Rscript db/script/metaphlan_hclust_heatmap.R -h +# # 默认读取result/metaphlan2/taxonomy.spf，按指定列合并、排序并取Top25种绘制热图，输出至输入目录 +# # 完整参数：-i输入Metaphlan2文件； +# # -t 分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain，界门纲目科属种株，推荐门，目，属 +# # -n 输出物种数量，默认为25，最大为合并后的数量 +# # -o输出图表前缀，默认根据输入文件、物种级别和数量自动生成； +# Rscript db/script/metaphlan_hclust_heatmap.R \ +# -i result/metaphlan2/taxonomy.spf \ +# -t Species \ +# -n 25 \ +# -o result/metaphlan2/heatmap_Species + + +# 1.2 解析命令行 +# 设置清华源加速下载 +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +# 判断命令行解析是否安装，安装并加载 +if (!suppressWarnings(suppressMessages(require("optparse", character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))) { + install.packages("optparse", repos=site) + require("optparse",character.only=T) +} + + +# 解析参数-h显示帮助信息 +# 此版本参数在windows中报错显示3，无法正常运行 +if (TRUE){ + option_list = list( + make_option(c("-i", "--input"), type="character", default="metaphlan4/taxonomy.spf", help="Metaphlan2 [default %default]"), + make_option(c("-t", "--taxonomy"), type="character", default="Genus", help="Taxonomy level [default %default]"), + make_option(c("-n", "--TopN"), type="numeric", default="25", help="Number of taxonomy showing [default %default]"), + make_option(c("-o", "--output"), type="character", default="", help="Output heatmap filename [default %default]"), + make_option(c("-w", "--width"), type="numeric", default=183, + help="Figure width in mm [default %default]"), + make_option(c("-e", "--height"), type="numeric", default=118, + help="Figure heigth in mm [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) + + # 调置如果无调设置输出，根据其它参数设置默认输出 + prefix = gsub("taxonomy.spf$", "", opts$input, perl = T) + if (opts$output==""){opts$output=paste0(prefix, "Heatmap", opts$taxonomy, opts$TopN)} + + # 显示输入输出确认是否正确 + # Metaphlan2物种组成表 + print(paste("The input file: ", opts$input, sep = "")) + # 绘制的分类级别, 默认为种 + print(paste("Taxonomy level: ", opts$taxonomy, ". Default if Species", sep = "")) + # 选择绘制高丰度物种数量，默认30，最大为物种级别非冗余条目数量 + print(paste("Number of taxonomy showing: ", opts$TopN, sep = "")) + # 输出文件名，不填则为输入目录+heatmap+taxonomy + print(paste("Output heatmap filename: ", opts$output, sep = "")) +} + + +# 1.3 安装CRAN来源常用包 +# 依赖包列表：参数解析、数据变换、绘图和开发包安装、安装依赖、ggplot主题 +package_list = c("pheatmap","ggplot2","dplyr","stringr") +# 判断R包加载是否成功来决定是否安装后再加载 +for(p in package_list){ + if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){ + install.packages(p, repos=site) + suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))) + } +} + + +# 1.4 读取输入文件 + +# 读取metaphlan2文件 +# 默认的quote会跳过2/3的数据，导致行减少产生NA，改默认值为空 +taxonomy = read.table(opts$input, header=T, sep="\t", quote = "", row.names=NULL, comment.char="") +print(paste0("All taxonomy annotations are ", dim(taxonomy)[1], " lines!")) +# 去除NA，否则无法计算 +taxonomy = na.omit(taxonomy) +# 显示样本总数据，有冗余 +# colSums(taxonomy) + + + +# 2. 计算过程 + +## 2.1 按指定组合并 + +grp = taxonomy[, opts$taxonomy, drop=F] + +# 判断metaphlan 8列 或kraken 7列，从第9或8列取 +n=9 +if (str_detect(opts$input,"kraken")) { + n=8 +} +abu = taxonomy[,n:dim(taxonomy)[2]] + +merge = cbind(abu, grp) +# group_by传变量，前面加".dots=" +# mergeTax = merge %>% group_by(.dots=opts$taxonomy) %>% summarise_all(sum) +mergeTax = merge %>% + group_by(across(all_of(opts$taxonomy))) %>% + summarise(across(everything(), sum), .groups = "drop") +# 合并后表格转换为数据框 +mergeTax = as.data.frame(mergeTax) +# 按丰度排序 +idx = order(rowMeans(mergeTax[,2:dim(mergeTax)[2]]), decreasing = T) +mergeTax = mergeTax[idx,] +# 添加行名 +rownames(mergeTax)=mergeTax[,1] + +## 2.1 筛选TopN绘图 + +# remove rownames line +Top = mergeTax[,-1] +# normalization to percentage 100 +Top = as.data.frame(t(t(Top)/colSums(Top,na=T) * 100)) +# Select TopN line for plotting +Top = head(Top, n=opts$TopN) + +pheatmap(Top) + +# 3. 结果输出 + +# 保存表格 +write.table(mergeTax, file=paste(opts$output, ".txt", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=T) + + +# 保存图片 +pheatmap(Top, + filename = paste0(opts$output, ".pdf"),fontsize=7, + width = opts$width/25.4, height = opts$height/25.4, + main = paste("Top", opts$TopN, opts$taxonomy, sep=" ")) + + +# width=dim(Top)[2], height=dim(Top)[1]/4, +# scale = "row", +# cutree_rows=2,cutree_cols = 2, +# annotation_col = anno_col, annotation_row = anno_row, +# annotation_names_row= T,annotation_names_col=T, +# show_rownames=F,show_colnames=T, # display_numbers=F \ No newline at end of file From faca6b0ec894e4f106ba7c74b265dace23d147b0 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:38:52 +0800 Subject: [PATCH 17/19] Add files via upload --- script/alpha_boxplot_new.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/script/alpha_boxplot_new.R b/script/alpha_boxplot_new.R index aa2dc41..f8c6eeb 100644 --- a/script/alpha_boxplot_new.R +++ b/script/alpha_boxplot_new.R @@ -143,7 +143,8 @@ alpha_boxplot2 <- function(alpha_div, metadata, index = "shannon", groupID = "gr max = max(df[, c(index)]) min = min(df[, index]) x = df[, c("group", index)] - y = x %>% group_by(group) %>% summarise_(Max = paste("max(", index, ")", sep = "")) + # y = x %>% group_by(group) %>% summarise_(Max = paste("max(", index, ")", sep = "")) + y = x %>% group_by(group) %>% summarise(Max = max(.data[[index]], na.rm = TRUE)) y = as.data.frame(y) rownames(y) = y$group df$y = y[as.character(df$group), ]$Max + (max - min) * 0.05 From ce540f3e6f9943f7c1377fb216c74da5a5d840aa Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 15:39:20 +0800 Subject: [PATCH 18/19] Add files via upload --- script/alpha_boxplot.R | 402 +++++++++++++++++++++++++++-------------- 1 file changed, 263 insertions(+), 139 deletions(-) diff --git a/script/alpha_boxplot.R b/script/alpha_boxplot.R index 4992c64..dc383e0 100644 --- a/script/alpha_boxplot.R +++ b/script/alpha_boxplot.R @@ -1,139 +1,263 @@ -#!/usr/bin/env Rscript - -# Copyright 2016-2026 Yong-Xin Liu - -# If used this script, please cited: -# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 - -# 手动运行脚本请，需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 - -# 更新 -# 2021/5/31: 更新引文，添加数据表转置和标准化的选项 -# 2025/11/25: 更新引文，更新模板 - -# 1.1 程序功能描述和主要步骤 - -# 程序功能：Alpha多样性箱线图+统计 -# Functions: Alpha boxplot - -options(warn = -1) # Turn off warning - - -## 设置输入输出文件和参数 - -# 修改下面`default=`后面的文件和参数。 -# -# 输入文件为原始alpha多样性(vegan.txt)+分组信息(metadata.txt) -# -# 输入文件"-i", "--input"，result/alpha/vegan.txt; alpha多样性表格 -# -# 实验设计"-d", "--design"，默认`metadata.txt`，可手动修改文件位置； -# -# 分组列名"-n", "--group"，默认将metadata.txt中的Group列作为分组信息，可修改为任意列名； -# -# 分组列名"-o", "--output"，默认为输出目录，图片文件名为alpha_boxplot_+多样性指数名+.pdf；统计文本位于代码运行目录中alpha_boxplot_TukeyHSD.txt； -# -# 图片宽"-w", "--width"，默认89 mm，根据图像布局可适当增大或缩小 -# -# 图片高"-e", "--height"，默认59 mm，根据图像布局可适当增大或缩小 - - -# 1.2 依赖包安装 - -site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" -a = rownames(installed.packages()) - -# install CRAN -install_CRAN <- c("ggplot2", "BiocManager", "optparse") -for (i in install_CRAN) { - if (!i %in% a) - install.packages(i, repos = site) - require(i,character.only=T) - a = rownames(installed.packages()) -} - -# install bioconductor -install_bioc <- c("ggplot2", "multcompView") -for (i in install_bioc) { - if (!i %in% a) - BiocManager::install(i, update = F) - a = rownames(installed.packages()) -} - - - -install_bioc <- c("ggplot2", "multcompView") - -for (i in install_bioc) { - if (!i %in% a) - BiocManager::install(i, update = F, site_repository=site) - a = rownames(installed.packages()) -} - -# install github -if (!"amplicon" %in% a){ - devtools::install_github("microbiota/amplicon") -} - - -# 1.3 解析命令行 -# 解析参数-h显示帮助信息 -if (TRUE){ - option_list = list( - make_option(c("-i", "--input"), type="character", default="result/alpha/vegan.txt", - help="Alpha diversity matrix [default %default]"), - make_option(c("-a", "--alpha_index"), type="character", default="richness", - help="Group name [default %default]"), - make_option(c("-d", "--design"), type="character", default="result/metadata.txt", - help="Design file or metadata [default %default]"), - make_option(c("-t", "--transpose"), type="logical", default=FALSE, - help="Design file or metadata [default %default]"), - make_option(c("-n", "--group"), type="character", default="Group", - help="Group name [default %default]"), - make_option(c("-s", "--scale"), type="logical", default=FALSE, - help="Normalize to 100 [default %default]"), - make_option(c("-o", "--output"), type="character", default="result/alpha/", - help="Output pdf directory, with prefix alpha_boxplot_; Stat in alpha_boxplot_TukeyHSD.txt [default %default]"), - make_option(c("-x", "--xlabAngle"), type="logical", default=FALSE, - help="X lab set in angle [default %default]"), - make_option(c("-w", "--width"), type="numeric", default=89, - help="Figure width in mm [default %default]"), - make_option(c("-e", "--height"), type="numeric", default=59, - help="Figure heigth in mm [default %default]") - ) - opts = parse_args(OptionParser(option_list=option_list)) -} -suppressWarnings(dir.create(dirname(opts$output), showWarnings = F)) - - -# 2. 依赖关系检查、安装和加载 - -suppressWarnings(suppressMessages(library(amplicon))) - - -# 3. 读取输入文件 - -# 读取OTU表 -alpha_div = read.table(opts$input, header=T, row.names=1, sep="\t", comment.char="") - - -# 条件判断是否转置 -if (opts$transpose){ - alpha_div = as.data.frame(t(alpha_div)) -} - -# 条件判断是否标准化 -if (opts$scale){ - alpha_div = alpha_div/rowSums(alpha_div,na=T)*100 -} - -# 读取实验设计 -metadata = read.table(opts$design, header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F) - -p = alpha_boxplot(alpha_div, index = opts$alpha_index, metadata, groupID = opts$group) -if (opts$xlabAngle){ - p = p + theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1)) -} -# Saving figure -# 保存图片，大家可以修改图片名称和位置，长宽单位为毫米 -ggsave(paste0(opts$output,"boxplot_",opts$alpha_index,".pdf"), p, width = opts$width, height = opts$height, units = "mm") +#!/usr/bin/env Rscript + +# Copyright 2016-2026 Yong-Xin Liu + +# If used this script, please cited: +# Bai, et al. 2025. EasyMetagenome: A User‐Friendly and Flexible Pipeline for Shotgun Metagenomic Analysis in Microbiome Research. iMeta 4: e70001. https://doi.org/10.1002/imt2.70001 + +# 手动运行脚本请，需要设置工作目录，使用 Ctrl+Shift+H 或 Session - Set Work Directory - Choose Directory / To Source File Location 设置工作目录 + +# 更新 +# 2021/5/31: 更新引文，添加数据表转置和标准化的选项 +# 2025/11/25: 更新引文，更新模板 + +# 1.1 程序功能描述和主要步骤 + +# 程序功能：Alpha多样性箱线图+统计 +# Functions: Alpha boxplot + +options(warn = -1) # Turn off warning + + +## 设置输入输出文件和参数 + +# 修改下面`default=`后面的文件和参数。 +# +# 输入文件为原始alpha多样性(vegan.txt)+分组信息(metadata.txt) +# +# 输入文件"-i", "--input"，result/alpha/vegan.txt; alpha多样性表格 +# +# 实验设计"-d", "--design"，默认`metadata.txt`，可手动修改文件位置； +# +# 分组列名"-n", "--group"，默认将metadata.txt中的Group列作为分组信息，可修改为任意列名； +# +# 分组列名"-o", "--output"，默认为输出目录，图片文件名为alpha_boxplot_+多样性指数名+.pdf；统计文本位于代码运行目录中alpha_boxplot_TukeyHSD.txt； +# +# 图片宽"-w", "--width"，默认89 mm，根据图像布局可适当增大或缩小 +# +# 图片高"-e", "--height"，默认59 mm，根据图像布局可适当增大或缩小 + + +# 1.2 依赖包安装 + +site="https://mirrors.tuna.tsinghua.edu.cn/CRAN" +a = rownames(installed.packages()) + +# install CRAN +install_CRAN <- c("ggplot2", "BiocManager", "optparse") +for (i in install_CRAN) { + if (!i %in% a) + install.packages(i, repos = site) + require(i,character.only=T) + a = rownames(installed.packages()) +} + +# install bioconductor +install_bioc <- c("ggplot2", "multcompView") +for (i in install_bioc) { + if (!i %in% a) + BiocManager::install(i, update = F) + a = rownames(installed.packages()) +} + + + +install_bioc <- c("ggplot2", "multcompView") + +for (i in install_bioc) { + if (!i %in% a) + BiocManager::install(i, update = F, site_repository=site) + a = rownames(installed.packages()) +} + +# install github +if (!"amplicon" %in% a){ + devtools::install_github("microbiota/amplicon") +} + + +# 1.3 解析命令行 +# 解析参数-h显示帮助信息 +if (TRUE){ + option_list = list( + make_option(c("-i", "--input"), type="character", default="result/alpha/vegan.txt", + help="Alpha diversity matrix [default %default]"), + make_option(c("-a", "--alpha_index"), type="character", default="richness", + help="Group name [default %default]"), + make_option(c("-d", "--design"), type="character", default="result/metadata.txt", + help="Design file or metadata [default %default]"), + make_option(c("-t", "--transpose"), type="logical", default=FALSE, + help="Design file or metadata [default %default]"), + make_option(c("-n", "--group"), type="character", default="Group", + help="Group name [default %default]"), + make_option(c("-s", "--scale"), type="logical", default=FALSE, + help="Normalize to 100 [default %default]"), + make_option(c("-o", "--output"), type="character", default="result/alpha/", + help="Output pdf directory, with prefix alpha_boxplot_; Stat in alpha_boxplot_TukeyHSD.txt [default %default]"), + make_option(c("-x", "--xlabAngle"), type="logical", default=FALSE, + help="X lab set in angle [default %default]"), + make_option(c("-w", "--width"), type="numeric", default=89, + help="Figure width in mm [default %default]"), + make_option(c("-e", "--height"), type="numeric", default=59, + help="Figure heigth in mm [default %default]") + ) + opts = parse_args(OptionParser(option_list=option_list)) +} +suppressWarnings(dir.create(dirname(opts$output), showWarnings = F)) + + +# 2. 依赖关系检查、安装和加载 +#install.packages("remotes") +#remotes::install_github("microbiota/amplicon") +suppressWarnings(suppressMessages(library(amplicon))) + + +# 3. 读取输入文件 + +# 读取OTU表 +alpha_div = read.table(opts$input, header=T, row.names=1, sep="\t", comment.char="") + + +# 条件判断是否转置 +if (opts$transpose){ + alpha_div = as.data.frame(t(alpha_div)) +} + +# 条件判断是否标准化 +if (opts$scale){ + alpha_div = alpha_div/rowSums(alpha_div,na=T)*100 +} + +# 读取实验设计 +metadata = read.table(opts$design, header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F) + +alpha_boxplot <- function(alpha_div, metadata, index = "richness", groupID = "Group", + outlier = TRUE +) { + # 依赖关系检测与安装 + p_list = c("ggplot2", "dplyr", "multcompView") # "agricolae" + for(p in p_list){ + if (!requireNamespace(p)){ + install.packages(p)} + suppressPackageStartupMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)) + } + + # 测试默认参数 + # library(amplicon) + # index = "richness" + # groupID = "Group" + # metadata = read.table("metadata2.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F) + # metadata = subset(metadata, Group %in% c("KO","OE")) + + # 交叉筛选 + idx = rownames(metadata) %in% rownames(alpha_div) + metadata = metadata[idx,,drop=F] + alpha_div = alpha_div[rownames(metadata),] + + # 提取样品组信息,默认为group可指定 + sampFile = as.data.frame(metadata[, groupID],row.names = row.names(metadata)) + # colnames(sampFile)[1] = "group" + + # 合并alpha_div和metadata + df = cbind(alpha_div[rownames(sampFile),index], sampFile) + colnames(df) = c(index,"group") + + # 统计各种显著性 + model = aov(df[[index]] ~ group, data=df) + # 计算Tukey显著性差异检验 + Tukey_HSD = TukeyHSD(model, ordered = TRUE, conf.level = 0.95) + # 提取比较结果 + Tukey_HSD_table = as.data.frame(Tukey_HSD$group) + + # 保存统计结果 + # 保存一个制表符，解决存在行名时，列名无法对齐的问题 + write.table(paste(date(), "\nGroup\t", groupID, "\n\t", sep=""), file=paste("alpha_boxplot_TukeyHSD.txt",sep=""),append = T, quote = F, eol = "", row.names = F, col.names = F) + # 保存统计结果，有warning正常 + suppressWarnings(write.table(Tukey_HSD_table, file=paste("alpha_boxplot_TukeyHSD.txt",sep=""), append = T, quote = F, sep="\t", eol = "\n", na = "NA", dec = ".", row.names = T, col.names = T)) + + # 函数：将Tukey检验结果P值转换为显著字母分组 + # 输入文件为图基检验结果和分组 + generate_label_df = function(TUKEY, variable){ + # library(multcompView) + # 转换P值为字母分组 + ## 提取图基检验中分组子表的第4列P adjust值 + Tukey.levels = TUKEY[[variable]][,4] + # 方法1.multcompLetters函数将两两p值转换为字母，data.frame并生成列名为Letters的数据框 + Tukey.labels = data.frame(multcompLetters(Tukey.levels)['Letters']) + # 方法2. 解决字母顺序相反的问题 + # library(multcomp) + # tuk <- cld(glht(model, alternative = 'two.sided', linfct = mcp(group = 'Tukey')), sig = p, decreasing = TRUE) + # Tukey.labels <- data.frame(Letters=tuk$mcletters$Letters, stringsAsFactors = FALSE) + + # 按分组名字母顺序 + ## 提取字母分组行名为group组名 + Tukey.labels$group = rownames(Tukey.labels) + # 按组名的字母顺序排列，默认的Levels + Tukey.labels=Tukey.labels[order(Tukey.labels$group), ] + return(Tukey.labels) + } + + # 当只有两组时，用LSD标注字母 + if (length(unique(df$group)) == 2){ + # LSD检验，添加差异组字母 + library(agricolae) + out = LSD.test(model, "group", p.adj="none") + stat = out$groups + # 分组结果添入Index + df$stat=stat[as.character(df$group),]$groups + # 当大于两组时，用multcompView标注字母 + }else{ + # library(multcompView) + LABELS = generate_label_df(Tukey_HSD , "group") + df$stat=LABELS[as.character(df$group),]$Letters + } + + # 设置分组位置为各组y最大值+高的5% + max=max(df[,c(index)]) + min=min(df[,index]) + x = df[,c("group",index)] + #y = x %>% group_by(group) %>% summarise(Max=paste('max(',index,')',sep="")) + y = x %>% group_by(group) %>% summarise(Max = max(.data[[index]], na.rm = TRUE)) + y=as.data.frame(y) + rownames(y)=y$group + df$y=y[as.character(df$group),]$Max + (max-min)*0.05 + + + if (outlier) { + # 绘图 plotting + p = ggplot(df, aes(x=group, y=.data[[index]], color=group)) + + geom_boxplot(alpha=1, + # outlier.shape = NA, + # outlier.size=0, + size=0.7, + width=0.5, fill="transparent") + + labs(x="Groups", y=paste(index, "index"), color=groupID) + theme_classic() + + geom_text(data=df, aes(x=group, y=y, color=group, label=stat)) + + geom_jitter(position=position_jitter(0.17), size=1, alpha=0.7)+ + theme(text=element_text(family="sans", size=7)) + p + } else{ + # 绘图 plotting + p = ggplot(df, aes(x=group, y=.data[[index]], color=group)) + + geom_boxplot(alpha=1, + outlier.shape = NA, + outlier.size=0, + size=0.7, + width=0.5, fill="transparent") + + labs(x="Groups", y=paste(index, "index"), color=groupID) + theme_classic() + + geom_text(data=df, aes(x=group, y=y, color=group, label=stat)) + + geom_jitter(position=position_jitter(0.17), size=1, alpha=0.7)+ + theme(text=element_text(family="sans", size=7)) + p + } + + +} + +p = alpha_boxplot(alpha_div, index = opts$alpha_index, metadata, groupID = opts$group) +if (opts$xlabAngle){ + p = p + theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1)) +} +# Saving figure +# 保存图片，大家可以修改图片名称和位置，长宽单位为毫米 +ggsave(paste0(opts$output,"boxplot_",opts$alpha_index,".pdf"), p, width = opts$width, height = opts$height, units = "mm") From 6b377fffc3331f7e32bb8448f978275a8381d7a7 Mon Sep 17 00:00:00 2001 From: baidefeng <64759199+baidefeng@users.noreply.github.com> Date: Tue, 4 Aug 2026 16:39:39 +0800 Subject: [PATCH 19/19] Bump version to v1.25 and update date Updated version and last update date in README. --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 75c5f9a..0adfeaa 100644 --- a/README.md +++ b/README.md @@ -4,9 +4,9 @@ Popular software, scripts and database annotation for EasyAmplicon and EasyMetagenome -版本(Version)：EasyMicrobiome v1.24 +版本(Version)：EasyMicrobiome v1.25 -更新时间(Update)：2025/11/25 +更新时间(Update)：2026/8/4 项目主页(Project homepage): https://github.com/yongxinliu/EasyMicrobiome