Calling SNPs用Samtools

基于reads_mapping-samtools-bcftools的流程常用于重reads中获取snp，基因组水平的长序列snp calling另有其他适配方法。

samtools获取

下载kneaddata后也有了samtools，但是版本过低，这里重新建环境下载最新版。

conda create -n snp
conda activate snp
conda install samtools
samtools --version

bcftools获取

conda install bcftools
bcftools --version

1 Bowtie比对 - 获取sam文件

建索引
-f fasta

bowtie2-build -f AF04-12_ref.fna AF04-12_ref

比对
-f fasta -x db index

bowtie2 -f -p 8 \
-U AF04-12_bgi.fna \
-x AF04-12_ref \
-S AF04-12_bgi.sam

(ERR): bowtie2-align died with signal 9 (KILL)
bowtie适合reads map to reference，基因组大片段比对不合适

2 ALFALFA比对 - 获取sam文件

ALFALFA：长序列对齐工具

建索引

# index
route="/hutongyuan/software/alfalfa"
mkdir index
$route/alfalfa index \
--reference AF04-12_ref.fna \
--prefix ./index/AF04-12_ref

比对

# mapping
$route/alfalfa align \
-r AF04-12_ref.fna \
-i ./index/AF04-12_ref \
-0 AF04-12_bgi.fna \
-o AF04-12_bgi.sam

*** Error in alfalfa': free(): invalid pointer: 0x00007f2ea089e1c0 ***
Aborted

3 Blast - 获取sam文件

安装blast

conda create -n blast
conda install -c bioconda blast
#blastn: 2.9.0+
#Package: blast 2.9.0, build May 31 2019 20:53:30

建索引

makeblastdb -in AF04-12_ref.fna \
-input_type fasta -dbtype nucl -title AF04-12_ref -parse_seqids \
-out ./Blastn/index/AF04-12_ref

比对

blastn -db ./index/AF04-12_ref \
-query ../AF04-12_bgi.fna \
-outfmt 17 \
-out ./AF04-12_bgi.sam \
-num_threads 8
#-perc_identity 80 -qcov_hsp_perc 80

sam -> bam -> sorted.bam -> bcf -> vcf

conda activate snp
# index
samtools faidx AF04-12_ref.fna
# sam -> bam
samtools view -@ 8 -b -S AF04-12_bgi.sam -o AF04-12_bgi.bam
# bam -> sorted.bam
samtools sort -@ 8 -l 9 -O BAM AF04-12_bgi.bam \
-o AF04-12_bgi.sorted.bam
# sorted.bam -> bcf
samtools mpileup -f ../AF04-12_ref.fna AF04-12_bgi.sorted.bam \
> test.bcf
# bcf -> vcf
bcftools view test.bcf > test.vcf

更多：
CGView Server:环形基因组可视化比较工具: blast比对基因组
Calling SNPs with Samtools & Calling SNPs with GATK
Samtools+bcftools Call SNP
Calling SNPs/INDELs with SAMtools/BCFtools
A tale of three next generation sequencing platforms: comparison of Ion Torrent, Pacific Biosciences and Illumina MiSeq sequencers. BMC Genomics 2012
突变注释工具SnpEff,Annovar,VEP,oncotator比较分析--转载
Bioinformatics at COMAV
A simple SNP calling pipeline - EMBL-EBI

使用Snippy构建细菌基因组SNP、core SNP和系统进化分析
细菌基因组 snps 、indels和core genome分析的软件