Kraken2 物种注释
#!/bin/bash
# 单样本 Kraken2 分析脚本
# kraken2 --db /public/database/kraken/kraken2/PlusPF/20241228 --output /data/khanh/hurmicro/5.annotation/kraken2/2_kraken2.output --report /data/khanh/hurmicro/5.annotation/kraken2/2_kraken2.report --threads 30 --paired /data/khanh/hurmicro/3.alignment/all/2.un.1.fq.gz /data/khanh/hurmicro/3.alignment/all/2.un.2.fq.gz --use-names --report-zero-counts --confidence 0.1 --memory-mapping --gzip-compressed --minimum-hit-groups 3
# nohup ./run_all_parallel.sh 2>&1 &
# 方案1: 为每个组创建单独的脚本，可以手动同时启动多个任务
# 这个脚本生成每个组的独立脚本文件
# 设置路径变量
BASE_DIR="/data/khanh/hurmicro"
INPUT_BASE_DIR="${BASE_DIR}/5.annotation/sgroup"
OUTPUT_DIR="${BASE_DIR}/5.annotation/kraken2"
DATABASE_DIR="/public/database/kraken/kraken2/PlusPF/20241228"
SCRIPTS_DIR="${INPUT_BASE_DIR}/scripts"
# 创建脚本目录
mkdir -p "$SCRIPTS_DIR"
# 查找所有组目录
GROUP_DIRS=$(find "$INPUT_BASE_DIR" -type d -name "group*" | sort)
# 为每个组创建单独的脚本
for GROUP_DIR in $GROUP_DIRS; do
GROUP_NAME=$(basename "$GROUP_DIR")
# 创建该组的脚本文件
SCRIPT_FILE="${SCRIPTS_DIR}/process_${GROUP_NAME}.sh"
# 写入脚本头部
cat > "$SCRIPT_FILE" << EOF
#!/bin/bash
# 自动生成的处理 ${GROUP_NAME} 的脚本
# 生成时间: $(date)
# 设置路径变量
BASE_DIR="/data/khanh/hurmicro"
INPUT_DIR="${INPUT_BASE_DIR}/${GROUP_NAME}"
OUTPUT_DIR="${BASE_DIR}/5.annotation/kraken2"
DATABASE_DIR="/public/database/kraken/kraken2/PlusPF/20241228"
LOG_DIR="${OUTPUT_DIR}/logs/${GROUP_NAME}"
# 设置参数
THREADS=15
READ_LEN=150
# 创建输出和日志目录
mkdir -p "\$OUTPUT_DIR"
mkdir -p "\$LOG_DIR"
# 日志文件
LOG_FILE="\$LOG_DIR/${GROUP_NAME}_process.log"
echo "===== 开始处理 ${GROUP_NAME} =====" | tee -a "\$LOG_FILE"
echo "开始时间: \$(date)" | tee -a "\$LOG_FILE"
# 查找所有 .un.1.fq.gz 文件
FQ1_FILES=\$(find "\$INPUT_DIR" -name "*.un.1.fq.gz")
# 处理每个样本
for FQ1 in \$FQ1_FILES; do
# 获取样本ID
SAMPLE_ID=\$(basename "\$FQ1" .un.1.fq.gz)
# 获取第二个配对文件
FQ2="\${INPUT_DIR}/\${SAMPLE_ID}.un.2.fq.gz"
# 检查第二个文件是否存在
if [ ! -f "\$FQ2" ]; then
echo "警告: 未找到配对文件 \$FQ2, 跳过 \$SAMPLE_ID" | tee -a "\$LOG_FILE"
continue
fi
# 设置输出文件路径
OUTPUT_FILE="\${OUTPUT_DIR}/\${SAMPLE_ID}_kraken2.output"
REPORT_FILE="\${OUTPUT_DIR}/\${SAMPLE_ID}_kraken2.report"
SAMPLE_LOG="\$LOG_DIR/\${SAMPLE_ID}_kraken2.log"
echo "运行 Kraken2: \$SAMPLE_ID" | tee -a "\$LOG_FILE"
# 检查是否已经完成 Kraken2 分析
if [ -f "\$REPORT_FILE" ] && [ -f "\$OUTPUT_FILE" ]; then
echo "Kraken2 结果已存在: \$SAMPLE_ID, 跳过" | tee -a "\$LOG_FILE"
continue
fi
# 运行 Kraken2
kraken2 --db "\$DATABASE_DIR" \\
--output "\$OUTPUT_FILE" \\
--report "\$REPORT_FILE" \\
--threads "\$THREADS" \\
--paired "\$FQ1" "\$FQ2" \\
--use-names \\
--report-zero-counts \\
--confidence 0.1 \\
--memory-mapping \\
--gzip-compressed \\
--minimum-hit-groups 3 2>&1 | tee -a "\$SAMPLE_LOG"
# 检查 Kraken2 是否成功
if [ \$? -eq 0 ]; then
echo "Kraken2 完成: \$SAMPLE_ID" | tee -a "\$LOG_FILE"
else
echo "错误: Kraken2 失败: \$SAMPLE_ID" | tee -a "\$LOG_FILE"
fi
done
echo "${GROUP_NAME} 处理完成时间: \$(date)" | tee -a "\$LOG_FILE"
EOF
# 添加执行权限
chmod +x "$SCRIPT_FILE"
echo "已创建脚本: $SCRIPT_FILE"
done
# 创建一个汇总脚本，用于 Bracken 分析(所有组完成后运行)
BRACKEN_SCRIPT="${SCRIPTS_DIR}/run_bracken.sh"
# 写入 Bracken 脚本
cat > "$BRACKEN_SCRIPT" << 'EOF'
#!/bin/bash
# Bracken 分析脚本 - 在所有 Kraken2 分析完成后运行
# 设置路径变量
BASE_DIR="/data/khanh/hurmiciro"
OUTPUT_DIR="${BASE_DIR}/5.annotation/bracken"
DATABASE_DIR="/public/database/kraken/kraken2/PlusPF/20241228"
LOG_DIR="${OUTPUT_DIR}/logs"
# 设置参数
READ_LEN=150
# 日志文件
LOG_FILE="${LOG_DIR}/bracken_process.log"
echo "===== 开始 Bracken 分析 =====" | tee -a "$LOG_FILE"
echo "开始时间: $(date)" | tee -a "$LOG_FILE"
# 创建分类级别输出目录
LEVELS=("S" "G" "F" "O" "C" "P" "D")
LEVEL_NAMES=("species" "genus" "family" "order" "class" "phylum" "domain")
# 创建各分类级别的目录
for LEVEL_NAME in "${LEVEL_NAMES[@]}"; do
mkdir -p "${OUTPUT_DIR}/${LEVEL_NAME}"
done
# 处理所有 Kraken2 报告文件
KRAKEN_REPORTS=$(find "$OUTPUT_DIR" -name "*_kraken2.report")
for REPORT_FILE in $KRAKEN_REPORTS; do
# 获取样本ID
SAMPLE_ID=$(basename "$REPORT_FILE" _kraken2.report)
echo "运行 Bracken: $SAMPLE_ID" | tee -a "$LOG_FILE"
# 对每个分类级别运行 Bracken
for i in "${!LEVELS[@]}"; do
LEVEL="${LEVELS[$i]}"
LEVEL_NAME="${LEVEL_NAMES[$i]}"
BRACKEN_OUT="${OUTPUT_DIR}/${LEVEL_NAME}/${SAMPLE_ID}_bracken_${LEVEL_NAME}.txt"
# 检查是否已经完成 Bracken 分析
if [ -f "$BRACKEN_OUT" ]; then
echo "Bracken $LEVEL_NAME 结果已存在: $SAMPLE_ID, 跳过" | tee -a "$LOG_FILE"
continue
fi
# 运行 Bracken
bracken -d "$DATABASE_DIR" \
-i "$REPORT_FILE" \
-o "$BRACKEN_OUT" \
-r "$READ_LEN" \
-l "$LEVEL" \
-t 0 2>&1 | tee -a "${LOG_DIR}/${SAMPLE_ID}_bracken_${LEVEL_NAME}.log"
# 检查 Bracken 是否成功
if [ $? -eq 0 ]; then
echo "Bracken $LEVEL_NAME 完成: $SAMPLE_ID" | tee -a "$LOG_FILE"
else
echo "错误: Bracken $LEVEL_NAME 失败: $SAMPLE_ID" | tee -a "$LOG_FILE"
fi
done
done
# 尝试合并结果
COMBINE_SCRIPT=$(which combine_bracken_outputs.py 2>/dev/null)
if [ -z "$COMBINE_SCRIPT" ]; then
echo "警告: 未找到combine_bracken_outputs.py脚本，跳过合并阶段" | tee -a "$LOG_FILE"
else
# 对每个分类级别合并结果
for LEVEL_NAME in "${LEVEL_NAMES[@]}"; do
LEVEL_DIR="${OUTPUT_DIR}/${LEVEL_NAME}"
COMBINED_OUT="${OUTPUT_DIR}/combined_bracken_${LEVEL_NAME}.txt"
echo "合并$LEVEL_NAME级别的Bracken结果" | tee -a "$LOG_FILE"
# 查找所有匹配的文件并保存到临时文件中
FILE_LIST="${LOG_DIR}/bracken_${LEVEL_NAME}_files.txt"
find "$LEVEL_DIR" -name "*_bracken_${LEVEL_NAME}.txt" > "$FILE_LIST"
# 检查是否找到了文件
if [ ! -s "$FILE_LIST" ]; then
echo "警告: 未找到${LEVEL_NAME}级别的Bracken结果文件，跳过合并" | tee -a "$LOG_FILE"
continue
fi
# 显示找到的文件
FILE_COUNT=$(wc -l < "$FILE_LIST")
echo "找到$FILE_COUNT个${LEVEL_NAME}级别的文件" | tee -a "$LOG_FILE"
# 构建文件参数列表
FILES_ARG=""
while IFS= read -r file; do
FILES_ARG+="\"$file\" "
done < "$FILE_LIST"
# 运行合并脚本
eval python "$COMBINE_SCRIPT" --files $FILES_ARG --output "$COMBINED_OUT" 2>&1 | tee -a "${LOG_DIR}/combine_${LEVEL_NAME}.log"
if [ $? -eq 0 ]; then
echo "合并$LEVEL_NAME完成: $COMBINED_OUT" | tee -a "$LOG_FILE"
else
echo "错误: 合并$LEVEL_NAME失败" | tee -a "$LOG_FILE"
fi
done
fi
echo "Bracken 分析完成时间: $(date)" | tee -a "$LOG_FILE"
EOF
# 添加执行权限
chmod +x "$BRACKEN_SCRIPT"
echo "已创建 Bracken 脚本: $BRACKEN_SCRIPT"
# 创建一个运行所有任务的脚本
MAIN_SCRIPT="${SCRIPTS_DIR}/run_all_parallel.sh"
# 写入主脚本
cat > "$MAIN_SCRIPT" << 'EOF'
#!/bin/bash
# 并行运行所有组处理脚本
# 注意: 这将同时启动多个任务，请确保您的系统有足够的资源
# 设置最大并行任务数 (根据服务器资源调整)
MAX_PARALLEL=9
# 脚本目录
SCRIPTS_DIR="$(dirname "$0")"
LOG_DIR="/data/khanh/hurmicro/5.annotation/kraken2/logs"
mkdir -p "$LOG_DIR"
MAIN_LOG="$LOG_DIR/parallel_execution.log"
echo "===== 开始并行处理所有组 =====" | tee -a "$MAIN_LOG"
echo "最大并行任务数: $MAX_PARALLEL" | tee -a "$MAIN_LOG"
echo "开始时间: $(date)" | tee -a "$MAIN_LOG"
# 查找所有组处理脚本
GROUP_SCRIPTS=$(find "$SCRIPTS_DIR" -name "process_group*.sh" | sort)
# 并行运行脚本
for SCRIPT in $GROUP_SCRIPTS; do
# 检查当前运行的任务数
RUNNING_JOBS=$(jobs -p | wc -l)
# 如果运行的任务达到最大值，等待任意一个任务完成
while [ $RUNNING_JOBS -ge $MAX_PARALLEL ]; do
sleep 10
RUNNING_JOBS=$(jobs -p | wc -l)
done
# 运行脚本
GROUP_NAME=$(basename "$SCRIPT" .sh | sed 's/process_//')
echo "启动 $GROUP_NAME 处理 ($(date))" | tee -a "$MAIN_LOG"
bash "$SCRIPT" &
done
# 等待所有组处理完成
echo "等待所有组处理完成..." | tee -a "$MAIN_LOG"
wait
# 所有组处理完成后，运行 Bracken 分析
echo "所有组处理完成，开始 Bracken 分析 ($(date))" | tee -a "$MAIN_LOG"
bash "${SCRIPTS_DIR}/run_bracken.sh"
echo "全部处理完成时间: $(date)" | tee -a "$MAIN_LOG"
EOF
# 添加执行权限
chmod +x "$MAIN_SCRIPT"
echo "已创建主并行脚本: $MAIN_SCRIPT"
echo ""
echo "使用方法:"
echo "1. 独立运行各组: 分别执行每个 process_groupX.sh 脚本"
echo "2. 并行运行所有组: 执行 $MAIN_SCRIPT 脚本 (同时处理多个组)"
echo "3. 在所有 Kraken2 分析完成后运行 Bracken: 执行 $BRACKEN_SCRIPT 脚本"
本文来自博客园，作者：Hahntoh，转载请注明原文链接：https://www.cnblogs.com/hahn/p/18803003
浙公网安备 33010602011771号