今天最新发表的土豆泛基因组Nature论文中基于基因组比对检测结构变异的分析代码

论文

Leveraging a phased pangenome for haplotype design of hybrid potato

https://www.nature.com/articles/s41586-024-08476-9

论文中的代码链接

https://github.com/Chenglin20170390/Haplotype-diversity/tree/main

论文中关于基于基因组比对检测结构变异的方法描述

Assembly-based calling was used to detect SVs with a minimum length of 50 bp. Genome sequences were aligned to the DMv6.1 reference genome to produce alignment BAM files using Minimap2 (v.2.17)⁹⁶ with “-ax asm5”. SyRI (v.1.5)⁹⁷ was then utilized to call genome-wide variants and generate assemblies-based VCFs. We kept variants of Ins, Del, Inv and Dup from the output of SyRI as the individual SV dataset. The accuracy of the SV dataset was validated by randomly selecting and manually verifying 100 SVs longer than 100 kb using Hi-C maps. Finally, we merged the individual SVs based on 80% overlap using SURVIVOR with the parameters “0.8 1 1 -1 -1 -1” to produce the final population SVs file. For SNPs and indels, we merged the output of assembly-based variation calling from SyRI using BCFtools (v.1.13)⁹

今天的推文里用拟南芥的数据集来复现一下这个流程

第一步是minimap2比对

minimap2 -ax asm5 -t 12 --eqx ../00.ref/at.col0.chr.fna ../00.qrys/An1.fa | samtools sort -O BAM - > An1.bam
samtools index An1.bam

minimap2 -ax asm5 -t 12 --eqx ../00.ref/at.col0.chr.fna ../00.qrys/C24.fa | samtools sort -O BAM - > C24.bam
samtools index C24.bam

minimap2 -ax asm5 -t 12 --eqx ../00.ref/at.col0.chr.fna ../00.qrys/Kyo.fa | samtools sort -O BAM - > Kyo.bam
samtools index Kyo.bam

第二步是syri检测变异

变异是 插入 删除 易位 和 重复

syri -c An1.bam -r ../00.ref/at.col0.chr.fna -q ../00.qrys/An1.fa -F B --prefix An1_
syri -c Kyo.bam -r ../00.ref/at.col0.chr.fna -q ../00.qrys/Kyo.fa -F B --prefix Kyo_
syri -c C24.bam -r ../00.ref/at.col0.chr.fna -q ../00.qrys/C24.fa -F B --prefix C24_

第三步是根据变异长度进行过滤

他这里的做法是根据输出的out文件进行过滤

mkdir 01_bed
awk '$11=="DEL" && $3-$2>49' An1_syri.out > An1_DEL.50.bed

awk '$11=="INS" && $8-$7>49' An1_syri.out > An1_INS.50.bed

awk '$11=="INV" && $3-$2>49' An1_syri.out > An1_INV.50.bed

awk '$11=="DUP" && $3-$2>49' An1_syri.out > An1_DUP.50.bed

###
awk '$11=="DEL" && $3-$2>49' C24_syri.out > C24_DEL.50.bed
awk '$11=="INS" && $8-$7>49' C24_syri.out > C24_INS.50.bed
awk '$11=="INV" && $3-$2>49' C24_syri.out > C24_INV.50.bed
awk '$11=="DUP" && $3-$2>49' C24_syri.out > C24_DUP.50.bed

###

awk '$11=="DEL" && $3-$2>49' Kyo_syri.out > Kyo_DEL.50.bed
awk '$11=="INS" && $8-$7>49' Kyo_syri.out > Kyo_INS.50.bed
awk '$11=="INV" && $3-$2>49' Kyo_syri.out > Kyo_INV.50.bed
awk '$11=="DUP" && $3-$2>49' Kyo_syri.out > Kyo_DUP.50.bed

mv *.bed 01_bed/

第四步是对变异进行过滤

这里的做法是根据hifi reads的覆盖度进行过滤

这里先跳过

第五步是bed文件转vcf文件

SURVIVOR bedtovcf 01_bed/An1_DEL.50.bed DEL 02_bed2vcf/An1_DEL.50.survivor.vcf
sed -i 's/Sample/An1/g' 02_bed2vcf/An1_DEL.50.survivor.vcf
## 这里的sed是替换默认的样本名Sample为自己想要的样本名

SURVIVOR bedtovcf 01_bed/C24_DEL.50.bed DEL 02_bed2vcf/C24_DEL.50.survivor.vcf
sed -i 's/Sample/C24/g' 02_bed2vcf/C24_DEL.50.survivor.vcf

SURVIVOR bedtovcf 01_bed/Kyo_DEL.50.bed DEL 02_bed2vcf/Kyo_DEL.50.survivor.vcf
sed -i 's/Sample/Kyo/g' 02_bed2vcf/Kyo_DEL.50.survivor.vcf

另外三种变异类型也是一样的操作，这里先不写了

这里有一个疑问是为什么不直接操作vcf文件，而是操作out文件再转换为vcf文件呢？

还需要将基因型替换，默认的都是./.

论文里的处理方式是把 Hap1 替换为 1|0

把Hap 2替换为 0|1

我这里直接替换成 1|1

sed -i 's/.\/./1|1/g' 02_bed2vcf/An1_DEL.50.survivor.vcf
sed -i 's/.\/./1|1/g' 02_bed2vcf/C24_DEL.50.survivor.vcf
sed -i 's/.\/./1|1/g' 02_bed2vcf/Kyo_DEL.50.survivor.vcf

第六步是合并多样本vcf

这里的合并标准是根据80%的overlap

ls 02_bed2vcf/*.vcf > merge_vcf_list_DEL
SURVIVOR merge merge_vcf_list_DEL 0.8 1 1 1 -1 50 merged_DEL.vcf
sed -i 's/.\/./0|0/g' merged_DEL.vcf

这里对合并后的vcf文件里的 ./.基因型的处理方式是替换成0/0

第七步是合并多个变异类型的vcf文件

使用的是 bcftools 这里我就不重复这个步骤了

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！