生物信息分析 mibiogen 菌群分析 源码

PBS作业管理系统
在Linux系统中，PBS（Portable Batch System）是常用的作业管理工具，用于提交任务到服务器队列执行。用户通过编写.pbs后缀的脚本文件（如task.pbs）来定义任务参数，并通过命令qsub提交执行。 ‌
1

#!/bin/bash
# 【指定脚本解释器】这一行是Bash脚本的标识，告诉系统用Bash执行后续代码，是所有Bash脚本的必备开头。

#===========================================================
#  RAISS imputation (MiBioGen - 1000G) — PBS array job
#  【脚本核心用途】针对MiBioGen菌群GWAS数据（来自fmicb-15-1356437.pdf的MiBioGen数据集），
#  结合1000G参考面板进行RAISS基因型插补；采用PBS数组作业（array job）批量处理多性状+多染色体任务。
#  
#  One array index = one (trait, chromosome) job.
#  【数组作业逻辑】每个PBS数组索引对应1个“菌群性状+染色体”的插补任务（避免重复手动提交）。
#  
#  Array size: 22 chromosomes  *  190 traits = 4180 indices.
#  【数组规模说明】共处理22条常染色体（人类常染色体范围，符合fmicb-15-1356437.pdf 2.1节“常染色体遗传变异”要求），
#  190个MiBioGen菌群性状（来自设计方案.docx中的“菌群分类单元”，如属、科水平性状），总计4180个任务索引。
#  
#  Note: one trait has only 21 z-score files; that index will SKIP.
#  【特殊情况说明】某1个MiBioGen菌群性状（如MiBioGen.id.961）仅提供21条染色体的z-score文件（缺失1条），
#  对应缺失文件的数组索引会自动跳过，避免脚本报错。
#
#  Input z-score files already have this 6-column header:
#      ID    rsID    pos    A0    A1    Z
#  【输入文件格式】MiBioGen菌群GWAS的z-score文件（预处理产物，用于RAISS插补）需包含6列固定表头：
#      ID：SNP内部标识（与rsID一致，避免后续步骤丢失信息）；
#      rsID：SNP公共标识（用于匹配1000G参考面板，符合设计方案.docx“工具变量匹配”需求）；
#      pos：SNP基因组位置（碱基坐标，用于染色体定位）；
#      A0：非效应等位基因（参考等位基因）；
#      A1：效应等位基因（与fmicb-15-1356437.pdf 2.2节“等位基因方向”一致）；
#      Z：SNP与菌群性状的关联Z值（GWAS核心统计量）。
#  
#  Why both ID and rsID?
#    RAISS expects both. I prepared inputs with ID=rsID; otherwise
#    rsID can get dropped and downstream steps return error!!
#  【双ID字段原因】RAISS工具强制要求同时提供ID和rsID；若仅保留rsID，后续插补后可能丢失SNP标识，
#  导致下游MR分析（如设计方案.docx中的工具变量筛选）报错，因此设置ID=rsID确保一致性。
#===========================================================

# PBS SETTINGS 
# 【PBS作业调度配置区】以下行均为PBS调度系统的专用指令，定义作业在集群中的运行属性，仅PBS系统识别。
#PBS -N RAISS.array.MiBioGen.1000G.R2_0.3
#  【作业名称】指定作业名为“RAISS.array.MiBioGen.1000G.R2_0.3”，便于在集群中识别（包含核心信息：工具+数据+参数）：
#  RAISS=插补工具，MiBioGen=数据来源（fmicb-15-1356437.pdf），1000G=参考面板，R2_0.3=LD阈值参数。

#PBS -l select=1:ncpus=1:mem=16gb
#  【资源需求】向集群请求计算资源：
#  select=1：请求1个计算节点；
#  ncpus=1：每个节点分配1个CPU核心（避免多核心抢占，确保数组任务并行效率）；
#  mem=16gb：每个节点分配16GB内存（满足RAISS处理单个性状+染色体的内存需求）。

#PBS -l walltime=48:00:00
#  【运行时长限制】设置作业最大运行时间为48小时（2天），超时未完成的任务会被集群自动终止，
#  预留充足时间处理MiBioGen大数据量的插补（避免因数据量大导致任务中断）。

#PBS -J 1-4180
#  【数组作业范围】定义PBS数组的索引范围为1到4180，对应“22染色体×190性状”的所有任务，
#  集群会自动并行调度每个索引对应的任务（无需手动逐个提交）。

# SHELL SAFETY 
set -euo pipefail
#  【Bash脚本安全配置】开启3个关键安全选项，避免MiBioGen数据插补过程中出现“静默失败”（错误未察觉但结果异常）：
#  -e：若任一简单命令返回非0状态（执行失败），立即终止脚本（防止错误累积）；
#  -u：将未定义的变量视为错误（如拼写错误的变量名，避免用空值导致逻辑错误）；
#  pipefail：管道命令中任一环节失败，整个管道返回失败（而非仅最后一步，如z-score文件读取失败时立即报错）。

# I found it useful to avoid silent failures!
#  【注释说明】作者备注：该安全配置可有效避免“静默失败”——即脚本看似运行完成，但实际因隐藏错误导致插补结果无效，
#  对MiBioGen菌群数据的后续MR分析（设计方案.docx核心步骤）至关重要（错误数据会导致因果推断偏差）。

# -e  : makes the shell exit immediately if any simple command returns a non-zero status. It’s a guardrail so the  script doesn’t keep running after a failure.
#  【-e选项解释】简单命令执行失败（返回非0）时立即退出脚本，防止错误继续传播（如z-score文件缺失仍继续插补）。

# -u  : Treat unset variables as errors (catches typos/missing env).
#  【-u选项解释】未定义的变量被视为错误（如把“ZS”错写为“ZSS”时，脚本会报错而非用空值），避免环境变量缺失导致路径错误。

# pipefail: Make a pipeline fail if *any* element fails (not just the last command).
#  【pipefail选项解释】管道命令（如“cat file.txt | grep xxx”）中任一环节失败，整个管道返回失败，
#  例如“读取z-score文件失败”但“grep过滤成功”时，仍会报错（而非仅忽略读取错误）。


# THREADING 
export OMP_NUM_THREADS=1   
#  【线程控制】设置OpenMP（多线程库）的线程数为1，强制每个PBS数组任务使用单线程运行：
#  原因1：PBS数组已通过“多索引并行”实现任务并行，单任务单线程可避免CPU资源抢占（如多个线程争抢1个核心导致变慢）；
#  原因2：防止RAISS依赖的库（如线性代数库）自动启用多线程，超出请求的ncpus=1资源限制（集群可能终止违规任务）。

# stops libraries from secretly using multiple CPU threads. 
#  【注释说明】防止依赖库“秘密启用多线程”——部分科学计算库（如BLAS）会默认使用多线程，若不限制，
#  单个数组任务可能占用多个CPU核心，导致其他数组任务资源不足（拖慢整体插补效率）。

# That way each PBS array task really uses one CPU, like I requested, avoiding slowdowns. 
#  【注释说明】确保每个数组任务仅使用1个CPU核心（与PBS -l ncpus=1请求一致），避免因资源超用导致的集群调度延迟。

# espicially as I have many many subjobs!
#  【注释说明】强调：因数组任务总数达4180个（数量极多），单任务单线程是保障整体调度效率的关键。


# PATHS 
RAISS_HOME="/home/n11385332/RAISS"
#  【路径变量】定义RAISS工具的安装根目录，便于后续引用（避免重复写长路径，减少拼写错误）。

RAISS_BIN="${RAISS_HOME}/raiss-env/bin/raiss"
#  【路径变量】定义RAISS工具的可执行文件路径：
#  ${RAISS_HOME}：引用上述根目录；
#  raiss-env/bin/raiss：RAISS虚拟环境中的可执行文件（确保依赖库版本正确）。

REF="/home/n11385332/1000G_EUR_Phase3_plink"    
#  【路径变量】定义1000G欧洲人群参考面板的根目录：
#  参考面板用于MiBioGen基因型插补（fmicb-15-1356437.pdf 2.1节提到MiBioGen数据以欧洲人群为主，匹配参考面板种族）；
#  面板文件格式为PLINK（.bim/.bed/.fam），后续通过前缀/后缀拼接完整路径。

LD="${RAISS_HOME}/1000G/ld"                     
#  【路径变量】定义LD（连锁不平衡）文件的存储目录：
#  LD文件用于RAISS插补时的基因型相位推断，文件格式为PLINK生成的.ld文件（如chr1_*.ld）。

ZS="${RAISS_HOME}/zscores"                     
#  【路径变量】定义MiBioGen菌群GWAS的z-score文件目录：
#  z-score文件是RAISS的输入数据，命名格式为“z_<TRAIT>_chrN.txt”（TRAIT=菌群性状，N=染色体号），
#  对应`fmicb-15-1356437.pdf`中MiBioGen菌群与SNP的关联统计量。

OUT="${RAISS_HOME}/out.MiBioGen.1000G.R2_0.3"
#  【路径变量】定义RAISS插补结果的输出目录：
#  结果按“菌群性状”分文件夹存储（如${OUT}/${TAG}/），便于后续整理用于MR分析（设计方案.docx的工具变量筛选）。

TAGS_FILE="${RAISS_HOME}/tags.list"
#  【路径变量】定义菌群性状列表文件路径：
#  tags.list是纯文本文件，每行对应1个MiBioGen菌群性状名（如genus.RuminococcaceaeUCG004，来自设计方案的菌群分类单元），
#  用于通过数组索引匹配对应的菌群性状。

# Create output folder
mkdir -p "${OUT}"
#  【创建输出目录】使用mkdir -p命令创建输出根目录：
#  -p选项：若目录已存在，不报错；若父目录不存在，自动创建（避免因目录缺失导致插补结果无法写入）。


# ARRAY INDEX > WORK 
NCHR=22
#  【变量定义】设置染色体总数为22（人类常染色体，排除性染色体，符合fmicb-15-1356437.pdf 2.1节“常染色体遗传变异”要求）。

NTRAITS=$(wc -l < "${TAGS_FILE}")
#  【变量定义】计算菌群性状总数：
#  wc -l < "${TAGS_FILE}"：统计tags.list文件的行数（即性状数量，此处为190）；
#  用变量存储避免硬编码，若tags.list更新（如新增性状），脚本无需修改即可适配。

# Read the array index from PBS 
INDEX="$PBS_ARRAY_INDEX"
#  【读取数组索引】从PBS调度系统获取当前运行的数组索引（如1、2、...、4180），
#  每个索引对应1个“菌群性状+染色体”的插补任务。

# Map 1D index to (trait, chromosome)
TR_IDX=$(( (INDEX - 1) / NCHR + 1 ))    
#  【索引映射】将1维数组索引转为“菌群性状索引”（TR_IDX）：
#  计算逻辑：(索引-1) ÷ 染色体数（22），取整数商后+1，结果范围1~NTRAITS（1~190）；
#  例：INDEX=1 → (1-1)/22 +1=1（第1个性状），INDEX=23 → (23-1)/22 +1=2（第2个性状）。

CHR_IDX=$(( (INDEX - 1) % NCHR + 1 ))   
#  【索引映射】将1维数组索引转为“染色体索引”（CHR_IDX）：
#  计算逻辑：(索引-1) 对染色体数（22）取余，结果+1，范围1~22（对应chr1~chr22）；
#  例：INDEX=1 → (1-1)%22 +1=1（chr1），INDEX=22 → (22-1)%22 +1=22（chr22）。

CHR="chr${CHR_IDX}"
#  【变量定义】将染色体索引转为标准染色体名（如CHR_IDX=1 → chr1），匹配z-score文件和参考面板的命名格式。

#  examples
#   INDEX | TR_IDX (trait) | CHR_IDX (chr) | CHR   |
# | ---- | ------------- | ------------ | ----- |    
# |     1 |              1 |             1 | chr1  |
# |     2 |              1 |             2 | chr2  |
# |    22 |              1 |            22 | chr22 |
# |    23 |              2 |             1 | chr1  |
# |    44 |              2 |            22 | chr22 |
# |  4180 |            190 |            22 | chr22 |
#  【示例说明】通过表格直观展示索引映射逻辑，帮助理解“1个索引对应1个性状+1条染色体”的分配规则。


# Get the trait name from tags.list
TAG="$(sed -n "${TR_IDX}p" "${TAGS_FILE}")"
#  【获取性状名】从tags.list文件中读取第TR_IDX行的内容，作为当前任务的菌群性状名（TAG）：
#  sed -n "${TR_IDX}p"：sed命令的“安静模式”（-n），仅打印第TR_IDX行（p=print）；
#  例：TR_IDX=3 → 读取tags.list第3行，得到菌群性状名（如genus.Oscillibacter，来自fmicb-15-1356437.pdf 3.2节的关键菌属）。

# Point to the expected z-score file for this (trait, chr)
ZFILE="${ZS}/z_${TAG}_${CHR}.txt"
#  【构建z-score文件路径】拼接当前“性状+染色体”对应的z-score文件完整路径：
#  格式：${ZS}/z_<性状名>_<染色体名>.txt（如z_genus.Oscillibacter_chr1.txt），
#  对应MiBioGen菌群GWAS中“某菌属+某染色体”的关联z-score数据（fmicb-15-1356437.pdf 2.1节的SNP-菌群关联数据）。

# If that one z-score is missing (the trait MiBioGen.id.961 has only 21 files), skip
if ! test -s "$ZFILE"; then
  echo "[SKIP] Missing z-score: $ZFILE"
  exit 0
fi
#  【处理缺失文件】检查z-score文件是否存在且非空，若缺失则跳过当前任务：
#  test -s "$ZFILE"：test命令的-s选项，判断文件“存在且大小非0”；
#  ! test -s ...：否定判断（文件缺失或为空）；
#  若满足条件：打印跳过提示（说明缺失的文件路径），执行exit 0（正常退出脚本，不视为失败）；
#  适配特殊情况：如MiBioGen.id.961性状仅21个染色体文件，缺失的1个文件会触发此逻辑，避免脚本报错终止。

# test is a built-in Unix command (also available as a Bash builtin) that checks a condition and returns an exit status:
# 0 = true (condition passed)
# non-zero = false (condition failed)
#  【test命令说明】test是Unix/Bash内置命令，用于判断条件，返回退出状态码：0=条件成立（真），非0=条件不成立（假）。

# test -s "$ZFILE" is "file exists and non-empty”.
#  【test -s选项说明】test -s "$ZFILE"的含义是“文件存在且大小大于0”（非空），用于验证z-score文件有效性。

# ! negates it  >  "missing or empty”.
#  【否定逻辑说明】! 符号否定test命令的结果，即“文件缺失或为空”。


# Create log folder 
LOG_DIR="${RAISS_HOME}/logs.RAISS.array_1000G_R2_0.3"
mkdir -p "${LOG_DIR}"
#  【创建日志目录】创建用于存储任务日志的文件夹：
#  日志按任务单独存储，便于后续排查失败原因（如某个性状插补失败，可查看对应日志）。

# Per-task logs (simple)
exec > >(tee -a "${LOG_DIR}/${TAG}.${CHR}.idx${INDEX}.${PBS_JOBID:-NA}.out")
#  【重定向标准输出（stdout）】将脚本的标准输出（如echo的信息、RAISS的运行日志）重定向到日志文件：
#  exec > ...：重新定义标准输出的流向；
#  >(tee -a ...)：进程替换（Bash特性），将输出同时写入日志文件（tee -a，-a=追加模式）和终端；
#  日志文件名格式：${性状名}.${染色体名}.数组索引.作业ID.out（如genus.Oscillibacter.chr1.idx1.12345.out）；
#  ${PBS_JOBID:-NA}：若PBS_JOBID（作业ID）未定义，用NA填充（避免变量未定义报错，符合set -u安全配置）。

exec 2> >(tee -a "${LOG_DIR}/${TAG}.${CHR}.idx${INDEX}.${PBS_JOBID:-NA}.err" >&2)
#  【重定向标准错误（stderr）】将脚本的标准错误（如命令执行失败信息、RAISS报错）重定向到错误日志文件：
#  exec 2> ...：重新定义标准错误（文件描述符2）的流向；
#  tee -a ... >&2：将错误信息同时写入.err日志文件（追加模式）和终端（>&2=将tee的输出重定向回标准错误，保持终端报错提示）；
#  错误日志与输出日志同名，后缀为.err，便于对应查看“正常输出”和“错误信息”。

echo "[INFO] TAG=${TAG}  CHR=${CHR}  ZFILE=${ZFILE}"
#  【打印任务信息】输出当前任务的核心参数（性状名、染色体名、z-score文件路径），
#  同时写入日志文件，便于后续确认任务对应的MiBioGen数据信息（如是否为fmicb-15-1356437.pdf中的关键菌属）。


# RUN RAISS 
# RAISS reads z-scores for TAG * CHR, uses 1000G PLINK + LD,
# and writes outputs under ${OUT}/${TAG}/
mkdir -p "${OUT}/${TAG}"
#  【创建性状输出目录】为当前菌群性状创建单独的输出文件夹（如${OUT}/genus.Oscillibacter/）：
#  避免不同性状的插补结果混杂，便于后续按性状整理数据（设计方案.docx中“菌群数据按分类单元分组”的需求）。

"${RAISS_BIN}" \
  --chrom "${CHR}" \
  #  【RAISS参数】指定当前处理的染色体（如chr1），用于匹配参考面板和LD文件。
  --gwas "${TAG}" \
  #  【RAISS参数】指定当前处理的GWAS性状名（即MiBioGen菌群性状名，如genus.Oscillibacter），用于输出文件命名。
  --ref-folder "${REF}" \
  #  【RAISS参数】指定1000G参考面板的根目录（${REF}），RAISS从该目录读取PLINK格式的参考基因型数据。
  --ref-panel-prefix "1000G.EUR.QC." \
  #  【RAISS参数】指定参考面板文件的前缀：
  #  完整参考文件名为“1000G.EUR.QC.chr1.bim”（前缀+染色体+后缀），符合1000G欧洲人群参考面板的标准命名。
  --ref-panel-suffix ".bim" \
  #  【RAISS参数】指定参考面板的后缀（.bim），RAISS通过“前缀+染色体+后缀”拼接完整的.bim文件路径，
  #  进而找到对应的.bed/.fam文件（PLINK三文件格式）。
  --ld-folder "${LD}" \
  #  【RAISS参数】指定LD文件的存储目录（${LD}），RAISS从该目录读取对应染色体的LD文件，用于基因型相位推断。
  --zscore-folder "${ZS}" \
  #  【RAISS参数】指定z-score文件的根目录（${ZS}），RAISS从该目录读取当前性状+染色体的z-score文件。
  --output-folder "${OUT}" \
  #  【RAISS参数】指定插补结果的输出根目录（${OUT}），RAISS会在该目录下创建以${TAG}命名的子文件夹，存储结果文件。
  --ld-type plink \
  #  【RAISS参数】指定LD文件的格式为PLINK生成的.ld格式（RAISS支持的LD格式之一，确保兼容性）。
  --R2-threshold 0.3
#  【RAISS参数】设置LD的R²阈值为0.3：
#  仅保留与目标SNP的LD R²≥0.3的标记用于插补，过滤弱连锁的标记（符合GWAS数据预处理的标准，设计方案.docx中的“连锁不平衡控制”）。

echo "[INFO] Done: ${TAG} ${CHR}"
#  【打印完成信息】当前“性状+染色体”的RAISS插补任务执行完成后，输出提示信息，
#  同时写入日志文件，便于后续确认任务完成状态（如是否所有MiBioGen关键菌属的插补都已完成）。

这个另一个版本

#!/bin/bash
#===========================================================
#  PRED-LD imputation — PBS array job (1 task = 1 trait)
#  Avoids collisions by running each task in its own WORK dir
#===========================================================

# PBS
#PBS -N pred-ld_array_MiBioGen
#PBS -o pred-ld_array_MiBioGen.out
#PBS -e pred-ld_array_MiBioGen.err
#PBS -l select=1:ncpus=1:mem=128gb
#PBS -l walltime=24:00:00
#PBS -J 1-190

set -euo pipefail
export OMP_NUM_THREADS=1

#  set the paths
LIST="/home/n11385332/Pred_LD/MiBioGen.predld.inputs.txt"  # one input file per line (*.PREDLD.txt)
OUTDIR="/home/n11385332/Pred_LD/out.MiBioGen.Pred-LD"      # final collected outputs
LOGDIR="${OUTDIR}/logs"                                    # per-task logs
PREDLD_SRC="/home/n11385332/Pred_LD/PRED-LD"               # folder containing pred_ld.py
VENV="$HOME/predld-venv/bin/activate"                      # the Python venv to run pred_ld.py

# Create output/log dirs
mkdir -p "$OUTDIR" "$LOGDIR"

# Which input for this array index?
INDEX="$PBS_ARRAY_INDEX"
INFILE="$(sed -n "${INDEX}p" "$LIST")"
BASE="$(basename "$INFILE" .PREDLD.txt)"   # trait basename, used to prefix outputs

# Per-task logs (simple)
exec > "${LOGDIR}/${BASE}.predld.out"
exec 2> "${LOGDIR}/${BASE}.predld.err"


# These two lines will tell us exactly which taxa (i.e., BASE) is running for which array index:
echo "[INFO] $(date) INDEX=${INDEX} BASE=${BASE}"
echo "[INFO] INFILE=${INFILE}"

#  Isolated working directory 
# Each array task uses its own folder; prevents filename collisions.
WORK="${OUTDIR}/${BASE}"
mkdir -p "$WORK"
echo "[INFO] OUT=${WORK}"
cd "$WORK"


# NOTE: pred_ld.py writes fixed filenames like:
#   imputation_results_chr_all.txt
#   imputation_results_chr1.txt 
# Because we run in a unique WORK dir, nothing collides.

# But we still need to do something regarding --ref; otherwise I’ll get an error like:
# FileNotFoundError: [Errno 2] No such file or directory: 
# '/mnt/hpccs01/home/n11385332/Pred_LD/out.MiBioGen.Pred-LD/genus.Ruminiclostridium6.id.11356/ref/TOP_LD/EUR/SNV/EUR_chr1_no_filter_0.2_1000000_info_annotation.csv.gz' 
# To fix that, I add this line:
ln -sfn /home/n11385332/Pred_LD/PRED-LD/ref ref

# ln = the Unix command to make a link to a file or folder.
# Symbolic link: If the target moves or is deleted, the symlink becomes broken.
# Hard link: If you delete one name, the other still has the data.

# more about symlink:
# A symlink (symbolic link) is a tiny file that points to another file or folder by path. 
# Think of it like a desktop shortcut:
# When you open/read the symlink, the OS transparently follows it to the target.
# If the target changes, you see the change through the symlink too.
# If the target is moved/deleted, the symlink becomes broken (it points nowhere).
# -s = make a symlink
# -f: replace an existing file or symlink at the destination name.
# -n: if the destination is a symlink to a directory, treat it as a file to replace (don’t descend into it). Handy when recreating ref links.
# -f and -n are handy especially when: 
# rerun this script for the same trait (e.g., resubmit the array, retry after a failure),


# Run pred_ld.py inside WORK
source "$VENV"
# stay in WORK; call pred_ld.py by full path
python "${PREDLD_SRC}/pred_ld.py" \
  --file-path "$INFILE" \
  --r2threshold 0.8 \
  --pop EUR \
  --maf 0.01 \
  --ref TOP_LD

# Files like imputation_results_chr10.txt now live in:
# $OUTDIR/$BASE/imputation_results_chr10.txt

echo "[INFO] Done ${BASE}  (outputs in $WORK)"