ART 使用笔记 - Lihuax's Blog

1. ART 是什么#

ART，全称通常指 Artificial Read Transcripter，是一个用于模拟二代或三代测序 reads 的工具包。

它的核心作用是：

1
输入参考序列 FASTA
2
  ↓
3
按照指定测序平台、read length、覆盖度、fragment size、错误模型
4
  ↓
5
生成模拟 FASTQ reads

ART 常用于：

模拟 Illumina / 454 / SOLiD reads
构建算法 benchmark 数据
测试 read mapping 流程
测试 SNP / indel / SV caller
测试自定义变异检测算法
模拟不同测序深度下的检测能力
模拟不同 VAF 条件下的突变检出能力

需要注意：

ART 不是突变生成器。

它不会自动生成 SNP、indel、SV、DSB 修复结果、易位、倒位、L1 插入等事件。
这些结构变化需要提前构造到输入 FASTA 中，ART 只负责从这些序列上模拟 reads。

2. ART 在模拟流程中的位置#

一个标准 benchmark 流程通常是：

1
reference genome
2
  ↓
3
人为构造变异 / 使用 simulator 构造 alt genome
4
  ↓
5
得到 mutated genome / alt genome / junction sequence
6
  ↓
7
ART 从这些序列上生成 FASTQ reads
8
  ↓
9
reads 比对回 reference genome
10
  ↓
11
运行 caller / 自己的算法
12
  ↓
13
与 truth set 比较

例如：

1
原始参考基因组
2
  ↓
3
引入 NHEJ insertion / MMEJ deletion / BND translocation
4
  ↓
5
生成 alt.fa
6
  ↓
7
ART 生成 alt reads
8
  ↓
9
bowtie2 / bwa 比对回原始 reference
10
  ↓
11
Manta / Delly / Lumpy / 自己的算法检测

因此：

1
ART 负责：模拟测序 reads
2
不负责：生成生物学变异事件

3. ART 支持的平台#

ART 包含多个程序：

1
art_illumina
2
art_454
3
art_SOLiD

其中最常用的是：

1
art_illumina

因为大多数短读长测序数据都是 Illumina paired-end reads。

4. 安装 ART#

4.1 使用 conda 安装#

推荐使用 conda / mamba 安装：

1
conda install -c conda-forge -c bioconda art

或者：

1
mamba install -c conda-forge -c bioconda art

检查是否安装成功：

1
which art_illumina
2
art_illumina

如果能看到 help / usage 信息，说明安装成功。

5. ART 输入和输出#

5.1 输入文件#

ART 的主要输入是 FASTA 文件：

1
reference.fa
2
alt.fa
3
target_regions.fa
4
junction.fa

FASTA 示例：

1
>chr_mock
2
ACTGACTGACTGACTGACTGACTGACTGACTG

可以输入：

1
1. 完整基因组 FASTA
2
2. 某几条染色体 FASTA
3
3. 小型 mock genome FASTA
4
4. target regions FASTA
5
5. junction sequence FASTA
6
6. mutated genome FASTA

5.2 输出文件#

以 Illumina paired-end 为例：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i reference.fa \
7
  -l 150 \
8
  -f 30 \
9
  -m 350 \
10
  -s 50 \
11
  -o sim/sample_

可能生成：

1
sample_1.fq
2
sample_2.fq
3
sample_.sam

含义：

1
sample_1.fq    R1 reads
2
sample_2.fq    R2 reads
3
sample_.sam    ART 生成 reads 时相对于输入 FASTA 的模拟比对结果

注意：

如果输入的是 alt.fa，则 ART 输出的 SAM 坐标是相对于 alt.fa 的，不是相对于原始 reference 的。

实际做 caller benchmark 时，通常应该把 ART 生成的 FASTQ 重新比对回原始 reference。

6. `art_illumina` 常用参数#

6.1 核心参数#

参数	含义
`-i`	输入 FASTA 文件
`-o`	输出文件前缀
`-l`	read length
`-f`	fold coverage，模拟覆盖度
`-c`	指定 reads 或 read pairs 数量
`-p`	paired-end 模式
`-m`	paired-end fragment length 均值
`-s`	paired-end fragment length 标准差
`-ss`	Illumina 测序平台 profile
`-sam`	输出 SAM 文件
`-na`	不输出 ALN 文件
`-rs`	random seed，固定随机种子
`-mp`	mate-pair 模式
`-amp`	amplicon 模式

7. `-ss` 测序平台 profile#

-ss 用来指定 Illumina 平台的质量分布模型。

常见值：

参数值	平台
`GA1`	Genome Analyzer I
`GA2`	Genome Analyzer II
`HS10`	HiSeq 1000
`HS20`	HiSeq 2000
`HS25`	HiSeq 2500
`MS`	MiSeq

常用选择：

1
-ss HS25

对于 150 bp paired-end reads，通常可以先用 HS25。

8. Single-end reads 模拟#

8.1 基本命令#

1
art_illumina \
2
  -ss HS25 \
3
  -i reference.fa \
4
  -l 150 \
5
  -f 30 \
6
  -rs 42 \
7
  -o sim/single_end_

参数解释：

1
-ss HS25        使用 HiSeq 2500 profile
2
-i reference.fa 输入 FASTA
3
-l 150          read length = 150 bp
4
-f 30           覆盖度 = 30x
5
-rs 42          固定随机种子
6
-o              输出前缀

输出：

1
single_end_.fq

如果加上 -sam：

1
art_illumina \
2
  -ss HS25 \
3
  -sam \
4
  -na \
5
  -i reference.fa \
6
  -l 150 \
7
  -f 30 \
8
  -rs 42 \
9
  -o sim/single_end_

输出：

1
single_end_.fq
2
single_end_.sam

9. Paired-end reads 模拟#

9.1 基本命令#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i reference.fa \
7
  -l 150 \
8
  -f 30 \
9
  -m 350 \
10
  -s 50 \
11
  -rs 42 \
12
  -o sim/paired_end_

参数解释：

1
-p        paired-end mode
2
-l 150    read length = 150 bp
3
-f 30     coverage = 30x
4
-m 350    mean fragment length = 350 bp
5
-s 50     fragment length standard deviation = 50 bp
6
-rs 42    random seed

输出：

1
paired_end_1.fq
2
paired_end_2.fq
3
paired_end_.sam

10. 覆盖度 `-f` 的含义#

-f 表示 fold coverage。

例如：

1
-f 30

表示模拟约 30x 覆盖度。

paired-end 覆盖度近似计算公式：

1
coverage ≈ read_pairs × 2 × read_length / genome_length

例如：

1
genome length = 100,000 bp
2
read length = 150 bp
3
read pairs = 10,000
4

5
coverage ≈ 10,000 × 2 × 150 / 100,000
6
         = 30x

所以对于 100 kb genome，150 bp paired-end，10,000 对 reads 大约就是 30x。

11. 使用 `-c` 指定 read pairs 数量#

-c 可以直接指定 reads 或 read pairs 的数量。

paired-end 示例：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i reference.fa \
7
  -l 150 \
8
  -c 10000 \
9
  -m 350 \
10
  -s 50 \
11
  -rs 42 \
12
  -o sim/sample_

注意：

1
-f 和 -c 不应该同时使用

通常：

1
-f 适合按覆盖度模拟
2
-c 适合精确控制 reads 数量

12. `-m` 和 `-s` 的含义#

-m 是 fragment length 的均值。
-s 是 fragment length 的标准差。

例如：

1
-m 350 -s 50

表示 DNA fragment size 大多数集中在 350 bp 左右。

注意：

这里的 fragment length 不是插入突变长度，而是测序文库中 DNA 片段的长度。

对于 150 bp paired-end WGS：

1
read length: 150 bp
2
mean fragment length: 300–500 bp
3
standard deviation: 30–100 bp

常用：

1
-m 350 -s 50

对于 CUT&Tag-like 或 target-enriched 数据，可以考虑：

1
-m 250 -s 80

具体最好参考真实数据的 fragment size distribution。

13. 固定随机种子 `-rs`#

-rs 用于设置 random seed。

例如：

1
-rs 42

优点：

1
1. 结果可复现
2
2. debug 更方便
3
3. 多次运行不会因为随机性产生不同结果

建议 benchmark 时一定加上 -rs。

14. 输出 SAM 的作用#

加上：

1
-sam

ART 会输出 SAM 文件。

例如：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i alt.fa \
7
  -l 150 \
8
  -f 30 \
9
  -m 350 \
10
  -s 50 \
11
  -rs 42 \
12
  -o sim/alt_

输出：

1
alt_1.fq
2
alt_2.fq
3
alt_.sam

SAM 用途：

1
1. 检查 reads 来源
2
2. 检查模拟覆盖度
3
3. 检查 fragment size
4
4. 检查 reads 是否覆盖目标 junction

但是要注意：

1
ART 的 SAM 坐标是相对于输入 FASTA 的

如果输入是 alt.fa，SAM 坐标就是 alt.fa 坐标。

实际变异检测时，应该重新比对 FASTQ 到原始 reference。

15. 常规 WGS-like 数据模拟#

15.1 模拟 control 样本#

1
mkdir -p sim/control
2

3
art_illumina \
4
  -ss HS25 \
5
  -p \
6
  -sam \
7
  -na \
8
  -i ref/control.fa \
9
  -l 150 \
10
  -f 30 \
11
  -m 350 \
12
  -s 50 \
13
  -rs 1001 \
14
  -o sim/control/control_

15.2 模拟 treated / tumor / alt 样本#

1
mkdir -p sim/treated
2

3
art_illumina \
4
  -ss HS25 \
5
  -p \
6
  -sam \
7
  -na \
8
  -i ref/treated_alt.fa \
9
  -l 150 \
10
  -f 30 \
11
  -m 350 \
12
  -s 50 \
13
  -rs 1002 \
14
  -o sim/treated/treated_

15.3 比对回原始 reference#

1
mkdir -p bam logs
2

3
bowtie2 \
4
  -x ref/bowtie2/reference \
5
  -1 sim/control/control_1.fq \
6
  -2 sim/control/control_2.fq \
7
  -p 8 \
8
  2> logs/control.bowtie2.log \
9
| samtools sort -@ 8 -o bam/control.sorted.bam
10

11
samtools index bam/control.sorted.bam

1
bowtie2 \
2
  -x ref/bowtie2/reference \
3
  -1 sim/treated/treated_1.fq \
4
  -2 sim/treated/treated_2.fq \
5
  -p 8 \
6
  2> logs/treated.bowtie2.log \
7
| samtools sort -@ 8 -o bam/treated.sorted.bam
8

9
samtools index bam/treated.sorted.bam

16. 模拟不同 VAF#

ART 本身没有直接的 --vaf 参数。
如果想模拟 VAF，需要分别从 reference allele 和 alt allele 生成 reads，然后混合。

16.1 模拟 50% VAF#

假设总深度 30x，可以设置：

1
reference allele: 15x
2
alternate allele: 15x

命令：

1
mkdir -p sim/vaf50

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i ref/reference.fa \
6
  -l 150 \
7
  -f 15 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 1 \
11
  -o sim/vaf50/ref_

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i ref/alt.fa \
6
  -l 150 \
7
  -f 15 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 2 \
11
  -o sim/vaf50/alt_

混合：

1
cat sim/vaf50/ref_1.fq sim/vaf50/alt_1.fq > sim/vaf50/sample_1.fq
2
cat sim/vaf50/ref_2.fq sim/vaf50/alt_2.fq > sim/vaf50/sample_2.fq

近似 VAF：

1
15 / (15 + 15) = 50%

16.2 模拟 30% VAF#

总深度 30x：

1
reference allele: 21x
2
alternate allele: 9x

1
9 / (21 + 9) = 30%

16.3 模拟 15% VAF#

总深度 30x：

1
reference allele: 25.5x
2
alternate allele: 4.5x

实际命令中可以近似为：

1
reference allele: 25x
2
alternate allele: 5x

近似 VAF：

1
5 / (25 + 5) = 16.7%

如果想更精确，可以用 -c 控制 read pairs 数量。

17. 使用 `-c` 精确模拟 VAF#

假设：

1
genome length = 100,000 bp
2
read length = 150 bp
3
total depth = 30x

总 read pairs：

1
read_pairs = coverage × genome_length / (2 × read_length)
2
           = 30 × 100000 / 300
3
           = 10000

如果 VAF = 30%：

1
alt read pairs = 10000 × 0.3 = 3000
2
ref read pairs = 10000 × 0.7 = 7000

命令：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i ref/reference.fa \
6
  -l 150 \
7
  -c 7000 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 1 \
11
  -o sim/vaf30/ref_

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i ref/alt.fa \
6
  -l 150 \
7
  -c 3000 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 2 \
11
  -o sim/vaf30/alt_

混合：

1
cat sim/vaf30/ref_1.fq sim/vaf30/alt_1.fq > sim/vaf30/sample_1.fq
2
cat sim/vaf30/ref_2.fq sim/vaf30/alt_2.fq > sim/vaf30/sample_2.fq

18. 模拟 target-enriched / CUT&Tag-like 数据#

ART 默认模拟的是 WGS-like reads。
如果输入完整 genome，输出就类似 WGS。

如果要模拟富集型数据，可以只从目标区域生成 reads。

18.1 准备 target BED#

示例：

1
chr1    1000000    1002000
2
chr1    3000000    3002500
3
chr2    5000000    5003000

18.2 给 target region 加 flank#

1
bedtools slop \
2
  -i targets.bed \
3
  -g genome.sizes \
4
  -b 1000 \
5
> targets.flank1k.bed

18.3 提取 target FASTA#

1
bedtools getfasta \
2
  -fi reference.fa \
3
  -bed targets.flank1k.bed \
4
  -fo targets.flank1k.fa

18.4 从 target FASTA 生成 reads#

1
mkdir -p sim/enriched
2

3
art_illumina \
4
  -ss HS25 \
5
  -p \
6
  -sam \
7
  -na \
8
  -i targets.flank1k.fa \
9
  -l 150 \
10
  -f 100 \
11
  -m 250 \
12
  -s 80 \
13
  -rs 42 \
14
  -o sim/enriched/target_

这样生成的 reads 主要来自目标区域。
将这些 reads 比对回完整 reference 后，就会表现为 target-enriched 信号。

19. 模拟富集区域加背景噪音#

真实 CUT&Tag / capture sequencing 数据通常不是只有 target reads，还会有 genome-wide background。

可以混合：

1
target-enriched reads
2
+
3
low-depth WGS background reads

19.1 生成 target reads#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i targets.flank1k.fa \
6
  -l 150 \
7
  -f 100 \
8
  -m 250 \
9
  -s 80 \
10
  -rs 1 \
11
  -o sim/enriched/target_

19.2 生成 low-depth background reads#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i reference.fa \
6
  -l 150 \
7
  -f 1 \
8
  -m 350 \
9
  -s 100 \
10
  -rs 2 \
11
  -o sim/enriched/background_

19.3 混合 FASTQ#

1
cat sim/enriched/target_1.fq sim/enriched/background_1.fq > sim/enriched/sample_1.fq
2
cat sim/enriched/target_2.fq sim/enriched/background_2.fq > sim/enriched/sample_2.fq

这种数据更接近：

1
局部高覆盖
2
全基因组低背景

20. 模拟不同 peak 强度#

如果想模拟不同强度的 peak，可以给不同区域设置不同覆盖度。

例如：

1
strong peaks: 200x
2
medium peaks: 80x
3
weak peaks: 30x
4
background: 1x

命令示例：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i peaks_strong.fa \
6
  -l 150 \
7
  -f 200 \
8
  -m 250 \
9
  -s 80 \
10
  -rs 11 \
11
  -o sim/peaks/strong_

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i peaks_weak.fa \
6
  -l 150 \
7
  -f 30 \
8
  -m 250 \
9
  -s 80 \
10
  -rs 12 \
11
  -o sim/peaks/weak_

混合：

1
cat sim/peaks/strong_1.fq sim/peaks/weak_1.fq > sim/peaks/sample_1.fq
2
cat sim/peaks/strong_2.fq sim/peaks/weak_2.fq > sim/peaks/sample_2.fq

21. 模拟 DSB 修复结果的推荐策略#

ART 不直接生成 DSB 修复结果。
需要先自己构造修复后的序列。

常见设计：

1
1. NULL
2
   无变化
3

4
2. NHEJ_INS
5
   DSB 位点插入少量随机碱基
6

7
3. NHEJ_DEL
8
   DSB 位点附近发生小缺失
9

10
4. MMEJ_DEL
11
   LmXmR -> LmR
12
   利用 microhomology 修复，删除中间 X 区域
13

14
5. BND_INS
15
   两个远距离断点异常连接，中间可能有短插入
16

17
6. INTER_TRANSLOCATION
18
   不同染色体之间异常连接
19

20
7. INTRA_REARRANGEMENT
21
   同一染色体内远距离异常连接

21.1 MMEJ 示例#

修复前：

1
Left flank + microhomology + deleted segment + microhomology + Right flank

修复后：

1
Left flank + one copy of microhomology + Right flank

例如：

1
Before:
2
AAAAAA CCGT TTTTTTTTTT CCGT GGGGGG
3

4
After:
5
AAAAAA CCGT GGGGGG

然后将修复后的序列保存为：

1
mmej_alt.fa

再用 ART：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i mmej_alt.fa \
6
  -l 150 \
7
  -f 100 \
8
  -m 250 \
9
  -s 50 \
10
  -rs 101 \
11
  -o sim/mmej/mmej_

21.2 NHEJ insertion 示例#

修复前：

1
AAAAAA | DSB | GGGGGG

修复后：

1
AAAAAA TCA GGGGGG

其中 TCA 是随机插入序列。

保存为：

1
nhej_ins_alt.fa

再用 ART 生成 reads。

21.3 BND / translocation 示例#

假设有两个断点：

1
chrA:posA
2
chrB:posB

可以构造 junction sequence：

1
chrA_left_flank + inserted_sequence + chrB_right_flank

例如：

1
>junction_chrA_chrB
2
AAAAAAATTTTTTCCCCCCGGGGGG

然后：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i junction_chrA_chrB.fa \
6
  -l 150 \
7
  -f 100 \
8
  -m 250 \
9
  -s 50 \
10
  -rs 201 \
11
  -o sim/bnd/junction_

这种方式适合快速测试 split read / soft clipping / discordant pair 识别逻辑。

如果要更接近真实 WGS，应构造完整 derivative chromosome：

1
derivative_A = chrA_left + insert + chrB_right
2
derivative_B = chrB_left + insert + chrA_right

再从完整 alt genome 上用 ART 生成 reads。

22. 推荐的 benchmark 数据层级#

22.1 Level 1：最小 debug 数据#

目标：

1
确认算法逻辑能跑通

特点：

1
1. 基因组很短
2
2. 只有一个事件
3
3. 覆盖度很高
4
4. 没有复杂背景

示例：

1
one_mmej.fa
2
one_nhej_ins.fa
3
one_bnd_junction.fa

ART：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i one_event.fa \
6
  -l 150 \
7
  -f 100 \
8
  -m 250 \
9
  -s 50 \
10
  -rs 1 \
11
  -o sim/debug/one_event_

22.2 Level 2：WGS-like benchmark#

目标：

1
测试传统 SV caller 和自定义算法

特点：

1
1. 使用小型 genome
2
2. 包含多个 truth events
3
3. 模拟 30x WGS
4
4. 比对回原始 reference

示例：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i control.fa \
6
  -l 150 \
7
  -f 30 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 1 \
11
  -o sim/wgs/control_

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i alt.fa \
6
  -l 150 \
7
  -f 30 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 2 \
11
  -o sim/wgs/alt_

22.3 Level 3：target-enriched benchmark#

目标：

1
模拟 CUT&Tag-like 或 capture-like 数据

特点：

1
1. 目标区域高覆盖
2
2. 全基因组低背景
3
3. 可设置 control / treated 差异
4
4. 可加入 alt junction reads

结构：

1
control = target_ref_reads + low_background_reads
2

3
treated = target_ref_reads
4
        + target_alt_reads
5
        + low_background_reads

23. ART 与比对工具结合#

23.1 使用 bowtie2#

1
bowtie2 \
2
  -x ref/bowtie2/reference \
3
  -1 sim/sample_1.fq \
4
  -2 sim/sample_2.fq \
5
  -p 8 \
6
  2> logs/sample.bowtie2.log \
7
| samtools sort -@ 8 -o bam/sample.sorted.bam
8

9
samtools index bam/sample.sorted.bam

23.2 使用 bwa mem#

1
bwa mem \
2
  -t 8 \
3
  reference.fa \
4
  sim/sample_1.fq \
5
  sim/sample_2.fq \
6
| samtools sort -@ 8 -o bam/sample.sorted.bam
7

8
samtools index bam/sample.sorted.bam

24. ART 与 SV caller benchmark#

常规流程：

1
1. 准备 reference.fa
2
2. 构造 alt.fa
3
3. 保存 truth.vcf / truth.bed / truth.tsv
4
4. ART 从 reference.fa 生成 control reads
5
5. ART 从 alt.fa 生成 treated/tumor reads
6
6. reads 比对回 reference.fa
7
7. 运行 SV caller
8
8. caller result 与 truth set 比较

示例：

1
reference.fa
2
alt.fa
3
truth.tsv
4

5
sim/control_1.fq
6
sim/control_2.fq
7
sim/alt_1.fq
8
sim/alt_2.fq
9

10
bam/control.sorted.bam
11
bam/alt.sorted.bam
12

13
sv/manta/
14
sv/delly/
15
sv/lumpy/

25. ART 常见使用模板#

25.1 150 bp PE，30x WGS#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i reference.fa \
7
  -l 150 \
8
  -f 30 \
9
  -m 350 \
10
  -s 50 \
11
  -rs 42 \
12
  -o sim/sample_

25.2 150 bp PE，100x target sequencing#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i targets.fa \
7
  -l 150 \
8
  -f 100 \
9
  -m 250 \
10
  -s 80 \
11
  -rs 42 \
12
  -o sim/target_

25.3 指定 read pairs 数量#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -sam \
5
  -na \
6
  -i reference.fa \
7
  -l 150 \
8
  -c 10000 \
9
  -m 350 \
10
  -s 50 \
11
  -rs 42 \
12
  -o sim/sample_

25.4 不输出 SAM#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i reference.fa \
6
  -l 150 \
7
  -f 30 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 42 \
11
  -o sim/sample_

25.5 模拟 MiSeq#

1
art_illumina \
2
  -ss MS \
3
  -p \
4
  -na \
5
  -i reference.fa \
6
  -l 250 \
7
  -f 30 \
8
  -m 500 \
9
  -s 100 \
10
  -rs 42 \
11
  -o sim/miseq_

26. 454 reads 模拟#

454 平台现在较少使用，但 ART 仍支持。

26.1 454 single-end#

1
art_454 reference.fa sim/454_single_ 20

含义：

1
reference.fa      输入 FASTA
2
sim/454_single_   输出前缀
3
20                覆盖度 20x

26.2 454 paired-end#

1
art_454 reference.fa sim/454_pair_ 20 500 20

含义：

1
20      覆盖度
2
500     mean fragment length
3
20      fragment length standard deviation

454 的错误模型更偏向 homopolymer 相关 indel error。

27. SOLiD reads 模拟#

SOLiD 现在也较少使用。

27.1 SOLiD single-end#

1
art_SOLiD reference.fa sim/solid_single_ 32 10

含义：

1
32      read length
2
10      coverage

27.2 SOLiD paired-end#

1
art_SOLiD reference.fa sim/solid_pair_ 25 10 500 20

含义：

1
25      read length
2
10      coverage
3
500     mean fragment length
4
20      fragment length standard deviation

28. 常见问题#

28.1 ART 可以直接模拟 SV 吗？#

不可以。

ART 只从输入 FASTA 上生成 reads。
如果要模拟 SV，需要先生成包含 SV 的 FASTA。

例如：

1
reference.fa
2
  ↓
3
手动修改 / 使用 RSVSim / 自己写脚本
4
  ↓
5
alt.fa
6
  ↓
7
ART

28.2 ART 可以直接模拟 DSB 修复吗？#

不可以。

需要先构造 DSB 修复后的序列，例如：

1
NHEJ insertion
2
MMEJ deletion
3
BND junction
4
translocation derivative chromosome

然后再用 ART 生成 reads。

28.3 ART 可以模拟 CUT&Tag 吗？#

不能完整模拟。

ART 可以模拟 reads，但不能真实模拟：

1
1. 抗体富集
2
2. pA-Tn5 插入偏好
3
3. peak shape
4
4. chromatin accessibility bias
5
5. ligase binding signal

但是可以通过以下方式近似：

1
1. 只从 target regions 生成 reads
2
2. 设置不同 target coverage
3
3. 混入低深度 WGS background
4
4. 加入特定 junction reads

28.4 ART 输出的 SAM 能直接用于 caller 吗？#

一般不建议。

28.5 为什么比对回 reference 后出现 soft clipping？#

如果 reads 来自 alt genome，而 alt genome 和 reference 有差异，那么比对回 reference 时，差异附近可能出现：

1
soft clipping
2
split alignment
3
discordant pair
4
indel CIGAR
5
low MAPQ
6
supplementary alignment

这正是测试 SV / DSB 检测算法时需要观察的信号。

29. 常见错误和排查#

29.1 输出目录不存在#

错误原因：

1
ART 不一定会自动创建输出目录

解决：

1
mkdir -p sim/output

再运行：

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -i reference.fa \
5
  -l 150 \
6
  -f 30 \
7
  -m 350 \
8
  -s 50 \
9
  -o sim/output/sample_

29.2 `-f` 和 `-c` 同时使用#

不要同时使用：

1
-f 30 -c 10000

二选一：

1
按覆盖度：使用 -f
2
按 read 数量：使用 -c

29.3 paired-end 忘记加 `-p`#

如果不加：

1
-p

ART 会按 single-end 模式生成 reads。

paired-end 必须加：

1
-p -m 350 -s 50

29.4 prefix 命名混乱#

建议输出 prefix 总是以 _ 结尾：

1
-o sim/sample_

这样输出更清晰：

1
sample_1.fq
2
sample_2.fq
3
sample_.sam

29.5 随机结果不可复现#

加上：

1
-rs 42

30. 推荐目录结构#

1
project/
2
├── ref/
3
│   ├── reference.fa
4
│   ├── reference.fa.fai
5
│   ├── alt.fa
6
│   ├── targets.bed
7
│   ├── targets.fa
8
│   └── bowtie2/
9
├── truth/
10
│   ├── truth.tsv
11
│   ├── truth.bed
12
│   └── truth.vcf
13
├── sim/
14
│   ├── control/
15
│   ├── treated/
16
│   ├── enriched/
17
│   └── vaf/
18
├── bam/
19
│   ├── control.sorted.bam
20
│   └── treated.sorted.bam
21
├── logs/
22
│   ├── control.bowtie2.log
23
│   └── treated.bowtie2.log
24
└── sv/
25
    ├── manta/
26
    ├── delly/
27
    └── lumpy/

31. 推荐实际工作流#

31.1 Step 1：准备 reference#

1
samtools faidx ref/reference.fa
2
cut -f1,2 ref/reference.fa.fai > ref/genome.sizes

31.2 Step 2：构造 alt genome#

可以使用：

1
1. 自己写脚本
2
2. RSVSim
3
3. SURVIVOR simSV
4
4. 自定义 sequence generator
5
5. 手动构造 junction FASTA

输出：

1
ref/alt.fa
2
truth/truth.tsv

31.3 Step 3：ART 生成 reads#

1
art_illumina \
2
  -ss HS25 \
3
  -p \
4
  -na \
5
  -i ref/alt.fa \
6
  -l 150 \
7
  -f 30 \
8
  -m 350 \
9
  -s 50 \
10
  -rs 42 \
11
  -o sim/alt_

31.4 Step 4：比对回原始 reference#

1
bowtie2 \
2
  -x ref/bowtie2/reference \
3
  -1 sim/alt_1.fq \
4
  -2 sim/alt_2.fq \
5
  -p 8 \
6
  2> logs/alt.bowtie2.log \
7
| samtools sort -@ 8 -o bam/alt.sorted.bam
8

9
samtools index bam/alt.sorted.bam

31.5 Step 5：运行 caller#

1
Manta
2
Delly
3
Lumpy
4
自定义算法

31.6 Step 6：与 truth 比较#

比较内容：

1
1. 是否检出事件
2
2. breakpoint 是否接近 truth
3
3. SVTYPE 是否正确
4
4. ALT junction 是否正确
5
5. VAF 估计是否合理
6
6. false positive 数量
7
7. false negative 数量

32. 对 DSB / SV benchmark 的建议#

如果目标是测试 DSB 修复或 SV caller，不建议一开始直接做复杂模拟。
推荐分阶段：

1
第一阶段：
2
单事件，高覆盖，无背景
3

4
第二阶段：
5
多个事件，WGS-like，30x
6

7
第三阶段：
8
加入不同 VAF
9

10
第四阶段：
11
加入 target enrichment
12

13
第五阶段：
14
加入 background noise
15

16
第六阶段：
17
加入真实数据中的 fragment size 和 coverage distribution

这样 debug 更容易。

33. 简要总结#

ART 的核心定位：

1
ART = reads simulator

它适合：

1
1. 从 reference / alt genome 生成 FASTQ
2
2. 模拟不同测序深度
3
3. 模拟不同 read length
4
4. 模拟 paired-end reads
5
5. 模拟不同 fragment size
6
6. 模拟不同 VAF 的混合样本
7
7. 构建 benchmark 数据

它不适合直接完成：

1
1. SV 生成
2
2. DSB 修复模拟
3
3. L1 插入模拟
4
4. CUT&Tag 生物学过程模拟
5
5. peak calling 模拟

最重要的原则：

1
先构造 truth sequence，再用 ART 生成 reads。

对于 DSB / SV / CUT&Tag-like benchmark，推荐逻辑是：

1
reference.fa
2
  ↓
3
alt.fa / junction.fa / targets.fa
4
  ↓
5
ART
6
  ↓
7
FASTQ
8
  ↓
9
比对回 reference.fa
10
  ↓
11
BAM
12
  ↓
13
caller / custom algorithm
14
  ↓
15
compare with truth