大阪大学医学部 Python会 (情報医科学研究会)

Now is better than never.

Mockinbirdを用いたPAR-CLIP解析

2019-01-18(Fri) - Posted by 平岡 in 技術ブログ    tag:Bioinformatics

Contents

    今回は PAR-CLIP 解析のAll-in-oneパイプラインソフトウェアである、Mockinbirdを紹介します。日本語での情報がほぼなく、説明が長くなってしまいそうなので、2回に分けて書きます。

    PAR-CLIPって何?

    PAR-CLIPは photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation の略でRNAやmicroRNA結合タンパクの結合サイトを特定するために使います。

    実験法としては、まず4SU(4-チオウリジン)を細胞培地中に添加し、4SUによるラベル標識を行います。次にUVの照射を行い、RNA結合タンパクとRNAの間で架橋反応を起こします。その後、RBPと結合したリードの抽出を行い、サンプルとして用います。

    サンプルRNAを逆転写するときに、4SUがシトシン( C )に置換されます(misread)。そのため、リードがピークを形成し、かつT -> C置換が入っていた場合、RBPと結合していた可能性が高いということになります。これによりあるタンパクが、RNAのどの領域に結合するのか?結合領域間でどのような相互作用があるのかを知ることができます。特に一塩基という高い解像度での解析ができることが特徴です。

    1

    Mockinbirdとは?

    公式ドキュメントGitHubに詳細は記載されていますが、日本語での情報が少ない、というか全くない。のでまとめていきたいと思います。

    PAR-CLIPの解析は生のFASTQデータから、Quality Check -> Trimming -> Mapping -> Downstream Analysisという流れで行いますが、それらすべての解析をMockinbirdでは一気に行ってくれます。また、MockinbirdModule(詳細は後述)を使うことでmock PAR-CLIP experimentを考慮に入れた解析を行うことができるのも大きな特徴です。

    このソフトでは、2つにパイプラインが分かれており、preprocessing phaseとpostprocessing phaseと呼ばれています。タンパク結合サイトの情報が書かれたテーブルを出力するまでが、preprocessing phaseでそれ以降のDownstream解析をpostprocessing phaseと呼んでいます。それぞれにYAMLファイルが用意されており、ユーザーはYAMLファイルで解析したいサンプルデータ、使いたいリファレンスデータ、使いたいモジュール(Mappingの時にSTARやBowtieを選べたりする。)閾値などのパラメータを指定します。最後にコマンドを一行実行するだけで、解析がすべて自動で流れていきます。あるモジュールで出力されたファイルのパスが次のモジュールに自動で設定されていくシステムになっているので、慣れると非常に便利です。

    ただ、ドキュメントでも、既存の解析手法に完全に取って代わるものではなく、PAR-CLIP実験のコンディションのトラブルシューティングを再現可能性が高い状態で行なったり、研究の最初の仮説づくりを短時間で行うために使ってもらうことが目的と書かれていました。

    2

    環境構築

    前置きが長くなりましたが、実際に環境構築から進めていきます。

    conda create -n mockinbird -c bioconda -c conda-forge python=3.6 mockinbird
    

    Anacondaでmockinbirdという名前の仮想環境を作ります。mockinbirdが仮想環境にインストールされることになり、これにより、下記のライブラリとツールが環境内で使えるようになります。

    source activate mockinbird
    

    でmockinbirdの仮想環境を立ち上げ、

    source deactivate mockinbird
    

    で仮想環境を閉じます。

    conda remove --all -n mockinbird
    

    で環境を削除します。

    git clone https://github.com/soedinglab/mockinbird.git
    

    でgithubからクローンします。

    cd mockinbird/mockinbird/data
    

    とディレクトリを進めていくと、preprocess.yaml, postprocess.yamlというファイルがあります。基本はこの二つのYAMLファイルを調整していくことになります。terminalでコマンドを出すときも mockinbird/mockinbird/dataディレクトリで実行するのが個人的にはおすすめです(ファイルの構成が綺麗になる。パスの調整も楽)。

    Tutorialデータがドキュメントの方に書かれていましたが、私が見たときはリンクがNotFoundになっていました。ので、私は公式ドキュメントGitHubを参考にYAMLファイルの調整などを行いました。

    preprocess, postprocessの実行のコードは以下のようになっています。

    # mockinbird preprocess [parclip_fastq] [output_dir] [prefix] [config_file]
    $ mockinbird preprocess nrd1.fastq nrd1 nrd1 preprocess.yaml
    
    # mockinbird postprocess {{genomefasta}} {{output_dir}} {{output_dir}} {{script_dir}}
    $ mockinbird postprocess nrd1 nrd1_pp postprocess.yaml
    

    ですが、コードの実行はYAMLの設定ができてからになります。ということでYAMLの中身について説明していきます。

    preprocess.yaml

    preprocess.yamlの4つの区画から構成されています。

    • ① 変数設定。
    • ② general. 必須情報の設定。
    • ③ reads. リードについての情報。
    • ④ pipeline. 使うモジュールと各引数の設定。

    ① 変数設定 下記のコードはドキュメントから引用していますが、モジュールで使う、ディレクトリの変数化と、mock_processing をFalseに設定しています。(ここ重要!)

    {% set data_dir = "data" %}
    
    {% set genome_fasta = data_dir + "/genome.fa" %}
    {% set mock_pileup = data_dir + "/mock.mpileup" %}
    {% set mock_statistics = data_dir + "/mock_stat.json" %}
    {% set norm_pileup = data_dir + "/normalization.mpileup" %}
    {% set bowtie_index = data_dir + "/bowtie_index/genome" %}
    
    # setting mock_processing to True will only process the mock. Setting to `False` will run the full
    # pipeline
    {% set mock_processing = False %}
    

    ② general. 必須情報の設定 アダプター配列の指定、リファレンスデータの指定、UMI(Unique molecular identifiers)の有無、スレッド数の指定。

    ③ reads. リードについての情報 リードの最短長、T -> C mutationの指定

    general:
        adapter5prime: GTTCAGAGTTCTACAGTCCGACGATC
        adapter3prime: TGGAATTCTCGGGTGCCAAGG
        genomefasta: {{ genome_fasta }}
        normalization_pileup: {{ norm_pileup }}
        rmTemp: yes
        n_threads: 4
    
    reads:
        bc_5prime: 5
        bc_3prime: 0
        min_len: 20
        reference_nucleotide: T
        mutation_nucleotide: C
    

    ④ pipeline. 使うモジュールと各引数の設定 クオリティチェック -> トリミング -> マッピング

    pipeline:
      - FastQCModule:
          outdir_name: fastQC_raw
    
      - UmiToolsExtractModule
    
      - SkewerAdapterClippingModule
    
      - ClippyAdapterClippingModule:
          clipped_5prime_bc: True
    
      - FastQCModule:
          outdir_name: fastQC_clipped
          genome_index: {{ bowtie_index }}
    
      - BamPPModule:
          remove_n_edge_mut: 2
          max_mut_per_read: 1
          min_mismatch_quality: 20
    
      - SortIndexModule:
         keep_all: yes
    
      - UmiToolsDedupModule
    
      - SortIndexModule:
         keep_all: yes
    
      - PileupModule:
         keep_all: yes
    
      - BamStatisticsModule
    

    MockinbirdModuleを使うことでmock PAR-CLIP experimentを考慮に入れた解析ができますと、すでに書きましたが、Mock Experimentの結果をバックグラウンドとして処理して、PAR-CLIP Experimentの結果を出すことができます。

    具体的な方法としては、まず mock_processing = Trueとし、 mockinbird preprocess nrd1_mock.fastq nrd1_mock nrd1_mock preprocess.yamlを実行します。これで{% if not mock_processing %}までで解析がストップします。

    PileupModule, BamStatisticsModuleによって、Mock Experimentの.mpileupファイルと_stat.jsonファイルがそれぞれ生成されます。(この二つのファイルを次に使うことになります。)

    次にmock_processing = Falseに戻し、 $ mockinbird preprocess nrd1.fastq nrd1 nrd1 preprocess.yamlを実行します。すると下記のコードも一緒に実行されます。

    {% if not mock_processing %}
    
        - PredictionSitesModule:
            sites_file: {{ data_dir }}/genome.sites
            fasta_file: {{ genome_fasta }}
            transition_nucleotide: T
    
        - MockTableModule:
            mock_table: {{ data_dir }}/mock.table
            mock_pileup: {{ mock_pileup }}
    
        - TransitionTableModule
    
        - LearnMockModule:
            mock_model: mock_model/model.pkl
            mock_statistics: {{ mock_statistics }}
            n_mixture_components: 5
            em_iterations: 250
    
        - MockinbirdModule
    
        - NormalizationModule
        - QuantileCapModule
    {% endif %}
    

    これによって、Mock Experimentの.mpileupファイル_stat.jsonファイルを受け取りつつ、最終的なCLIP-seqのテーブルを作ることができ、preprocessは終了となります。下記のようなテーブルが出力され、これをもとにpostprocessに進んでいくことになります。

    seqid   position        transitions     coverage        score   strand  occupancy       posterior
    chrI    32601   4       4       6.352710158543762       -       0.15384615384615385     0.772109761967
    chrI    35562   4       5       5.292152965717071       +       0.0223463687150838      0.539841709238
    chrI    35805   5       5       9.376254703871467       +       0.03184713375796178     0.98585024515
    

    長くなってしまったので、モジュールとpostprocessの説明は次回に回したいと思います。長い間お付き合いいただきありがとうございました。