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的文件做直方图统计展示,比较简单。