Rose Pangenome Analysis Scripts 详细说明
目录
- Genome Assembly — 基因组组装
- Genome Annotation — 基因组注释
- RNA-seq — 转录组分析
- Pangenome — 泛基因组构建
- Phylogenetic — 系统发育与群体遗传
1. Genome Assembly — 基因组组装
脚本: Genome_Assembly/assembly.sh
功能概述: 使用 PacBio HiFi 读长进行初始组装，再结合 Hi-C 数据进行挂载（scaffolding），最后评估组装质量。
1.1 初始组装 — Hifiasm
1
| hifiasm -o hifiasm_out -s 0.3 --h1 HiC_R1.fq.gz --h2 HiC_R2.fq.gz HiFi.fq.gz
|
| 参数 |
说明 |
-o hifiasm_out |
输出文件前缀 |
-s 0.3 |
设置去冗余 (purge duplicated) 时的覆盖率阈值比例为 0.3 |
--h1 HiC_R1.fq.gz |
Hi-C 数据 read 1 |
--h2 HiC_R2.fq.gz |
Hi-C 数据 read 2 |
HiFi.fq.gz |
PacBio HiFi 输入读长 |
工具: Hifiasm — 基于 HiFi 读长的基因组组装软件，可同时利用 Hi-C 数据进行组装。
1.2 Hi-C 挂载
Step 1: 构建 BWA 索引
工具: BWA — 为 Hi-C reads 比对构建参考基因组索引。
Step 2: 生成酶切位点文件
1
2
| python ~/software/juicer-1.6/misc/generate_site_positions.py DpnII genome genome.fa
awk 'BEGIN{OFS="\t"}{print $1, $NF}' genome_DpnII.txt > draft.genome.chrom.sizes
|
| 参数 |
说明 |
DpnII |
限制性内切酶类型 |
genome |
项目名称 |
genome.fa |
基因组 fasta 文件 |
说明: 生成 DpnII 酶切位点位置文件，并提取染色体大小信息。
Step 3: 运行 Juicer
1
| ~/software/juicer-1.6/CPU/juicer.sh -g hap1_hic -z genome.fa -s DpnII -y genome_DpnII.txt -p draft.genome.chrom.sizes -t 56 -S early
|
| 参数 |
说明 |
-g hap1_hic |
组名称 |
-z genome.fa |
参考基因组 |
-s DpnII |
限制性内切酶类型 |
-y genome_DpnII.txt |
酶切位点文件 |
-p draft.genome.chrom.sizes |
染色体大小文件 |
-t 56 |
使用 56 线程 |
-S early |
在 early 阶段停止（生成 merged_nodups.txt） |
工具: Juicer — Hi-C 数据处理流程。
Step 4: 3D-DNA 挂载与手动纠错
1
2
| ~/software/3d-dna/run-asm-pipeline.sh -r 0 reference/genome.fa aligned/merged_nodups.txt
~/software/3d-dna/run-asm-pipeline-post-review.sh -r genome.review.assembly ./reference/genome.fa aligned/merged_nodups.txt
|
| 参数 |
说明 |
-r 0 |
不自动校正 |
reference/genome.fa |
参考基因组路径 |
aligned/merged_nodups.txt |
Juicer 生成的比对文件 |
-r genome.review.assembly |
使用手动 review 后的纠错文件 |
工具: 3D-DNA — 利用 Hi-C 数据对基因组进行染色体级别挂载和纠错。
1.3 组装质量评估
BUSCO 评估
1
| busco -l ~/database/busco/eudicots_odb10 --offline -m genome -c 40 --out_path ./BUSCO_out -o genome.fa -i genome.fa
|
| 参数 |
说明 |
-l ~/database/busco/eudicots_odb10 |
使用真双子叶植物 (eudicots) BUSCO 数据库 |
--offline |
离线模式 |
-m genome |
评估模式为基因组模式 |
-c 40 |
使用 40 线程 |
--out_path ./BUSCO_out |
输出目录 |
-o genome.fa |
输出名称 |
-i genome.fa |
输入基因组 fasta |
工具: BUSCO — 基于单拷贝直系同源基因评估基因组组装/注释完整度。
KAT K-mer 分析
1
| kat comp -t 40 -o example_kat ~/raw_data/HiFi.fastq genome.fa
|
| 参数 |
说明 |
-t 40 |
40 线程 |
-o example_kat |
输出前缀 |
~/raw_data/HiFi.fastq |
原始 HiFi reads |
genome.fa |
组装后的基因组 |
工具: KAT — K-mer 分析工具，用于比较 reads 和基因组之间的 K-mer 一致性，评估组装完整度和污染。
2. Genome Annotation — 基因组注释
脚本: Genome_Annotation/annotation.sh
功能概述: 完整的基因组注释流程，包括重复序列注释、编码区域注释（HISAT2 + StringTie + TransDecoder + PASA）、同源蛋白预测（Exonerate）、de novo 基因预测（Augustus）、多证据整合（EVidenceModeler）以及 BUSCO 评估。
2.1 重复序列注释 — EDTA
1
| EDTA.pl --genome genome.fa --anno 1 --force 1 --debug 1 -sensitive 1 -t 24
|
| 参数 |
说明 |
--genome genome.fa |
输入基因组 fasta 文件 |
--anno 1 |
运行完整注释流程 |
--force 1 |
强制覆盖已有结果 |
--debug 1 |
开启调试模式，保留中间文件 |
-sensitive 1 |
使用敏感模式检测重复序列 |
-t 24 |
使用 24 线程 |
工具: EDTA — Extensive de novo Transposable Element Annotation，全基因组重复序列注释工具。
2.2 编码区域注释 — HISAT2 + StringTie + TransDecoder
Step 1: 构建 HISAT2 索引
1
| hisat2-build genome.fa genome.fa -p 4
|
| 参数 |
说明 |
genome.fa |
输入基因组（也作为索引前缀） |
-p 4 |
使用 4 线程 |
Step 2: 比对转录组 reads
1
| hisat2 -x genome.fa -1 trans_1.fq.gz -2 trans_2.fq.gz --dta -S trans.sam
|
| 参数 |
说明 |
-x genome.fa |
HISAT2 索引前缀 |
-1 trans_1.fq.gz |
双端测序 read 1 |
-2 trans_2.fq.gz |
双端测序 read 2 |
--dta |
为 StringTie/Assemble 转录本优化输出 |
-S trans.sam |
输出 SAM 格式比对文件 |
Step 3: SAM 转 BAM 并排序
1
2
| samtools view -b -S trans.sam > trans.bam
samtools sort -@ 16 trans.bam -o trans.sort.bam
|
| 参数 |
说明 |
-b -S |
以 BAM 格式输出 |
-@ 16 |
使用 16 线程排序 |
Step 4: StringTie 转录本组装
1
| stringtie -o genome.fa.gtf -l genome.fa trans.sort.bam
|
| 参数 |
说明 |
-o genome.fa.gtf |
输出 GTF 文件 |
-l genome.fa |
转录本前缀标签 |
工具: StringTie — 转录本组装与定量工具。
Step 5: TransDecoder 预测编码区
1
2
3
4
5
| perl gtf_genome_to_cdna_fasta.pl genome.fa.gtf genome.fa > genome.fa.trans.fa
TransDecoder.LongOrfs -t genome.fa.trans.fa
TransDecoder.Predict -t genome.fa.trans.fa
gtf_to_alignment_gff3.pl genome.fa.gtf > genome.fa.gff3
cdna_alignment_orf_to_genome_orf.pl genome.fa.trans.fa.transdecoder.gff3 genome.fa.gff3 genome.fa.trans.fa > transdecoder.genome.gff3
|
| 步骤 |
说明 |
gtf_genome_to_cdna_fasta.pl |
从 GTF 和基因组中提取 cDNA FASTA |
TransDecoder.LongOrfs |
识别长开放阅读框 (ORF) |
TransDecoder.Predict |
预测最终编码区域 |
gtf_to_alignment_gff3.pl |
GTF 转 GFF3 格式 |
cdna_alignment_orf_to_genome_orf.pl |
将 cDNA ORF 坐标映射回基因组坐标 |
工具: TransDecoder — 转录本编码区域预测工具。
2.3 Trinity 从头转录本组装 + PASA 优化
Step 1: Trinity 组装
1
| Trinity --seqType fq --max_memory 150G --left trans_1.fq --right trans_2.fq --SS_lib_type RF --CPU 20 --output trans_trinity --min_kmer_cov 2
|
| 参数 |
说明 |
--seqType fq |
输入为 fastq 格式 |
--max_memory 150G |
最大内存限制 150GB |
--left trans_1.fq |
左端 reads |
--right trans_2.fq |
右端 reads |
--SS_lib_type RF |
链特异性文库类型（R 表示 read1 反向，F 表示 read2 正向） |
--CPU 20 |
使用 20 线程 |
--output trans_trinity |
输出目录 |
--min_kmer_cov 2 |
K-mer 最低覆盖度为 2 |
工具: Trinity — 从头转录本组装软件。
Step 2: PASA 流程
1
2
3
4
5
6
7
8
9
| perl -pi -e "s#DATABASE.*#DATABASE=/program/rosa/genome.fa.database#g" alignAssembly.config
seqkit seq /program/rosa/5.anno/genome.fa -w 0 > genome.fa
samtools faidx genome.fa
seqkit seq trans_trinity/Trinity.fasta -w 0 > trinity.fa
samtools faidx trinity.fa
seqclean-x86_64/seqclean trinity.fa -v seqclean-x86_64/data/UniVec
samtools faidx trinity.fa.clean
perl -e 'while (<>) { print "$1\n"if />(\S+)/ }' trinity.fa > tdn.accs
Launch_PASA_pipeline.pl -c alignAssembly.config --ALIGNERS blat -R -C -g genome.fa -t trinity.fa.clean -T -u trinity.fa --TDN tdn.accs --trans_gtf genome.fa.gtf --CPU 10
|
| 步骤/参数 |
说明 |
alignAssembly.config |
PASA 配置文件，修改数据库路径 |
seqkit seq -w 0 |
将序列转为单行格式 |
seqclean |
清理 Trinity 组装结果中的载体/接头污染 |
-v seqclean-x86_64/data/UniVec |
使用 UniVec 数据库进行污染检测 |
Launch_PASA_pipeline.pl |
启动 PASA 流程 |
--ALIGNERS blat |
使用 BLAT 作为比对工具 |
-R |
保留旧注释 |
-C |
清理旧数据库 |
-g genome.fa |
基因组 fasta |
-t trinity.fa.clean |
转录本 fasta（已清理） |
-T |
使用 Trinity 模式 |
-u trinity.fa |
未清理的转录本（用于比对） |
--TDN tdn.accs |
转录本名称列表 |
--trans_gtf genome.fa.gtf |
已有转录本 GTF |
--CPU 10 |
使用 10 线程 |
工具: PASA — Program to Assemble Spliced Alignments，利用转录本比对提升基因注释质量。
Step 3: 提取 PASA 训练集
1
2
| pasa_asmbls_to_training_set.dbi --pasa_transcripts_fasta genome.fa.database.assemblies.fasta --pasa_transcripts_gff3 genome.fa.database.pasa_assemblies.gff3
pasa_asmbls_to_training_set.extract_reference_orfs.pl genome.fa.database.assemblies.fasta.transdecoder.genome.gff3 > best_candidates.gff3
|
| 说明 |
|
| 将 PASA 注释结果转换为 Augustus 训练集格式 |
|
2.4 同源蛋白预测 — Exonerate
1
| exonerate --model protein2genome ref.pep.fa genome.fa --showalignment false --showtargetgff true --percent 50 --bestn 1 --minintron 10 --maxintron 100000 > exonerate.ref.pep.fa.out
|
| 参数 |
说明 |
--model protein2genome |
蛋白质到基因组的比对模式 |
ref.pep.fa |
参考蛋白序列 |
genome.fa |
目标基因组 |
--showalignment false |
不显示比对详情 |
--showtargetgff true |
以 GFF 格式输出 |
--percent 50 |
最低 50% 蛋白覆盖度 |
--bestn 1 |
每个蛋白只保留最佳 1 个比对 |
--minintron 10 |
最短内含子 10 bp |
--maxintron 100000 |
最长内含子 100 kb |
工具: Exonerate — 蛋白质-基因组比对工具。
2.5 De novo 基因预测 — Augustus
训练 Augustus
1
| perl autoAugTrain.pl --genome=genome.fa --trainingset=best_candidates.gff3 --species=rosa_rna --useexisting 2> auto_training.log
|
| 参数 |
说明 |
--genome=genome.fa |
基因组文件 |
--trainingset=best_candidates.gff3 |
PASA 提取的训练集 |
--species=rosa_rna |
新建物种参数名为 rosa_rna |
--useexisting |
使用已有中间文件 |
运行 Augustus
1
| augustus --species=rosa_rna --AUGUSTUS_CONFIG_PATH=/usr/software/augustus-3.3.3/config --uniqueGeneId=true --noInFrameStop=true --gff3=on --strand=both genome.fa.mod.MAKER.masked > genome.fa.mod.MAKER.masked.augustus
|
| 参数 |
说明 |
--species=rosa_rna |
使用训练好的玫瑰物种参数 |
--AUGUSTUS_CONFIG_PATH |
Augustus 配置文件路径 |
--uniqueGeneId=true |
确保基因 ID 唯一 |
--noInFrameStop=true |
不允许框内终止密码子 |
--gff3=on |
输出 GFF3 格式 |
--strand=both |
正负链均预测 |
genome.fa.mod.MAKER.masked |
输入为 masking 后的基因组 |
工具: Augustus — De novo 基因预测工具。
2.6 多证据整合 — EVidenceModeler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
| perl EVidenceModeler/EvmUtils/partition_EVM_inputs.pl \
--genome genome.fa.mod.MAKER.masked \
--gene_predictions denovo.gff3 \
--protein_alignments gene.gff3 \
--transcript_alignments trans.gff3 \
--segmentSize 5000000 \
--overlapSize 10000 \
--partition_listing partitions_list.out
perl EVidenceModeler/EvmUtils/write_EVM_commands.pl \
--genome genome.fa.mod.MAKER.masked \
--gene_predictions denovo.gff3 \
--protein_alignments gene.gff3 \
--transcript_alignments trans.gff3 \
--weights /program/rosa/5.anno/1.lcfh1/2.geneanno/evm/weights.txt \
--output_file_name evm.out \
--partitions partitions_list.out > commands.list
EVidenceModeler/EvmUtils/recombine_EVM_partial_outputs.pl \
--partitions partitions_list.out \
--output_file_name evm.out
EVidenceModeler/EvmUtils/convert_EVM_outputs_to_GFF3.pl \
--partitions partitions_list.out \
--output_file_name evm.out \
--genome genome.fa.mod.MAKER.masked
cat */evm.out.gff3 > genome.fa.gene.gff3
|
| 参数 |
说明 |
--genome |
基因组文件 |
--gene_predictions denovo.gff3 |
Augustus de novo 预测结果 |
--protein_alignments gene.gff3 |
Exonerate 蛋白同源比对结果 |
--transcript_alignments trans.gff3 |
PASA 转录本比对结果 |
--segmentSize 5000000 |
分区大小 5 Mb |
--overlapSize 10000 |
分区间重叠 10 kb |
--weights weights.txt |
各类证据权重配置文件 |
--output_file_name evm.out |
输出文件名 |
--partitions |
分区列表文件 |
工具: EVidenceModeler (EVM) — 整合多种证据（基因预测、蛋白比对、转录本比对）生成一致的基因模型。
2.7 注释结果评估 — BUSCO
1
| busco -l ~/database/busco/eudicots_odb10 --offline -o prot_out -m prot -c 56 --out_path ./BUSCO_out -i pro.fa
|
| 参数 |
说明 |
-l ~/database/busco/eudicots_odb10 |
真双子叶植物 BUSCO 数据库 |
--offline |
离线模式 |
-o prot_out |
输出名称 |
-m prot |
评估蛋白序列 |
-c 56 |
使用 56 线程 |
--out_path ./BUSCO_out |
输出目录 |
-i pro.fa |
输入蛋白 fasta 文件 |
3. RNA-seq — 转录组分析
脚本: RNA_seq/RNA.sh
功能概述: RNA-seq 比对、转录本定量与表达量估计。
3.1 比对与定量
1
2
3
4
5
6
7
| hisat2 -x genome -1 RNA_R1.fastq.gz -2 RNA_R2.fastq.gz --rna-strandness RF --summary-file genome_summary --new-summary -p 24 \
| samtools view -@ 24 -Sb - \
| samtools sort -@ 24 -o genome.sorted.bam - >&2 2>genome_log
samtools index genome.sorted.bam
~/software/stringtie genome.sorted.bam -e -G genome.gff3 -b genome.Ballgown -o genome.gtf
|
HISAT2 参数
| 参数 |
说明 |
-x genome |
HISAT2 索引前缀 |
-1 RNA_R1.fastq.gz |
RNA-seq read 1 |
-2 RNA_R2.fastq.gz |
RNA-seq read 2 |
--rna-strandness RF |
链特异性文库（R=reverse, F=forward） |
--summary-file genome_summary |
输出比对统计摘要 |
--new-summary |
使用新版摘要格式 |
-p 24 |
24 线程 |
| 参数 |
说明 |
-Sb - |
从 stdin 读取 SAM，转为 BAM 输出 |
-@ 24 |
24 线程 |
StringTie 参数
| 参数 |
说明 |
-e |
只对参考注释中的转录本进行丰度估计 |
-G genome.gff3 |
参考基因注释 GFF3 |
-b genome.Ballgown |
输出 Ballgown 格式（用于后续差异表达分析） |
-o genome.gtf |
输出 GTF 文件 |
4. Pangenome — 泛基因组构建
4.1 基因组图构建 — PGGB
脚本: Pangenome/pggb.sh
功能概述: 将多个基因组 fasta 合并后，使用 PGGB 构建泛基因组图（variation graph）。
1
2
3
4
5
6
7
8
9
10
11
| id=$(pwd | cut -d "/" -f 9)
ls *.fa | while read f; do
sample_name=$(echo $f | cut -f 1 -d '.')
fastix -p "${sample_name}#" $f >> ${id}_all.fasta
done
bgzip -@ 15 ${id}_all.fasta
samtools faidx ${id}_all.fasta.gz
mash triangle ${id}_all.fasta.gz > ${id}_all.fasta.gz.mash_triangle.txt
sed 1,1d Chr1_all.fasta.gz.mash_triangle.txt | tr '\t' '\n' | grep chr -v | LC_ALL=C sort -g -k 1nr | uniq | less -S
pggb -i ${id}_all.fasta.gz -o ${id}_pggb_out -n 55 -t 20 -p 90 -S -m -V A01_1:
|
样本合并
| 步骤 |
说明 |
pwd | cut -d "/" -f 9 |
从路径中提取染色体号（如 Chr1） |
fastix -p "${sample_name}#" |
给每条序列添加样本名前缀 |
bgzip -@ 15 |
使用 15 线程压缩合并后的 fasta |
samtools faidx |
建立 fasta 索引 |
距离矩阵
| 参数 |
说明 |
mash triangle |
计算所有样本间的 Mash 距离矩阵 |
PGGB 运行
| 参数 |
说明 |
-i |
输入压缩 fasta 文件 |
-o |
输出目录 |
-n 55 |
样本数量为 55 |
-t 20 |
使用 20 线程 |
-p 90 |
最小序列相似度 90% |
-S |
将大型 segmental duplications 切分为 unitig |
-m |
采用精确模式 |
-V A01_1:#,A01_2:#,... |
55 个样本的路径名映射（用于图中的路径标识） |
工具: PGGB — PanGenome Graph Builder，基于 VG 工具包构建泛基因组图。
包含的 55 个样本:
A01_1, A01_2, A02_1, A02_2, A03_1, A03_2, A04_1, A04_2, A05_1, A05_2, A06_1, A06_2, A10_1, A10_2, A11_1, A11_2, A12_1, A12_2, A13_1, A13_2, A14_1, A14_2, A15_1, A15_2, A16_1, A16_2, A17_1, A17_2, A18_1, A18_2, A19_1, A19_2, A20_1, A20_2, A21_1, A21_2, A23_1, A23_2, A25_1, A25_2, A27_1, A27_2, A28_1, A28_2, A29_1, A29_2, A31_1, A31_2, DR21, HR20
4.2 直系同源基因分析 — OrthoFinder
脚本: Pangenome/orthofinder.sh
功能概述: 识别多个基因组间的直系同源基因群（orthogroups），用于泛基因组核心基因和可变基因的划分。
1
| orthofinder -f orthofinder_workdir -t 20 -S diamond -M msa
|
| 参数 |
说明 |
-f orthofinder_workdir |
输入目录（包含各物种蛋白 fasta） |
-t 20 |
使用 20 线程 |
-S diamond |
使用 DIAMOND 作为序列比对工具（速度更快） |
-M msa |
推断基因树时进行多序列比对 (MSA) |
工具: OrthoFinder — 直系同源基因推断、基因树构建和物种树推断工具。
4.3 结构变异检测 — SyRI
脚本: Pangenome/syri.sh
功能概述: 比对两个基因组，检测结构变异（包括倒位、易位、重复等）。
1
2
| minimap2 -ax asm5 -t 4 --eqx ref.fa quesry.fa | samtools sort -O BAM - > ref.query.bam
syri -c ref.query.bam -r ref.fa -q query.fa -F B --prefix ref_query
|
Minimap2 参数
| 参数 |
说明 |
-ax asm5 |
使用 asm5 预设（适合近缘物种全基因组比对） |
-t 4 |
4 线程 |
--eqx |
输出 =/=X CIGAR 操作符（区分匹配和错配） |
ref.fa |
参考基因组 |
quesry.fa |
查询基因组（注：脚本中拼写为 quesry，可能为笔误） |
SyRI 参数
| 参数 |
说明 |
-c ref.query.bam |
输入排序后的 BAM 比对文件 |
-r ref.fa |
参考基因组 |
-q query.fa |
查询基因组 |
-F B |
仅检测 between-scaffold 的重排 |
--prefix ref_query |
输出文件前缀 |
工具: SyRI — 基因组结构变异与重排检测工具。
5. Phylogenetic — 系统发育与群体遗传
5.1 SNP 变异检测 — DeepVariant + GLnexus
脚本: Phylogenetic/Mapping_calling_SNP.sh
功能概述: 将 HiFi reads 比对到参考基因组，使用 DeepVariant 进行变异检测，再使用 GLnexus 合并多样本变异。
Step 1: Minimap2 比对
1
| minimap2 -t 10 -ax map-hifi --eqx --secondary=no ref.fa ~/data/genome_hifi.fastq.gz 2> log | samtools view -Sb - > genome.bam
|
| 参数 |
说明 |
-t 10 |
10 线程 |
-ax map-hifi |
HiFi 读长比对模式 |
--eqx |
输出详细的 =/=X CIGAR 操作符 |
--secondary=no |
不输出 secondary alignments |
ref.fa |
参考基因组 |
~/data/genome_hifi.fastq.gz |
HiFi reads |
Step 2: 标记重复并排序
1
2
3
4
| samtools view -h genome.bam | samblaster -e --maxReadLength 100000 --ignoreUnmated | samtools view -Sb - > genome.mark_dup.bam
samtools sort -@ 10 genome.mark_dup.bam -o genome.mark_dup.sorted.bam
samtools index -@ 5 genome.mark_dup.sorted.bam
rm -rf genome.bam genome.mark_dup.bam
|
| Samblaster 参数 |
说明 |
-e |
输出重复标记的 reads（不输出 non-clonal 分割 reads） |
--maxReadLength 100000 |
最大读长 100 kb |
--ignoreUnmated |
忽略未配对的 reads |
工具: samblaster — 标记 PCR/optical 重复 reads。
Step 3: DeepVariant 变异检测
1
2
3
4
5
6
7
8
9
10
| singularity exec -B $PWD:/input -B $PWD:/output \
~/software/deepvariant_calling/deepvariant_1.4.0.sif \
/opt/deepvariant/bin/run_deepvariant \
--model_type PACBIO \
--ref=/input/outgroup.fa \
--reads=/input/genome.mark_dup.sorted.bam \
--output_gvcf=/output/genome.g.vcf.gz \
--output_vcf=/output/genome.vcf.gz \
--num_shards 20 \
--intermediate_results_dir=/output/deepvariant_tmp_output/genome
|
| 参数 |
说明 |
--model_type PACBIO |
使用 PacBio 模型 |
--ref |
参考基因组 |
--reads |
排序并去重的 BAM 文件 |
--output_gvcf |
输出 GVCF 文件（用于后续合并） |
--output_vcf |
输出 VCF 文件 |
--num_shards 20 |
并行分片数 20 |
--intermediate_results_dir |
中间结果目录 |
运行环境: Singularity 容器 (deepvariant_1.4.0.sif)
工具: DeepVariant — 基于深度学习的变异检测工具。
Step 4: GLnexus 合并多样本变异
1
2
| singularity exec -B $PWD:/input -B $PWD:/output ~/software/GLnexus/glnexus_cli.sif \
glnexus_cli --config DeepVariantWGS --dir /input/GLnexus.DB --bed /input/input.bed -t 10 --list /input/input.list.use_change > merge.bcf
|
| 参数 |
说明 |
--config DeepVariantWGS |
使用 DeepVariant WGS 配置 |
--dir |
GLnexus 数据库目录 |
--bed |
BED 文件（指定分析区间） |
-t 10 |
10 线程 |
--list |
输入 GVCF 文件列表 |
运行环境: Singularity 容器 (glnexus_cli.sif)
工具: GLnexus — 多样本变异联合基因型推断工具。
5.2 群体基因流分析 — Dsuite
脚本: Phylogenetic/Dsuite.sh
功能概述: 使用 Dsuite 的 Dtrios 模块进行 ABBA-BABA 检验，检测群体间的基因流/渗入（introgression）。
1
| ~/software/Dsuite/Build/Dsuite Dtrios merge.snp.filter.vcf_pass species_sets.txt -n gene_flow_with_tree -t tree_file
|
| 参数 |
说明 |
Dtrios |
Dsuite 的 trios 分析模块 |
merge.snp.filter.vcf_pass |
输入 VCF 文件（已过滤的 SNP） |
species_sets.txt |
物种/群体分组定义文件 |
-n gene_flow_with_tree |
输出前缀 |
-t tree_file |
输入物种树文件（Newick 格式） |
工具: Dsuite — 利用 D 统计量（ABBA-BABA test）检测群体间基因流和种系不一致。
5.3 群体历史大小推断 — PSMC
脚本: Phylogenetic/PSMC.sh
功能概述: 使用 PSMC (Pairwise Sequentially Markovian Coalescent) 方法推断群体有效大小历史变化。
1
2
3
4
| for id in {1..100} ; do
echo "~/software/psmc/psmc -N25 -t13 -r5 -b -p 4+25*2+4+6 -o genome.round-${id}.psmc genome.psmcfa" >> psmc_bs.sh
done
ParaFly -c psmc_bs.sh -CPU 5
|
生成 100 次 bootstrap 脚本
使用 for 循环生成 100 个 bootstrap 命令，每次使用不同的随机种子。
PSMC 参数
| 参数 |
说明 |
-N25 |
最多 25 次迭代 |
-t13 |
最大 coalescent 时间单位为 13 |
-r5 |
初始 θ/ρ 比值为 5 |
-b |
执行 bootstrap（由脚本在循环中每次生成） |
-p 4+25*2+4+6 |
时间间隔参数：自由突变率间隔设置 |
-o genome.round-${id}.psmc |
每轮 bootstrap 输出文件 |
genome.psmcfa |
PSMC 输入格式（consensus 序列） |
ParaFly 参数
| 参数 |
说明 |
-c psmc_bs.sh |
包含 100 个 bootstrap 命令的脚本 |
-CPU 5 |
使用 5 个 CPU 并行执行 |
工具: PSMC — 基于单个基因组的杂合度推断群体历史大小变化。
工具: ParaFly — 简单的并行执行工具。
附录：流程总览
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
| 原始数据 (HiFi reads, Hi-C reads, RNA-seq reads)
│
▼
┌───────────────────────┐
│ Genome Assembly │ Hifiasm + 3D-DNA + Hi-C
│ 基因组组装 │
└──────────┬────────────┘
│
▼
┌───────────────────────┐
│ Genome Annotation │ EDTA + HISAT2 + StringTie + TransDecoder
│ 基因组注释 │ + Trinity + PASA + Exonerate + Augustus + EVM
└──────────┬────────────┘
│
┌─────┴──────┐
▼ ▼
┌─────────┐ ┌────────────────────────────────────┐
│ RNA-seq │ │ Pangenome │
│ 定量 │ │ PGGB (图构建) + OrthoFinder (直系 │
│ │ │ 同源) + SyRI (结构变异) │
└─────────┘ └──────────────┬─────────────────────┘
│
▼
┌──────────────────────────┐
│ Phylogenetic │
│ 系统发育与群体遗传 │
│ DeepVariant (SNP) │
│ Dsuite (基因流) │
│ PSMC (历史群体大小) │
└──────────────────────────┘
|