Snp-calling流程（BWA+SAMTOOLS+BCFTOOLS）

来源：生信菜鸟团评论12,035

比对可以选择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

发表评论

匿名网友