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

<think>好的,我现在需要帮助用户了解ChIP-seq的息学分析流程。用户已经提供了一些参考引用,里面提到了几个关键步骤和工具,比如比对到参考基因组、峰值调用,以及使用的工具如FastQC、Bowtie2、MACS2等。首先,我应该先回忆一下ChIP-seq的标准流程,然后结合用户提供的引用内容,整理出一个清晰的步骤说明。 首先,用户给的三个引用中,引用[1]详细列出了处理SRR14879780数据的步骤,包括环境准备、质量控制、比对、峰值调用和基序分析。引用[2]提到了ENCODE的标准流程,可能涉及数据标准和处理的不同阶段。引用[3]则讨论了数据可视化工具的选择,如IGV和UCSC Genome Browser。我需要将这些息整合成一个连贯的流程。 接下来,我需要确定回答的结构是否清晰。通常,ChIP-seq流程包括数据质控、预处理、比对、后处理、峰值检测、注释与富集分析、可视化、基序分析等步骤。用户的问题是关于步骤和方法的,所以每个步骤需要简要说明使用的工具和方法,并适当引用提供的资料。 需要注意的是,用户提到必须使用中文回答,并且行内数学公式用$...$,独立公式用$$...$$。不过ChIP-seq分析流程可能较少涉及数学公式,但如果有相关统计方法,可能需要提到,比如峰值调用中的统计模型。例如,MACS2使用泊松分布或负二项分布来识别显著富集区域,这时候可能需要公式。 另外,用户要求成相关问题,我需要根据回答内容提出几个相关问题,比如各个步骤的详细方法、工具比较、数据标准等。 现在需要检查是否有遗漏的步骤。引用[1]中的步骤包括准备环境、质量控制、去接头、比对、过滤排序、峰值调用、基序分析。引用[3]补充了可视化的部分。可能还需要数据下载和预处理,比如使用SRA Toolkit下载数据,但用户提供的引用中没有提到,可能需要假设用户已有数据。或者根据引用[1]中的步骤1,准备环境可能包括软件安装和数据准备。 另外,是否需要解释每个工具的作用?比如FastQC用于质控,Trim Galore用于去接头和裁剪,Bowtie2用于比对,samtools处理SAM/BAM文件,MACS2用于峰值检测,MEME-ChIP用于基序分析。这些都需要在步骤中说明。 在回答时,需要按照逻辑顺序排列步骤,并引用对应的引用文献。例如,引用[1]中的步骤1到7,引用[3]中的可视化部分。同时,确保每个步骤的说明简洁明了,避免过于技术化,但足够详细让用户理解流程。 最后,成相关问题需要考虑用户可能的后续疑问,比如各个步骤的详细参数设置、工具的比较(如为什么用Bowtie2而不用其他比对工具)、数据质控的标准、峰值调用的算法差异、可视化工具的选择等。 需要确保回答中没有错误,比如工具的正确使用顺序,以及引用对应的部分。例如,MACS2用于峰值调用是引用[1]提到的,而可视化部分引用[3]。同时,注意格式要求,如行内代码用反引号,步骤分点列出,使用中文等。 总结,回答的结构应该是: 1. 数据质控(FastQC) 2. 数据预处理(Trim Galore) 3. 序列比对(Bowtie2) 4. 后处理(samtools) 5. 峰值检测(MACS2) 6. 注释与富集分析 7. 可视化(IGV, UCSC) 8. 基序分析(MEME-ChIP) 9. 高级分析(可选) 每个步骤简要说明工具和目的,并添加引用。然后列出相关问题,如性能分析、工具选择、参数优化等。</think>### ChIP-seq息学分析流程 ChIP-seq(染色质免疫沉淀测序)分析流程通常包含以下步骤,结合ENCODE标准及常用方法[^1][^2][^3]: --- #### **1. 数据质量控制** 使用`FastQC`评估原始测序数据质量,检查测序错误率、GC含量、重复序列等指标。低质量数据会导致后续分析误差,需特别关注Per base sequence quality和Adapter contamination模块。 ```bash fastqc SRR14879780_1.fastq SRR14879780_2.fastq ``` --- #### **2. 数据预处理** 通过`Trim Galore`自动去除接头序列并裁剪低质量碱基(通常质量值阈值设为$Q \geq 20$),处理后的数据保存为`.trimmed.fastq`: ```bash trim_galore --paired --quality 20 SRR14879780_1.fastq SRR14879780_2.fastq ``` --- #### **3. 序列比对** 使用`Bowtie2`将reads比对到参考基因组(如hg38),成SAM/BAM文件。比对率需满足ENCODE标准(通常>70%): ```bash bowtie2 -x hg38 -1 SRR14879780_1.trimmed.fastq -2 SRR14879780_2.trimmed.fastq -S aligned.sam ``` --- #### **4. 比对后处理** 利用`samtools`过滤低质量比对结果(如MAPQ < 30)并排序: ```bash samtools view -@ 8 -bS -q 30 aligned.sam | samtools sort -@ 8 -o sorted.bam ``` --- #### **5. 峰值检测** 使用`MACS2`识别显著富集区域。其核心算法基于局部泊松分布: $$ \lambda_{local} = \max(\lambda_{BG}, \lambda_{1kb} \times \frac{L}{1kb}) $$ 其中$\lambda_{BG}$为全局背景噪声,$L$为峰宽[^1]。 ```bash macs2 callpeak -t sorted.bam -c control.bam -f BAM -g hs -n output --broad ``` --- #### **6. 注释与富集分析** 通过`HOMER`或`ChIPseeker`对峰值进行基因注释,结合`DAVID`或`GREAT`分析功能富集。 --- #### **7. 可视化** - **IGV**:查看特定基因位点的reads分布(如TP53基因区域)[^3]。 - **UCSC Genome Browser**:成可共享的基因组轨迹图。 --- #### **8. 基序分析** 使用`MEME-ChIP`预测转录因子结合基序,典型输出包括位置权重矩阵: $$ PWM_{i,j} = \log_2 \frac{p_{i,j}}{b_j} $$ 其中$p_{i,j}$为碱基$j$在位置$i$的频率,$b_j$为背景频率。 --- #### **9. 高级分析(可选)** - 差异结合分析(如`DiffBind`) - 与ATAC-seq或RNA-seq数据整合 ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

生信与基因组学

每一份鼓励是我坚持下去动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值