生信软件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