SYRI 共线性分析

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/