跟着Nature Genetics学数据分析：nucmer+lastz+svum流程全基因组比对鉴定CNV

论文

Super-pangenome analyses highlight genomic diversity and structural variation across wild and cultivated tomato species

https://www.nature.com/articles/s41588-023-01340-y

西红柿NG_superPan正文.pdf

数据分析的代码

https://github.com/HongboDoll/TomatoSuperPanGenome

论文里提供了绝大部分的数据处理代码，很好的学习材料，今天的推文我们学习一下其中利用minimap2 和 syri两个软件最全基因组水平比对然后鉴定结构变异的代码

还有 nucmer+lastz+svum流程全基因组比对鉴定CNV的代码

首先是minimap2基因组水平比对的代码

minimap2 -t 20 -ax asm5 --eqx ref.fasta query.fasta | samtools view -b > query.bam

minimap2 -t 20 -ax asm10 --eqx ref.fasta query.fasta | samtools view -b > query.bam

这里asm5是同一个种不同物种个体之间的基因组比对用到的参数，asm10是跨种比对用到的参数

利用bam文件检测变异的代码

syri -c query.bam -r ref.fa -q query.fa -k -F B --prefix query --nc 12

输出文件里会有syri.vcf结尾的vcf文件，但是vcf文件中有的变异信息记录的不是很完整，比如INV和TRANS

INV是Inversion
TRANS是Translocation

这两种变异类型在vcf文件中ref和alt列记录的是

N

N

这个论文里提供了代码处理另外的输出文件invOut.txt获得新的vcf文件，需要重点学习一下这部分代码，搞定了再来记录

论文中还鉴定CNV这个是拷贝数变异，这个是用numer比对

用拟南芥的数据试试

nucmer -t 8 ../ref.fa ../only.chr.genome/Sha.fa --prefix sha_ref

seqkit split -i ../ref.fa
seqkit split -i ../only.chr.genome/Sha.fa


lastz ../ref.fa.split/ref.part_chr1.fa ../only.chr.genome/Sha.fa.split/Sha.part_chr1.fa --chain --format=general:name1,strand1,start1,end1,name2,strand2,start2,end2 --ambiguous=iupac > chr1.sha_ref.sam_lastz.txt

lastz这个工具是只能一条序列和另外一条序列比 （这个比对时间相对还是比较长的）

https://github.com/lastz/lastz

~/biotools/software.package/svmu-master/svmu sha_ref.delta ../ref.fa ../only.chr.genome/Sha.fa h chr1.sha_ref.sam_lastz.txt sha_ref_svmu

这一步生成4个文件

sv.sha_ref_svmu.txt
cords.sha_ref_svmu.txt
cm.sha_ref_svmu.txt
small.sha_ref_svmu.txt
cnv_all.sha_ref_svmu.txt

我这里cnv开头的文件是空的

lastz比对fasta文件里如果有多条序列也可以比，写法如下。NG文章里提供的代码应该是有问题的（当然不确定，很大可能是我自己理解有问题）

~/lastz-distrib/bin/lastz ../ref.fa[multiple] ../only.chr.genome/Sha.fa[multiple] --chain --format=general:name1,strand1,start1,end1,name2,strand2,start2,end2 --ambiguous=iupac > sha_ref.sam_lastz.txt

拟南芥的一对基因组比就用了将近400分钟6个多小时

real    378m26.295s
user    376m9.733s
sys    2m9.278s

不知道有没有参数可以改或者设置多线程之类的

 ~/biotools/software.package/svmu-master/svmu sha_ref.delta ../ref.fa ../only.chr.genome/Sha.fa h sha_ref.sam_lastz.txt full_genome_sha_ref_svmu

论文中提供的代码接下来用到了 filter_delta_based_svmu_cm.py这个脚本，我暂时没有找到这个脚本

把svum软件的输出结果保存成vcf格式，先用grep和awk对数据进行一个过滤

grep "CNV" sv.full_genome_sha_ref_svmu.txt | awk '$(NF-2)>50' | awk '$NF!=0&&$(NF-1)!=0&&$NF!="inf"&&$(NF-1)!="inf"' | awk '$NF/$(NF-1)>=2||$(NF-1)/$NF>=2' > sv.full_genome_sha_ref_svmu_CNV.xls

这里awk的 命令暂时搞不太懂是什么意思

转vcf用到的代码

cat sv.full_genome_sha_ref_svmu_CNV.xls | while read n
do
        chr=`echo $n | awk '{print $1}'`
        pos=`echo $n | awk '{print $2}'`
        ref_pos=`echo $n | awk '{if($2<=$3){print $1":"$2"-"$3}else{print $1":"$3"-"$2}}'`
        alt_pos=`echo $n | awk '{if($6<=$7){print $5":"$6"-"$7}else{print $5":"$7"-"$6}}'`
        ref_seq=`samtools faidx ../ref.fa $ref_pos | grep -v '>' | sed ':a;N;s/\n//g;ta'`
        alt_seq=`samtools faidx ../only.chr.genome/Sha.fa $alt_pos | grep -v '>' | sed ':a;N;s/\n//g;ta'`
        echo -e "${chr}\t${pos}\t.\t${ref_seq}\t${alt_seq}\t30\tPASS\t.\tGT\t1/1" >> sv.full_genome_sha_ref_svmu_CNV.vcf
done

bash cnv2vcf.sh

这个vcf文件是没有表头的那些内容的，需要再单独添加

svum 这个软件的github主页

https://github.com/mahulchak/svmu

下载下来解压，然后直接make就可以

对应的论文

https://www.nature.com/articles/s41467-019-12884-1

Structural variants exhibit widespread allelic heterogeneity and shape variation in complex traits

示例数据和代码可以给推文点赞，然后点击在看，最后留言获取

欢迎大家关注我的公众号

小明的数据分析笔记本

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