pbmm2+deepvariant识别变异位点

pbmm2支持pbcbio的bam文件输入，对hifi数据非常友好。

deepvariant据文献说是目前变异识别更为准确的一款软件。

参考链接

https://anaconda.org/bioconda/deepvariant
https://github.com/PacificBiosciences/pbmm2
https://www.jianshu.com/p/2b72fabca049
https://anaconda.org/bioconda/vcftools
https://vcftools.sourceforge.net/documentation.html#freq

一、安装

conda安装及使用比较麻烦，如果有root权限，可以使用docker安装。我这里直接使用的师姐服务器上的docker（80服务器）。

法一：
conda create -n deepvariant python=3.6 -c bioconda pbmm2 deepvariant -y
法二：
BIN_VERSION="1.4.0"
sudo apt -y update
sudo apt-get -y install docker.io
sudo docker pull google/deepvariant:"${BIN_VERSION}"

二、使用

STEP1 pbmm2比对

****************************
pbmm2 index ref.fa ref.mmi
pbmm2 align ref.mmi movie.subreads.bam ref.movie.bam
samtools sort -@8 ref.movie.sorted.bam ref.movie.bam

#创建deepvariant需要的索引文件
samtools index -@8 ref.movie.sorted.bam
samtools faidx ref.fa
***************************

STEP02 deepvariant识别变异

该步骤需要注意，在生成的quickstart-output目录下存放所有的输入文件，包括ref.fa ref.fai ref.movie.sorted.bam ref.movie.sorted.bam.bai。（其中-v参数应该是设置输入输出容器以input代替输入输入路径，output代替输出路径））

docker运行包括三个阶段： make_examples, call_variants, and postprocess_variants。

程序并行运行，若需要中止，只需kill掉其中一个进程便会全部中止。

**************************************************
OUTPUT_DIR="/home/liulu/liuxin/quickstart-output"
INPUT_DIR="/home/liulu/liuxin/quickstart-output"
mkdir -p "${OUTPUT_DIR}"
BIN_VERSION="1.4.0"

sudo docker run \
-v "${INPUT_DIR}":"/input" \
-v "${OUTPUT_DIR}":"/output" \
google/deepvariant:"${BIN_VERSION}" \
/opt/deepvariant/bin/run_deepvariant \
--model_type=PACBIO \ `
--ref=/input/mav4.fa \
--reads=/input/mav4.BXJ.sorted.bam \ 
--output_vcf=/output/output.vcf.gz \
--output_gvcf=/output/output.g.vcf.gz \
--intermediate_results_dir /output/intermediate_results_dir \ 
--num_shards=48 
**************************************************

STEP03 gatk分离SNP以及INDEL

从上一步输出的vcf文件中分离出SNP和INDEL

gatk  SelectVariants -R /input/YOUR_REF  \   
    -V /input/YOUR_VCF  \
    -O /output/YOUR_RAW_SNP  \
    -select-type SNP|INDEL

STEP04 vcftools

使用vcftools对vcf文件中的SNP进行等位基因频率统计

./vcftools --vcf input_data.vcf --freq --out output

STEP05 R数据展示

根据freq的文件做直方图统计展示，比较简单。