SYRI 共线性分析技能
用于执行基因组间的 SYRI 共线性和结构变异分析。
技能说明
本技能提供完整的 SYRI 共线性分析流程，包括：
- 染色体提取与重命名
- minimap2 比对
- SYRI 结构变异分析
- plotsr 可视化
软件安装
1. minimap2
# conda 安装
conda install -c bioconda minimap2
# 或从源码编译
git clone https://github.com/lh3/minimap2.git
cd minimap2 && make
2. SYRI
# conda 安装（推荐）
conda install -c bioconda syri
# pip 安装
pip install syri
3. plotsr
# conda 安装
conda install -c bioconda plotsr
# pip 安装
pip install plotsr
4. Biopython（用于染色体处理）
pip install biopython
关键发现
1. 线程数设置
重要：使用 64-120 线程加速分析
N_THREADS=64 # 根据 CPU 核心数调整
minimap2 -ax asm5 --eqx -t $N_THREADS ref.fa query.fa > output.sam
2. minimap2 输出缓冲问题
问题：minimap2 运行但 SAM 文件始终为 0 字节（进程正在运行）
原因：stdout/stderr 输出缓冲导致文件写入延迟
解决：使用 stdbuf 禁用输出缓冲
stdbuf -oL -eL minimap2 -ax asm5 --eqx -t $N_THREADS ref.fa query.fa > output.sam
3. SAM 文件清理
minimap2 可能输出日志行到文件，需要清理：
grep -v "^\[" A_vs_Vc.sam > A_vs_Vc_clean.sam
mv A_vs_Vc_clean.sam A_vs_Vc.sam
标准流程
步骤 1：染色体提取与重命名
from Bio.Seq import Seq
from Bio import SeqIO
# 染色体对应关系（根据实际项目调整）
# 格式：(新染色体名称, 原始染色体名称, 是否反向互补)
mapping = [
("chr01", "Vcev1_p0.Chr02", True), # 反向互补
("chr02", "Vcev1_p0.Chr11", False),
# ... 其他染色体
]
src_file = "V_caesariense_W85-20_P0_v2.fasta"
out_file = "Vc_for_SYRI.fa"
records = {rec.id: rec for rec in SeqIO.parse(src_file, "fasta")}
with open(out_file, "w") as f:
for new_name, old_name, rc in mapping:
if old_name in records:
seq = records[old_name].seq
if rc:
seq = seq.reverse_complement()
f.write(f">{new_name}\n")
f.write(str(seq) + "\n")
print(f"{new_name} <- {old_name} {'(RC)' if rc else '(direct)'}")
步骤 2：生成染色体长度文件
from Bio import SeqIO
# Vc.len
with open("Vc.len", "w") as f:
for rec in SeqIO.parse("Vc_for_SYRI.fa", "fasta"):
f.write(f"{rec.id}\t{len(rec.seq)}\n")
# A.len
with open("A.len", "w") as f:
for rec in SeqIO.parse("subgenome_a_renamed.fa", "fasta"):
f.write(f"{rec.id}\t{len(rec.seq)}\n")
步骤 3：创建基因组信息文件
cat > genomes.txt << 'EOF'
#file name tags
A.len A ft:cl
Vc.len Vc ft:cl
EOF
步骤 4：minimap2 比对
# 关键参数：
# -ax asm5: 用于长读长的基因组比对
# --eqx: 将 CIGAR 中的 =/X 转换为 M（SYRI 需要）
# -t: 线程数（推荐 64-120）
# stdbuf: 解决输出缓冲问题
N_THREADS=64
stdbuf -oL -eL minimap2 -ax asm5 --eqx -t $N_THREADS \
/path/to/ref.fa \
/path/to/query.fa > A_vs_Vc.sam 2>&1 &
echo "Started minimap2, PID: $!"
# 等待完成（约 10-15 分钟）
步骤 5：SAM 文件预处理
# 1. 移除 minimap2 日志行（以 [ 开头的行）
grep -v "^\[" A_vs_Vc.sam > A_vs_Vc_clean.sam
mv A_vs_Vc_clean.sam A_vs_Vc.sam
# 2. 可选：过滤低质量比对（MAPQ >= 20）
awk 'BEGIN {FS="\t"} /^@/ {print; next} $5 >= 20 {print}' A_vs_Vc.sam > A_vs_Vc_filtered.sam
步骤 6：SYRI 分析
mkdir -p syri_out
syri -c A_vs_Vc_filtered.sam \
-r /path/to/ref.fa \
-q /path/to/query.fa \
-F S -k --dir syri_out
# 参数说明：
# -c: 比对文件 (SAM/BAM/PAF)
# -r: 参考基因组
# -q: 查询基因组
# -F S: 输入格式为 SAM
# -k: 保留中间文件
步骤 7：plotsr 可视化
# 单染色体
plotsr --sr syri_out/syri.out --genomes genomes.txt \
-o output.png --chr chr01 -H 8 -W 10
# 全染色体
plotsr --sr syri_out/syri.out --genomes genomes.txt \
-o output_all.png -H 8 -W 15
注意事项
1. 染色体命名一致性
- SYRI 要求 ref 和 query 的染色体名称和数量完全匹配
- 建议使用统一的命名格式（如 chr01-chr12）
2. CIGAR 格式
- SYRI 需要标准化的 CIGAR 字符串（使用 =/X 而非 M）
- 必须使用
--eqx 参数
3. 输入格式
- SAM 文件必须有正确的格式
- 使用
-F S 参数明确指定 SAM 格式
4. 内存和存储
- 大型基因组比对会产生很大的 SAM 文件（10-20GB）
- 确保有足够的磁盘空间
- 考虑过滤低质量比对以加速分析
5. 线程数选择
- 推荐使用 CPU 核心数的 50%-100%
- 过高线程数可能因内存带宽受限而无法提速
错误排查
错误：SAM 文件格式问题
Error: invalid literal for int() with base 10: '*'
解决：确保使用 -a 参数生成标准 SAM 格式
错误：CIGAR 格式问题
Error: Incorrect CIGAR string found
解决：添加 --eqx 参数
错误：染色体长度不匹配
Error: length in genome fasta is less than the maximum coordinate
解决：检查并确保 genomes.txt 中的长度与实际染色体匹配
错误：minimap2 进程运行但文件为 0
SAM file size = 0 bytes
解决：使用 stdbuf -oL -eL 禁用输出缓冲
输出文件说明
| 文件 | 说明 |
|---|
| syri.out | 结构变异注释（主要结果） |
| syri.vcf | VCF 格式结果 |
| synOut.txt | 共线性区域 |
| invOut.txt | 倒位区域 |
| TLOut.txt | 易位区域 |
| dupOut.txt | 重复区域 |
| snps.txt | SNP 位置 |
| sv.txt | 结构变异位置 |
文件位置（示例）
- 工作目录：
/path/to/syri/AVc_v2/
- SYRI 结果：
syri_out_v2/
- 图片输出：
/path/to/syri/plots/
