跟着Nature学数据分析：minimap2+bcftools基于全基因组比对鉴定SNP

用户7010445

发布于 2023-08-23 10:45:20

7570

发布于 2023-08-23 10:45:20

文章被收录于专栏：小明的数据分析笔记本

论文

Cycles of satellite and transposon evolution in Arabidopsis centromeres

https://www.nature.com/articles/s41586-023-06062-z

代码链接

https://github.com/vlothec/pancentromere/tree/main

论文中有一部分内容是基于基因组比对检测SNP和InDel，论文中的方法部分写到

Genome assemblies of each *A.* *thaliana* accession were aligned to the HiFi assembly of Col-0 using minimap2 (v.2.24) with the parameters ‘minimap2 -a -x asm5 --cs -r2k -t 16’. SNPs and short indels were called from the alignment using bcftools

论文中提供了这部分的代码，今天的推文我们来学习一下这部分的代码

拟南芥基因组的数据还是用之前提到的NC论文里的数据

这里有一个问题暂时没搞明白：是否需要将contig水平的数据去掉，只保留染色体水平的序列？

首先是比对

minimap2 -a -x asm5 --cs -r2k -t 16 C24.fa Cvi.fa > Cvi_C24.sam
minimap2 -a -x asm5 --cs -r2k -t 16 C24.fa Kyo.fa > Kyo_C24.sam

这里 --cs参数和-r2k参数是什么意思暂时还没搞清楚

这个比对很快 84秒

sam文件转bam

samtools sort -@ 16 -O bam -o Kyo_C24.bam Kyo_C24.sam
samtools sort -@ 16 -O bam -o Cvi_C24.bam Cvi_C24.sam

samtools index Kyo_C24.bam
samtools index Cvi_C24.bam

bcftools变异检测

这里是分染色体的，我就只尝试一号染色体了

~/biotools/software.package/bcftools-1.17/bcftools mpileup -f C24.fa --threads 16 -r chr1 -o Cvi.chr1.vcf -A -O v Cvi_C24.bam
~/biotools/software.package/bcftools-1.17/bcftools mpileup -f C24.fa --threads 16 -r chr1 -o Kyo.chr1.vcf -A -O v Kyo_C24.bam

A参数和O参数v参数是什么意思

~/biotools/software.package/bcftools-1.17/bcftools call -c Kyo.chr1.vcf -o Kyo.chr1.called.vcf

~/biotools/software.package/bcftools-1.17/bcftools call -c Cvi.chr1.vcf -o Cvi.chr1.called.vcf

这里还有一个问题是如果是用基因组来做比对然后检测变异，时候需要包染色体倍性的参数设置为1

多样本vcf文件合并

bgzip -f -@ 16 Cvi.chr1.called.vcf
bgzip -f -@ 16 Kyo.chr1.called.vcf

tabix Kyo.chr1.called.vcf.gz
tabix Cvi.chr1.called.vcf.gz

~/biotools/software.package/bcftools-1.17/bcftools merge -m all --threads 16 -o chr1.vcf *.chr1.called.vcf.gz

bgzip -f -@ 16 chr1.vcf 
tabix chr1.vcf.gz

对数据进行过滤

~/biotools/software.package/bcftools-1.17/bcftools view -e 'GT[*]="mis"' chr1.vcf.gz > chr1.nomissing.vcf

~/biotools/software.package/bcftools-1.17/bcftools view -e 'GT[*]="het"' chr1.nomissing.vcf > chr1.nomissing.nohet.vcf

~/biotools/software.package/bcftools-1.17/bcftools view -M2 -m2 -v snps chr1.nomissing.nohet.vcf > chr1.nomissing.nohet.justSNPs.vcf

 ~/biotools/software.package/bcftools-1.17/bcftools view -C0 chr1.nomissing.nohet.vcf > chr1.nomissing.nohet.justNV.vcf

bgzip -f -@ 16 chr1.nomissing.nohet.justSNPs.vcf
tabix  chr1.nomissing.nohet.justSNPs.vcf.gz

去除重复序列区间内的snp位点

~/biotools/software.package/bcftools-1.17/bcftools view -R repeat.bed -o chr1.nomissing.nohet.justSNPs.norepeats.vcf chr1.nomissing.nohet.justSNPs.vcf.gz

提供一个bed文件

我这里随便构造一个