B10K生信入门教程 · Guangji Chen

这是一个持续会更新的教程，专为B10K入门的新人们编写，旨在快速理解并使用一些常用脚本、软件等……

1. Bash

1.1 环境变量

当我们使用mobaxterm或iterm2以ssh usrname@ipaddress登录到集群上时，文件/home/usrname/.bash_profile会被加载。

我们常说的环境变量记录在了这一文件中。手动加载该文件的方式为source ~/.bash_profile。此处 ~ 即为 /home/usrname/。

为了更好的知道，我们上一行命令是否执行完成，我们有两种做法：

1. 执行命令并输出一个tag文件，需要用到&&命令。

   • 例如sleep 60 && touch sleep_60.done 或者 sleep && echo sleep_60 done。
   • 两者的差异是执行完sleep 60后touch命令生成了一个新文件 和 echo命令在屏幕端输出了一行字符。

2. 修改环境变量中PS1参数，通过^_^ 与O_O两个符号判断命令行是否执行成功，可以减少屏幕端的输出或一个新文件的生成。

   # 良渚端PS1参数示例
   #### BASIC ENV ####
   export PS1="\
   \[\033[33;1m\]\h:\u \033[35;1m\t \
   \$(if [[ \$? -gt 0 ]]; then printf \"\\[\\033[01;31m\\]O_O\"; else printf \"\\[\\033[36;1m\\]^_^\"; fi)
   \[\033[36;1m\]\w\
   \[\033[0m\] LiangZhu
   \[\e[32;1m\]$ \
   \[\e[0m\]"
   export LANG="en_US.UTF-8"
   export LC_ALL="en_US.UTF-8"
   export HISTSIZE=9999
   export HISTFILESIZE=9999

   将上述代码添加至文件/home/usrname/.bash_profile后，通过source手动加载该文件，输入sleep 60 后使用ctrl+c的方式感受一下^_^ 与O_O。

1.2 常用缩写

仔细观察文件/home/usrname/.bash_profile后，会发现除了加载当前文件以外还额外加载了文件/home/usrname/.bashrc，我会把常用的一些缩写写在该文件中。此处，我直接放出我的常用示例：

alias rm='rm -i'
alias mv='mv -i'
alias le='less -S'
alias ls="ls -v --color=auto"
alias ll="ls -lh -v"
alias lt="ls -lh -t -r"
alias la="ls -a -lh -v"
alias wl='ll | wc -l'
alias c='clear'  # 快速清屏
alias p='pwd'

alias vi="vim"
alias ps="/bin/ps ux --sort=pcpu,pmem,user,pid | cat"
alias h.my="htop -u chenguangji"

## Colorize the grep command output for ease of use (good for log files)##
alias grep='grep --color=auto'
alias egrep='egrep --color=auto'
alias fgrep='fgrep --color=auto'

alias tf='tail -f ' 
alias j='jobs -l'
alias tf='tail -f '
alias vi='vim'
alias which='alias | /usr/bin/which --tty-only --read-alias --show-dot --show-tilde'
alias ta='tmux attach'

# 进入目录并列出文件
alias ..="cdl .."
alias .2="cd ../.."   # 快速进入上上层目录
alias .3="cd ../../.."
alias cd..='cdl ..'

#Start calculator with math support
alias bc='bc -l'

# sbatch
alias sb='sbatch '
alias sq='squeue -u chenguangji'
alias sj='scontrol show job '
alias sdir='echo "Running" && sq | awk '\''$5=="R"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "WorkDir"|sort|uniq -c && echo "PENDING" && sq | awk '\''$5=="PD"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "WorkDir"|sort|uniq -c'
alias sdirR='echo "Running" && sq | awk '\''$5=="R"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "WorkDir"|sort|uniq -c'
alias sdirP='echo "PENDING" && sq | awk '\''$5=="PD"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "WorkDir"|sort|uniq -c'
alias scom='echo "Running" && sq | awk '\''$5=="R"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "Command"|sort|uniq -c && echo "PENDING" && sq | awk '\''$5=="PD"'\''|awk '\''{print "scontrol show job "$1}'\''|sh|grep "Command"|sort|uniq -c'

# quota
alias quota='mmlsquota -u chenguangji --block-size auto;mmlsquota -g zhanglab --block-size auto;echo "";df -h ./'

我最为常用的命令有：ta、sq、sdir、quota以及..与.3。

2. B10K 数据与 Codex_bin 工作流（新增）

下面的内容是在原教程基础上的补充，介绍 B10K Family Phase 数据来源，以及 ChenGuangji/Codex_bin 中 FASTA、phylogeny 和 trait 相关脚本的基本用法。命令默认在 Codex_bin 仓库根目录运行。

一个常见的基础流程是：

FASTA ID 整理与序列检查
        ↓
CDS 翻译、蛋白比对和 codon projection
        ↓
系统树裁剪、格式整理和重根
        ↓
B10K 类群分组与绘图
        ↓
性状表整理、缺失值填补和比较分析

脚本成功运行只说明输入格式和程序接口基本正确，不等于相关生物学结论已经完成验证。

2.1 数据和分析流程从哪里获得？

本部分主要对应B10K Family Phase 的数据和分析流程来源于这篇文章：

• Complexity of avian evolution revealed by family-level genomes

论文的支撑数据和分析材料，包括相关数据文件与脚本，可以在下面的长期存档中找到：

• Family Phase supporting data and scripts

良渚集群内部也下载了一份工作镜像：

/share/home/project/zhanglab/B10K/B10K_PMO/02.Family_Phase/Nature2024/B10K/data_upload/

该路径仅供具有相应集群和项目权限的成员访问。正式分析时应记录 DOI 存档中的数据版本和说明，不要只记录内部复制路径；也不要把受限数据或未经确认的大文件重新公开上传。

2.2 Codex_bin 的代码结构

仓库位置:

/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/

仓库把通用函数和具体工作流分开保存：

/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/src/R/                         # 可复用 R 函数
/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/src/python/for_codex_utils/    # 轻依赖 Python 函数
/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/                     # 文件输入、参数和结果输出

B10K 树图相关文件为：

/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/src/R/phylogeny/clade_grouping.R
/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/B10K/363_clade_tree.R

前者负责读取和校验系统树及分组表，后者保存 B10K 的类群颜色、扇形布局和输出参数。

2.3 准备环境 [可选]

基础 R 依赖：

install.packages(c("ape", "phangorn", "ggplot2"))

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("ggtree")

性状分析可能还需要：

install.packages(c("missForest", "phytools"))

FASTA 整理脚本只使用 Python 标准库。MAFFT 或 MUSCLE 需要单独安装，并确保程序位于 PATH。

可先查看帮助：

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/fasta_tool.py --help
python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/normalize_fasta.py --help
python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/translate_cds.py --help

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/trait/impute_trait_table.R
Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/extract_subtree.R
Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/reroot_tree.R

2.4 FASTA：先统一 ID

同一个物种在 CDS、蛋白、alignment、tree 和 trait table 中使用不同 ID，是比较基因组分析中常见的问题。建议先统一 ID，并保留旧 ID 到新 ID 的映射。

2.4.1 规范 FASTA ID

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/normalize_fasta.py \
  data/raw_sequences.fa \
  results/sequences.normalized.fa \
  --input-key first \
  --replacement-character _ \
  --max-id-length 60 \
  --mapping-file results/sequences.id_mapping.tsv

工作流会写出规范化 FASTA，以及包含 source_id 和 normalized_id 的映射表。需要显式替换残基字符时，可以重复使用 --mask：

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/normalize_fasta.py \
  data/raw_sequences.fa \
  results/sequences.masked.fa \
  --mask '?=N' \
  --mask 'X=N'

替换前应确认序列类型和缺失字符的真实含义。

2.4.2 查看序列统计

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/fasta_tool.py stats \
  results/sequences.normalized.fa \
  results/sequences.stats.tsv

默认统计 ID、序列长度、去除 gap 后长度、GC 比例、N 数量和 X 数量，可用于快速寻找异常序列。

2.4.3 提取区间和反向互补序列

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/fasta_tool.py subseq \
  data/genome.fa \
  results/region.revcomp.fa \
  --record chromosome_1 \
  --start 1001 \
  --end 2000 \
  --reverse-complement

这里的坐标是 1-based closed。和 BED 等 0-based half-open 文件联用时，必须先转换或确认坐标定义。

2.5 从 CDS 到 codon alignment

2.5.1 翻译 CDS

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/translate_cds.py \
  data/ortholog.cds.fa \
  results/ortholog.aa.fa \
  --frame 0 \
  --incomplete error \
  --internal-stop error \
  --terminal-stop trim

--frame 可取 0、1 或 2。正式分析建议先使用严格策略检查不完整密码子和内部终止密码子，而不是为了继续运行而直接忽略异常。

2.5.2 蛋白序列比对

使用 MAFFT：

bash /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/run_sequence_alignment.sh \
  results/ortholog.aa.fa \
  results/ortholog.aa.aln.fa \
  mafft \
  16 \
  --auto

同一脚本也支持 muscle3 和 muscle5。脚本会检查输入文件、线程数和外部程序，但 alignment 方法及参数仍需要根据数据特点选择。

2.5.3 映射回密码子 alignment

python3 /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/sequences/protein_alignment_to_codon.py \
  results/ortholog.aa.aln.fa \
  data/ortholog.cds.fa \
  results/ortholog.codon.aln.fa \
  --missing error \
  --trailing-cds error

蛋白 alignment 和 CDS 的 ID 必须一致。进入 PAML 或其他 codon 模型前，应检查序列长度、reading frame、内部终止密码子和异常 gap 区域。

2.6 系统树裁剪、整理和重根

树中的 tip 名称应与 FASTA、分组表和性状表完全一致，包括大小写、下划线和版本后缀。

2.6.1 提取目标物种子树

准备每行一个 tip 的 tips.txt：

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/extract_subtree.R \
  data/b10k_tree.nwk \
  data/tips.txt \
  results/b10k_subset.nwk

缺失 tip 的策略可以选择 error、warn 或 ignore。正式分析建议先使用默认的 error。

2.6.2 处理多分叉和节点标签

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/resolve_polytomies.R \
  results/b10k_subset.nwk \
  results/b10k_subset.resolved.nwk \
  auto \
  FALSE \
  - \
  - \
  _ \
  TRUE \
  Node \
  FALSE \
  FALSE \
  FALSE \
  results/b10k_subset.tips.txt

如果使用随机多分叉解析，应提供固定 seed，并明确说明随机解析不代表有数据支持的拓扑关系。

2.6.3 使用外群重根

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/reroot_tree.R \
  results/b10k_subset.resolved.nwk \
  results/b10k_subset.rooted.nwk \
  outgroup \
  data/outgroups.txt \
  FALSE \
  FALSE \
  -

使用 midpoint rooting：

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/reroot_tree.R \
  results/b10k_subset.resolved.nwk \
  results/b10k_subset.midpoint.nwk \
  midpoint \
  - \
  FALSE \
  FALSE \
  -

midpoint 模式需要 phangorn，但 midpoint rooting 不能替代真实的外群证据。

2.7 绘制 B10K 363 类群系统树

分组表默认包含两列：

Sample,Clade
species_A,Palaeognathae
species_B,Galloanseres
species_C,Australaves
species_D,Afroaves

新建 run_b10k_tree.R：

source("/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/src/R/phylogeny/clade_grouping.R")
source("/share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/B10K/363_clade_tree.R")

run_b10k_363_clade_tree(
  tree_file = "results/b10k_subset.rooted.nwk",
  group_file = "data/b10k_groups.csv",
  output_file = "results/b10k_clades.pdf",
  width = 12,
  height = 12,
  show_legend = FALSE
)

运行：

Rscript run_b10k_tree.R

当前预设类群包括：

Palaeognathae
Galloanseres
Mirandornithes
Columbimorphae
Otidimorphae
Opisthocomiformes
Cursorimorphae
Strisores
Phaethontimorphae
Aequornithes
Afroaves
Australaves

CSV 中的 Clade 拼写应与预设名称一致；未分组枝条默认显示为黑色。

2.8 Trait：整理和填补混合型性状表

示例性状表：

Species	BodyMass	WingLength	Habitat
species_A	125.3	88.1	Forest
species_B	NA	91.7	Grassland
species_C	84.6	NA	Forest

2.8.1 简单填补用于流程检查

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/trait/impute_trait_table.R \
  data/traits.tsv \
  results/traits_simple \
  Species \
  BodyMass,WingLength,Habitat \
  simple

简单填补默认对数值列使用 median，对分类列使用最常见水平。它适合检查接口，但不能自动替代正式研究中的缺失值模型。

主要输出包括：

results/traits_simple.imputed.csv
results/traits_simple.missingness.csv
results/traits_simple.replacements.csv
results/traits_simple.parameters.csv

2.8.2 使用 missForest

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/trait/impute_trait_table.R \
  data/traits.tsv \
  results/traits_missforest \
  Species \
  BodyMass,WingLength,Habitat \
  missforest \
  data/factor_levels.csv \
  BodyMass,WingLength \
  median \
  error \
  error \
  error \
  1234

factor_levels.csv 示例：

column,level,order
Habitat,Forest,1
Habitat,Grassland,2
Habitat,Wetland,3

应记录参与填补的列、factor level、随机 seed、缺失比例、OOB error，以及是否进行了 log10 转换。

2.9 连续性状与系统树分析

2.9.1 祖先性状重建

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/analyze_continuous_trait.R \
  results/b10k_subset.rooted.nwk \
  results/traits_missforest.imputed.csv \
  results/body_mass \
  Species \
  BodyMass \
  - \
  asr \
  error \
  error \
  3

脚本会先按物种 ID 对齐树和性状表，并输出对齐数据、tip 和 internal node 的重建结果及参数记录。

2.9.2 Phylogenetic ANOVA

性状表需额外包含分组列，例如 Clade：

Rscript /share/home/zhanglab/user/chenguangji/Usrbin/local/Codex_bin/workflows/phylogeny/analyze_continuous_trait.R \
  results/b10k_subset.rooted.nwk \
  results/traits_missforest.imputed.csv \
  results/body_mass_by_clade \
  Species \
  BodyMass \
  Clade \
  both \
  error \
  error \
  3 \
  1000 \
  TRUE \
  holm \
  1234

这个例子同时运行 ancestral reconstruction 和 phylogenetic ANOVA，使用 1,000 次模拟、post-hoc comparison、Holm 校正和固定 seed。

模型运行成功后仍需检查各组样本量、性状分布、异常值、树的枝长与根、缺失值处理以及分组的生物学合理性。

3. 并行化处理：paralleltask-ng

B10K 数据处理中经常需要执行几十到上千个彼此独立的任务，例如逐基因比对、逐基因树构建、逐物种统计或批量格式转换。paralleltask-ng 可以把一组命令拆分成独立子任务，并在本地或 SGE、PBS、Slurm、LSF 等调度系统中维持指定数量的并行任务。

项目地址：

• ChenGuangji/ParallelTask

3.1 安装

可以直接从 GitHub 安装当前版本：

python3 -m pip install "paralleltask-ng @ git+https://github.com/ChenGuangji/ParallelTask.git"

安装后先确认命令可用：

paralleltask-ng -h

如果集群不允许写系统 Python 环境，可以在虚拟环境或 Conda 环境中安装。安装命令需要系统中存在 git，并能够访问 GitHub。

3.2 准备任务文件

paralleltask-ng 的输入是一个命令集合文件。默认情况下，每一行形成一个独立子任务。例如创建 test_all.sh：

python3 workflow_A.py sample_001 > logs/sample_001.log 2>&1
python3 workflow_A.py sample_002 > logs/sample_002.log 2>&1
python3 workflow_A.py sample_003 > logs/sample_003.log 2>&1

这些命令应当彼此独立，并写入不同的输出文件。不要让多个任务同时覆盖同一个临时文件或结果文件。正式提交前建议先检查：

mkdir -p logs
bash -n test_all.sh
head -n 3 test_all.sh

还可以先复制少量命令到一个测试文件中运行，确认软件环境、输入路径和输出权限都正确，再提交完整任务列表。

当一个子任务需要由多行命令组成时，可以通过 -l 或 --lines 指定每个子任务包含的行数。默认值为 1。

3.3 在 Slurm 上运行

示例命令：

/share/home/zhanglab/user/chenguangji/Usrbin/anaconda3/bin/paralleltask-ng \
  -t slurm \
  -m 10 \
  -i 3 \
  --global-maxjob 10 \
  --submit-batch-size 5 \
  --submit-interval 3 \
  --submit-retry \
  test_all.sh

等价的单行写法是：

paralleltask-ng -t slurm -m 10 -i 3 --global-maxjob 10 --submit-batch-size 5 --submit-interval 3 --submit-retry test_all.sh

参数含义如下：

• -t slurm：使用 Slurm 的 sbatch、squeue 和相关控制命令；
• -m 10：本次 paralleltask-ng 进程最多同时维护 10 个任务；
• -i 3：每隔 3 秒检查任务状态或进行后续提交；
• --global-maxjob 10：根据 squeue 结果，将当前用户在 Slurm 中的全部任务数限制在 10 个以内，而不只计算本次运行提交的任务；
• --submit-batch-size 5：每成功提交 5 个任务后暂停一次，避免短时间向调度器发送过多请求；
• --submit-interval 3：每批提交之间暂停 3 秒；若不设置，则使用 -i/--interval 的值；
• --submit-retry：当 sbatch 提交失败或无法解析任务 ID 时，间隔一段时间后持续重试，而不是立即退出；
• test_all.sh：需要执行的命令集合文件。

这里同时使用 -m 10 和 --global-maxjob 10。前者限制当前这次运行的并发数，后者限制当前用户在整个 Slurm 队列中的任务总数。例如用户已经有 6 个其他任务时，这次运行最多只能再占用约 4 个名额，直到已有任务结束。

3.4 CPU 和内存

若不显式设置，当前 CLI 默认每个子任务申请 1 个 CPU 和 3G 内存。可以使用下面的参数覆盖：

paralleltask-ng \
  -t slurm \
  -p 4 \
  -M 16G \
  -m 10 \
  test_all.sh

其中：

• -p 4：每个子任务申请 4 个 CPU；
• -M 16G：每个子任务申请 16G 内存。

资源参数应与命令实际使用方式一致。为单线程程序盲目申请多个 CPU 不会自动提升速度；内存申请过小可能导致任务被调度器终止，申请过大则可能延长排队时间。

3.5 运行目录、日志和失败任务

运行 test_all.sh 后，程序会创建相应的工作目录：

test_all.sh.work/

其中保存拆分后的子任务脚本及相关运行文件。程序自身也会写日志；可以通过 --log FILE 指定日志文件。建议同时让每条任务命令把标准输出和标准错误写到明确的文件中，这样更容易定位具体失败的样本或基因。

程序支持跳过已经成功的子任务并重新执行未完成任务。-r/--rerun 控制失败子任务的重新运行轮数，默认最多重新运行 3 轮。重启 Slurm 流程时，程序还会查询活动任务，并尝试根据子任务脚本路径恢复仍在排队或运行的任务，避免重复提交。

3.6 使用 --submit-retry 的注意事项

--submit-retry 会对提交失败持续重试，适合调度器暂时繁忙、网络瞬时异常或提交服务短暂不可用的情况。但如果失败源于长期存在的错误，例如：

• Slurm 参数或集群配置不正确；
• 分区、账户或资源申请无效；
• 输入脚本不存在；
• 用户达到无法恢复的配额或权限限制；

程序可能一直重复尝试。看到持续失败时，应检查日志和 sbatch 错误信息，并在必要时使用 Ctrl+C 停止主程序，修正问题后再重新运行。

主程序收到终止信号时会尝试停止已经提交的子任务，但在大型任务批次中仍建议随后使用 squeue -u "$USER" 确认队列状态。

3.7 推荐的使用步骤

一个相对稳妥的提交顺序是：

1. 检查任务文件语法与输入路径；
2. 用 2—3 个任务进行小规模测试；
3. 确认单个任务的 CPU、内存和运行时间；
4. 设置合理的 -m 与 --global-maxjob；
5. 提交完整任务文件；
6. 查看主日志、子任务日志和 squeue；
7. 对失败任务先判断原因，再决定是否重跑。

并行化只负责调度彼此独立的命令，不会自动判断分析步骤之间的依赖关系。存在前后依赖的命令应合并为同一个子任务，或拆成多个阶段，在上一阶段全部完成并验证后再启动下一阶段。

Next to update

1. Pixi related pipeline