GATK流程学习

GATK pipeline 学习。本文结合gatk流程应用于检查GapFilling区域的正确性。

下载及安装细节参见GATK官网
流程理解请参考https://www.jianshu.com/p/859c0345624c

一、WGS分析流程
主要包括三个部分：原始数据质控、数据预处理以及变异检测。

二、软件安装
BWA: NGS数据比对
Samtools: 处理比对数据工具，如格式转换sam->bam
Picard: NGS数据处理工具，可用与去除重复比对，避免假阳性
GATK: call snp

PS: bwa、samtools克隆下载需编译，picard采用conda安装，gatk克隆下载后开箱即用，此处具体过程省略，非conda安装别忘了配置环境变量。

三、流程构建

#原始数据质量控制
待补充

1、建立索引库（以便快速比对）

$ bwa index hifi_fill.scaff_seqs.fasta

2、序列比对

$ bwa mem -t 48 -R '@RG\tID:unknown\tPL:illumina\tLB:wgs\tSM:cbp' hifi_fill.scaff_seqs.fasta D2124163A_122_1.fq D2124163A_122_1.fq | samtools view -@ 48 -bS > cbp.bam 

3、序列排序

$ samtools sort -@ 48 -o cbp.sorted.bam cbp.bam

4、去除PCR重复

$ picard MarkDuplicates -Xmx64g I=cbp.sorted.bam O=cbp.sorted.markdup.bam M=cbp.sorted.markdup_metics.txt REMOVE_DUPLICATES=true        
$ picard BuildBamIndex -Xmx64g I=cbp.sorted.markdup.bam

5、创建索引

$ samtools index cbp.sorted.markdup.bam

6、识别变异

$ gatk HaplotypeCaller -R hifi_fill.scaff_seqs.fasta -I cbp.sorted.markdup.bam -O cbp.vcf

PS: 此处报错，原因应该是在第一步创建索引库时缺少两条命令，补充完成即可:  

$ gatk CreateSequenceDictionary -R hifi_fill.scaff_seqs.fasta 
$ samtools faidx hifi_fill.scaff_seqs.fasta

PS: 此外还可使用带Spark模块，可利用所有可用线程运行.

7、gatk过滤

Filter Variants by Hard-filter (for SNPs only)

$ gatk SelectVariants -R hifi_fill.scaff_seqs.fasta -V raw.cbp.vcf -select-type SNP --restrict-alleles-to BIALLELIC -O /home/liuxin/output/filtering/snps.vcf
$ gatk VariantFiltration -R hifi_fill.scaff_seqs.fasta -V /home/liuxin/output/filtering/snps.vcf -O /home/liuxin/output/filtering/snps.filtered.vcf --filter-name "QD_filter" -filter "QD < 2.0" --filter-name "FS_filter" -filter "FS > 60.0" --filter-name "MQ_filter" -filter "MQ < 40.0" --filter-name "SOR_filter" -filter "SOR > 3.0" --filter-name "MQRankSum_filter" -filter "MQRankSum < -12.5" --filter-name "ReadPosRankSum_filter" -filter "ReadPosRankSum < -8.0"

 Exclude filtered variants

$ gatk SelectVariants --exclude-filtered -V /home/liuxin/output/filtering/snps.filtered.vcf -O /home/liuxin/output/filtering/snps.filtered.pass.vcf

四、查看统计

1、查看vcf统计信息

$ vcftools --vcf snps.filtered.pass.vcf --counts
#After filtering,kept 1 out of 1 Individuals
#After filtering,kept 12114 out of a possible 12114 Sites
$ vcftools --vcf snps.filtered.pass.vcf --depth -c > depth_summary.txt
#INDV N_SITES MEAN_DEPTH
#cbp 12114 124.526

2、检查GapFilling区域的snp分布情况

3、统计（cbp.sorted.markdup.bam/snps.filtered.pass.vcf作为输入）

###平均深度统计初步查看
############################################################################################
vcftools --vcf snps.filtered.pass.vcf --site-depth --out all
cut -f3 all.ldepth > all.site.depth
mawk '!/D/' all.site.depth | mawk -v x=1 '{print $1/x}' > all.meandepthpersite
#1为个体数
gnuplot << \EOF 
set terminal dumb size 120, 30
set autoscale
set xrange [10:150] 
unset label
set title "Histogram of mean depth per site"
set ylabel "Number of Occurrences"
set xlabel "Mean Depth"
binwidth=1
bin(x,width)=width*floor(x/width) + binwidth/2.0
set xtics 5
plot 'all.meandepthpersite' using (bin($1,binwidth)):(1.0) smooth freq with boxes
pause -1
EOF
##############################################################################################

#统计Mapping

$ samtools flagstats -@ 48 cbp.sorted.markdup.bam > bamstats.txt

#统计Coverage

$ samtools coverage cbp.sorted.markdup.bam
#16条染色体的整体mapping率均>99%，但igvviewer结果显示补洞区域的覆盖度并不好。这可能与二代区域在gap区域本身的覆盖度较差有关。但HiFi补洞的区域显示了一定的覆盖度>3x。

五、尝试用hifi reads再走一遍mapping流程。

pacbio片段比对，不存在pcr偏好性无需去重，只需排序，随后将bam结果文件、bam索引以及genome序列输入igv-viewer查看mapping结果，重点关注GapFill区域的覆盖度。

bwa index hifi_fill.scaff_seqs.fasta
bwa mem -t 48 -x pacbio /home/liuxin/input/hifi_fill.scaff_seqs.fasta /home/liuxin/input/cbp.ccs.fasta | samtools view -@ 48 -bS > hifi_aln.bam 
samtools sort -@ 48 -o hifi_aln.sorted.bam hifi_aln.bam
samtools index hifi_aln.sorted.bam 

比对结果使用IGV-Viewer查看GapFilling区域的覆盖度。