论文
Cycles of satellite and transposon evolution in Arabidopsis centromeres
代码链接
论文中有一部分内容是基于基因组比对检测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.bambcftools变异检测
这里是分染色体的，我就只尝试一号染色体了
~/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.bamA参数和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文件
我这里随便构造一个
以上的命令是只保留这个bed文件区间内的位点
bcftools好多参数都不知道是什么意思，还需要仔细看看帮助文档
推文记录的是自己的学习笔记，内容可能会存在错误，请大家批判着看，欢迎大家指出其中的错误
声明：文中观点不代表本站立场。本文传送门：https://eyangzhen.com/44901.html
