Archives

庖丁解牛,微生物基因组重测序:一条命令完成从sra到vcf

一、引子

前面在 Twitter上看到一个针对细菌基因组重测序流程(鉴定SNP)的使用工具调查,这里我们就庖丁解牛般解读下这些流程工具是如何工作的,流程又是如何设计的。

SNP Pipeline

该调查数据来自 google.docs

从图中可以看出snippy还是比较受欢迎的,Torsten Seemann的流程一向比较稳定,接口一向比较简单。

这次解读的主角是snippy,使用起来确实很方便,比如:

fastq-dump --defline-seq '@$sn/$ri' --defline-qual '+' ERR114970.sra -Z --split-3 |\
seqtk trimfq - | \
snippy --cpus 16 --outdir ERR114970 --ref GCF_000069185.1_ASM6918v1_genomic.gbff  \
--prefix ERR114970 -peil /dev/stdin &>log.txt

一条命令可以解决我们从sra到vcf到的问题,当然其中使用的工具很多了,但是都被封装在脚本里了。

二、流程细节

2.1. sra生成interleave格式fastq文件

fastq-dump --defline-seq '@$sn/$ri' --defline-qual '+' ERR114970.sra -Z --split-3

该结果使用-Z 选项将结果输出至标准输出流,格式是Interleave格式,这种模式有很多好处, 比如质量值修剪,序列比对等,后续很多程序都支持,这样比较方便以steaming方式工作。

内部链接: ENA和SRA数据预处理: SRA Toolkit 和 seqtk_sra

2.2 去除低质量序列

seqtk trimfq - 

steaming 模式切除低质量序列。

内部链接: NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

2.3. 序列比对

至此,snippy 预处理工作已经做完,现在可以进入snippy的工作模式,这里只记录主要步骤。

gbk转gff和fasta:

使用`BioPerl`从GBK文件中提取fasta和gff文件用于构建序列比对索引和snpEff库。

构建基因组文件序列

samtools faidx reference/ref.fa

构建基因组索引:

bwa index reference/ref.fa

序列比对:

bwa mem  -v 2 -M -R '@RG    ID:ERR114970    SM:ERR114970' \
    -p -t 16 reference/ref.fa /dev/stdin | \
    samtools view -@ 16 -F 3844 -q 60 -S -b -u -T reference/ref.fa - |\
    samtools sort -O bam -o ERR114970.bam -@ 16 -

构建bam文件索引:

samtools index ERR114970.bam

计算比对深度:

samtools depth -q 20 ERR114970.bam | bgzip > ERR114970.depth.gz

这一步完成了序列比对工作,使用的都是比较常规的操作。

2.4. Variation calling

freebayes-parallel reference/ref.txt 16 -p 1 -q 20 -m 60 \
     --min-coverage 10 -V -f \
     reference/ref.fa ERR114970.bam > ERR114970.raw.vcf

使用freebayes (Bayesian haplotype-based polymorphism discovery)鉴定 Variation, Freebayes的单倍体模型比较适合细菌基因组SNP/INDEL鉴定。

-q –min-base-quality, 最低碱基质量; -q 20; -m –min-mapping-quality,最低比对质量; -m 60; -p –ploidy 设置单倍体模型, -p 1; –min-coverage: 最小覆盖度。 –min-coverage 10;

snippy使用的是freebayes并行模式 freebayes-parallel

需要使用分割regions文件, 使用fasta_generate_regions.py 生成。

fasta_generate_regions.py reference/ref.fa.fai 80866 > reference/ref.txt

Variation 过滤:

snippy-vcf_filter --minqual 10 --mincov 10 --minfrac 0.9  ERR114970.raw.vcf > ERR114970.filt.vcf

使用下面参数控制variation 质量:

–mincov :位点最低覆盖度, 10 –minfrac :alt碱基最低比率,0.9

2.5. SNP/INDE 注释

创建配置文件和snpEff build库文件,然后使用snpEff

snpEff build -c reference/snpeff.config -dataDir . -gff3 ref
snpEff ann -no-downstream -no-upstream -no-intergenic -no-utr -c reference/snpeff.config -dataDir . -noStats ref ERR114970.filt.vcf > ERR114970.vcf

使用-no-downstream -no-upstream -no-intergenic -no-utr 过滤掉基因上游、下游、间区以及UTR区域的变异位点, 生成注释的VCF文件。

至此我们已经可以获得从SRA文件到VCF的细节过程。

2.6. reference guided assembly

snippy 提供了另外一个我们称之为”reference guided”的拼装,实现方式为:

vcf-consensus ERR114970.vcf.gz < reference/ref.fa > ERR114970.consensus.fa

三、总结

3.1. 设计程序时一定要注意区分:

标准输出流: stdout 标准错误流: stderr

准则: 如果不是结果文件,一律使用标准错误流打印。

3.2 Pipeline 设计

很多工作流,其实细节都很好实现,工作流主要起的作用就是封装,让接口更加简单, snippy 集成了 bioperl、bgzip、bwa、freebayes、samtools、tabix、vcfallelicprimitives、vcfstreamsort、vcfuniq等外部程序,使用这些程序可以使用相对路径调用,也可以设置PATH路径,但是最佳实践是环境变量只有程序能够访问,这样就不会造成程序版本污染。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-14 7:46 PM

Comments are closed.