Run GATK Best Practices for data pre-processing by Cromwell/WDL

与GATK4正式发布的还有WDL（workflow description langaue，https://software.broadinstitute.org/wdl/），WDL将工作流程分为了workflow, task, call, command 和 output。

与以往GATK提供Best practice的PPT介绍不同，现在Broad提供的是Best practice（https://software.broadinstitute.org/gatk/best-practices/）的WDL文件。WDL文件运行通过cromwell运行，并且有json格式的参数输入文件指定WDL文件中流程所需要的参数。比如

sudo java -jar cromwell.jar run workflow.wdl --inputs workflow.inputs.json

我们只需要修改json文件中的参数就可以运行gatk4 Best Practices，而不需要自己去搭建流程，简化了工作量，也遵循了Broad提供的推荐设置和流程。本文只介绍突变检测前的序列比对和recalibrate这部分的GATK best practices,该流程生成了用于variant calling的bam文件。

1，文件准备

WDL文件和json文件 Broad在github上提供了进行突变检测call variant之前的数据处理data proceesing流程，见https://github.com/gatk-workflows/gatk4-data-processing 从github上，我们需要下载两个文件 processing-for-variant-discovery-gatk4.wdl （用于data pre-processing 的 pipeline） processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json（指定WDL的参数文件）

ubam文件： 要求ubam文件中要有RG tag，经过排序sort之后，该文件可以通过picard将fastq文件转换得到

GATK resoure bundle，从中下载GATK需要的dbsnp文件，known site等文件 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/ https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0/

有幸跟同事去蹭了中国遗传学协会遗传咨询分会CBGC在苏州召开的CBGC年会，也第一次来苏州。

关于会议

遗传学协会遗传咨询分会一致推动这中国的遗传咨询事业的发展，从标准、体系、人才培养、临床应用，各位专家都倾注心血做事。我是做数据分析的，会接触后续的解读方面的工作，想谈谈解读和咨询这方面的看法。

1，解读不完全等于咨询：解读是遗传咨询过程中的一部分，利用遗传学和临床医学的知识，提供疾病和突变的关系，遗传咨询还涉及诊断前和诊断后与医生和患者之间的沟通，帮助医生确定疾病类型，指导患者优生优育，诊断和预后等。

2，表型很重要：包括遗传咨询和解读，尽可能详细准确的记录患者的表型，搜集患者的家族史信息等其他信息，可以使工作更加有目的性，成功概率大大提高。

3，HPO很好用：基因组内的基因和突变有千千万，如果采用普通的突变筛选方式，比如频率，风险预测等方法，依然会有很多候选突变。根据HPO，可以让我们重点关注特定基因panel，大大减少了候选突变的数目。

4，生信和解读：生信为解读提供线索，解读需要了解生信的线索。生信要尽可能的利用各种信息，缩小候选突变的范围，传统的根据特定指标就将突变过滤的方式太生硬，现在机器学习、大数据技术的发展，为生信的分析提供了一种思路，那就是可以不再依靠单一维度进行hard filter，而且通过多维度建模进行过滤分析。解读应了解生信的基本知识，这样才能更好的读懂生信提供的文件，理解生信的思路，帮助解读更加精准。或许有更好的模型出现，可以预测疾病风险。

5，基因检测不是万能药：现阶段，遗传咨询更多的是涉及罕见病、遗传病的诊断、优生优育等。很多疾病在医院都已经确定疾病类型了，并且有对应的指南进行操作，这个时候再做基因检测，个人觉得是浪费金钱。就好比基因检测和生化检测同样是一种技术和手段，是为了服务医生和患者的，当用其他手段确诊之后，再做基因检测有点浪费。如果疾病在某个人的家族中成家族性发生，又没有找到致病原因，那么可以做基因检测试下。如果医生对某个罕见病不能确定类型，可以从基因的角度试下。但患者也应了解，并不是做了基因检测，就会有针对应的指导就能找到致病原因，精准的指导和现阶段科学的发展有关。 6，再谈关联性位点：当大家都唾弃儿童天赋的时候，难道通过几个关联性风险位点预测癌症风险的项目就是好的？那么请告诉我，在那么多癌症关联位点中，就选选5个点的癌症风险预测和选10个点的癌症风险预测谁更准确？乳腺癌的BRAC1和2的基因检测难道测的是这里面的关联性分析位点？呵呵了，遗传咨询中前提是遗传，撇开遗传谈关联性分析得到的风险位点，都是耍流氓。不要拿几个GWAS位点，就吹的跟算命一样。 7，中国人群数据库：生信分析和数据解读很大程度上依赖数据库的支持，建立中国人群数据库，可以降低VUS意义不明的突变。政府、高校、社会应有意识的建立这些数据库，尽可能的设计完善，搜集多维度的数据，后续才能挖掘出更有价值的东西。

移动硬盘是ntfs格式的，服务器不能mount，报错 mount: unknown filesystem type ‘ntfs’

解决方法：安装 NTFS-3G，官网 https://www.tuxera.com/community/open-source-ntfs-3g

安装

1
2
3
4
5
6

wget -c https://tuxera.com/opensource/ntfs-3g_ntfsprogs-2017.3.23.tgz
tar -xvzf ntfs-3g_ntfsprogs-2017.3.23.tgz 
cd ntfs-3g_ntfsprogs-2017.3.23
./configure 
make 
sudo make install

挂载

1

mount -t ntfs-3g /dev/sd**  /target

Should I switch to a newer reference？

GRCh38 consists of several components: chromosomal assembly, unlocalized contigs (chromosome known but location unknown), unplaced contigs (chromosome unknown) and ALT contigs (long clustered variations). The combination of the first three components is called the primary assembly. It is recommended to use the complete primary assembly for all analyses.

参考基因组包括chromosomal assembly, unlocalized contigs (chromosome known but location unknown), unplaced contigs (chromosome unknown)和ALT contigs (long clustered variations)，前三个是primary assembly,alt contigs代表的是部分区域单体型的多样性，这些区域过于复杂不能用一条序列表示。

In addition to adding many alternate contigs, GRCh38 corrects thousands of SNPs and indels in the GRCh37 assembly that are absent in the population and are likely sequencing artifacts. It also includes synthetic centromeric sequence and updates non-nuclear genomic sequence.

除了添加了很多alternative contigs，GRCh38更正了数以千计的GRCh37版本中的SNPs和indels，这些SNPs和indels在人群中没有出现过、可能因为测序错误导致。GRCh38版本也包含了人工的着丝粒序列和更新了非核基因组序列。

The ability to recognize alternate haplotypes for loci is a drastic improvement that GRCh38 makes possible. Going forward, expanding genomics data will help identify variants for alternate haplotypes, improve existing and add additional alternate haplotypes and give us a better accounting of alternate haplotypes within populations. We are already seeing improvements and additions in the patch releases to reference genomes, e.g. the seven minor releases of GRCh38 available at the time of this writing. GRCh38大幅提高了识别alternate haplotype的能力，进一步提高了识别alternate haplotype的突变的能力。

BWA的作者Li Heng推荐GRCh37 primary assembly+ALT+decoy组成的参考基因组hs38DH，可以通过bwa下载，见https://github.com/lh3/bwa/blob/master/README-alt.md ：

bwa.kit/run-gen-ref hs38DH

作者的用NA12878做测试的比对结果如下，

1
2
3

Assembly	hs37	hs38	hs38DH
FP	255706	168068	142516
TP	2142260	2163113	2150844

在做基因组坐标转换的时候，用crossmap和liftover的时候，会用到chain file。大家都在讲坐标转换会用到这个文件，却没有讲过chain文件的具体内容（百度和谷歌搜索结果都没有中文介绍，本文应该是第一个）。本文的内容都翻译自UCSC网站，原文https://genome.ucsc.edu/goldenpath/help/chain.html，希望能帮到大家了解这个文件。

UCSC chain文件 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/liftOver/ Ensembl chain文件 https://sourceforge.net/projects/crossmap/files/Ensembl_chain_files/

chain file里面包含许多块alignment的信息（个人觉得可以理解为同源的地方，chain），其中每一块有一个header，记录alignment在两个版本中坐标，以及许多行alignment data line记录具体比对情况。

即每一块有 Header Line 和 Alignment Data Lines组成。形式如下，有两个chain

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

    chain 4900 chrY 58368225 + 25985403 25985638 chr5 151006098 - 43257292 43257528 1
     9       1       0
     10      0       5
     61      4       0
     16      0       4
     42      3       0
     16      0       8
     14      1       0
     3       7       0
     48

     chain 4900 chrY 58368225 + 25985406 25985566 chr5 151006098 - 43549808 43549970 2
     16      0       2
     60      4       0
     10      0       4