搜索
查看: 59|回复: 0

使用Gatk Germline spns-indels Pipeline分析遗传病(耳聋)

[复制链接]

13

主题

17

帖子

317

积分

中级会员

Rank: 3Rank: 3

积分
317
发表于 7 天前 | 显示全部楼层 |阅读模式
这是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  54MC261_R2.fastq.gz  a5f6fcac16cbbe2b2a0c98982243a949  54M
概览:
用到的分析系统及分析流程文件
名称 (点击下载)
备注

SliverWorkspace 2.0.290364提供运行控制平台/社区版
GATK Germline SNP/Indel V1.0分析流程文件,可以一键导入分析平台(点击查看操作) <br />当然可以参照图片中运行脚本,shell里运行,效果也是一样
最终结果过滤脚本(python2.7 )及编译版本hg19.exon.bed   用到的bed,intelval文件,因为没有原作者的150个基因的bed,这里用一个全外的替代,所以最后的QC结果计算也是不准确的 <br />GermlineSnvAnnotationFilter.py  过滤脚本 <br />GermlineSnvAnnotationFilter 使用pyinstaller编译的直接可执行文件<br />GermlineQcProcessor.py  获取整体QC数据的脚本 <br />GermlineQCProcessor 使用pyinstaller编译的直接可执行文件<br />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.threads32数值
envis.memory32G字符分析流程
  • 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推荐的数值,
  • 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


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-9-18 19:34 , Processed in 0.034874 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.