Snp-calling流程(BWA+SAMTOOLS+BCFTOOLS)

来源:生信菜鸟团评论11,774

比对可以选择BWA或者bowtie,测序数据可以是单端也可以是双端,我这里简单讲一个,但是脚本都列出来了。而且我选择的是bowtie比对,然后单端数据。

首先进入hg19的目录,对它进行两个索引

samtools faidx hg19.fa
Bowtie2-build hg19.fa hg19

我这里随便从26G的测序数据里面选取了前1000行做了一个tmp.fa文件,进入tmp.fa这个文件的目录进行操作

Bowtie的使用方法详解见https://www.plob.org/article/932.html

bowtie2 -x ../../../ref-database/hg19 -U tmp1.fa -S tmp1.sam
samtools view -bS tmp1.sam > tmp1.bam
samtools sort tmp1.bam tmp1.sorted
samtools index tmp1.sorted.bam 
samtools mpileup -d 1000 -gSDf ../../../ref-database/hg19.fa tmp1.sorted.bam |bcftools view -cvNg – >tmp1.vcf

然后就能看到我们产生的vcf变异格式文件啦!

当然,我们可能还需要对VCF文件进行再注释!

要看懂以上流程及命令,需要掌握BWA,bowtie,samtools,bcftools,

数据格式fasta,fastq,sam,vcf,pileup

如果是bwa把参考基因组索引化,然后aln得到后缀树,然后sampe对双端数据进行比对

首先bwa index 然后选择算法,进行索引。

然后aln脚本批量处理

==> bwa_aln.sh <==

while read id
do
echo $id
bwa aln hg19.fa $id >$id.sai
done <$1

然后sampe脚本批量处理

==> bwa_sampe.sh <==

while read id
do
echo $id
bwa sampe hg19.fa $id*sai $id*single >$id.sam
done <$1

然后是samtools的脚本

==> samtools.sh <==

while read id
do
echo $id
samtools view -bS $id.sam > $id.bam
samtools sort $id.bam $id.sorted
samtools index $id.sorted.bam
done <$1

然后是bcftools的脚本

==> bcftools.sh <==

while read id
do
echo $id
samtools mpileup -d 1000 -gSDf ref.fa $id*sorted.bam |bcftools view -cvNg – >$id.vcf
done <$1

==> mpileup.sh <==

while read id
do
echo $id
samtools mpileup -d 100000 -f hg19.fa $id*sorted.bam >$id.mpileup
done <$1

原文来自:http://www.bio-info-trainee.com/439.html

发表评论

匿名网友