生信分析进阶5 - 全外显子组变异检测和ANNOVAR注释Snakemake分析流程

最新推荐文章于 2024-10-25 20:44:17 发布
最新推荐文章于 2024-10-25 20:44:17 发布
阅读量742 收藏 1
点赞数
CC 4.0 BY-SA版权
分类专栏： 生信分析项目进阶技能合集 文章标签： 数据分析 数据挖掘 linux
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/LittleComputerRobot/article/details/139725847

基于yaml或ini配置文件，配置文件包含例如样本名称、参考基因组版本、exon capture bed文件路径、参考基因组路径和ANNOVAR注释文件等信息。

基于该流程可以实现全外显测序的fastq文件输入到得到最终变异VCF文件。

1. Snakemake分析流程基础软件安装

# conda安装
conda install -c bioconda snakemake -y
conda install -c bioconda fastqc -y
conda install -c bioconda trim-galore -y
conda install -c bioconda bwa -y
conda install -c bioconda samtools -y

# GATK 4.2.1版本
conda install -c bioconda gatk -y

ANNOVAR注释软件安装和使用参考文章：
https://www.jianshu.com/p/461f2cd47564

ANNOVAR注释全外显子有些坑和经常会报错的问题，
本人后续有空会写下对全外显子注释的ANNOVAR软件安装和使用：

2. AgilentSureSelect 外显子序列捕获流程图

外显子捕获基本原理参考下图。
AgilentSureSelect

3. GATK SNP和INDEL分析基本流程

GATK

4. Snakemake全外显子分析流程

其中regions_bed配置为外显子捕获所采用试剂盒对应的bed文件路径。

##################### parser configfile #########################
import os
configfile: "config.yaml"
# configfile: "config.ini"

sample = config["sample"]
genome = config["genome"]
#################################################################
PROJECT_DATA_DIR = config["project_data_dir"]
PROJECT_RUN_DIR = config["project_run_dir"]

os.system("mkdir -p {0}".format(PROJECT_RUN_DIR))

######################### Output #################################
ALL_FASTQC_RAW_ZIP1 = PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_1_fastqc.zip"
ALL_FASTQC_RAW_HTML1 = PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_1_fastqc.html"
ALL_FASTQC_RAW_ZIP2 = PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_2_fastqc.zip"
ALL_FASTQC_RAW_HTML2 = PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_2_fastqc.html"

ALL_CELAN_FASTQ1 = PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_1_val_1.fq.gz"
ALL_CELAN_FASTQ2 = PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_2_val_2.fq.gz"

ALL_SORTED_BAM = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam"
ALL_SORTED_BAM_BAI1 = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam.bai"
ALL_BAM_FLAGSTAT = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam.flagstat"

ALL_EXON_INTERVAL_BED = PROJECT_RUN_DIR + "/gatk/" + sample + ".Exon.Interval.bed"
ALL_STAT_TXT = PROJECT_RUN_DIR + "/gatk/" + sample + ".stat.txt"
ALL_MKDUP_BAM = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam"
ALL_METRICS = PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.bam.metrics"
ALL_SORTED_BAM_BAI2 = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam.bai"

ALL_INSERT_SIZE_HISTOGRAM_pdf = PROJECT_RUN_DIR + "/gatk/" + sample + ".insert_size_histogram.pdf"

ALL_RECAL_DATA_TABLE = PROJECT_RUN_DIR + "/gatk/" + sample + ".recal_data.table"
ALL_BQSR_SORTED_BAM = PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.mkdup.BQSR.bam"
ALL_depthOfcoverage = PROJECT_RUN_DIR + "/gatk/depthOfcoverage"

ALL_RAW_VCF = PROJECT_RUN_DIR + "/vcf/" + sample + ".raw.vcf"
All_SNP_SELECT = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.vcf"
ALL_SNP_FILTER_VCF = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.vcf"
ALL_SNP_FILTER_VCF_STAT = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.filter.vcf_stat.csv"

All_INDEL_SELECT = PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.vcf"
ALL_INDEL_FILTER_VCF = PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.vcf"

ALL_SNP_AVINPUT = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.avinput"
ALL_INDEL_AVINPUT = PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.avinput"

ALL_MERG_FILTERE_VCF = PROJECT_RUN_DIR + "/vcf/" + sample + ".merge.filter.vcf"

ALL_SNPANNO = PROJECT_RUN_DIR + "/vcf/" + sample + "_snpanno." + genome + "_multianno.csv"
ALL_INDELANNO = PROJECT_RUN_DIR + "/vcf/" + sample + "_indelanno." + genome + "_multianno.csv"
ALL_ALLANNO = PROJECT_RUN_DIR + "/vcf/" + sample + "_allanno." + genome + "_multianno.csv"

QUALIMAP_REPORT = PROJECT_RUN_DIR + "/gatk/bam_QC/" + sample + "_bamQC_report.pdf"


######################## Workflow #################################
rule all:
    input:
        ALL_FASTQC_RAW_ZIP1,
        ALL_FASTQC_RAW_HTML1,
        ALL_FASTQC_RAW_ZIP2,
        ALL_FASTQC_RAW_HTML2,
        ALL_SORTED_BAM,
        ALL_SORTED_BAM_BAI1,
        ALL_BAM_FLAGSTAT,
        ALL_EXON_INTERVAL_BED,
        ALL_STAT_TXT,
        ALL_MKDUP_BAM,
        ALL_METRICS,
        ALL_SORTED_BAM_BAI2,
        ALL_RECAL_DATA_TABLE,
        ALL_BQSR_SORTED_BAM,
        ALL_RAW_VCF,
        All_SNP_SELECT,
        All_INDEL_SELECT,
        ALL_SNP_FILTER_VCF,
        ALL_INDEL_FILTER_VCF,
        ALL_MERG_FILTERE_VCF,
        ALL_SNPANNO,
        ALL_INDELANNO,
        ALL_ALLANNO


rule fastqc_raw:
    input:
        r1 = PROJECT_DATA_DIR + "/" + sample + "_1.fq.gz",
        r2 = PROJECT_DATA_DIR + "/" + sample + "_2.fq.gz"
    output:
        PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_1_fastqc.zip",
        PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_2_fastqc.zip",
        PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_1_fastqc.html",
        PROJECT_RUN_DIR + "/fastqc/raw_data/" + sample + "_2_fastqc.html"
    threads: 8
    params:
        raw_data_dir = PROJECT_RUN_DIR + "/fastqc/raw_data",
        clean_data_dir = PROJECT_RUN_DIR + "/fastqc/clean_data",	
        align_dir = PROJECT_RUN_DIR + "/align", 
        gatk_dir = PROJECT_RUN_DIR + "/gatk",
        vcf_dir = PROJECT_RUN_DIR + "/vcf",
        tmp_dir = PROJECT_RUN_DIR + "/tmp",
	fastqc_report_dir = PROJECT_RUN_DIR + "/report/fastqc_report",
	align_report_dir = PROJECT_RUN_DIR + "/report/align_report"
    shell:
        """
        mkdir -p {params.raw_data_dir}
        mkdir -p {params.clean_data_dir}
        mkdir -p {params.align_dir}
        mkdir -p {params.gatk_dir}
        mkdir -p {params.vcf_dir}
        mkdir -p {params.tmp_dir}
        fastqc -f fastq --extract -t {threads} -o {params.raw_data_dir} {input.r1} {input.r2}        
	"""

rule fastqc_clean:
    input:
        r1 = PROJECT_DATA_DIR + "/" + sample + "_1.fq.gz",
        r2 = PROJECT_DATA_DIR + "/" + sample + "_2.fq.gz"
    output:
        temp(PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_1_val_1.fq.gz"),
        temp(PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_2_val_2.fq.gz")
    threads: 8
    params:
        quality = 20,
        length = 20,
	    clean_data_dir = PROJECT_RUN_DIR + "/fastqc/clean_data"
    shell:
        """
        trim_galore --paired --quality {params.quality} --length {params.length} -o {params.clean_data_dir} {input.r1} {input.r2}
	    """

# BWA mem比对
rule bwa_mem:
    input:
        reads1 = PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_1_val_1.fq.gz",
        reads2 = PROJECT_RUN_DIR + "/fastqc/clean_data/" + sample + "_2_val_2.fq.gz"
    output:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam"
    threads:8
    params:
        hg = config["hg_fa"]
    shell:
        """
        bwa mem -t {threads} -M -R "@RG\\tID:{sample}\\tPL:ILLUMINA\\tLB:{sample}\\tSM:{sample}" {params.hg} {input.reads1} {input.reads2}|samtools sort -@ {threads} -o {output}
        """

rule bam_index_1:
    input:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam"
    output:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam.bai"
    threads: 2
    shell:
        "samtools index {input} > {output}"

rule flagstat:
    input:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam"
    output:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam.flagstat"
    threads: 1
    shell:
        "samtools flagstat {input} > {output}"

# 创建reference fasta 字典，GATK需要用到
# rule CreateSequenceDictionary:
#     input:
#         "/reference/hg19/ucsc.hg19.fa"
#     output:
#         PROJECT_RUN_DIR + "/" + sample + "/gatk/" + sample + ".hg19.dict"
#         # /GATK_bundle/hg19/ucsc.hg19.dict
#     threads: 1
#     params:
#         gatk4 = "/public/software/gatk-4.0.6.0/gatk",
#         gatk_dir = PROJECT_RUN_DIR + "/" + sample + "/gatk"
#     shell:
#         """
#         mkdir -p {params.gatk_dir}
#         gatk CreateSequenceDictionary -R {input} -O {ouptut}
#         """ 

rule CreateExonIntervalBed:
    input:
        regions_bed = config["regions_bed"],
        hg_dict = config["hg_dict"]
    output:
        Exon_Interval_bed = PROJECT_RUN_DIR + "/gatk/" + sample + ".Exon.Interval.bed"
    threads: 1
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk"
    shell:
        """
        gatk BedToIntervalList -I {input.regions_bed} -O {output.Exon_Interval_bed} -SD {input.hg_dict}
        """ 

# 外显子区域覆盖度
rule exon_region:
    input:
        bam = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam",
        Exon_Interval_bed = PROJECT_RUN_DIR + "/gatk/" + sample + ".Exon.Interval.bed"
    output:
        stat_txt = PROJECT_RUN_DIR + "/gatk/" + sample + ".stat.txt"
    threads: 2
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk"
    shell:
        """
        gatk CollectHsMetrics -BI {input.Exon_Interval_bed} -TI {input.Exon_Interval_bed} -I {input.bam} -O {output.stat_txt}
        """ 

# 标记PCR重复序列
rule mark_duplicates:
    input:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.bam"
    output:
        mkdup_bam = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam",
        metrics = PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.bam.metrics"
    threads: 4
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk"
    shell:
        """
        gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" MarkDuplicates -I {input} -O {output.mkdup_bam} -M {output.metrics}
        """

rule bam_index_2:
    input:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam"
    output:
        PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam.bai"
    threads: 2
    shell:
        "samtools index {input} > {output}"

# 插入片段分布统计
#rule InsertSizeStat:
#    input:
#        PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.mkdup.bam"
#    output:
#        insert_size_metrics_txt = temp(PROJECT_RUN_DIR + "/gatk/" + sample + ".insert_size_metrics.txt"),
#        insert_size_histogram_pdf = PROJECT_RUN_DIR + "/gatk/" + sample + ".insert_size_histogram.pdf"
#    threads:1
#    params:
#        picard = "/public/software/picard.jar",
#        M = 0.5
#    shell:
#        """
#       java -jar {params.picard} CollectInsertSizeMetrics \
#        I={input} \
#        O={output.insert_size_metrics_txt} \
#        H={output.insert_size_histogram_pdf} \
#        M=0.5
#        """

# 重新校正碱基质量值
# 计算出所有需要进行重校正的read和特征值，输出为一份校准表文件.recal_data.table
rule BaseRecalibrator:
    input:
        mkdup_bam = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam",
        mkdup_bam_index = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam.bai",
        vcf_1000G_phase1_indels_sites = config["vcf_1000G_phase1_indels_sites"],
        vcf_Mills_and_1000G_gold_standard_indels_sites = config["vcf_Mills_and_1000G_gold_standard_indels_sites"],
        vcf_dbsnp = config["vcf_dbsnp"]
    output:
        PROJECT_RUN_DIR + "/gatk/" + sample + ".recal_data.table"
    threads: 8
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk",
        hg_fa = config["hg_fa"],
        regions_bed = config["regions_bed"]
    shell:
        """
        gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" BaseRecalibrator \
        -R {params.hg_fa} \
        -I {input.mkdup_bam} \
        --known-sites {input.vcf_1000G_phase1_indels_sites} \
        --known-sites {input.vcf_Mills_and_1000G_gold_standard_indels_sites} \
        --known-sites {input.vcf_dbsnp} \
        -L {params.regions_bed} \
        -O {output}
        """
        
# gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" BaseRecalibrator -R genome/hg19.fa -I ./TKQX.sorted.mkdup.bam --known-sites /GATK_bundle/hg19/1000G_phase1.indels.hg19.sites.vcf --known-sites /GATK_bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf --known-sites /GATK_bundle/hg19/dbsnp_138.hg19.vcf -L /GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed -O ./TKQX.recal_data.table
# java -jar CreateSequenceDictionary.jar R= Homo_sapiens_assembly18.fasta O= Homo_sapiens_assembly18.dict

# 利用校准表文件重新调整原来BAM文件中的碱基质量值,使用这个新的质量值重新输出新的BAM文件
rule ApplyBQSR:
    input:
        mkdup_bam = PROJECT_RUN_DIR + "/align/" + sample + ".sorted.mkdup.bam",
        recal_data = PROJECT_RUN_DIR + "/gatk/" + sample + ".recal_data.table",
        regions_bed = config["regions_bed"]
        #regions_bed = "/GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed"
    output:
        PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.mkdup.BQSR.bam"
    threads: 4
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk",
        hg_fa = config["hg_fa"]
        #hg_fa = "genome/hg19.fa"
    shell:
        """
        gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" ApplyBQSR \
        -R {params.hg_fa} \
        -I {input.mkdup_bam} \
        -bqsr {input.recal_data} \
        -L {input.regions_bed} \
        -O {output}
        """
# gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" ApplyBQSR -R genome/hg19.fa -I ./test_sample.sorted.mkdup.bam -bqsr ./test_sample.recal_data.table -L /GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed -O ./test_sample.sorted.mkdup.BQSR.bam

# # 校正碱基图
# rule recal_data_plot:
#     input:
#         PROJECT_RUN_DIR + "/" + sample + "/gatk/" + sample + ".recal_data.table"
#     output:
#         PROJECT_RUN_DIR + "/" + sample + "/gatk/" + sample + ".recal_data.table.plot"
#     threads: 1
#     params:
#         gatk4 = "/public/software/gatk-4.0.6.0/gatk",
#         gatk = "/public/software/anaconda3/envs/PGS_1M/bin/gatk"
#     shell:
#         """
#         gatk AnalyzeCovariates -bqsr {input} -plots {output}
#         """
# gatk AnalyzeCovariates -bqsr ./test_sample.recal_data.table -plots ./test_sample.recal_data.table.plot

# 统计测序深度和覆盖度
#rule depthOfcoverage:
#    input:
#        BQSR_mkdup_bam = PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.mkdup.BQSR.bam"
#    output:
#        PROJECT_RUN_DIR + "/gatk/depthOfcoverage"
#    threads:1
#    params:
#        #hg19_fa = "genome/hg19.fa",
#        hg_fa = config["hg_fa"],
#        bamdst = "/public/software/bamdst/bamdst",
#        regions_bed = config["regions_bed"]
#        #regions_bed = "/GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed"
#    shell:
#        """
#        {params.bamdst} -p {params.regions_bed} -o {output} {input.BQSR_mkdup_bam}
#        """
# """
# gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" DepthOfCoverage  \
# -R {params.hg19} \
# -o {output} \
# -I {input.BQSR_mkdup_bam} \
# -L {input.regions_bed} \
# --omitDepthOutputAtEachBase --omitIntervalStatistics \
# -ct 1 -ct 5 -ct 10 -ct 20 -ct 30 -ct 50 -ct 100
# """
#gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" DepthOfCoverage -R genome/hg19.fa -o ./TKQX -I ./TKQX.sorted.mkdup.BQSR.bam -L /GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed --omitDepthOutputAtEachBase --omitIntervalStatistics -ct 1 -ct 5 -ct 10 -ct 20 -ct 30 -ct 50 -ct 100

# /public/software/bamdst/bamdst -p /GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed -o ./depth_coverage ./TKQX.sorted.mkdup.BQSR.bam

# # 单样本VCF calling
rule vcf_calling:
    input:
        BQSR_mkdup_bam = PROJECT_RUN_DIR + "/gatk/" + sample + ".sorted.mkdup.BQSR.bam"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".raw.vcf"
    threads: 4
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk",
        hg_fa = config["hg_fa"],
        vcf_dbsnp = config["vcf_dbsnp"],
        regions_bed = config["regions_bed"]
    shell:
        """
        gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" HaplotypeCaller \
        -R {params.hg_fa} \
        -I {input.BQSR_mkdup_bam} \
        -D {params.vcf_dbsnp} \
        -L {params.regions_bed} \
        -O {output}
        """
# gatk --java-options "-Xmx30G -Djava.io.tmpdir=./" HaplotypeCaller -R genome/hg19.fa -I ./TKQX.sorted.mkdup.BQSR.bam -D /GATK_bundle/hg19/dbsnp_138.hg19.vcf -L /GATK_bundle/bed/SureSelectAllExonV6/S07604514_Regions.bed -O ./TKQX.raw.vcf

# 变异质控和过滤，对raw data进行质控，剔除假阳性的标记
rule vcf_select_SNP:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".raw.vcf"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.vcf"
    threads: 2
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk"
    shell:
        """
        gatk SelectVariants -select-type SNP -V {input} -O {output}
        """

rule vcf_filter_SNP:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.vcf"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.filter.vcf"
    threads: 2
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk",
        QD = 2.0,
        MQ = 40.0,
        FS = 60.0,
        SOR = 3.0,
        MQRankSum = -12.5,
        ReadPosRankSum = -8.0
    shell:
        """
        gatk VariantFiltration \
        -V {input} \
        --filter-expression "QD < {params.QD} || MQ < {params.MQ} || FS > {params.FS} || SOR > {params.SOR} || MQRankSum < {params.MQRankSum} || ReadPosRankSum < {params.ReadPosRankSum}" --filter-name "PASS" \
        -O {output}
        """
# --filter-expression "QD < 2.0|| MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS"‘
# gatk VariantFiltration -V TKQX.snp.vcf --filter-expression "QD < 2.0|| MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O TKQX.snp.filter.vcf

# rule snp_frequency_stat:
#     input:
#         PROJECT_RUN_DIR + "/" + sample + "/vcf/" + sample + ".snp.filter.vcf"
#     output:
#         PROJECT_RUN_DIR + "/" + sample + "/vcf/" + sample + ".snp.filter.vcf_stat.csv"
#     threads: 1
#     params:
#         gatk4 = "/public/software/gatk-4.0.6.0/gatk",
#         snp_frequency = "/public/analysis/pipeline/WES/scripts/snp_frequency.py"
#     shell:
#         """
#         python {params.snp_frequency} {input}
#         """

rule vcf_select_InDel:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".raw.vcf"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.vcf"
    threads: 2
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk"
    shell:
        """
        gatk SelectVariants -select-type INDEL -V {input} -O {output}
        """

rule vcf_filter_InDel:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.vcf"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.filter.vcf"
    threads: 2
    params:
        gatk4 = "/public/software/gatk-4.0.6.0/gatk",
        QD = 2.0,
        FS = 200.0,
        SOR = 10.0,
        MQRankSum = -12.5,
        ReadPosRankSum = -8.0
    shell:
        """
        gatk VariantFiltration \
        -V {input} \
        --filter-expression "QD < {params.QD} || FS > {params.FS} || SOR > {params.SOR} || MQRankSum < {params.MQRankSum} || ReadPosRankSum < {params.ReadPosRankSum}" --filter-name "PASS" \
        -O {output}
        """
# --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0|| MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS"

rule merge_vcf:
    input:
        snp_filter_vcf = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.filter.vcf",
        indel_filter_vcf = PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.filter.vcf"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".merge.filter.vcf"
    threads: 1
    params:
        vcf_dir = PROJECT_RUN_DIR + "/vcf"
    shell:
        """
        cd {params.vcf_dir}
        gatk MergeVcfs \
        -I {input.snp_filter_vcf} \
        -I {input.indel_filter_vcf} \
        -O {output}
        """

rule transfer_avinput:
    input:
        snp = PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.filter.vcf",
        indel= PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.filter.vcf",
        merge = PROJECT_RUN_DIR + "/vcf/" + sample + ".merge.filter.vcf"
    output:
        snp = temp(PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.avinput"),
        indel = temp(PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.avinput"),
        merge = temp(PROJECT_RUN_DIR + "/vcf/" + sample + ".merge.avinput")
    threads: 4
    params:
        convert2annovar = "/public/software/annovar/convert2annovar.pl"
    shell:
        """
        perl {params.convert2annovar} -format vcf4 {input.snp} > {output.snp}
        perl {params.convert2annovar} -format vcf4 {input.indel} > {output.indel}
        perl {params.convert2annovar} -format vcf4 {input.merge} > {output.merge}
        """
# SNP注释
rule snp_annotate:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".snp.avinput"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + "_snpanno." + genome + "_multianno.csv"
    threads: 2
    params:
        table_annovar = "/public/software/annovar/table_annovar.pl",
        humandb = "/public/software/annovar/humandb/",
        prefix = sample + "_snpanno",
        genome = config["genome"],
        vcf_dir = PROJECT_RUN_DIR + "/vcf"
    shell:
        """
        cd {params.vcf_dir}
        perl {params.table_annovar} {input} {params.humandb} -buildver {params.genome} -out {params.prefix} \
        -remove -protocol refGene,cytoBand,clinvar_20221231,avsnp147,exac03,gnomad211_exome,SAS.sites.2015_08,esp6500siv2_all,intervar_20180118,dbnsfp42c \
        -operation g,r,f,f,f,f,f,f,f,f -nastring . -csvout
        """
        
# INDEL注释
rule indel_annotate:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".indel.avinput"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + "_indelanno." + genome + "_multianno.csv"
    threads: 2
    params:
        table_annovar = "/public/software/annovar/table_annovar.pl",
        humandb = "/public/software/annovar/humandb/",
        prefix = sample + "_indelanno",
        genome = config["genome"],
        vcf_dir = PROJECT_RUN_DIR + "/vcf",
        perl = "/home/miniconda3/pkgs/perl-5.26.2-h36c2ea0_1008/bin/perl"
    shell:
        """
        cd {params.vcf_dir}
        perl {params.table_annovar} {input} {params.humandb} -buildver {params.genome} -out {params.prefix} \
        -remove -protocol refGene,cytoBand,clinvar_20221231,gnomad211_exome,exac03,esp6500siv2_all,ALL.sites.2015_08 \
        -operation g,r,f,f,f,f,f -nastring . -csvout
        """

#
rule merge_annotate:
    input:
        PROJECT_RUN_DIR + "/vcf/" + sample + ".merge.avinput"
    output:
        PROJECT_RUN_DIR + "/vcf/" + sample + "_allanno." + genome + "_multianno.csv"
    threads: 2
    params:
        table_annovar = "/public/software/annovar/table_annovar.pl",
        humandb = "/public/software/annovar/humandb/",
        prefix = sample + "_allanno",
        genome = config["genome"],
        vcf_dir = PROJECT_RUN_DIR + "/vcf"
    shell:
        """
        cd {params.vcf_dir}
        perl {params.table_annovar} {input} {params.humandb} -buildver {params.genome} -out {params.prefix} \
        -remove -protocol refGene,cytoBand,clinvar_20221231,avsnp147,exac03,gnomad211_exome,SAS.sites.2015_08,esp6500siv2_all,intervar_20180118,dbnsfp42c \
        -operation g,r,f,f,f,f,f,f,f,f -nastring . -csvout
        """

生信分析进阶文章推荐

生信分析进阶1 - HLA分析的HLA区域reads提取及bam转换fastq

生信分析进阶2 - 利用GC含量的Loess回归矫正reads数量

生信分析进阶3 - pysam操作bam文件统计unique reads和mapped reads高级技巧合辑

生信分析进阶4 - 比对结果的FLAG和CIGAR信息含义与BAM文件指定区域提取

全外显子测序分析流程5 - ANNOVAR注释VCF变异位点
LittleComputerRobot的博客
11-04 514
ANNOVAR注释VCF变异位点
生信分析进阶6 - 基于karyoploteR包进行SNP芯片CNV分析及可视化
LittleComputerRobot的博客
06-19 693
基于karyoploteR包进行SNP芯片CNV分析及可视化
0084-【生信软件】-ANNOVAR-使用
leadingsci的博客
10-07 6536
文章目录ANNOVAR的程序模块ANNOVAR的输入文件ANNOVAR输入文件的格式转换ANNOVAR注释功能下载数据库Annotate_variation.plAnnotate_variation.pl 实例Region-based annotationFilter-based annotation下面给大家介绍常用的两种过滤注释dbSNP annotations ANNOVAR的程序模块 (A...
基因组变异注释 — ANNOVAR（一）
long_1998的博客
12-05 4758
识别特定基因组区域的变异，例如 44 个物种的保守区域，预测转录因子结合位点，片段重复区域，GWAS hits，基因组变异数据库，DNAse I超敏位点，ENCODE H3K4Me1/H3K4Me3/H3K27Ac/CTCF 位点，ChIP-Seq 峰，RNA-Seq 峰，或许多其他基因组区间的注释。可以看到，总共的48万位点，其中有13万是在千人基因组计划出现的，有17万是在EXAC数据库出现，但是只有区区7527个位点是在COSMIC数据库出现，，在clinvar数据库的，有15575位点。
全外显子数据分析流程
bio_meimei的博客
08-28 9741
文章目录WES全外显子数据分析流程WES相关软件下载安装下载GATK4.X版本配置环境变量下载GATK需要的参考文件下载好后用bwa构建索引比对查看bam文件一步法找变异流程先查看bam文件最简单的找变异流程GATK完整找变异流程定义变量标记PCR重复reads WES全外显子数据分析流程 WES相关软件下载安装 #conda安装前要先search一下 conda search snpeff conda install sra-tools conda install samtools conda insta
资料总结分享：《全外显子测序数据的流程和原理》
weixin_69558614的博客
04-17 4148
所以这里大家需要记住一个重点，PCR扩增原本的目的是为了增大微弱DNA序列片段的密度，但由于整个反应都在一个试管中进行，因此其他一些密度并不低的DNA片段也会被同步放大，那么这时在取样去上机测序的时候，这些DNA片段就很可能会被重复取到相同的几条去进行测序，但是由同一个模板分子扩增出来的重复子文库只对应单一模板，在分析过程中应将重复片段予以去除。因此，我们需要先把这一大堆的短序列捋顺，一个个去跟该物种的参考基因组比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。
一文理解全外显子家系和病例-对照参考基因组比对和变异检测全流程
weixin_42192188的博客
08-14 2487
这篇文章是我的生信社群的部分内容，现在开个遗传学专栏，分享一些遗传学数据分析教程，更多教程可以加入我的生信社群查看。看完这篇文章，完全理解参考基因组比对和变异检测全流程！
全外显子测序分析流程4 - GATK碱基校正 + Calling变异
LittleComputerRobot的博客
09-10 422
全外数据分析：GATK碱基校正 + Calling变异
生信技能59 - 基于GATK CallingSNP变异检测及注释流程
LittleComputerRobot的博客
09-23 912
基于GATK CallingSNP变异检测及注释流程
对annovar注释的总结
Cassiel60的博客
01-20 7400
新的一年到了，到了该更新数据库的时候了，首先就从常用的annovar软件以及所使用的数据库开始更新吧。 1.下载最新版的annovar软件，annovar.latest版本，需要一个edu的邮箱。 2.下载注释所用的数据库，官网上太多，可以根据需要自行下载，比如下载hg19版本的, clinvar数据库更新到20190305版本了，但是NCBI数据库已经更新到20191223了，可以自己下载后...
ANNOVAR region-based annotation-下篇
庐州月光的博客
06-12 705
欢迎关注"生信修炼手册"！在UCSC网站上，对于不同的参考基因组，提供了许多数据，这些数据中只要是提供了基因组区域的，理论上都可以用于annovar region-bas...
9月14日-全外显子组测序分析流程
qq_39777567的博客
09-14 9576
全外显子组测序分析及应用 Whole Exome Sequencing，WES 全外显子组测序分析及应用 一 背景 二 技术应用 外显子捕获试剂盒有哪些？ 外显子捕获效率是什么？ 外显子测序一般建议做多少倍的覆盖？ 外显子测序能够测出多大的片段缺失？ 外显子捕获可以做CNV吗？ 外显子测序的优点是什么？ 外显子测序样品准备： 一 背景 外显子(exon)是真核...
ANNOVAR gene-based annotation
庐州月光的博客
06-05 1240
欢迎关注"生信修炼手册"！通过基因相关注释，可以知道变异位点在基因组上的位置和对蛋白质编码的影响。在进行注释之前，首先需要下载物种对应的数据库，以human为例，命令如下...
外显子bed文件获取
baidu_38679951的博客
04-18 2126
去年打通了WES流程，但是中间一个环节一直不放心，就是bed文件的获取。 开始是根据CCDS导出的。今天发现安捷伦官网上有靶向的区域。附上网站。 Agilent SureDesign 需要注册，然后找到数据位置，因为我的数据是SureSelect V6版本，所以选了划线的bed。 选取“region”这个文件。 ...
生信软件40 - bedtools经典使用方法合集
LittleComputerRobot的博客
10-25 2805
bedtools经典使用方法合集
变异分析
lifezzu的博客
07-17 7149
变异检测（BWA+SAMtools+picard+GATK） 原文： 一、 简介 首先说说GATK可以做什么。它主要用于从sequencing 数据中进行variant calling，包括SNP、INDEL。比如现在风行的exome sequencing找variant，一般通过BWA+GATK的pipeline进行数据分析。 要run GATK，首先得了解它的网站(http
Annovar各程序的功能，自建avdb，构建索引
yuyuena的专栏
09-20 3184
ANNOVAR的程序模块（本人目录） ├── annotate_variation.pl //annovar主程序，功能包括下载数据库，三种不同的注释 ├── annovar_index.pl //index构建的程序，对于染色体是第一列的文件有效，修改网上程序 ├── coding_change.pl //用来推断蛋白质序列的程序 ├── convert2annovar.pl /...
全外显子组生信分析流程-4-数据质控
kevin_zhang2018的博客
03-21 3714
数据质控 测序数据的产生经过了DNA提取、建库、测序等多个步骤，这些步骤会产生低质量或者无效的数据，需要对下机的原始数据进行质控。 1.原始序列数据解读 高通量测序得到的原始图像数据经过碱基识别（base calling）分析转化为原始测序序列（reads），我们称之为raw data，结果以fastq文件格式存储，该文件包含序列信息和序列的质量信息。一条read由4行描述： 第1行：以@开头，...
转： Annovar 软件注释流程介绍
dibiao4034的博客
12-08 1236
第一步：下载Annovar 上Annovar官网下载（http://annovar.openbioinformatics.org/en/latest/user-guide/download/），现在要邮件注册后才能下载。邮件注册后会给你最新版软件下载地址， 下载后文件为annovar.latest.tar.gz。 第二步：安装Annovar linux系统下用该命令解压 tar...
请给我一个详细的基于ubuntu系统，完整的全外显子测序的生信分析流程。包括整个分析流程的每个步骤，所用的软件及代码参数。
最新发布
09-09
<think> 我们正在构建一个完整的全外显子测序（WES）生物信息学分析流程，基于Ubuntu系统。流程通常包括以下步骤： 1. 数据质控（Quality Control, QC） 2. 序列比对（Alignment） 3. 标记重复序列（Marking Duplicates） 4. 碱基质量重校正（Base Quality Score Recalibration, BQSR） 5. 变异检测（Variant Calling） 6. 变异质控与过滤（Variant Filtering） 7. 变异注释（Variant Annotation） 注意：我们将使用GRCh38参考基因组，并假设数据为Illumina平台产生的双端测序数据（paired-end）。 下面逐步说明每个步骤，包括使用的软件和典型命令参数。由于实际分析中可能需要根据具体情况调整参数，这里给出的是常用参数。 ### 1. 数据质控（QC） 使用FastQC进行原始测序数据的质量评估，然后使用Trimmomatic进行数据过滤。 **软件安装：** ```bash sudo apt-get update sudo apt-get install fastqc sudo apt-get install trimmomatic ``` **FastQC运行：** ```bash fastqc -o ./qc_results sample_R1.fastq.gz sample_R2.fastq.gz ``` 生成的质量报告需要人工检查，然后决定过滤参数。 **Trimmomatic运行：** ```bash java -jar /usr/share/trimmomatic/trimmomatic.jar PE -phred33 \ sample_R1.fastq.gz sample_R2.fastq.gz \ sample_R1_trimmed_paired.fastq.gz sample_R1_trimmed_unpaired.fastq.gz \ sample_R2_trimmed_paired.fastq.gz sample_R2_trimmed_unpaired.fastq.gz \ ILLUMINACLIP:/usr/share/trimmomatic/adapters/TruSeq3-PE-2.fa:2:30:10 \ LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 ``` 参数说明： - PE: 双端数据 - ILLUMINACLIP: 去除接头，后面接接头文件、允许的最大错配数等 - LEADING: 切除首端碱基质量低于3的碱基 - TRAILING: 切除末端碱基质量低于3的碱基 - SLIDINGWINDOW: 滑动窗口大小4，平均质量低于15则切除 - MINLEN: 最小长度36 ### 2. 序列比对 使用BWA-MEM将过滤后的reads比对到参考基因组。 **安装BWA：** ```bash sudo apt-get install bwa ``` **建立参考基因组索引（如果尚未建立）：** ```bash bwa index /path/to/GRCh38.fa ``` **比对：** ```bash bwa mem -t 8 -R '@RG\tID:sample\tSM:sample\tPL:ILLUMINA' \ /path/to/GRCh38.fa \ sample_R1_trimmed_paired.fastq.gz sample_R2_trimmed_paired.fastq.gz \ > sample.sam ``` 参数说明： - -t: 线程数 - -R: 设置reads组信息，包括样本ID（SM）等，用于后续分析 将SAM转为BAM并排序（使用samtools）： ```bash samtools view -@ 8 -Sb sample.sam | samtools sort -@ 8 -o sample_sorted.bam ``` ### 3. 标记重复序列 使用GATK的MarkDuplicates来标记PCR重复。 **安装GATK：** 从官网下载GATK（需要注册）并安装。假设GATK安装在`/path/to/gatk`。 **运行：** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar MarkDuplicates \ -I sample_sorted.bam \ -O sample_marked_duplicates.bam \ -M sample_metrics.txt ``` ### 4. 碱基质量重校正（BQSR） 使用GATK的BaseRecalibrator和ApplyBQSR。 **第一步：生成重校准表** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar BaseRecalibrator \ -R /path/to/GRCh38.fa \ -I sample_marked_duplicates.bam \ --known-sites /path/to/dbsnp.vcf \ # 已知的SNP位点，用于训练模型 -O sample_recal_data.table ``` **第二步：应用重校准** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar ApplyBQSR \ -R /path/to/GRCh38.fa \ -I sample_marked_duplicates.bam \ --bqsr-recal-file sample_recal_data.table \ -O sample_recalibrated.bam ``` ### 5. 变异检测 使用GATK HaplotypeCaller进行变异检测。 **运行HaplotypeCaller：** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar HaplotypeCaller \ -R /path/to/GRCh38.fa \ -I sample_recalibrated.bam \ -O sample_raw_variants.vcf ``` ### 6. 变异质控与过滤 使用GATK的VariantFiltration对原始变异进行过滤。 **对于SNP：** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar VariantFiltration \ -R /path/to/GRCh38.fa \ -V sample_raw_variants.vcf \ -O sample_filtered_snps.vcf \ --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \ --filter-name "my_snp_filter" ``` **对于INDEL：** ```bash java -jar /path/to/gatk/gatk-package-4.x.x.x-local.jar VariantFiltration \ -R /path/to/GRCh38.fa \ -V sample_raw_indels.vcf \ # 注意：需要先通过SelectVariants提取INDEL -O sample_filtered_indels.vcf \ --filter-expression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0" \ --filter-name "my_indel_filter" ``` 注意：在实际流程中，可能需要先使用SelectVariants将原始VCF分成SNP和INDEL两部分，然后分别过滤，再合并。 ### 7. 变异注释 使用ANNOVAR或VEP（Variant Effect Predictor）进行注释。 **安装VEP：** 可以按照VEP官方文档安装（使用cpan或conda）。 **运行VEP：** ```bash vep -i sample_filtered.vcf \ # 假设已经合并了SNP和INDEL --species homo_sapiens \ --assembly GRCh38 \ --offline \ # 离线模式，需提前下载缓存数据库 --cache \ --dir_cache /path/to/vep_cache \ --everything \ --output_file sample_annotated.vcf ``` 以上是一个全外显子测序分析的标准流程概览。实际项目中，可能还需要进行批次效应校正、样本间比较（如肿瘤-正常配对分析）等。 注意：每一步骤的输出都应检查质量，确保上一步成功再进行下一步。 §§相关问题§§ 1. 全外显子测序分析中，如何评估数据质控的结果？ 2. 除了BWA，还有哪些常用的比对工具？它们各有什么优缺点？ 3. 变异注释中，VEP和ANNOVAR有何异同？如何选择？ 4. 在肿瘤样本的全外显子测序分析中，有哪些特别的步骤？ </think> 以下是基于Ubuntu系统的完整全外显子测序（WES）生信分析流程，包含各步骤、软件及代码参数。流程基于行业标准实践设计[^1]，适用于人类基因组分析（GRCh38参考基因组）： --- ### **1. 原始数据质控（FastQC + MultiQC）** **软件**：FastQC + MultiQC **作用**：评估原始FASTQ文件质量 ```bash # 单个样本质控 fastqc -t 8 -o qc_report/ SRR28391647_1.fastq.gz SRR28391647_2.fastq.gz # 整合报告 multiqc qc_report/ -o multiqc_report/ ``` **关键参数**： - `-t 8`：使用8线程 - 需检查：Per base sequence quality > Q30，Adapter contamination --- ### **2. 数据过滤（Trimmomatic）** **软件**：Trimmomatic **作用**：去除低质量碱基和接头 ```bash java -jar trimmomatic PE -phred33 \ SRR28391647_1.fastq.gz SRR28391647_2.fastq.gz \ SRR28391647_1_trimmed_paired.fq SRR28391647_1_trimmed_unpaired.fq \ SRR28391647_2_trimmed_paired.fq SRR28391647_2_trimmed_unpaired.fq \ ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 \ LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:36 ``` **参数说明**： - `ILLUMINACLIP`：切除Illumina接头 - `SLIDINGWINDOW:4:20`：窗口大小4bp，平均质量≥Q20 - `MINLEN:36`：保留长度≥36bp的reads --- ### **3. 序列比对（BWA-MEM）** **软件**：BWA + SAMtools **作用**：将reads比对到参考基因组 ```bash # 生成索引（只需一次） bwa index GRCh38.fa # 比对 bwa mem -t 16 -R '@RG\tID:16177_CCPM_1300019\tSM:BB5\tLB:Illumina\tPL:ILLUMINA' \ GRCh38.fa \ SRR28391647_1_trimmed_paired.fq SRR28391647_2_trimmed_paired.fq \ | samtools sort -@ 8 -o SRR28391647.sorted.bam # 建立索引 samtools index SRR28391647.sorted.bam ``` **关键参数**： - `-R`：设置Read Group信息（样本ID/SM必填） - `-t 16`：16线程加速 --- ### **4. 标记重复序列（GATK MarkDuplicates）** **软件**：GATK4 ```bash gatk MarkDuplicates \ -I SRR28391647.sorted.bam \ -O SRR28391647.marked_dups.bam \ -M metrics.txt \ --CREATE_INDEX true ``` --- ### **5. 碱基质量校正（GATK BQSR）** **软件**：GATK4 **作用**：校正系统偏差 ```bash # 生成重校准表 gatk BaseRecalibrator \ -R GRCh38.fa \ -I SRR28391647.marked_dups.bam \ --known-sites dbsnp_146.hg38.vcf.gz \ -O recal_data.table # 应用校正 gatk ApplyBQSR \ -R GRCh38.fa \ -I SRR28391647.marked_dups.bam \ --bqsr-recal-file recal_data.table \ -O SRR28391647.recal.bam ``` --- ### **6. 变异检测（GATK HaplotypeCaller）** **软件**：GATK4 ```bash gatk HaplotypeCaller \ -R GRCh38.fa \ -I SRR28391647.recal.bam \ -O sample_raw_variants.vcf.gz \ -ERC GVCF \ --dbsnp dbsnp_146.hg38.vcf.gz ``` **参数说明**： - `-ERC GVCF`：输出gVCF便于多样本联合分析 - 肿瘤样本需添加`--af-of-alleles-not-in-resource`参数 --- ### **7. 变异质控与过滤（GATK VariantFiltration）** ```bash # SNP过滤 gatk VariantFiltration \ -R GRCh38.fa \ -V sample_raw_snps.vcf.gz \ -O filtered_snps.vcf \ --filter-name "lowQD" --filter-expression "QD < 2.0" \ --filter-name "highFS" --filter-expression "FS > 60.0" # Indel过滤 gatk VariantFiltration \ -V sample_raw_indels.vcf.gz \ --filter-name "lowQD" --filter-expression "QD < 2.0" \ --filter-name "highFS" --filter-expression "FS > 200.0" ``` --- ### **8. 变异注释（VEP）** **软件**：Ensembl VEP ```bash vep -i filtered_variants.vcf \ --species homo_sapiens \ --assembly GRCh38 \ --cache --dir_cache ${VEP_CACHE} \ --offline \ --fork 8 \ --everything \ --output_file annotated_variants.vcf ``` **关键参数**： - `--everything`：启用所有注释模块 - `--fork 8`：并行加速 --- ### **9. 结果可视化（IGV）** **软件**：Integrative Genomics Viewer 手动载入BAM和VCF文件检查变异位点可靠性 --- ### **环境配置建议** ```bash # 通过Bioconda安装所有软件 conda create -n wes_analysis \ fastqc multiqc trimmomatic bwa samtools gatk4 vep ``` ---
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
表情包 代码片
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

打赏作者

生信与基因组学

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

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

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

打赏作者

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

抵扣说明：

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

余额充值