使用Gatk Germline spns-indels Pipeline分析遗传病(耳聋)
这是GATK Best Practice系列学习文章中的一篇，本文尝试使用Gatk Germline spns-indels Pipeline来分析遗传病（耳聋）
数据
这次没有拿到遗传病的室间质评的数据，直接从NCBI上找一些数据来分析。NCBI上搜索deaf，点击第一条搜索结果，最后几经跳转找到数据下载页面：https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP218677
可以看到：Targeted next generation sequencing of 139 known deafness-related genes in 44 deaf patients with mono-allelic GJB2 mutations
这里一共44个样本，有兴趣的可以全部下载来分析。
奇怪的是，这里不是sra下载，是直接下载fastq文件，不管了。使用下载工具，最后得到的文件是C261 R1.fastq.1这种文件名，这里首先改名为C261_R1.fastq，然后使用gzip，bgzip压缩得到 C261_R1.fastq.gz；C261_R2.fastq.gz也是同样.
也可以直接从一下链接下载：
C261_R1.fastq.gz bc72a0b02876273e2457c753a87241aa 54M
C261_R2.fastq.gz a5f6fcac16cbbe2b2a0c98982243a949 54M
概览：
用到的分析系统及分析流程文件
| 名称 （点击下载） | 备注 |
|---|---|
| SliverWorkspace 2.0.290364 | 提供运行控制平台/社区版 |
| GATK Germline SNP/Indel V1.0 |
分析流程文件，可以一键导入分析平台(点击查看操作) 当然可以参照图片中运行脚本，shell里运行，效果也是一样 |
| 最终结果过滤脚本（python2.7 )及编译版本 | hg19.exon.bed 用到的bed，intelval文件，因为没有原作者的150个基因的bed，这里用一个全外的替代，所以最后的QC结果计算也是不准确的 GermlineSnvAnnotationFilter.py 过滤脚本 GermlineSnvAnnotationFilter 使用pyinstaller编译的直接可执行文件 GermlineQcProcessor.py 获取整体QC数据的脚本 GermlineQCProcessor 使用pyinstaller编译的直接可执行文件 blacklist.csv 记录了耳聋相关的基因和部分突变位点，用来过滤耳聋相关结果 |
| 分析结果 | result.zip分析结果 |
流程用到的变量
| 变量名 | 变量值 | 类型 |
|---|---|---|
| refs.blacklist | /opt/ref/deaf/blacklist.csv | 文件 |
| tools.fastp | /opt/ref/fastp | 程序 |
| tools.samtools | /opt/samtools/samtools | 程序 |
| tools.java | /opt/jdk1.8.0_162/bin/java | 程序 |
| tools.bwa | /opt/bwa/bwa | 程序 |
| tools.sambamba | /opt/ref/sambamba-0.7.0-linux-static | 程序 |
| tools.gatk | /opt/ref/gatk-4.1.3.0/gatk | 程序 |
| refs.gene | /opt/ref/hg19/hg19_refGene.txt | 文件 |
| refs.dict | /opt/ref/hg19/ucsc.hg19.dict | 文件 |
| refs.hum | /opt/ref/hg19/ucsc.hg19.fa | 文件 |
| refs.interval | /opt/ref/projects/hg19.exon.interval_list | 文件 |
| refs.dbsnp | /opt/ref/hg19/dbsnp_138.hg19.vcf | 文件 |
| refs.mills | /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.vcf | 文件 |
| refs.1000G | /opt/ref/hg19/1000G_phase1.indels.hg19.vcf | 文件 |
| refs.af_only | /opt/ref/hg19/af-only-gnomad.raw.sites.hg19.vcf.gz | 文件 |
| refs.bed | /opt/ref/projects/hg19.exon.bed | 文件 |
| refs.small.exac | /opt/ref/hg19/small_exac_common_3_b37.vcf | 文件 |
| envis.threads | 32 | 数值 |
| envis.memory | 32G | 字符 |
分析流程
-
01 INPUT 输入文件
-
02 fastp,因为数据没有去接头，这里用fastp直接去接头，默认参数处理
-
03 Map Sort Dedup 一套操作
-
04-Recalibrator
-
05 Process Dup Bam
-
06-Call SNV INDEL
-
07-Hard Filter
因为call出来的突变数量不够运行VQSR，所以这里选择使用硬过滤（Hard Filter），
过滤参数使用gatk推荐的数值，
见链接#https://gatk.broadinstitute.org/hc/en-us/articles/360035532412?id=11097
-
08-Annotation，使用Annovar注释结果
-
09-Annotation Filter对注释结果过滤，过滤工具使用python开发，见文章头部链接
- 初步过滤：过滤掉突变类型为NAN和synonymous SNV的结果
- 初步过滤：过滤掉突变位点深度低于10，或者reads数少于10的结果，因为遗传突变要么是50%要么100%，低于一定reads数该结果的可靠性不高。这里暂时定为10.
- 匹配blacklist.csv文件中记录的和耳聋相关的基因，过滤掉无关结果
-
10-QC
获取整体数据的QC状态，平均覆盖深度，总的reads数，等等。
-
11-Output 最后的输出文件：
-
12-分析结果匹配：
最后我们得到一个结果表格，根据Clinvar数据库的注释来判断的话，只有图中黄色一行的状态为：
Pathogenic/Likely_pathgenic(致病的/可能致病的)；
其余行记录均为：NAN或者Benign/Like_benign（良性的/可能良性的）
这与下载数据的内容相匹配：Targeted next generation sequencing of deaf patients with mono-allelic GJB2 mutations
说明我们的分析结果准确度还可以。
不足
- 因为未能下载到原数据用到150个和耳聋相关基因的bed文件，这里使用的是hg19全外的bed文件，所以QC工具最后计算的整体数据的QC数值不准确，平均测序深度、测序深度中位数值、1X，4X等深度下的覆盖度也很低，insert size更是不准确。这只有获得精确的bed文件之后可以得到改进。
- 因为流程分析中使用的是一个hg19全外的bed文件，未包含耳聋中部分线粒体突变的区域。这里最终结果不包含线粒体上的耳聋致病突变，
- 分析流程中使用的一个数据文件，blacklist.csv是搜集到的一些和耳聋相关的基因和突变位点，过滤的时候只匹配了基因。数据还相对比较粗糙，如果改进的话改进数据结构和提高数据的质量：如使用这个网站的数据https://hereditaryhearingloss.org/ 或者搜集更高质量的数据，注释时候关联更多的数据库如OMIM