HUMAnN3 处理脚本
#!/bin/bash
# HUMAnN3 处理脚本
# 作者: Hahn
# 创建日期: 2025-03-31
# 设置路径变量
BASE_DIR="/data/khanh/hurmicro"
INPUT_DIR="${BASE_DIR}/3.alignment/allx" # 输入质控后的reads所在路径
OUTPUT_DIR="${BASE_DIR}/5.annotation/humann3" # HUMAnN3输出路径
LOG_DIR="${OUTPUT_DIR}/logs" # 日志路径
SCRIPTS_DIR="${OUTPUT_DIR}/scripts" # 脚本目录
TEMP_DIR="${OUTPUT_DIR}/temp" # 临时文件目录
MERGED_DIR="${OUTPUT_DIR}/merged_results" # 合并结果目录
# 设置HUMAnN3数据库路径 - 请根据实际安装位置修改
NUCLEOTIDE_DB="/public/database/humann/humann3_databases/chocophlan" # ChocoPhlAn核苷酸数据库
PROTEIN_DB="/public/database/humann/humann3_databases/uniref" # UniRef蛋白数据库
UTILITY_MAPPING_DB="/public/database/humann/humann3_databases/utility_mapping" # 功能映射数据库
METAPHLAN_OPTIONS="--bowtie2db /home/khanh/.conda/envs/biobakery/lib/python3.7/site-packages/metaphlan/metaphlan_databases" # MetaPhlAn选项，包括数据库路径
# 设置HUMAnN3参数
THREADS=30 # 每个任务使用的CPU核心数
MEMORY="500G" # 最大内存使用量
# 设置并行参数
MAX_PARALLEL=3 # 最大并行样本数
# 创建必要的目录
mkdir -p "${OUTPUT_DIR}"
mkdir -p "${LOG_DIR}"
mkdir -p "${SCRIPTS_DIR}"
mkdir -p "${TEMP_DIR}"
mkdir -p "${MERGED_DIR}"
# 主日志文件
MAIN_LOG="${LOG_DIR}/humann3_parallel_$(date +%Y%m%d).log"
echo "===== 开始HUMAnN3并行分析 =====" | tee -a "${MAIN_LOG}"
echo "最大并行样本数: ${MAX_PARALLEL}" | tee -a "${MAIN_LOG}"
echo "开始时间: $(date)" | tee -a "${MAIN_LOG}"
# 首先，需要将双端reads合并成单个文件，因为HUMAnN3处理单端文件
echo "开始处理样本..." | tee -a "${MAIN_LOG}"
# 查找所有配对的双端reads文件
FQ1_FILES=$(find "${INPUT_DIR}" -name "*.un.1.fq.gz" | sort)
# 为每个样本创建单独的脚本
for FQ1 in ${FQ1_FILES}; do
# 获取样本名称
SAMPLE_NAME=$(basename "${FQ1}" .un.1.fq.gz)
# 构建第二个配对文件的路径
FQ2="${INPUT_DIR}/${SAMPLE_NAME}.un.2.fq.gz"
# 检查第二个文件是否存在
if [ ! -f "${FQ2}" ]; then
echo "警告: 未找到配对文件 ${FQ2}, 跳过 ${SAMPLE_NAME}" | tee -a "${MAIN_LOG}"
continue
fi
# 创建样本输出目录
SAMPLE_OUTPUT_DIR="${OUTPUT_DIR}/${SAMPLE_NAME}"
mkdir -p "${SAMPLE_OUTPUT_DIR}"
# 检查是否已完成分析
if [ -f "${SAMPLE_OUTPUT_DIR}/${SAMPLE_NAME}_genefamilies.tsv" ] && [ -f "${SAMPLE_OUTPUT_DIR}/${SAMPLE_NAME}_pathabundance.tsv" ]; then
echo "HUMAnN3结果已存在: ${SAMPLE_NAME}, 跳过" | tee -a "${MAIN_LOG}"
continue
fi
# 为这个样本创建单独的脚本
SAMPLE_SCRIPT="${SCRIPTS_DIR}/humann3_${SAMPLE_NAME}.sh"
# 写入脚本内容
cat > "${SAMPLE_SCRIPT}" << EOF
#!/bin/bash
# HUMAnN3分析脚本 - ${SAMPLE_NAME}
# 自动生成时间: $(date)
# 激活conda环境（如果需要）
source /path/to/conda/etc/profile.d/conda.sh
conda activate biobakery # 修改为包含HUMAnN3的环境名称
# 设置输出目录
OUTPUT_DIR="${SAMPLE_OUTPUT_DIR}"
LOG_FILE="${LOG_DIR}/${SAMPLE_NAME}_humann3.log"
TEMP_DIR="${TEMP_DIR}/${SAMPLE_NAME}"
# 创建临时目录
mkdir -p "\${TEMP_DIR}"
# 设置临时目录环境变量
export TMPDIR="\${TEMP_DIR}"
echo "===== 开始处理样本: ${SAMPLE_NAME} =====" | tee -a "\${LOG_FILE}"
echo "开始时间: \$(date)" | tee -a "\${LOG_FILE}"
# 步骤1: 合并双端reads
echo "合并双端reads..." | tee -a "\${LOG_FILE}"
MERGED_FASTQ="\${TEMP_DIR}/${SAMPLE_NAME}_merged.fastq.gz"
# 使用cat直接合并（HUMAnN3推荐的方法）
cat "${FQ1}" "${FQ2}" > "\${MERGED_FASTQ}"
# 检查合并是否成功
if [ ! -f "\${MERGED_FASTQ}" ] || [ ! -s "\${MERGED_FASTQ}" ]; then
echo "错误: reads合并失败: ${SAMPLE_NAME}" | tee -a "\${LOG_FILE}"
exit 1
fi
# 步骤2: 运行HUMAnN3
echo "运行HUMAnN3..." | tee -a "\${LOG_FILE}"
# 注意：HUMAnN3实际命令是humann而不是humann3
# 另外，HUMAnN不支持--tmp-dir参数，它使用TMPDIR环境变量
humann --input "\${MERGED_FASTQ}" \\
--output "\${OUTPUT_DIR}" \\
--threads ${THREADS} \\
--memory-use maximum \\
--nucleotide-database ${NUCLEOTIDE_DB} \\
--protein-database ${PROTEIN_DB} \\
--metaphlan-options "${METAPHLAN_OPTIONS}" \\
--output-basename "${SAMPLE_NAME}" \\
--remove-temp-output \\
--o-log "\${LOG_FILE}.humann3.log" \\
2>&1 | tee -a "\${LOG_FILE}"
# 检查HUMAnN3是否成功
if [ \$? -eq 0 ]; then
echo "HUMAnN3完成: ${SAMPLE_NAME}" | tee -a "\${LOG_FILE}"
# 标准化结果
echo "标准化结果..." | tee -a "\${LOG_FILE}"
# 标准化基因家族丰度
if [ -f "\${OUTPUT_DIR}/${SAMPLE_NAME}_genefamilies.tsv" ]; then
humann_renorm_table --input "\${OUTPUT_DIR}/${SAMPLE_NAME}_genefamilies.tsv" \\
--output "\${OUTPUT_DIR}/${SAMPLE_NAME}_genefamilies_cpm.tsv" \\
--units cpm \\
2>&1 | tee -a "\${LOG_FILE}"
else
echo "警告: 未找到基因家族文件: \${OUTPUT_DIR}/${SAMPLE_NAME}_genefamilies.tsv" | tee -a "\${LOG_FILE}"
fi
# 按相对丰度标准化通路丰度
if [ -f "\${OUTPUT_DIR}/${SAMPLE_NAME}_pathabundance.tsv" ]; then
humann_renorm_table --input "\${OUTPUT_DIR}/${SAMPLE_NAME}_pathabundance.tsv" \\
--output "\${OUTPUT_DIR}/${SAMPLE_NAME}_pathabundance_relab.tsv" \\
--units relab \\
2>&1 | tee -a "\${LOG_FILE}"
else
echo "警告: 未找到通路丰度文件: \${OUTPUT_DIR}/${SAMPLE_NAME}_pathabundance.tsv" | tee -a "\${LOG_FILE}"
fi
# 记录成功状态
touch "\${OUTPUT_DIR}/.analysis_complete"
# 清理临时文件
rm -rf "\${TEMP_DIR}"
else
echo "错误: HUMAnN3失败: ${SAMPLE_NAME}" | tee -a "\${LOG_FILE}"
exit 1
fi
echo "样本处理完成: ${SAMPLE_NAME}, 时间: \$(date)" | tee -a "\${LOG_FILE}"
exit 0
EOF
# 添加执行权限
chmod +x "${SAMPLE_SCRIPT}"
echo "已创建样本脚本: ${SAMPLE_SCRIPT}" | tee -a "${MAIN_LOG}"
done
# 创建并行执行脚本
PARALLEL_SCRIPT="${SCRIPTS_DIR}/run_humann3_parallel.sh"
cat > "${PARALLEL_SCRIPT}" << 'EOF'
#!/bin/bash
# 并行运行HUMAnN3分析任务
# 生成时间: $(date)
# 设置最大并行任务数
MAX_PARALLEL=3
# 脚本目录
SCRIPTS_DIR="$(dirname "$0")"
LOG_DIR="/data/khanh/hurmicro/5.annotation/humann3/logs"
MERGED_DIR="/data/khanh/hurmicro/5.annotation/humann3/merged_results"
UTILITY_MAPPING_DB="/path/to/utility_mapping" # 功能映射数据库
# 创建日志目录
mkdir -p "${LOG_DIR}"
mkdir -p "${MERGED_DIR}"
PARALLEL_LOG="${LOG_DIR}/parallel_execution.log"
echo "===== 开始并行HUMAnN3分析 =====" | tee -a "${PARALLEL_LOG}"
echo "最大并行任务数: ${MAX_PARALLEL}" | tee -a "${PARALLEL_LOG}"
echo "开始时间: $(date)" | tee -a "${PARALLEL_LOG}"
# 查找所有样本脚本
SAMPLE_SCRIPTS=$(find "${SCRIPTS_DIR}" -name "humann3_*.sh" | sort)
# 并行运行脚本
for SCRIPT in ${SAMPLE_SCRIPTS}; do
# 检查当前运行的任务数
RUNNING_JOBS=$(jobs -p | wc -l)
# 如果运行的任务达到最大值，等待任意一个任务完成
while [ ${RUNNING_JOBS} -ge ${MAX_PARALLEL} ]; do
sleep 30
RUNNING_JOBS=$(jobs -p | wc -l)
done
# 运行脚本
SAMPLE_NAME=$(basename "${SCRIPT}" .sh | sed 's/humann3_//')
echo "启动 ${SAMPLE_NAME} 处理 ($(date))" | tee -a "${PARALLEL_LOG}"
bash "${SCRIPT}" &
done
# 等待所有样本处理完成
echo "等待所有样本处理完成..." | tee -a "${PARALLEL_LOG}"
wait
# 合并所有样本结果
echo "所有样本处理完成，合并结果 ($(date))" | tee -a "${PARALLEL_LOG}"
# 合并基因家族结果
echo "合并基因家族丰度结果..." | tee -a "${PARALLEL_LOG}"
# 使用-s参数搜索子目录
humann_join_tables -i "/data/khanh/hurmicro/5.annotation/humann3" \
-o "${MERGED_DIR}/all_samples_genefamilies.tsv" \
--file_name "genefamilies.tsv" \
-s \
2>&1 | tee -a "${PARALLEL_LOG}"
# 合并通路丰度结果
echo "合并通路丰度结果..." | tee -a "${PARALLEL_LOG}"
# 使用-s参数搜索子目录
humann_join_tables -i "/data/khanh/hurmicro/5.annotation/humann3" \
-o "${MERGED_DIR}/all_samples_pathabundance.tsv" \
--file_name "pathabundance.tsv" \
-s \
2>&1 | tee -a "${PARALLEL_LOG}"
# 合并通路覆盖度结果
echo "合并通路覆盖度结果..." | tee -a "${PARALLEL_LOG}"
# 使用-s参数搜索子目录
humann_join_tables -i "/data/khanh/hurmicro/5.annotation/humann3" \
-o "${MERGED_DIR}/all_samples_pathcoverage.tsv" \
--file_name "pathcoverage.tsv" \
-s \
2>&1 | tee -a "${PARALLEL_LOG}"
# 标准化合并结果
echo "标准化合并结果..." | tee -a "${PARALLEL_LOG}"
# 标准化基因家族丰度(CPM)
if [ -f "${MERGED_DIR}/all_samples_genefamilies.tsv" ]; then
humann_renorm_table --input "${MERGED_DIR}/all_samples_genefamilies.tsv" \
--output "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" \
--units cpm \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 合并的基因家族文件不存在，跳过标准化" | tee -a "${PARALLEL_LOG}"
fi
# 标准化通路丰度(相对丰度)
if [ -f "${MERGED_DIR}/all_samples_pathabundance.tsv" ]; then
humann_renorm_table --input "${MERGED_DIR}/all_samples_pathabundance.tsv" \
--output "${MERGED_DIR}/all_samples_pathabundance_relab.tsv" \
--units relab \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 合并的通路丰度文件不存在，跳过标准化" | tee -a "${PARALLEL_LOG}"
fi
# 分组基因家族到GO功能
echo "分组基因家族到GO功能..." | tee -a "${PARALLEL_LOG}"
if [ -f "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" ]; then
humann_regroup_table --input "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" \
--output "${MERGED_DIR}/all_samples_GO_cpm.tsv" \
--groups uniref90_go \
--custom "${UTILITY_MAPPING_DB}" \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 标准化的基因家族文件不存在，跳过GO分组" | tee -a "${PARALLEL_LOG}"
fi
# 分组基因家族到KO功能
echo "分组基因家族到KO功能..." | tee -a "${PARALLEL_LOG}"
if [ -f "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" ]; then
humann_regroup_table --input "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" \
--output "${MERGED_DIR}/all_samples_KO_cpm.tsv" \
--groups uniref90_ko \
--custom "${UTILITY_MAPPING_DB}" \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 标准化的基因家族文件不存在，跳过KO分组" | tee -a "${PARALLEL_LOG}"
fi
# 分组基因家族到EC功能
echo "分组基因家族到EC功能..." | tee -a "${PARALLEL_LOG}"
if [ -f "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" ]; then
humann_regroup_table --input "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" \
--output "${MERGED_DIR}/all_samples_EC_cpm.tsv" \
--groups uniref90_level4ec \
--custom "${UTILITY_MAPPING_DB}" \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 标准化的基因家族文件不存在，跳过EC分组" | tee -a "${PARALLEL_LOG}"
fi
# 分组基因家族到MetaCyc反应
echo "分组基因家族到MetaCyc反应..." | tee -a "${PARALLEL_LOG}"
if [ -f "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" ]; then
humann_regroup_table --input "${MERGED_DIR}/all_samples_genefamilies_cpm.tsv" \
--output "${MERGED_DIR}/all_samples_RXN_cpm.tsv" \
--groups uniref90_rxn \
--custom "${UTILITY_MAPPING_DB}" \
2>&1 | tee -a "${PARALLEL_LOG}"
else
echo "警告: 标准化的基因家族文件不存在，跳过MetaCyc反应分组" | tee -a "${PARALLEL_LOG}"
fi
echo "HUMAnN3分析完成，结果位于: ${MERGED_DIR}" | tee -a "${PARALLEL_LOG}"
echo "全部处理完成时间: $(date)" | tee -a "${PARALLEL_LOG}"
EOF
# 添加执行权限
chmod +x "${PARALLEL_SCRIPT}"
echo "已创建并行执行脚本: ${PARALLEL_SCRIPT}" | tee -a "${MAIN_LOG}"
echo ""
echo "使用方法:" | tee -a "${MAIN_LOG}"
echo "运行并行HUMAnN3分析: bash ${PARALLEL_SCRIPT}" | tee -a "${MAIN_LOG}"
echo ""
echo "准备工作完成时间: $(date)" | tee -a "${MAIN_LOG}"
# 自动启动并行执行
echo "是否立即启动并行执行 [y/n]?"
read ANSWER
if [ "${ANSWER}" = "y" ] || [ "${ANSWER}" = "Y" ]; then
echo "启动并行执行..." | tee -a "${MAIN_LOG}"
bash "${PARALLEL_SCRIPT}"
else
echo "准备完成。请手动执行: bash ${PARALLEL_SCRIPT}" | tee -a "${MAIN_LOG}"
fi
本文来自博客园，作者：Hahntoh，转载请注明原文链接：https://www.cnblogs.com/hahn/p/18802947
浙公网安备 33010602011771号