olavurmortensen/linkseq

LinkSeq — 适配linked-reads的GATK最佳实践 pipeline

https://img.shields.io/badge/Docker%20build-Available-informational](https://hub.docker.com/repository/docker/olavurmortensen/linkseq)

目录

• 概述
• 工作流
• 安装
• 使用示例
  • 运行LinkSeq
  • 调试说明
• 参考资源
  • Barcode白名单
  • GATK资源
  • 外显子测序目标区域
  • SnpEff数据

概述

本 pipeline 对 10x Genomics 的 linked-reads 进行比对，并使用GATK进行变异检测。具体而言，根据 https://software.broadinstitute.org/gatk/best-practices/workflow?id=*** 执行生殖系短变异（SNP和插入缺失）发现。

该 pipeline 基于 https://www.nextflow.io/ 编写，主要步骤如下：

• 使用 https://github.com/arshajii/ema/ 比对 reads 并校准BAM文件（BQSR）
• 使用GATK的HaplotypeCaller生成GVCF（可用于联合基因分型）
• 使用GATK的GenotypeGVCFs对GVCF进行基因分型，生成单样本VCF
• 使用SnpEff注释变异效应并过滤变异
• 使用 https://github.com/vibansal/HapCUT2 对VCF进行相位分析
• 使用 https://whatshap.readthedocs.io/en/latest/index.html 将相位信息附加到BAM文件
• 使用GATK的VariantEval进行变异质量控制
• 使用 Qualimap 进行BAM质量控制
• 使用 MultiQC 生成综合质量控制报告

工作流

LinkSeq有多个关联 pipeline，以下是它们的配合关系：

• https://github.com/olavurmortensen/linkseq-demux%EF%BC%9A%E5%AF%B9%E5%8E%9F%E5%A7%8B%E6%B5%8B%E5%BA%8F%E6%95%B0%E6%8D%AE%E8%BF%9B%E8%A1%8C%E7%A2%B1%E5%9F%BA%E8%AF%86%E5%88%AB%E3%80%81%E8%A7%A3%E5%A4%8D%E7%94%A8%E5%92%8Creads%E4%BF%AE%E5%89%AA
• LinkSeq：对每个样本进行序列比对、变异检测、VCF和BAM相位分析
• https://github.com/olavurmortensen/gatk-joint-genotyping%EF%BC%9A%E5%AF%B9%E6%89%80%E6%9C%89%E6%A0%B7%E6%9C%AC%E8%BF%9B%E8%A1%8C%E8%81%94%E5%90%88%E5%9F%BA%E5%9B%A0%E5%88%86%E5%9E%8B
• https://github.com/olavurmortensen/linkseq-phase%EF%BC%9A%E5%AF%B9%E5%A4%9A%E6%A0%B7%E6%9C%AC%E8%81%94%E5%90%88%E5%9F%BA%E5%9B%A0%E5%88%86%E5%9E%8B%E7%9A%84VCF%E8%BF%9B%E8%A1%8C%E7%9B%B8%E4%BD%8D%E5%88%86%E6%9E%90

安装

您有两种选择：(1) 使用conda安装依赖（如下所述）；(2) 使用 https://hub.docker.com/repository/docker/olavurmortensen/linkseq Docker镜像。

使用environment.yml文件通过conda安装依赖：

bash
conda env create -f environment.yml

激活环境（环境名称在environment.yml中定义，通常为linkseq）：

bash
conda activate linkseq

使用nextflow拉取项目：

bash
nextflow pull https://github.com/olavurmortensen/linkseq

使用示例

我们将在10x Genomics的"tiny-bcl"示例数据集上运行linkseq。运行前，请先执行linkseq-demux pipeline（https://github.com/olavurmortensen/linkseq-demux%EF%BC%89%E7%9A%84%E7%A4%BA%E4%BE%8B%EF%BC%8C%E5%AE%8C%E6%88%90%E5%8E%9F%E5%A7%8B%E5%BA%8F%E5%88%97%E7%9A%84%E7%A2%B1%E5%9F%BA%E8%AF%86%E5%88%AB%E3%80%81%E8%A7%A3%E5%A4%8D%E7%94%A8%E5%92%8C%E4%BF%AE%E5%89%AA%E3%80%82

关于 pipeline 使用的参考数据，请参见下文的参考资源部分。

运行LinkSeq

最简便的运行方式是为nextflow定义配置文件。以下配置文件定义了 pipeline 的输入文件和运行时设置：

groovy
params {
    sample = "Sample1"
    
    // 输入的FASTQ文件
    fastq_r1 = "tiny-fastq/Sample1/fastqs/Sample1_L005_R1*fastq.gz"
    fastq_r2 = "tiny-fastq/Sample1/fastqs/Sample1_L005_R2*fastq.gz"
    
    // 参考数据
    whitelist = "resources/whitelist/4M-with-alts-february-2016.txt"
    reference = "resources/gatk_bundle/Homo_sapiens_assembly38/Homo_sapiens_assembly38.fasta"
    dbsnp = "resources/gatk_bundle/Homo_sapiens_assembly38.dbsnp138/Homo_sapiens_assembly38.dbsnp138.vcf"
    targets = "resources/sureselect_human_all_exon_v6_utr_grch38/S07604624_Padded_modified.bed"
    snpeff_datadir = "resources/snpeff_data"
    
    // EMA比对使用的bin数量
    bcbins = 10
    
    outdir = "outs"
}

// 每个进程可用的资源
process {
    executor = 'local'
    memory = '10GB'
    cpus = 1 
}

// pipeline可用的总资源
executor {
    name = 'local'
    cpus = 10
    memory = '100GB'
    queueSize = 100 
}

// 管道执行时捕获上游进程的退出码
process.shell = ['/bin/bash', '-euo', 'pipefail']

注意：EMA比对的数量设置为10仅适用于小型示例数据集。对于全基因组数据，EMA开发者建议使用500个。

使用以下命令运行 pipeline：

bash
nextflow run olavurmortensen/linkseq -resume -with-trace

其中-with-trace用于生成进程进度日志，-resume用于在 pipeline 失败后从断点继续运行。

运行完成后，可通过tree -L 3 outs/查看输出结构：

bash
outs/
├── Sample1
│   ├── bam
│   │   ├── Sample1.bam
│   │   ├── Sample1.bam.bai
│   │   └── bx_stats.csv
│   ├── gvcf
│   │   ├── gvcf.g.vcf
│   │   └── gvcf.g.vcf.idx
│   └── vcf
│       ├── Sample1.vcf.gz
│       ├── Sample1.vcf.gz.tbi
│       └── phasing
├── multiqc
│   ├── multiqc_data
│   │   ├── multiqc.log
│   │   ├── multiqc_data.json
│   │   ├── multiqc_fastqc.txt
│   │   ├── multiqc_gatk_varianteval.txt
│   │   ├── multiqc_general_stats.txt
│   │   ├── multiqc_qualimap_bamqc_genome_results.txt
│   │   ├── multiqc_snpeff.txt
│   │   └── multiqc_sources.txt
│   └── multiqc_report.html
└── multiqc_logs
    ├── AnalyzeCovariates
    │   └── Sample1
    ├── SnpEff
    │   └── Sample1
    ├── VariantEval
    │   └── Sample1
    ├── WhatsHap
    │   └── Sample1
    ├── bqsr_after
    │   └── Sample1
    ├── bqsr_before
    │   └── Sample1
    ├── bx_stats
    ├── fastqc
    │   └── Sample1_fastqc.zip
    └── qualimap
        └── Sample1

调试说明

若 pipeline 运行失败需调试，可查看进程跟踪日志：

bash
less trace.txt

或检查work/目录下的文件：

bash
tree work/ | less

参考资源

Barcode白名单

10X Genomics GemCode技术的barcode白名单可从LongRanger软件包获取： longranger-2.y.z/longranger-cs/2.y.z/tenkit/lib/python/tenkit/barcodes/4M-with-alts-february-2016.txt

LongRanger可从10X Genomics官网下载： [***]

也可直接从以下链接下载：

bash
http://cb.csail.mit.edu/cb/ema/data/4M-with-alts-february-2016.txt

GATK资源

运行 pipeline 需要 https://software.broadinstitute.org/gatk/download/bundle 中的参考数据库及外显子目标区域文件。

reference/gatk_bundle.sh脚本可下载所需资源；注意若GATK团队更新资源，该脚本可能失效。我们于2019年3月22日从以下站点下载资源：

GATK资源包
https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0

目标区域BED文件

运行 pipeline 时应始终使用目标区域BED文件。若进行全基因组测序，可使用GATK资源包中的WGS调用区域文件：

GATK资源包
https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle
WGS调用区域hg38
https://storage.cloud.google.com/genomics-public-data/resources/broad/hg38/v0/wgs_calling_regions.hg38.interval_list

本示例使用Agilent SureSelect Human All Exon V6 UTR试剂盒捕获外显子，reference/sureselect_human_all_exon_v6_utr_grch38目录包含示例BED文件及说明文档，可用于调试BED文件兼容性问题。

SnpEff数据

使用reference/snpeff_data.sh脚本下载hg38的SnpEff数据库。在LinkSeq中，该数据库路径与SnpEff的dataDir参数配合使用。

也可通过snpeff databases | grep hg38查看数据库名称。

在Docker容器中运行 pipeline 时，预下载该数据库可节省时间并避免下载失败导致的 pipeline 崩溃。