一個(gè)完整的WGS變異檢測(cè)snakemake流程詳解

寫在前邊的話

??在此要感謝堿基礦工，第一次做全外顯子測(cè)序的變異檢測(cè)時(shí)便是參考了他寫的教程，當(dāng)時(shí)的GATK還是3.8版本，而現(xiàn)在已經(jīng)是4.1了。使用GATK來(lái)做變異檢測(cè)非常好用，但是在有些步驟是非常耗時(shí)的，如果不能用多線程去處理的話，尤其是有多個(gè)樣本時(shí)，就會(huì)效率很差。所以才打算寫一個(gè)snakemake流程化處理，最后會(huì)附上一個(gè)完整腳本供大家參考。

第一步：建立index以及下載數(shù)據(jù)庫(kù)


??Index 使用UCSC數(shù)據(jù)庫(kù)的hg38作為reference，使用bwa進(jìn)行比對(duì)。下載hg38數(shù)據(jù)存為hg38.chrom.fasta, 建立bwa比對(duì)index：

$ bwa index -a bwtsw -p hg38.chrom.fasta hg38.chrom.fasta

參數(shù)說(shuō)明：
-p 指定輸出index文件名字
第二個(gè)hg38.chrom.fasta為當(dāng)前文件夾下的reference

??另外，還需要一個(gè)hg38.chrom.dict文件和一個(gè)hg38.chrom.fasta.fai文件, 命令如下：

$ gatk CreateSequenceDictionary -R hg38.chrom.fasta -O hg38.chrom.dict
$ samtools faidx hg38.chrom.fasta

參數(shù)說(shuō)明：
-R 基因組reference文件
-O 輸出文件路徑
生成文件和輸入文件都是在當(dāng)前目錄下

??index基本建立完成，然后是已知突變數(shù)據(jù)庫(kù)下載，gatk提供了相關(guān)網(wǎng)站bundle可下載相關(guān)數(shù)據(jù)。下載hg38的所有數(shù)據(jù)到本地：

wget -r ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/*

第二步：數(shù)據(jù)處理

1. 去接頭

??使用trim_galore進(jìn)行去接頭，自動(dòng)檢測(cè)接頭序列：

input:
    R1 = '/Volumes/RawData1/WGS/Fastq/{sample}_R1.fastq.gz',
    R2 = '/Volumes/RawData1/WGS/Fastq/{sample}_R2.fastq.gz'
output:
    'Trim/{sample}_R1_val_1.fq.gz',
    'Trim/{sample}_R2_val_2.fq.gz'
shell:
    '''
        trim_galore \
        -q 20 \         
        --length 25 \ 
        --e 0.1 \
        --paired \
        {input.R1} \
        {input.R2} \
        --gzip \
        -o Trim
    '''

2. 序列比對(duì)

??BWA 的比對(duì)速度實(shí)在是太慢了，而hisat2速度相對(duì)快一些，而且一樣可以用于dna測(cè)序比對(duì)；比對(duì)之后利用samtools進(jìn)行轉(zhuǎn)換為bam并排序：

input:
    r1 = 'Trim/{sample}_R1_val_1.fq.gz',
    r2 = 'Trim/{sample}_R2_val_2.fq.gz',
    index = index
output:
    bam = 'Mapping/{sample}.sorted.bam',
    sum = 'Mapping/{sample}_aln_sum.txt'
shell:
    '''
        hisat2 \
        -p {threads} \
        -x {input.index}/genome \
        -1 {input.r1} \
        -2 {input.r2} \
        --summary-file {output.sum} | \
        samtools view -Sb -q 30 - | \
        samtools sort -@ {threads} -m 2G -O bam \
        -T Mapping/{wildcards.sample}.tmp -o {output.bam}
    '''

3. 添加頭信息

??由于不是用bwa比對(duì)，所以沒(méi)辦法在一開(kāi)始就為reads加上頭信息，所以這一步使用gatk去添加：

input:
    'Mapping/{sample}.sorted.bam'
output:
    'Mapping/{sample}.addhead.bam'
shell:
    '''
        gatk AddOrReplaceReadGroups \
        -I {input} \
        -O {output} \
        -LB {wildcards.sample} \
        -PL illumina \ # 測(cè)序平臺(tái)不能亂寫，其他隨意
        -PU {wildcards.sample} \
        -SM {wildcards.sample} \
        -SO coordinate
    '''

4.去掉pcr重復(fù)

??這一步只是把重復(fù)reads標(biāo)記了出來(lái)，并沒(méi)有刪除，可以通過(guò)更改參數(shù)去刪除reads：

input:
    'Mapping/{sample}.addhead.bam'
output:
    bam = 'MarkDup/{sample}.markdup.bam',
    met = 'MarkDup/{sample}.metrics.txt'
shell:
    '''
        gatk MarkDuplicates \
        -I {input} \
        -O {output.bam} \
        -M {output.met}
    '''

??做完這一步需要建立一個(gè)index，方便后續(xù)調(diào)用bam文件：

input:
    'MarkDup/{sample}.markdup.bam'
output:
    'MarkDup/{sample}.markdup.bam.bai'
shell:
    'samtools index {input}'

5. 堿基質(zhì)量校正

??由于比對(duì)到SNP或INDEL上的reads附近會(huì)有很多錯(cuò)配，為了避免出現(xiàn)過(guò)多假陽(yáng)性，需要對(duì)這部分reads進(jìn)行局部重新比對(duì);這個(gè)過(guò)程用到了很多已知變異集，即已知的可靠的變異位點(diǎn)，重比對(duì)將主要圍繞這些位點(diǎn)進(jìn)行：

# SNP and INDEL datasets
omni='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/1000G_omni2.5.hg38.vcf.gz'
thg='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/1000G_phase1.snps.high_confidence.hg38.vcf.gz'
mill='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz'
db146='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/dbsnp_146.hg38.vcf.gz'
hap='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/hapmap_3.3.hg38.vcf.gz'
dbsnp='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/beta/Homo_sapiens_assembly38.dbsnp.vcf.gz'
kn_indel='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/beta/Homo_sapiens_assembly38.known_indels.vcf.gz'
gold='/Volumes/RawData1/WGS/Reference/hg38_snp_datasets/beta/Homo_sapiens_assembly38.variantEvalGoldStandard.vcf.gz'
  • 第一步是重比對(duì)的過(guò)程：
input:
    bam = 'MarkDup/{sample}.markdup.bam',
    ref = ref,
    db146 = db146,
    mill = mill,
    thg = thg,
    hap = hap,
    omni = omni,
    kn_indel = kn_indel,
    gold = gold,
    dbsnp = dbsnp
output:
    'BQSR/{sample}.markdup.recal.table'
shell:
    '''
        gatk BaseRecalibrator \
        -R {input.ref} \
        -I {input.bam} \
        --known-sites {input.db146} \
        --known-sites {input.mill} \
        --known-sites {input.thg} \
        --known-sites {input.hap} \
        --known-sites {input.omni} \
        --known-sites {input.kn_indel} \
        --known-sites {input.gold} \
        --known-sites {input.dbsnp} \
        -O {output}
    '''
  • 第二步是把重比對(duì)的結(jié)果再寫入到bam文件中：
input:
    bam = 'MarkDup/{sample}.markdup.bam',
    table = 'BQSR/{sample}.markdup.recal.table',
    ref = ref
output:
    'BQSR/{sample}.BQSR.bam'
shell:
    '''
        gatk ApplyBQSR \
        --bqsr-recal-file {input.table} \
        -R {input.ref} \
        -I {input.bam} \
        -O {output}
    '''
  • 第三步依然是為bam文件建立index：
input:
    'BQSR/{sample}.BQSR.bam'
output:
    'BQSR/{sample}.BQSR.bam.bai'
shell:
    'samtools index {input}''

6. 開(kāi)始真正的變異calling

??在對(duì)reads做了校正之后，就可以拿來(lái)做Variant calling 了，這一步只是拿到初始的變異位點(diǎn)，因?yàn)楹罄m(xù)還有對(duì)這些變異位點(diǎn)進(jìn)行過(guò)濾和校正:

input:
    bam = 'BQSR/{sample}.BQSR.bam',
    ref = ref
output:
    'HC/{sample}.HC.vcf.gz'
shell:
    '''
        gatk HaplotypeCaller \
        -R {input.ref} \
        -I {input.bam} \
        -O {output}
    '''

7. 對(duì)上一步得到的變異進(jìn)行過(guò)濾

??該過(guò)程也分為兩步，第一步是根據(jù)已知變異集的變異位點(diǎn)信息，利用自己的測(cè)序數(shù)據(jù)建立一個(gè)高斯模型，用來(lái)區(qū)分好的變異位點(diǎn)和壞的變異位點(diǎn)；好的變異位點(diǎn)即為已知變異集相同或相似的位點(diǎn)，壞的變異位點(diǎn)則相反：

  • 首先是SNP篩選

  • 第一步建立模型
input:
    vcf = 'HC/{sample}.HC.vcf.gz',
    ref = ref,
    hap = hap,
    omi = omni,
    thg = thg,
    dbs = db146,
    dbsnp = dbsnp,
    gold = gold
output:
    R = 'VQSR/{sample}.snps.plots.R',
    tr = 'VQSR/{sample}.snps.tranches',
    recal = 'VQSR/${sample}.snps.recal'
shell:
    '''
        gatk VariantRecalibrator \
        -R {input.ref} \
        -V {input.vcf} \
        --resource hapmap,known=false,training=true,truth=true,prior=15.0:{input.hap} \
        --resource omini,known=false,training=true,truth=true,prior=12.0:{input.omi} \
        --resource 1000G,known=false,training=true,truth=false,prior=10.0:{input.thg} \
        --resource dbsnp,known=true,training=false,truth=false,prior=2.0:{input.dbs} \
        --resource dbsnp38,known=true,training=false,truth=false,prior=2.0:{input.dbsnp} \
        --resource gold,known=true,training=false,truth=false,prior=2.0:{input.gold} \
        -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \ # 全基因組測(cè)序還有一個(gè)參數(shù)是 -an QD ，外顯子測(cè)序沒(méi)有
        -mode SNP \
        -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
        --rscript-file {output.R} \
        --tranches-file {output.tr} \
        -O {output.recal}
    '''
  • 第二步篩選變異位點(diǎn)：
input:
    vcf = 'HC/{sample}.HC.vcf.gz',
    ref = ref,
    tr = 'VQSR/{sample}.snps.tranches',
    recal = 'VQSR/${sample}.snps.recal'
output:
    'VQSR/{sample}.snp.vcf'
shell:
    '''
        gatk ApplyVQSR \
        -R {input.ref} \
        -V {input.vcf} \
        --ts-filter-level 99.0 \
        --tranches-file {input.tr} \
        --recal-file {input.recal} \
        -mode SNP \
        -O {output}
    '''

  • 然后是INDEL篩選

??INDEL的篩選是建立在snp篩選基礎(chǔ)上的，所以snp篩選用過(guò)的變異集在這一步就不再用了。

  • 第一步建立模型
input:
    vcf = 'VQSR/{sample}.snp.vcf',
    ref = ref,
    mill = mill,
    kn_indel = kn_indel
output:
    tr = 'VQSR/{sample}.indel.tranches',
    R = 'VQSR/{sample}.indel.plots.R',
    recal = 'VQSR/${sample}.indel.recal'
shell:
    '''
        gatk VariantRecalibrator \
        -R {input.ref} \
        -V {input.vcf} \
        --resource mills,known=true,training=true,truth=true,prior=12.0:{input.mill} \
        --resource kn_indel,known=true,training=true,truth=true,prior=10.0:{input.kn_indel} \
        -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \ # 全基因組測(cè)序還有一個(gè)參數(shù)是 -an QD ，外顯子測(cè)序沒(méi)有
        -mode INDEL \
        --max-gaussians 6 \
        --rscript-file {output.R} \
        --tranches-file {output.tr} \
        -O {output.recal}
    '''
  • 第二步篩選變異位點(diǎn)
input:
    vcf = 'VQSR/{sample}.snp.vcf',
    ref = ref,
    tr = 'VQSR/{sample}.indel.tranches',
    recal = 'VQSR/${sample}.indel.recal'
output:
    'VQSR/{sample}.snp.indel.vcf'
shell:
    '''
        gatk ApplyVQSR \
        -R {input.ref} \
        -V {input.vcf} \
        --ts-filter-level 99.0 \
        --tranches-file {input.tr} \
        --recal-file {input.recal} \
        -mode INDEL \
        -O {output}
    '''

8.拆分SNP和INDEL結(jié)果

??由于前兩步產(chǎn)生的結(jié)果都在同一個(gè)文件里，所以需要從中把SNP和INDEL分別拆分出來(lái)：

input:
    vcf = 'VQSR/{sample}.snp.indel.vcf',
    ref = ref
output:
    snp = 'SNP_INDEL/{sample}.snp.vcf.gz',
    indel = 'SNP_INDEL/{sample}.indel.vcf.gz'
shell:
    '''
        gatk SelectVariants -R {input.ref} -select-type SNP --variant {input.vcf} -O {output.snp}
        gatk SelectVariants -R {input.ref} -select-type INDEL --variant {input.vcf} -O {output.indel}
    '''

9. 選擇通過(guò)篩選的變異

??在前兩步的結(jié)果中，通過(guò)篩選的變異會(huì)加上一個(gè)'PASS'的標(biāo)簽，我們通過(guò)這個(gè)標(biāo)簽可以選擇出對(duì)應(yīng)的變異位點(diǎn)：

input:
    snp = 'SNP_INDEL/{sample}.snp.vcf.gz',
    indel = 'SNP_INDEL/{sample}.indel.vcf.gz'
output:
    snp = 'PASS/{sample}.filtered.snp.vcf',
    indel = 'PASS/{sample}.filtered.indel.vcf'
shell:
    '''
        zgrep 'PASS' {input.snp} > {output.snp}
        zgrep 'PASS' {input.indel} > {output.indel}
    '''

10. 對(duì)變異位點(diǎn)進(jìn)行注釋

??最后得到的變異位點(diǎn)我們并不能直接看出它的功能等信息，因此需要利用已知的變異庫(kù)去注釋，然后去觀察它是否和一些疾病等有關(guān)。注釋軟件我們選擇ANNOVAR,因?yàn)樗螺d即可使用，而且提供常用的變異庫(kù)，使用起來(lái)非常方便：

input:
    snp = 'PASS/{sample}.filtered.snp.vcf',
    indel = 'PASS/{sample}.filtered.indel.vcf',
    humandb = humandb,
    xref = xref
output:
    snp = 'ANNOVAR/{sample}.filtered.snp',
    indel = 'ANNOVAR/{sample}.filtered.indel'
shell:
    '''
        table_annovar.pl {input.snp} {input.humandb} -buildver hg38 -out {output.snp} \
        -remove -protocol refGene,cytoBand,exac03,dbnsfp33a,avsnp150,cosmic70,dbscsnv11 \
        -operation gx,r,f,f,f,f,f -nastring . \
        -polish -xref {input.xref}
    

        table_annovar.pl {input.indel} {input.humandb} -buildver hg38 -out {output.indel} \
        -remove -protocol refGene,cytoBand,exac03,dbnsfp33a,avsnp150,cosmic70,dbscsnv11 \
        -operation gx,r,f,f,f,f,f -nastring . \
        -polish -xref {input.xref}
    '''

最后，我們得到的是Tab分割的表格，分別儲(chǔ)存著變異位點(diǎn)的位置、序列、相關(guān)基因等信息，可以利用這些信息去進(jìn)行可視化，例如利用IGV去觀察或者用circos或cirlize去作圖。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容