跟着PNAS学数据分析：MUM&Co软件基于基因组做结构变异检测

MUM&Co软件对应的论文

MUM&Co: accurate detection of all SV types through whole-genome alignment

https://academic.oup.com/bioinformatics/article/36/10/3242/5756209?login=false

工具的github主页 https://github.com/SAMtoBAM/MUMandCo

工具是一个shell脚本，需要安装MUMmer4

能够检测的变异类型有 Deletions, insertions, tandem duplications and tandem contractions (>=50bp & <=150kb) Inversions (>=1kb) , translocations (>=10kb)

可可树基因组结构变异对应的论文

Genomic structural variants constrain and facilitate adaptation in natural populations of Theobroma cacao, the chocolate tree

https://www.pnas.org/doi/10.1073/pnas.2102914118

巧克力树结构变异PNAS.pdf

这篇论文里有对应的分析代码，github主页是

https://github.com/thamala/cacaoSV/tree/main

31个基因组，二倍体，做的是单倍型的组装，每个单倍型分别与参考基因组比对做结构变异检测，论文里提供了工具可以报两个单倍型分别检测sv的结果合并成一个二倍体的结果。之前看人类的结构变异分析的论文里也提供了类似的工具

用拟南芥的数据尝试一下这个工具

首先是比对做结构变异检测

time bash ../MUMandCo-master/mumandco_v3.8.sh -r An1.fa -q C24.fa -g 120000000 -o C24 -t 8
time bash ../MUMandCo-master/mumandco_v3.8.sh -r An1.fa -q Eri.fa -g 120000000 -o Eri -t 8
time bash ../MUMandCo-master/mumandco_v3.8.sh -r An1.fa -q Kyo.fa -g 120000000 -o Kyo -t 8
time bash ../MUMandCo-master/mumandco_v3.8.sh -r An1.fa -q Sha.fa -g 120000000 -o Sha -t 8

用拟南芥的数据，基因组大小在120M左右，5分钟一个样本

输出的文件内容

image.png

这里有一个vcf文件，有一个tsv文件，tsv文件是下面合并多个样本需要用到的数据

合并多个样本的工具是mumco2vcf.c

首先是编译

gcc mumco2vcf.c -o mumco2vcf -lm

mumco2vcf 这个是可执行文件

将需要合并的样本名整理到一个文本文件里，一行一个样本，如果是一个样本的两个单倍型，放在同一行，用制表符分隔

比如我这里是四个样本，准备的示例文件

ls */*all.tsv > vcf.list


C24_output/C24.SVs_all.tsv
Eri_output/Eri.SVs_all.tsv
Kyo_output/Kyo.SVs_all.tsv
Sha_output/Sha.SVs_all.tsv

合并运行的命令

../cacaoSV-main/mumco2vcf -i vcf.list > out.vcf

还有好几个参数

Usage: -i [file] File listing samples to combine. Haplotypes are in one row, separated by a tab. -t [str] either INS, DEL, DUP, TRA, or INV. (optional) -l [int]-[int] minimum and maximum length of the SVs. (optional) -o [int] overlap in base pairs to combine variants. (optional)

如果是单倍型

vcf.list的文件内容

C24_output/C24.SVs_all.tsv      Eri_output/Eri.SVs_all.tsv
Kyo_output/Kyo.SVs_all.tsv      Sha_output/Sha.SVs_all.tsv

合并的命令是一样的

欢迎大家关注我的公众号

小明的数据分析笔记本

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

image.png