snakemakeでバリアントコールパイプライン構築

今回は snakemake でバリアントコールパイプラインを構築してみたいと思います。

bcftoolsによるバリアントコール

まずはbamファイルからバリアントコールを行うルールを書いていきます。

rule bcftools_mpileup_call:
    input:
        rg="read_groups.tsv",
        bam=expand("bam/{sample}_sorted.bam", sample = samples['sample'])
    output:
        "vcf/calls.bcf"
    params:
        ref=config['ref']
    shell:
        "bcftools mpileup -G {input.rg} -f {params.ref} {input.bam} | bcftools call -mv -Ob -o {output}"

inputは前回紹介した bcftools mpileup の -G オプションに指定するファイル (read_groups.tsv) とbamファイルです。リファレンスのfastaはconfigfileから設定できるようにしています。

read_groups.tsvを生成

バリアントコールに必要なinputの1つである read_groups.tsv を生成するルールを記述します。

rule make_read_groups:
    input:
    output:
        "read_groups.tsv"
    run:
        read_groups = samples.assign(v1='*', bam='bam/'+samples['sample']+'_sorted.bam')[['v1', 'bam', 'sample']]
        read_groups.to_csv(output[0], sep='\t', index=False, header=False)

python3によってファイルを生成しています。ここで samples はサンプル情報が格納されたpandasによる表形式のデータです。このデータからサンプル名とそれに対し '_sorted.bam' という文字列をつなげてbamファイルを指定しています。

bamファイルの生成

バリアントコールに必要なもう1つのinputであるbamファイルを生成するルールを記述します。

rule bwa_index:
    input:
        config['ref']
    output:
        expand("{ref}.{ext}", ref = config['ref'], ext = ["amb", "ann", "bwt", "pac", "sa"])
    shell:
        "bwa index {input}"

rule bwa_mem:
    input:
        expand("{ref}.{ext}", ref = config['ref'], ext = ["amb", "ann", "bwt", "pac", "sa"]),
        r1="fastq/{sample}_trimmed_1.fq.gz",
        r2="fastq/{sample}_trimmed_2.fq.gz"
    output:
        "bam/{sample}.bam"
    params:
        ref=config['ref']
    shell:
        "bwa mem {params.ref} {input.r1} {input.r2} | samtools view -bS - > {output}"

rule samtools_sort:
    input:
        "bam/{sample}.bam"
    output:
        "bam/{sample}_sorted.bam"
    shell:
        "samtools sort -o {output} {input}"

bwa mem によるマッピングを行ったあと samtools でソートをかけています。

fastqのトリミング

bwa mem によるマッピングでinputとなっていたトリミング済みfastqを生成するルールを記述します。

rule trimmomatic:
    input:
        lambda wildcards: samples.loc[wildcards.sample].read1,
        lambda wildcards: samples.loc[wildcards.sample].read2
    output:
        "fastq/{sample}_trimmed_1.fq.gz",
        "fastq/{sample}_trimmed_unpaired_1.fq.gz",
        "fastq/{sample}_trimmed_2.fq.gz",
        "fastq/{sample}_trimmed_unpaired_2.fq.gz"
    shell:
        "java -jar trimmomatic-0.39.jar PE {input} {output} "
        "ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 MINLEN:36"

trimmomatic によるトリミングを行います。

その他細かい設定部分

サンプル名とそれらに対応するリードファイル名は samples.tsv に保存しておき、pandas で読み込みます。また、設定ファイル config.yaml も読み込んでおきます。また、最終成果物の指定のために rule all を記述します。

import pandas as pd
from snakemake.utils import validate

configfile: "config.yaml"

samples = pd.read_table("samples.tsv").set_index("sample", drop=False)

rule all:
    input:
        "vcf/calls.bcf"

まとめ

最後に上記ルールを一つのファイルにまとめて完成です。

import pandas as pd
from snakemake.utils import validate

configfile: "config.yaml"

samples = pd.read_table("samples.tsv").set_index("sample", drop=False)

rule all:
    input:
        "vcf/calls.bcf"

rule trimmomatic:
    input:
        lambda wildcards: samples.loc[wildcards.sample].read1,
        lambda wildcards: samples.loc[wildcards.sample].read2
    output:
        "fastq/{sample}_trimmed_1.fq.gz",
        "fastq/{sample}_trimmed_unpaired_1.fq.gz",
        "fastq/{sample}_trimmed_2.fq.gz",
        "fastq/{sample}_trimmed_unpaired_2.fq.gz"
    shell:
        "java -jar trimmomatic-0.39.jar PE {input} {output} "
        "ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 MINLEN:36"

rule bwa_index:
    input:
        config['ref']
    output:
        expand("{ref}.{ext}", ref = config['ref'], ext = ["amb", "ann", "bwt", "pac", "sa"])
    shell:
        "bwa index {input}"

rule bwa_mem:
    input:
        expand("{ref}.{ext}", ref = config['ref'], ext = ["amb", "ann", "bwt", "pac", "sa"]),
        r1="fastq/{sample}_trimmed_1.fq.gz",
        r2="fastq/{sample}_trimmed_2.fq.gz"
    output:
        "bam/{sample}.bam"
    params:
        ref=config['ref']
    shell:
        "bwa mem {params.ref} {input.r1} {input.r2} | samtools view -bS - > {output}"

rule samtools_sort:
    input:
        "bam/{sample}.bam"
    output:
        "bam/{sample}_sorted.bam"
    shell:
        "samtools sort -o {output} {input}"

rule make_read_groups:
    input:
    output:
        "read_groups.tsv"
    run:
        read_groups = samples.assign(v1='*', bam='bam/'+samples['sample']+'_sorted.bam')[['v1', 'bam', 'sample']]
        read_groups.to_csv(output[0], sep='\t', index=False, header=False)
        
rule bcftools_mpileup_call:
    input:
        rg="read_groups.tsv",
        bam=expand("bam/{sample}_sorted.bam", sample = samples['sample'])
    output:
        "vcf/calls.bcf"
    params:
        ref=config['ref']
    shell:
        "bcftools mpileup -G {input.rg} -f {params.ref} {input.bam} | bcftools call -mv -Ob -o {output}"