使用refs获取节点_使用Gatk Germline spns-indels Pipeline分析遗传病(耳聋)

最新推荐文章于 2024-12-10 15:42:57 发布

原创最新推荐文章于 2024-12-10 15:42:57 发布 · 145 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#使用refs获取节点

本文是GATK Best Practice系列文章，尝试用Gatk Germline spns - indels Pipeline分析耳聋遗传病。从NCBI获取数据，介绍了分析系统、流程文件、变量，阐述分析流程，最后得到结果表格，分析结果准确度尚可，但存在因bed文件不精确导致QC数值不准等不足。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

这是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