生信软件39 - GATK最佳实践流程重构，提高17倍分析速度的LUSH流程

1. LUSH流程简介

基因组测序通常用于分子诊断、分期和预后，而大量测序数据在分析时间方面提出了挑战。

对于从FASTQ到VCF的整个流程，LUSH流程在非GVCF和GVCF模式下都大大降低了运行时间，30 X WGS数据耗时不到2 h，从BAM到VCF约需12分钟，比GATK流程快约17倍，比GATK-Spark流程快4.5倍；在准确性方面，LUSH和GATK的结果同样准确且高度一致。这对于儿科重症监护室（PICU）和新生儿重症监护室（NICU）中的婴儿等急性病患者至关重要。

LUSH可以有效地加速所有类型的DNAseq数据（如WGS，WES，PANEL等）的分析。由LUSH生成的BAM文件（遵循原始的BWA算法）也可用于随后的结构变异（SV）和拷贝数变异（CNV）的检测，用作常见SV调用者如Manta、Delly 和Lumpy 等软件的输入。

论文: Fast and accurate DNASeq variant calling workflow composed of LUSH toolkit

https://humgenomics.biomedcentral.com/articles/10.1186/s40246-024-00666-w

github: https://github.com/Bgi-LUSH/LUSH-DNASeq-pipeline

2. LUSH DNASeq工作流程

LUSH DNASeq工作流程是基于GATK最佳实践的优化管道，由LUSH_Aligner、LUSH_BQSR、LUSH_HC和LUSH_GenotypeGVCF组成。

工作流程

LUSH_Aligner集成了多个功能模块，如SOAPnuke，Bwa MEM，Samtools sort和GATK-MarkDuplicates（Picard），同时基于原始算法进行了完全重新开发。

2.1 LUSH_Aligner

LUSH_Aligner由三个主要功能模块组成：“FqFilterAlignent”、“bwaMEM Alignent”和“SortDuplicateAlignent”，分别进行fastq过滤、比对齐以及BAM排序和标记重复任务。

2.2 LUSH_BQSR

测序仪产生的碱基质量分数受到各种系统性技术误差的影响，导致获得的质量分数过高或过低。碱基质量分数矫正（Base Quality Score Recalibration， BQSR）是使用机器学习方法对这些错误进行经验建模，并调整质量分数，主要包括2个步骤：Base Recalibration and applying BQSR。

LUSH_BQSR实现了生产者-消费者并行计算结构，以优化并行任务并提高CPU利用率，减少了冗余IO消耗，最终提高了处理速度。

2.3 LUSH_HC和LUSH_GenotypeGVCF

HaplotypeCaller采用区域内的局部从头组装方法，用于准确检测单核苷酸多态性（SNP）和小的插入缺失（Indels），主要包括4个连续的步骤：识别区域、区域内局部组装推断单倍型、估计可能值和利用隐马尔可夫模型（HMM）在贝叶斯推理的基础上确定基因型。LUSH_GenotypeGVCF（0.21小时）比GATK-GenotypeGVCF（1.13小时）快5倍。

LUSH_HC采用了进一步的任务细分策略，结合资源动态分配，以实现负载平衡并优化资源分配。

当以12个线程运行时，流程在约4.89小时内完成，当以56个线程运行时，流程在约1.6小时内完成，这表明LUSH流水线具有很大的线程可扩展性。

3. LUSH管道与GATK和GATK-Spark管道的准确性

LUSH的底层算法与GATK或GATK-Spark大致相同，因此它们预计会产生相同的结果。

LUSH和GATK管道的非GVCF模式和GVCF模式的结果比较表明，前者表现出更高的准确率和略低的召回率。F1分数显示非GVCF模式在两个流水线的准确性方面表现出更好的性能，与两个家系WGS数据完全一致。

4. 使用方法

# 克隆github项目
git clone https://github.com/Bgi-LUSH/LUSH-DNASeq-pipeline

4.1 lush_aligner构建参考基因组索引

# 构建hg19参考基因组序列索引
./bin/LUSH_toolkit-Aligner/lush_aligner index /path/hg19.fa

4.2 lush_aligner执行fastq过滤+比对+排序+标记重复

# 创建结果目录
mkdir -p ./outdir/clean_data

./bin/LUSH_toolkit-Aligner/lush_aligner filter4mem \
        -6 ./outdir/ \
        -n 0.1 -J 0.5 -l 12 -g 2 -b 2 -t 20 -M \
        -r /path/hg19.fa \
        # 输出排序+标记重复bam
        -o ./outdir/NA12878.sort.dup.bam \
        -Z ./outdir/clean_data \
        # fastq路径及-R参数配置文件
        -i ./example_data/lush.config

./example_data/lush.config内容：

./example_data/NA12878_l01_1.fq.gz  NA12878_l01_1       @RG\tID:NA12878.1\tLB:LibA\tSM:NA12878\tPL:COMPLETE\tCN:BGI
./example_data/NA12878_l01_2.fq.gz  NA12878_l01_2
./example_data/NA12878_l02_1.fq.gz  NA12878_l02_1       @RG\tID:NA12878.2\tLB:LibA\tSM:NA12878\tPL:COMPLETE\tCN:BGI
./example_data/NA12878_l02_2.fq.gz  NA12878_l02_2

4.3 lush_bqsr执行碱基质量矫正和ApplyBQSR

export LD_LIBRARY_PATH=./bin/LUSH_toolkit-BQSR:$LD_LIBRARY_PATH

./bin/LUSH_toolkit-BQSR/lush_bqsr \
 --bam_path /INPUT_PATH/NA12878.sort.dup.bam \
 --out_dir ./outdir/LUSH_BQSR  \
 --plugin_path ./bin/LUSH_toolkit-BQSR/libbqsr.so \
 --producer_number 2 \
 --worker_number 21 \
 --fasta /path/hg19.fa \
 # 金标准indels vcf文件
 --known_site Mills_and_1000G_gold_standard.indels.hg19.vcf \
 --writer_thread 5 \
 --pr_one_bam 1

4.4 lush_hc执行GenotypeGVCFs变异检测

export LD_LIBRARY_PATH=./bin/LUSH_toolkit-HC:$LD_LIBRARY_PATH

./bin/LUSH_toolkit-HC/lush_hc HaplotypeCaller \
        --pcr-indel-model NONE \
        -I /INPUT_PATH/NA12878.sort.dup.bam \
        -R hg19.fa \
        -O ./outdir/NA12878.vcf.gz

4.5 LUSH_GenotypeGVCFs

LUSH_GenotypeGVCFs是GATK GenotypeGVCFs功能使用C/C++的重新实现。

# UASGE: 
# LUSH_GenotypeGVCF inputGvcfFile outputVcfFile stand-call-conf

# inputGvcfFile   input VCF file
# outputGvcfFile  output file name:/file/NA12878_PCR.vcf.gz
# stand-call-conf The minimum phred-scaled confidence threshold at which variants # should be called:10.0

# 示例
export LD_LIBRARY_PATH=./bin/LUSH_toolkit-GenotypeGVCFs:$LD_LIBRARY_PATH

./bin/LUSH_toolkit-GenotypeGVCFs/lush_genotypegvcfs \
INPUT_PATH/NA12878.g.vcf.gz \
./outdir/NA12878.vcf.gz 10

5. GATK管道Shell脚本

Usage:
  GATK_pipeline.sh [-i FQFile] [-t THREAD] [-o OUTDIR] [-m MODEL] [-s PREFIX] [-p SPARK]

Description:
    FQFile, the path of INPUT fastq file, should be like '/path/fastq1,/path/fatq2'
    THREAD, the number of thread [10]
    OUTDIR, the path of outdir [./]
    PREFIX, the prefix of outputfile [GATKtest]
    MODE, GVCF or not [Y/N]
    SPARK, Spark or not [Y/N]

GATK_pipeline.sh \
-i /PATH/MGISEQ2000_PCR-free_NA12878_30X_1.fq.gz,/PATH/MGISEQ2000_PCR-free_NA12878_30X_2.fq.gz \
-t 40 \
-o ./  \
-m N \
-s samplename \
-p N

6. 比较LUSH和GATK管道获得的变异VCF文件

Haplotype_Comparison.sh \
-i LUSHtest.vcf.gz,GATKtest.vcf.gz \
-t 40 \
-o ./ -s sample

生信软件文章推荐

生信软件1 - 测序下机文件比对结果可视化工具 visNano

生信软件2 - 下游比对数据的统计工具 picard

生信软件3 - mapping比对bam文件质量评估工具 qualimap

生信软件4 - 拷贝数变异CNV分析软件 WisecondorX

生信软件5 - RIdeogram包绘制染色体密度图

生信软件6 - bcftools查找指定区域的变异位点信息

生信软件7 - 多线程并行运行Linux效率工具Parallel

生信软件8 - bedtools进行窗口划分、窗口GC含量、窗口测序深度和窗口SNP统计

生信软件9 - 多公共数据库数据下载软件Kingfisher

生信软件10 - DNA/RNA/蛋白多序列比对图R包ggmsa

生信软件11 - 基于ACMG的CNV注释工具ClassifyCNV

生信软件12 - 基于Symbol和ENTREZID查询基因注释的R包（easyConvert ）

生信软件13 - 基于sambamba 窗口reads计数和平均覆盖度统计

生信软件14 - bcftools提取和注释VCF文件关键信息

生信软件15 - 生信NGS数据分析强大的工具集ngs-bits

生信软件16 - 常规探针设计软件mrbait

生信软件17 - 基于fasta文件的捕获探针设计工具catch

生信软件18 - 基于docker部署Web版 Visual Studio Code

生信软件19 - vcftools高级用法技巧合辑

生信软件20 - seqkit+awk+sed+grep高级用法技巧合辑

生信软件21 - 多线程拆分NCBI-SRA文件工具pfastq-dump

生信软件22 - 测序数据5‘和3‘端reads修剪工具sickle

生信软件23 - Samtools和GATK去除PCR重复方法汇总

生信软件24 - 查询物种分类学信息和下载基因组TaxonKit和ncbi-genome-download

生信软件25 - 三代测序数据灵敏比对工具ngmlr

生信软件26 - BWA-MEM比对算法性能更好的bwa-mem2

生信软件27 - 基于python的基因注释数据查询/检索库mygene

生信软件28 - fastq与bam的reads数量计算与双端fastq配对检测工具fastq-pair

生信软件29 - 三代数据高效映射精确的长读段比对工具mapquik

生信软件30 - 快速单倍型分析工具merlin

生信软件31 - Bcftools操作VCF/BCF文件高级用法合集

生信软件32 - 变异位点危害性评估预测工具合集

生信软件33 - Wgsim生成双端(PE) fastq模拟数据

生信软件34 - 大幅提升Python程序执行效率的工具Pypy

生信软件35 - AI代码编辑器Cursor

生信软件36 - SAM/BAM/CRAM文件插入SNV/INDEL/SV工具Bamsurgeon

生信软件37 - 基于测序reads变异进行单倍型分型工具WhatsHap

生信软件38 - 基因型填充软件IMPUTE2