小马的生信笔记

引言

在基因组学研究中,结构变异(structural variation, SV)作为驱动基因组多样性和表型分化的重要遗传基础,已逐渐成为继SNP之后的研究热点。尤其是在多倍体基因组、重复序列富集区域以及复杂重排频发的基因组背景下,SV不仅广泛存在,而且往往以嵌套、组合或多类型叠加的形式出现。这类复杂结构变异(complex SVs),如插入-缺失复合事件、倒位伴随重复扩增、以及多断点重排等,显著增加了变异解析的难度。

然而,现有主流SV检测工具大多依赖于单一信号(如read pair、split read或depth)或针对简单SV类型(DEL、INS、INV等)进行建模,在面对复杂SV时往往存在断点解析不准确、事件拆分不合理以及跨样本一致性差等问题。这种方法学上的局限,直接制约了我们对复杂基因组结构变异全貌的认识,尤其是在群体水平和进化研究中,复杂SV的系统性解析仍然面临巨大挑战。

在这样的背景下,Swave 作为发表于 Nature Genetics 的新方法,提供了一种面向复杂SV解析的新思路。Swave并不是简单地在线性参考上寻找变异信号,而是直接基于由一个或多个基因组组装构建的泛基因组图(pangenome graph)进行SV检测,从而能够更自然地刻画复杂的变异模式。根据论文摘要与软件说明,Swave利用循环神经网络(RNN)识别图中的结构变异模式,不仅能够检测常规SV,也特别强调对复杂SV的识别能力。

Swave的另一个突出优势在于,它面向的不只是单个样本,而是能够在群体水平对SV进行统一表征。软件可对多样本、多单倍型数据进行分析,并输出多等位(multi-allelic)和拆分后的双等位(bi-allelic)结果,这使其更适合后续群体遗传学分析、关联分析以及复杂变异在不同材料中的比较研究。对于当前越来越常见的高质量组装数据和泛基因组场景而言,这种从“单样本变异发现”走向“群体水平复杂SV刻画”的设计,正是Swave最具吸引力的地方。

Swave

Swave 是一种面向 pangenome graph 的结构变异检测与分型工具,其核心思想是将参考与变异等位序列之间的局部比对关系由传统的线性差异(如长度变化)转化为 dotplot 表征,并进一步压缩为具有方向信息的 projection waves,从而捕捉复杂结构模式;在此基础上,通过 Bi-LSTM 模型对波形特征进行序列建模与分类,实现对插入、缺失、倒位、重复及倒位重复等简单 SV 的识别,并可将共享断点的多个成分整合为复杂结构变异(CSV)。功能上,Swave不仅支持从 pangenome graph 中系统解析 SV,还能直接输出群体层面的基因型与频率信息,适用于大规模群体变异图谱构建。

下载

				
					git clone https://github.com/songbowang125/Swave.git
cd Swave
mamba env create -f ./environment.yml 
mamba activate Swave-env
				
			

示例数据运行

下载示例数据

				
					git clone https://github.com/songbowang125/Swave.git
cd Swave
mamba env create -f ./environment.yml 
mamba activate Swave-env
				
			

运行

				
					python ../Swave.py call --input_path assemblies-trio-hg00514.tsv  --ref_path t2t.chr20.fa --output_path ./ --gfa_path pg.minigraph.gfa --gfa_source minigraph
				
			

主要的输出文件如下:

  • swave.sample_level.split.vcf:保留原始 SV 结构,**一个位点可以有多个 ALT(多等位)**
  • swave.sample_level.vcf:将多等位拆开,**每个记录只有一个 ALT(双等位)**

示例数据运行

需要准备的数据:

  • 图泛的gfa文件
  • vcf文件
  • 基因组文件
  • 配置文件

配置文件格式如下:

				
					HGSVC3_HG00512  HG00512.vrk-ps-sseq.asm-hap1.fasta      HG00512.vrk-ps-sseq.asm-hap2.fasta
HGSVC3_HG00513  HG00513.vrk-ps-sseq.asm-hap1.fasta      HG00513.vrk-ps-sseq.asm-hap2.fasta
HGSVC3_HG00514  HG00514.vrk-ps-sseq.asm-hap1.fasta      HG00514.vrk-ps-sseq.asm-hap2.fasta
				
			

第一列为sample name,第二列第三列分别为两个hap的基因组文件(实测下来,只有一个hap也是可以的)

运行

				
					python ../Swave.py call --input_path sample.tsv  \
--ref_path Rename_Chrysanthemum_rhombifolium.fasta \
--gfa_source pggb \
--gfa_path Combine.fasta.gz.b1c69d8.11fba48.4f40d20.smooth.final.gfa \
--decomposed_vcf --output_path ./
				
			

第一次运行的时候会报错,具体信息是因为加载*gfa2fa.fa文件的时候出错了,试了几次是因为这个软件会调用gfatools基于gfa文件生成对应的文件,但是官方的environment.yml文件里却没有gfatools,应该是是一个小bug,没关系手动安装即可。

				
					git clone https://github.com/lh3/gfatools
gfatools && make

# 编辑配置文件,设置环境变量,方便调用
export PATH=/home/majunpeng/Software/gfatools:$PATH
				
			

再次运行

				
					python ../Swave.py call --input_path sample.tsv  \
--ref_path Rename_Chrysanthemum_rhombifolium.fasta \
--gfa_source pggb \
--gfa_path Combine.fasta.gz.b1c69d8.11fba48.4f40d20.smooth.final.gfa \
--decomposed_vcf --output_path ./
				
			

和原vcf文件进行比较

				
					# Swave 输出的VCF
grep "191189867" swave.sample_level.vcf

Chrysanthemum_rhombifolium#1#Chr6       191189863       >40024932>40024960      >40024933>40024936>40024937>40024939>40024940>40024941>40024944 >40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024958>40024959,>40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024945>40024947>40024948>40024950>40024952>40024953>40024955>40024956>40024959,>40024933>40024934>40024943>40024946>40024947>40024948>40024949>40024951>40024952>40024954>40024955>40024957>40024959,>40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024957>40024959        .       PASS    END=191190266;AC=5;AF=1;AN=5;NS=6;LV=0  GT:TYPE:LENGTH:BKPS     1:INS+DEL:513:INS_131_Chrysanthemum-rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum-rhombifolium#1#Chr6_191189867_191190248     2:INS+DEL:513:INS_131_Chrysanthemum-rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum-rhombifolium#1#Chr6_191189867_191190248     3:INS+DEL:2142:INS_1761_Chrysanthemum-rhombifolium#1#Chr6_191189879_191189879,DEL_381_Chrysanthemum-rhombifolium#1#Chr6_191189879_191190259   4:INS+DEL:513:INS_131_Chrysanthemum-rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum-rhombifolium#1#Chr6_191189867_191190248     .:.:.:. 1:INS+DEL:513:INS_131_Chrysanthemum-rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum-rhombifolium#1#Chr6_191189867_191190248

# Swave 输出的split.vcf
grep "191189867" swave.sample_level.split.vcf 
Chrysanthemum_rhombifolium#1#Chr6       191189867       >40024932>40024960_>40024932>40024960   >40024933>40024936>40024937>40024939>40024940>40024941>40024944       >40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024945>40024947>40024948>40024950>40024952>40024953>40024955>40024956>40024959     .PASS    END=191190248;SVLEN=513;SVTYPE=INS+DEL;BKPS=INS_131_Chrysanthemum_rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum_rhombifolium#1#Chr6_191189867_191190248;AC=1;AF=0.2;AN=5;NS=6;LV=0    GT      0       1       0       0       .       0
Chrysanthemum_rhombifolium#1#Chr6       191189867       >40024932>40024960_>40024932>40024960   >40024933>40024936>40024937>40024939>40024940>40024941>40024944       >40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024958>40024959     .PASS    END=191190248;SVLEN=513;SVTYPE=INS+DEL;BKPS=INS_131_Chrysanthemum_rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum_rhombifolium#1#Chr6_191189867_191190248;AC=2;AF=0.4;AN=5;NS=6;LV=0    GT      1       0       0       0       .       1
Chrysanthemum_rhombifolium#1#Chr6       191189867       >40024932>40024960_>40024932>40024960   >40024933>40024936>40024937>40024939>40024940>40024941>40024944       >40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024957>40024959     .PASS    END=191190248;SVLEN=513;SVTYPE=INS+DEL;BKPS=INS_131_Chrysanthemum_rhombifolium#1#Chr6_191189867_191189867,DEL_382_Chrysanthemum_rhombifolium#1#Chr6_191189867_191190248;AC=1;AF=0.2;AN=5;NS=6;LV=0    GT      0       0       0       1       .       0

# 原VCF
Chrysanthemum_rhombifolium#1#Chr6       191189863       >40024932>40024960      ATGAAAATAAAATCCGAAAAAAAAGTATATATAGGTTGATGTACTTATTGCGGTGTTTAAAGAATTTAGGCTCATGTTTCATAGCTCGTCTCGGGGTGTGTTTGACGAAAATAACACCAGCTATATTTTTGTAGCAGAAAAAAGTAGTACTCCCTCCATCCCTAAATAAGTGTCCAGTTGGACTTTTCAAAGTCAAACTGGATCAACTTTGACTATATCTATTTTTGCGAGTGTTATATAATACTTGATGCAAAATATATGAATTGATTGGGTTTTAAGTGTGTTTTTTCATTGATATAGTTTTTCATCAAGTATTATATAACACAAAAAAAATATTTACGGTCAAAGTTGAACTAGTTTTACTTTGAAAATTCAAACCGGACACTTATTTAGGGACGGAGGG       AAGAAAATAAAAACCGAAAAAAAGTATATATAGGTGGATGTACTTATTGCGGTGTTTAAAGAATTTAGGCTCATGTTTCATAGCTCGCCACAAGGTGTGTTTGGCGAAAATAACACTAGCTATATTTCTGTAACAGAAAAAAGT,AAGAAAATAAAAACCGAAAAAAAGTATATATAGGTGGATGTACTTATTGCGGTGTTTAAAGAATTTAGGCTCTTGTTTCATAGCTCGCCACGAGGTGTGTTTGGCGAAAATAACACTAGCTATATTTCTGTATCAGAAAAAAGT,ATGAAAATAAAAACCGAAATTTTTTATAGAGTAAACTACACCATTCGTCCTTGAGTAATACGTTGACCATCAAGTTTCGTCCTTTGACCAAAAAAATTGCACCGTTGGTCCTTAACTGCGAAACGTTAGTCATCGTTGGTCCCTCACCCTAACTGGTTGGACTAACGGGCGTTAAGCCTTGTTAAGTGGAGGGTATATTTGTCTTTTCCCACCCAAGATCTTTATTTTAATGCCCCCATCAACCCTAAACAATCACCAGACGCCCCCATCATTTTTTTTCTGTAACATATCAAGAAACATAAAAAAAGGAAAAAGAATATGAATCACTTGTCCTTTCATTACGCCTCATGCGTAATTAACCACAAACATCAACATTCAATTATTCTCTTGTGATTACGAGACTCATAAGTACCACGAATTCAAGATCTTTGGAAAAAATATCAGGTAGGTGTTTTAGTAATTTGTTTTCGGTGTGAATGATCTTTTGAATGATTAGAGTGCACGAAAAAAATTCAACTAAGATCTTAAAAATAGCTACTTTTCAGTCCTTCTCAGTCCTTGTAGTCATAGATCTGTGATGGATTTATCTTTATGGTTTTTGTCTTGGACGATTTTTGAGTACCGGAAAGGATAAAACTGAGATATAATTTTGTCATTTAGATTCGGTTCAGTGATGATATCATTCGTGTTTATGACCTTTGAGTCATCATGTATTACTTGAAAAAAATAACCACATGAAGAATGAATCTGAGATCCATTTATCCTTGTTTCTGTTTGAGTCGAGAAAGGTTGAGATTAGTCTAGGATAGTGGAGTTGATTGGATTTGTTGAAATACGTAAGATTGGTGATGGGTTTAGTATTAAAGTTGTGTTGATGAGTTTTGTGCTTCTCGAATTCTTCTTCTTCTTCTAGGATAGTGGAGTTGATTGGATTTGTTGAAATACATAAGATTGGTGATAGGTTTAGTATTAAAGTTGTGTTGATGGGTTTAGTATTAAAGGGTATCCTTGAGGCTGAGATCGAGGTGGGTATGGACTTGTGGTCGTGGTTATTGCTAGGTCGAAGAGTGGAGGATGGTTGGTGTTGCTGTCGTGGGATCACCTAACCTATCAAGAATCACCTTTTGCTACCTGAAAAGCCTTGCTGCTACTGCAAAAAAATTTTTTTTTCTTGCTTACCCATTTGTTTTTGTTGTTGGTGTTGTGGTTATTGATCTATTTGATGGGTTTTTTGAACTTGTGGTAAAAAACAATTTAATGAATACGATTATCGGGCTTTGTGAAGGTGGCCGGAGGTGGTGAGAATGGTGGTGGGGGAAGATGAAGATTGTGGGTTTAAAAAGATGGTGGGGAAGATGAAGATGGTGGGTTTAAAAAGATGGTGGGGGGAGATGAAGATGGTGATGGGGGAAGATGATGATTCAAGAAGAAGAAAATGATGGATGGTGCTGGCAAATATAATTGAGGGTGTGGTGTGGGTGGGAAAAGACAAATATACCCTCCACTTAACAAGGCTTAACGCCCGTTAGTCCAACCAGTTAGGGTGAGAGACCAACGATGACTAACGTTTCGCAGTTAAGGACCAACGGTGCAATTTTTTTGGTCAAAGGACGAAACTTGATGGTCAACGTATTACTCAAGGACGAATGATGTAGTTTACTCTTTTTTATATATAGGTTGATGTACTTATTGCGGTGTTTAAAGAATTTAGGCTCATGTTTCATAGCTCGCCACGAGGGTGTGTTTGGTGAAAATAACACTAGCTATATTTCTGTAGCAGAAAAAAGT,AAGAAAATAAAAACCGAAAAAAAGTATATATAGGTGGATGTACTTATTGCGGTGTTTAAAGAATTTAGGCTCATGTTTCATAGCTCGCCACAAGGTGTGTTTGGCGAAAATAACACTAGCTATATTTCTGTAGCAGAAAAAAGT     60      .       AC=2,1,1,1;AF=0.4,0.2,0.2,0.2;AN=5;AT=>40024932>40024933>40024936>40024937>40024939>40024940>40024941>40024944>40024960,>40024932>40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024958>40024959>40024960,>40024932>40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024945>40024947>40024948>40024950>40024952>40024953>40024955>40024956>40024959>40024960,>40024932>40024933>40024934>40024943>40024946>40024947>40024948>40024949>40024951>40024952>40024954>40024955>40024957>40024959>40024960,>40024932>40024935>40024936>40024938>40024939>40024941>40024942>40024943>40024946>40024947>40024949>40024950>40024952>40024953>40024955>40024957>40024959>40024960;NS=5;LV=0    GT      1       2       3       4       1
				
			

和原始的vcf文件相比,swave把这个插入重新鉴定成了一个deletion和一个insertion,并且输出的vcf文件有了变异的类型和长度等信息。至于为什么split.vcf把一个SV拆成了3个,我能力有限,有点不理解,欢迎大家在评论区讨论!

结语

Swave确实在解析SV,特别是复杂的SV有很大的优势,但是对于输出文件的格式我还是没听懂,大家可以在评论区讨论一下

发表评论

基因组专栏