PAR-CLIP解析データ（Peak calling, Motif discovery, Peak annotation）

NGS PAR-CLIP CLIP-seq RNA結合タンパク質 RBP 次世代シーケンス

概要

CLIP-seqとは、RNA結合タンパク質（RBP）の結合サイトを網羅的に同定する手法です。CLIP-seqには主に、HITS-CLIP、iCLIP、PAR-CLIP、eCLIPの４種類が存在します。今回は、中でも幅広く利用されているPAR-CLIPのデータについて取り上げたいと思います。

ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、PAR-CLIPのデータ解析の詳細をStep by Stepで説明します。

ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。

このTutorialを通した到達目標

PAR-CLIPのデータ解析に対する理解

  • SRAファイルを解凍してFASTQファイルを生成することができる。
  • アダプター配列を各リードからトリミングすることができる。
  • リピート配列に該当するリードをFASTQファイルから除外することができる。
  • FASTQファイルをゲノムとトランスクリプトームにマッピングすることができる。
  • マッピングされたリードのStrandnessを確認することができる。
  • BAMファイルからBase Substitution frequency（T>C置換の割合）を計算することができる。
  • UCSC genome browser上でデータを可視化することができる。
  • PAR-CLIPのデータからPARalyzerを用いてPeak callを行うことができる。
  • MEMEでモチーフ検索を行うことができる。

NGSで扱う種々のデータに対する理解

  • FASTQ, BAM/SAM, BedGraph, BED, GTFがそれぞれどのようなデータか説明することができる。
  • （特に）SAMファイルからMutationの情報を取得することができる。
  • FastQCから得られたシークエンスリードのクオリティデータから、FASTQファイルの状態を確認できる。
  • STARのIndexファイルを作成することができる。
  • MEMEに用いるMarkov Backgroundモデルを作成することができる。

PAR-CLIP法の特徴

シークエンスのデータを解析を始める前に、PAR-CLIPの手法の特徴について概説します。

PAR-CLIPでは一般的に、ウリジンのアナログである4-チオウリジン（4SU）を細胞培地中に添加し細胞に取り込ませ、4SUを含むRNAを細胞内で転写させます（4SUによるラベル標識）。

次に、長波長の紫外線（UV 365 nm）を細胞に照射させることにより、RNA結合タンパク質とRNAとの間で架橋反応（共有結合）を起こします。

RNAase処理、特定のRNA結合タンパク質（RBP）をIP、Protease K処理した後、精製して得られたRNA（RBPと結合しているサイト）を次世代シーケンサーのサンプルとして使用します。

精製したRNAからDNAライブラリを作成する過程（逆転写・cDNA合成）で4SUがシトシン（C）に置換されます。そのため、マッピング後に検出される各ピークを形成するリードにこのT→C置換が入っていた場合、その領域でRBPと架橋していた可能性が高いといえます。

このことから、PARalyzerを用いたPeak callでは、このT→C置換を指標にNon-specificでないピークを検出しています。

使用するデータ

Genome-wide analysis of pre-mRNA 3' end processing reveals a decisive role of human cleavage factor I in the regulation of 3' UTR length: CLIP

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE37398

上記のデータのうち、3'endのcleavageに関与するRNA結合タンパク質である「CstF-64」のPAR-CLIPデータを解析する。詳細は下記に載っている。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM917676

使用するソフトウェア

解析ワークフロー

今回の解析の流れを以下にまとめました。 PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-01-19 14.53.30.png (190.0 kB) 上記の解析はENCODEプロジェクトのeCLIP解析を参考に作成しました。
https://www.encodeproject.org/documents/dde0b669-0909-4f8b-946d-3cb9f35a6c52/@@download/attachment/eCLIP_analysisSOP_v1.P.pdf

サンプルコード

Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

シェルスクリプトを使用するときの注意点

権限の問題

実行権限をシェルスクリプトに与えた後、実行すること。

$ chmod 755 <hoge.sh>

改行コードの問題

WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。

ここではAtomと呼ばれるテキストエディタでの設定方法を説明する。

  • [Ctrl] + [ , ]（カンマ）キーを押して、Settingsの画面を呼び出す。
  • 左のメニューから、Packagesをクリック。
  • Installed Packagesでline-ending-selectorと検索。
  • Core Packagesの項目から、line-ending-selectorのSettingsをクリック。
  • 設定画面で、Default line endingをLFに変更。

Settings — C__Users_Naoto_OneDrive_NGS解析_NGS-Tutorial_RNA-seq_Tutorial — Atom 2017-01-31 11.20.21.png (250.6 kB)

これで、新規に作成したファイルの改行コードがデフォルトでLF（Linuxの改行コード）になります。

各ステップの説明

0. 必要なソフトウェアのインストール

まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。

imamachi-n.hatenablog.com

$ conda install sra-tools
$ conda install fastqc
$ conda install cutadapt
$ conda install fastx_toolkit
$ conda install bowtie2
$ conda install star
$ conda install samtools
$ conda install bedtools
$ conda install paralyzer
$ conda install meme

1. SRAファイルのダウンロード

NGSデータはシークエンスされたリード（塩基配列）の情報とクオリティスコアのデータをまとめたFASTQ形式のファイルとして保存されています。

FASTQファイルの詳細についてはこちらを参照。
https://en.wikipedia.org/wiki/FASTQ_format

FASTQデータは非常に大きいファイルなので、NCBI GEOでは、FASTQファイルを圧縮したSRAファイルと呼ばれるファイル形式でデータを公開しています。

まず、下記のURLにアクセスすると、
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM917676

サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。

GEO Accession viewer - Google Chrome 2017-01-13 11.12.08.png (273.9 kB)

そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。

GEO Accession viewer - Google Chrome 2017-01-13 11.12.38.png (141.2 kB)

SRR488740をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。

_sra_sra-instant_reads_ByExp_sra_SRX_SRX143_SRX143276 のインデックス - Google Chrome 2017-01-14 23.22.04.png (44.8 kB)

SRR488740.sraのURLをコピー。

_sra_sra-instant_reads_ByExp_sra_SRX_SRX143_SRX143276_SRR488740_ のインデックス - Google Chrome 2017-01-14 23.22.50.png (50.9 kB)

wgetで該当するSRAファイルをダウンロードします。

$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX/SRX143/SRX143276/SRR488740/SRR488740.sra

2. SRAファイルからFASTQファイルに変換

先ほどインストールしたSRA Toolkitを使って、NCBI GEOからダウンロードしたSRAファイルをFASTQファイルに展開します。

$ fastq-dump SRR488740.sra
$ mv SRR488740.fastq SRR488740_PAR-CLIP_CstF-64.fastq

カレントディレクトリにSRR488740.fastqファイルが展開されているはずです。次に、mvコマンドを使ってファイル名をSRR488740.fastqからSRR488740_PAR-CLIP_CstF-64.fastqに変更します（SRRのIDだけだと何のファイルかわからなくなるので）。

3. Quality check

FastQCを使って、FASTQファイルに記載されている各リードのクオリティをチェックします。

FASTQファイルはPhred Quality Scoreと呼ばれる指標によって、各塩基のクオリティを記述しています。

詳しくはこちらのサイトを参照のこと。 https://en.wikipedia.org/wiki/Phred_quality_score

実行例

$ mkdir fastqc_SRR488740_PAR-CLIP_CstF-64
$ fastqc -t 8 -o ./fastqc_SRR488740_PAR-CLIP_CstF-64 ./SRR488740_PAR-CLIP_CstF-64.fastq -f fastq

mkdirコマンドでfastqc_SRR488740_PAR-CLIP_CstF-64という名前のディレクトリをカレントディレクトリに作ります。

次に、fastqcを使ってSRR488740_PAR-CLIP_CstF-64.fastqのデータのクオリティをチェックします。 -tオプションで使用するCPUのコア数を指定します。コア数を増やせば単純に計算速度が向上しますが、FastQCではメモリを1コアあたり250MB消費するので、要求できるメモリ量を超えるコア数を指定するとエラーが起こります。

-oオプションで出力先のディレクトリを指定し、-fオプションでInputファイルのファイル形式を指定しています。

FastQCから出力されるデータが、fastqc_SRR488740_PAR-CLIP_CstF-64ディレクトリの中のSRR488740_PAR-CLIP_CstF-64_fastqc.htmlになります。 ファイルをブラウザで開くと以下のようなグラフを確認できます。

Basic Statistics

基本的な統計量が示されています。

  • 全リード数
  • リード長
  • GC contents

などなど

Per base sequence quality

リードの各塩基のクオリティスコアを示しています。 Phred quality scoreがだいたいグリーンの領域（Scoreが28以上）に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。

ダウンロード.png (8.5 kB)

Per tile sequence quality

フローセルの各タイルごとのクオリティスコアを示しています。

Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。

シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。

特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。

詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf

ダウンロード (1).png (8.6 kB)

Per sequence quality scores

各リードの全長のクオリティスコアの分布を示しています。 ダウンロード (2).png (16.2 kB)

Per base sequence content

各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。

ダウンロード (3).png (33.0 kB)

Per sequence GC content

リードのGC contentsの分布を示しています。 ダウンロード (4).png (32.4 kB)

Per base N content

各塩基中に含まれるNの含有率（塩基を読めなかった箇所）を示しています。 ダウンロード (5).png (7.2 kB)

Sequence Length Distribution

リード長の分布を示しています。 ダウンロード (10).png (20.3 kB)

Sequence Duplication Levels

Duplidate readsの含まれている数を示しています。 ダウンロード (6).png (22.9 kB)

Overrepresented sequences

頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。

Adapter Content

各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。 ダウンロード (7).png (8.8 kB)

Kmer Content

特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。 ダウンロード (8).png (60.2 kB)

4. Adapter trimming

PAR-CLIPのリードはアダプター配列の一部が含まれているので、Cutadaptを利用してそれを除きます。 アダプター配列の設定はサンプルに応じて変更が必要です。

どのようなアダプター配列がリードに混入しているかどうか調べるには、ダウンロード時に参照したNCBI GEOの各サンプルのページから得るか、もしくは論文のMaterial & Methodsを参照して得ます。

今回は、NCBI GEOのサイトにアダプター配列の情報が載っているので、それを使います。

GEO Accession viewer - Google Chrome 2017-01-16 09.40.49.png (270.9 kB)

まず、TCGTATGCCGTCTTCTGCTTGTというアダプター配列がそもそもFASTQファイル内に保存されている各リードに含まれているかどうか、目で確認したいと思います。

$ less SRR488740_PAR-CLIP_CstF-64.fastq

ここでは、lessコマンドでFASTQファイルの中身を見ています。lessで表示したとき、十字キーでさらに下の情報も見ることができます。

また、/に続けて文字列を入力してEnterを押すと、その文字列を含む列を検索することができます。アダプター配列全長を検索にかけるとヒットしにくくなるので、5'end側から6-7bpぐらいの配列で検索にかけるとよいかもしれません（TCGTATGとか）。

slogin.hgc.jp_22 - Tera Term VT 2017-01-18 17.29.15.png (126.9 kB)

検索した文字列と一致する箇所に、このようにハイライトがつきます。

slogin.hgc.jp_22 - Tera Term VT 2017-01-18 17.30.28.png (165.7 kB)

上図では、今回のサンプルとは異なるサンプルを見ています。ちなみにCstF64に対するPAR-CLIPデータでは、20%程度しかアダプター配列が含まれていません。これは、そもそもシークエンスしたリード長が38bpと短いため、アダプター配列の混入率が低かったと考えられます。

実行例

$ cutadapt -f fastq --match-read-wildcards --times 1 -e 0.1 -O 5 --quality-cutoff 6 -m 18 \
-a TCGTATGCCGTCTTCTGCTTGT \
SRR488740_PAR-CLIP_CstF-64.fastq > SRR488740_PAR-CLIP_CstF-64_1_trimmed_adapter.fastq 2>> ./log_SRR488740_PAR-CLIP_CstF-64.txt
  • -f: フォーマットの指定。
  • --match-read-wildcards: IUPACのワイルドカードを許容。
  • --times: 各リードに対してトリミングを実行する回数。
  • -e: 指定したアダプター配列とのミスマッチ率の許容範囲。
  • -O: アダプター配列とOverlapする最低配列長。例えば、5と指定すると、最低5塩基分Overlapしている必要がある。
  • --quality-cutoff: クオリティの低い塩基を除く。
  • -m: トリミング後のリードの最低配列長。例えば、18と指定すると、トリミング後にリード長が18未満のリードは排除する。
  • -a: 3'側にライゲーションされている、指定したアダプター配列を除去する。

5. Quality filtering

FASTQファイル中にクオリティスコアの低いリードが含まれていることが確認されたので、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。

コード例

$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./SRR488740_PAR-CLIP_CstF-64_1_trimmed_adapter.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o SRR488740_PAR-CLIP_CstF-64_2_filtered.fastq
fastq_quality_trimmerのパラメータ。
  • -Q: Phred quality scoreのバージョン指定（Solexa, Illumina 1.3+, Illumina 1.5+といったかなり古いバージョンのFASTQファイルでは'64'を指定。最近のFASTQファイルは'33'と指定しておけばOK）。
  • -t: 指定したPhred quality scoreを下回る塩基をトリミングする。
  • -l: トリミング後の最低リード長。
fastq_quality_filterのパラメータ。
  • -q 20 -q 80: リードの80%でPhred quality scoreが20を下回る場合、そのリードを排除する。

6. リボソームRNAとRepetitive elementsを除く

次に、リボソームRNAとRepetitive elements (RepBase databaseから得た配列)由来のリードを除去します。

Repetitive elementsの用意

RepBaseと呼ばれるデータベースから、Repetitive elementsの配列情報を取得します。 下記のURLから最新のRepBaseのデータ (FASTAファイル)を入手することができます。
ただし、データの取得にはRegistrationが必要になるので注意してください。
http://www.girinst.org/server/RepBase/index.php

FASTA format (42.11 MB) 01-24-2017:と書いてあるほうをダウンロードします。
http://www.girinst.org/server/RepBase/protected/RepBase22.01.fasta.tar.gz

ファイルを展開して、

$ tar xzvf RepBase22.01.fasta.tar.gz

appendix/humapp.refとhumrep.refとhumsub.refをマージします。ファイル名はRepBase_human_v22_01.faとしておきます。

$ cat humapp.ref  humrep.ref  humsub.ref > RepBase_human_v22_01.fa

次に、リボソームRNAの配列を用意します。 NCBIの登録されているX12811.1| Human 5S DNAとU13369.1|HSU13369 Human ribosomal DNA complete repeating unitの配列情報を使います。 下記のURLからFASTA形式のデータをコピペしてFASTAファイルをつくります。

Human 5S DNA

https://www.ncbi.nlm.nih.gov/nuccore/23898

Human ribosomal DNA complete repeating unit

https://www.ncbi.nlm.nih.gov/nuccore/555853

上記のrRNAのFASTA形式のデータをcontam_Ribosomal_RNA.faとしてまとめ、先ほど取得したRepetitive elementsのFASTAファイルとマージさせます。名前をRepBase_v22_01_human_rRNA.faとしておきます。

cat RepBase_human_v22_01.fa contam_Ribosomal_RNA.fa > RepBase_v22_01_human_rRNA.fa

rRNAやRepetitive elementsを除去するために、FASTQファイルに格納されている各リードをそれらの配列にアライメント（マッピング）して、それらの配列にマッピングされないリードを抽出します。こうすることで、rRNAやRepetitive elementsに由来しないリードのみをもとのFASTQファイルから抽出することができます。

STARと呼ばれるソフトでマッピングさせるために、まず、rRNAやRepetitive elementsのIndexファイルを用意する必要があります。以下で、先ほど作成したFASTAファイルをもとに、STAR用のIndexファイルを作成しています。

実行例

$ mkdir ./STAR_Index_repBase_v22_01_rRNA_contam
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./STAR_Index_repBase_v22_01_rRNA_contam \
     --genomeFastaFiles ./RepBase_v22_01_human_rRNA.fa

まず、Indexファイルを保存するディレクトリをmkdirコマンドで作成します。ここでは、STAR_Index_repBase_v22_01_rRNA_contamという名前のディレクトリを用意しました。

STARのパラメータ。
  • --runThreadN: 使用するコア数を指定。
  • --runMode: genomeGenerateを指定して、Indexファイルを作成するモードに設定。
  • --genomeDir: 出力先のディレクトリを指定。
  • --genomeFastaFiles: Indexを作成したいFASTAファイルを指定。ここでは、先ほど用意したリピート配列とrRNAの配列データ（FASTAファイル）を指定する。

rRNAとRepetitive elementsを除く

Indexファイルを作成した後、rRNAとRepetitive elementsへのマッピングを行います。マッピングされなかったリードについては、FASTQファイルとして出力されるようにオプションを指定します。

実行例

$ STAR --runMode alignReads --runThreadN 8 --genomeDir ./STAR_Index_repBase_v22_01_rRNA_contam \
--readFilesIn ./SRR488740_PAR-CLIP_CstF-64_2_filtered.fastq \
--outSAMunmapped Within --outFilterMultimapNmax 30 --outFilterMultimapScoreRange 1 \
--outFileNamePrefix ./SRR488740_PAR-CLIP_CstF-64_3_rm_repbase_rrna.fastq \
--outSAMattributes All --outStd BAM_Unsorted --outSAMtype BAM Unsorted \
--outFilterType BySJout --outReadsUnmapped Fastx --outFilterScoreMin 10 \
--alignEndsType EndToEnd > ./SRR488740_PAR-CLIP_CstF-64_3_repbase_rrna_comtam.bam
$ mv ./SRR488740_PAR-CLIP_CstF-64_3_rm_repbase_rrna.fastqUnmapped.out.mate1 ./SRR488740_PAR-CLIP_CstF-64_3_rm_repbase_rrna.fastq
STARのパラメータ。
  • --runMode: alignReadsを指定して、アライメントを行うモードに設定。
  • --runThreadN: 使用するコア数を指定。
  • --genomeDir: Indexファイルの保存されているディレクトリを指定。
  • --readFilesIn: マッピングさせたいFASTQファイルを指定。
  • --outSAMunmapped Within: 出力されるBAMファイルにマッピングされなかったリード情報も含める。
  • --outFilterMultimapNmax: 複数箇所にマッピング（Multi-mapping）されることをどこまで許容するか設定する。指定した数を超える箇所にマッピングされたリードは排除する。
  • --outFilterMultimapScoreRange: Multi-mappingアライメントのスコアの範囲。複数箇所マッピングされてもアライメントスコアが高ければ残す。
  • --outFileNamePrefix: 出力ファイルのPrefixを指定。
  • -outSAMattributes: 基本情報以外に、BAMファイルに追加する情報を指定。
  • --outStd BAM_Unsorted: 出力されるBAMファイルをソートしない。
  • --outSAMtype BAM Unsorted: SAMファイルでなく、BAMファイルとしてマッピングされた情報を直接出力する。ただし、ソートはしない。
  • --outFilterType BySJout: SJ.out.tabファイルとして出力される信頼性の高いExon junctionをまたぐリードのみを残す。Filteringをクリアできないリードは排除する。
  • --outReadsUnmapped Fastx: マッピングされなかったリードをFASTQファイルとして出力する。
  • --outFilterScoreMin 10: アライメントスコアが10を下回るリードは排除する。
  • --alignEndsType EndToEnd: End-to-Endアライメントを行う（Localアライメントでない）。

7. Quality check

rRNAやRepetitive elementsに由来するリードを除いた後、フィルタリング後のFASTQファイルをFastQCを用いて再確認します。

実行例

$ mkdir fastqc_${file}_filtered
$ fastqc -o ./fastqc_SRR488740_PAR-CLIP_CstF-64_filtered ./SRR488740_PAR-CLIP_CstF-64_3_rm_repbase_rrna.fastq -f fastq

8. Mapping

PAR-CLIPから得られるリード長は短い (25 ntぐらい)ので、リードの複雑度によってはゲノム上に重複してマッピングされるおそれがあります。そのため、ここではゲノムとトランスクリプトームにユニークにマッピングされるリードしか扱わないよう、オプションを指定しています。

Indexファイルの作成

ここでも先ほどのマッピング時と同じように、GenomeとTranscriptomeのSTAR用のIndexファイルを作成する必要があります。このIndexファイルの作成には、50G程度のメモリが必要になるので注意してください。

まず、Human Genome（hg19）のFASTAファイルのダウンロードしてきます。 各染色体ごとにファイルが用意されているので、catコマンドでダウンロードしたすべてのファイルをマージして、hg19.faという名前で出力します。

$ for file in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X Y M
  do
    wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/chr${file}.fa.gz
    gunzip chr${file}.fa.gz
  done

$ cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chr20.fa chr21.fa chr22.fa chrX.fa chrY.fa chrM.fa > hg19.fa

Transcriptomeの情報として、Gencodeのアノテーション情報を利用します。
https://www.gencodegenes.org/releases/19.html

上記のサイトにアクセスすると、さまざまなデータが置いてあるが基本的にgencode.v19.annotation.gtfというメインのアノテーション情報を利用すればOKです。
ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

wgetコマンドで該当のファイルをダウンロードし解凍します。

$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
$ tar zxvf ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

カレントディレクトリにgencode.v19.annotation.gtfが保存されます。

最後に、STARを使ってhg19+Gencode_v19のIndexファイルを作成します。

実行例

$ mkdir hg19_Genome_Gencode_v19
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir hg19_Genome_Gencode_v19 \
     --genomeFastaFiles ./hg19.fa --sjdbGTFfile ./gencode.v19.annotation.gtf --sjdbOverhang 100
STARのパラメータ。
  • --runThreadN: 使用するコア数を指定。
  • --runMode: genomeGenerateを指定して、Indexファイルを作成するモードに設定。
  • --genomeDir: 出力先のディレクトリを指定。
  • --genomeFastaFiles: Indexを作成したいFASTAファイルを指定。ここでは、先ほど用意したヒトゲノムのデータ（hg19、FASTAファイル）を指定する。
  • --sjdbGTFfile: STARではアノテーション情報のGTFファイルからExon junctionの情報を抽出する。ここでは、先ほどダウンロードしたgencode.v19.annotation.gtfを使用する。
  • --sjdbOverhang 100: Exon junctionの周囲、100塩基のゲノム配列をIndex化する。マッピングに使用するFASTQファイルのリード長と設定すると良い。公式のおすすめは100。

ゲノムとトランスクリプトームへのマッピング

Indexファイルが用意できたので、ゲノム（hg19）とトランスクリプトーム（Gencode v19）にリードをマッピングします。

実行例

$ mkdir STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd
$ STAR --runMode alignReads --runThreadN 8 --genomeDir ./hg19_Genome_RefSeq \
--readFilesIn ./SRR488740_PAR-CLIP_CstF-64_3_rm_repbase_rrna.fastq \
--outFilterMultimapNmax 1 -- outFilterMultimapScoreRange 1 \
--outFileNamePrefix ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_ \
--outSAMattributes All --outSAMtype BAM SortedByCoordinate --limitBAMsortRAM 32000000000 \
--outFilterType BySJout --outReadsUnmapped Fastx --outFilterScoreMin 10 \
--alignEndsType EndToEnd
STARのパラメータ。
  • --outSAMtype BAM SortedByCoordinate: SAMファイルでなく、BAMファイルとしてマッピングされた情報を直接出力する。出力ファイルは染色体の座標でソートする。
  • --limitBAMsortRAM 32000000000: BAMファイルをソートするときに割り当てられるメモリの上限を設定。

マッピング後に、samtoolsを用いてBAMファイルをSAMファイルに変換します（この後の作業で、SAMファイルを使っていきます）。

実行例

$ samtools view -h ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.bam \
  > ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.sam

9. Base Substitution frequencyの確認

PAR-CLIPの場合、TからCの置換が見られるはずなので、SAMファイルからリファレンスゲノムと比較して、ミスマッチしている塩基配列を抽出し、各塩基のミスマッチの割合を調べます。

SAMファイルの中身について

詳しくは下記のURL・PDFファイルを参照のこと。
https://samtools.github.io/hts-specs/
https://samtools.github.io/hts-specs/SAMv1.pdf
https://samtools.github.io/hts-specs/SAMtags.pdf

Mandatory field

SAMファイルの各行には、以下のような必須項目が含まれています。 image.png (194.6 kB) ミスマッチの情報は、CIGAR stringとMD fieldに記載されています。 CIGAR stringは「6カラム目」に記載されており、MD fieldはオプションとして各行の後半に記載されています。

SAMファイルの中身の例

SRR488740_PAR-CLIP_CstF-64.13819410     256     chr1    11488   0       38M     *       0       0       ATGCTGTTTGGTCTCAGTAGACTCCTAAATATGGGACT  BCCCCCBCCCCBBBBACABB>BBBBCBAAB?AABBBBB  NH:i:7  HI:i:2  AS:i:35 nM:i:1  NM:i:1  MD:Z:36T1       jM:B:c,-1       jI:B:i,-1
SRR488740_PAR-CLIP_CstF-64.462991       16      chr1    12953   0       38M     *       0       0       GGTTCATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCA  ;='9@@;9<2@5@9@BB<5BB9<)B=+=9)&7B5BB@A  NH:i:6  HI:i:1  AS:i:35 nM:i:1  NM:i:1  MD:Z:2G35       jM:B:c,-1       jI:B:i,-1
SRR488740_PAR-CLIP_CstF-64.14573742     272     chr1    14320   0       10M1D28M        *       0       0       CAACAGGGGCGGAGGCAGTCACTGACCCCGAGACGTTT  >B;1=A@B5<@=2A<<=BAB:A6B;B2;BA>ABCBACB  NH:i:8  HI:i:4  AS:i:33 nM:i:0  NM:i:1  MD:Z:10^A28     jM:B:c,-1       jI:B:i,-1
SRR488740_PAR-CLIP_CstF-64.2153554      272     chr1    17024   0       32M177N6M       *       0       0       AGGTCTGGCACATAGAGGTAGTTCTCTGGGACCTGCTG  55A=AA@;@7B@;@B=??B>??>B9BA=?A9ABBABB@  NH:i:7  HI:i:7  AS:i:36 nM:i:1  NM:i:1  MD:Z:16A21      jM:B:c,22       jI:B:i,17056,17232
SRR488740_PAR-CLIP_CstF-64.15174491     256     chr1    360396  0       19M1I18M        *       0       0       TCAGGCACCGTGGAGAAAAAGTCGCCGTGTGTAGGCAG  BBB?ABBB@BABBABBBCBBB?ABAAB<B<>29=AA;A  NH:i:8  HI:i:4  AS:i:32 nM:i:0  NM:i:1  MD:Z:37 jM:B:c,-1       jI:B:i,-1

上記の表と照らし合わせると、各列の情報は以下の通りになります。 ここでは、上記のSAMファイルの１行目の情報を見ています。

  • QNAME: SRR488740_PAR-CLIP_CstF-64.13819410
  • FLAG: 256
  • RNAME: chr1
  • POS: 11488
  • MAPQ: 0
  • CIGAR: 38M
  • RNEXT: *
  • PNEXT: 0
  • TLEN: 0
  • SEQ: GGTTCATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCA
  • QUAL: ;='9@@;9<2@5@9@BB<5BB9<)B=+=9)&7B5BB@A
  • オプション領域: NH:i:6 HI:i:1 AS:i:35 nM:i:1 NM:i:1 MD:Z:2G35 jM:B:c,-1 jI:B:i,-1

STARでマッピングして得られたSAMファイルの中身を見てみると、CIGAR stringには38M, 10M1D28M, 32M177N6M, 19M1I18Mなどがあります。

CIGER stringには以下の種類があり、 image.png (100.1 kB) 38Mは、ミスマッチを含めて38塩基長のリードであることを示しています。

10M1D28Mは、最初の10塩基（10M）と後半の28塩基（28M）はミスマッチを含めてアライメントされており、その間にリード中にDeletionの箇所（1D）が存在することを意味しています。

An image of meaning of CIGAR    <---10--->D<--------28-------->
D: Deletion from the reference

32M177N6Mは、最初の32塩基（32M）と後半の6塩基（6M）はミスマッチを含めてアライメントされており、その間はExon junctionである（177N）ことを意味しています。RNA-seqで今回のようにExon junctionを考慮してリードをマッピングしている場合、このようなCIGAR stringが散見されます。

An image of meaning of CIGAR    <--------32-------->NNNNN...NNNNN<-6->
                                [------exon_1------]             [-exon_2-]
N: Skipped region from the reference

19M1I18Mは、最初の19塩基（19M）と後半の18塩基（18M）はミスマッチを含めてアライメントされており、その間にリファレンスゲノムに特定の塩基が挿入されている（1I）ことを意味しています。

An image of meaning of CIGAR    <-----19-----> I <-----18----->
I: Insertion region from the reference

今回例示したPythonスクリプトでは、データマイニングを簡単にするために、Deletion、Exon juction、Insertionを含むリードを解析対象から除外しています。つまりは、38Mのようなミスマッチ箇所のみを含むリードだけを抽出して解析しています。

CIGAR stringだけでは、いったいATGCのどの塩基からどの塩基へのMutationが起こっているかの情報がわかりません。そこで、さらにMD fieldを参照します。

先ほどのリードでCIGAR stringとMD fieldを比較すると以下のようになります。 38Mに対して、MD:Z:36T1というMD fieldが記述されています。

An image of meaning of CIGAR    <-------36------->N<-1->    :Sequence Read
An image of meaning of MD tag   <-------36------->T<-1->    :Reference genome
N: A, T, G, Cのいずれかの塩基

これは、リファレンスゲノムでTだった箇所で、リード中で塩基のミスマッチが見られるということを意味しています。

10M1D28Mに対して、MD:Z:10^A28という記述がされています。

An image of meaning of CIGAR    <--10-->D<-------28------->    :Sequence Read
An image of meaning of MD tag   <--10-->A<-------28------->    :Reference genome
D: Deletion from the reference

これは、リファレンスゲノムでAだった箇所で、リード中で塩基の欠失が見られるということを意味しています。

19M1I18Mに対して、MD:Z:37という記述がされています。

An image of meaning of CIGAR    <-------36------->N<-6->    :Sequence Read
An image of meaning of MD tag   <-------36------->D<-6->    :Reference genome
D: Deletion from the reference
N: A, T, G, Cのいずれかの塩基

これは、リファレンスゲノムの対して、リード中に特定の塩基の挿入が見られるということを意味しています。

必要なスクリプト

以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

$ python E_mismatch_call_from_bam.py ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.sam \
./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_indel.txt \
Base_substitution_frequency.txt
  • 1つ目の引数: 解析するSAMファイルを指定。
  • 2つ目の引数: ミスマッチを含むリードを抽出した中間ファイルを出力。
  • 3つ目の引数: 各塩基のミスマッチの数を出力。

結果の図

上記のスクリプトから得られた各塩基の置換の割合をもとに、Excelを用いて結果を図示した。STARでのマッピング条件によって結果が異なります。

STARでのマッピング条件（PARalyzer default）

--outFilterMultimapNmax 10

image.png (19.0 kB)

STARでのマッピング条件（ENCODE eCLIP）

--outFilterMultimapNmax 1 -- outFilterMultimapScoreRange 1

image.png (19.5 kB)

予想通り、T>C substitutionが起こっているリードがEnrichされていることがわかります。この後行うPARalyzerによるPeak callは、このT>C substitutionを指標にしてピークを検出するので、T>C substitutionが高いことが以降の解析で重要となってきます。

上記の結果を見てみると、ユニークにマッピングされたリードのみを使ったほうが、T>C substitutionの割合が大きいことがわかります。このことから、今回の解析では--outFilterMultimapNmax 1 -- outFilterMultimapScoreRange 1の条件でマッピングを行ったほうがよいと考えました。

10. マッピングされたリードのStrandnessの確認

PAR-CLIPではリードがStrand-specificにマッピングされると想定されることから、リードのStrandnessを確認します。ここでは、RSeQCのinfer_experiment.pyというスクリプトを用いました。

アノテーション情報としてBEDファイルを用意

下記のURLに置いてあるgtf2bed.plを使用します。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

$ perl gtf2bed.pl gencode.v19.annotation.gtf gencode.v19.annotation.bed
  • 1つ目の引数: 変換したいGTFファイルを指定
  • 2つ目の引数: 出力するBEDファイルを指定。

Strandnessの確認

$ infer_experiment.py -i ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.bam -r ./gencode.v19.annotation.bed -s 400000 >  ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_strand_stat.txt
  • -i: Strandnessを調べるBAMファイルを指定。
  • -r: 先ほど用意したBEDファイルのアノテーション情報を指定。
  • s: 解析に使用するリードを指定。ここでは、40万リードを抽出してStrandnessを調べています。指定する数値を大きくすれば計算精度も上がると思いますが、その分時間もかかります。

下記がその結果になります。"++“や”+-“のように2文字で表現されていますが、1つ目がアノテーション情報でのStrandの向き、2つ目がマッピングされたリードのStrandの向きを表しています。

unstranded RNA-seqの場合、Fraction of reads explained by “++,–"とFraction of reads explained by ”+-,-+“の割合が、50%に近くなりますが、今回のようにStrand-specificに読まれているサンプルの場合、下記のように偏りが見られます。

PAR-CLIPのデータでは、+鎖のRNA配列は+鎖として、-鎖のRNA配列は-鎖としてそれぞれ読まれていることがわかりました。

This is SingleEnd Data
Fraction of reads failed to determine: 0.0622
Fraction of reads explained by "++,--": 0.8506
Fraction of reads explained by "+-,-+": 0.0873

この情報をもとに、以下でUCSC genome browserで可視化させるデータを作成する際に、+と-鎖をそれぞれ分けてファイルを作成します。

11. Visualization (For UCSC genome browser)

Bedtoolsを利用して、STARから得られたBAMファイルを、UCSC genome browser上で可視化できるファイル（Bedgraphファイル）に変換します。

$ mkdir UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd

$ bedtools genomecov -ibam ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.bam -bg -split -strand + > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_Forward.bg
$ bedtools genomecov -ibam ./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.bam -bg -split -strand - > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_Reverse.bg

$ echo "track type=bedGraph name=SRR488740_PAR-CLIP_CstF-64_STAR_End_Fw description=SRR488740_PAR-CLIP_CstF-64_STAR_End_Fw visibility=2 maxHeightPixels=40:40:20 color=255,0,0" > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/tmp_SRR488740_PAR-CLIP_CstF-64_forward.txt
$ echo "track type=bedGraph name=SRR488740_PAR-CLIP_CstF-64_STAR_End_Re description=SRR488740_PAR-CLIP_CstF-64_STAR_End_Re visibility=2 maxHeightPixels=40:40:20 color=0,0,255" > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/tmp_SRR488740_PAR-CLIP_CstF-64_reverse.txt

$ cat ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/tmp_SRR488740_PAR-CLIP_CstF-64_forward.txt ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_Forward.bg > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_forward_for_UCSC.bg
$ cat ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/tmp_SRR488740_PAR-CLIP_CstF-64_reverse.txt ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_Reverse.bg > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_reverse_for_UCSC.bg

$ bzip2 -c ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_forward_for_UCSC.bg > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_forward_for_UCSC.bg.bz2
$ bzip2 -c ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_reverse_for_UCSC.bg > ./UCSC_visual_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out_reverse_for_UCSC.bg.bz2

bedtools genomecovで特定のゲノム領域でのCoverageを計算する。このとき、-bgオプションで、出力ファイルの形式をBedGraphファイルに指定しておく。また、-splitオプションで、Splicing juctionにマッピングされたリード（リードが分割されてマッピングされているリード）を考慮する。

+鎖と-鎖で別々のBedGraphファイルを作るために、-strandオプションを指定する。+鎖もしくは-鎖のみを抽出することができる。

先ほどのチェックで、+鎖のRNA配列は+鎖として、-鎖のRNA配列は-鎖としてそれぞれ読まれていることがわかっているので、ここでは-strand +と-strand -と指定して、それぞれを+鎖と-鎖のデータとして保存する。

UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要がある。

それらの情報をヘッダー行（ファイルの一行目）に記載する。 echoコマンドのところでヘッダー行の情報をtxtファイルとして出力している。

そのあとに、catコマンドを用いて、先程作ったBedGraphファイルと結合させる。

出力されたファイルをアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2ファイルに圧縮しています。

12. UCSC genome browserへのアップロード

下記のURLからUCSC genome browserにアクセスします。
http://genome.ucsc.edu/

[MyData] -> [Custom Tracks]をクリック。

image.png (1.0 MB)

Add Custom Tracksのページに飛ぶので、「ファイルを選択」をクリックし、アップロードしたいファイル（先ほど作成したbzip2に圧縮したBedGraph）を選択します。

Add Custom Tracks - Google Chrome 2017-01-18 17.11.59.png (91.5 kB)

「Submit」ボタンをクリックします。

Add Custom Tracks - Google Chrome 2017-01-18 17.12.47.png (91.3 kB)

「add custom tracks」をクリックすると、先ほどのページに飛び再び別のファイルをアップロードできます。また、「go」をクリックするとGenome browser上でアップロードしたデータを閲覧することができます。

Manage Custom Tracks - Google Chrome 2017-01-18 17.13.26.png (72.4 kB)

最初は、chr1の先頭に飛ばされるので、検索欄に興味のある遺伝子名・もしくはGencodeのIDを入力して検索することができます。

Human chr1_10611-10629 - UCSC Genome Browser v343 - Google Chrome 2017-01-18 17.16.46.png (312.5 kB)

実際に検索すると、マッピングされたCsfF64のPAR-CLIPのデータを見ることができ、3'endにピークが確認されます。

Human chr19_2475589-2478791 - UCSC Genome Browser v343 - Google Chrome 2017-01-18 17.17.50.png (629.1 kB)

13. PARalyzerの実行 (Peak calling)

必要なスクリプト

下記のURLに必要なスクリプト群が置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

  • PARalyzer_v1_5.sh
  • PARalyzer_ini_file_generator.py

内容

STARでマッピングした後、カレントディレクトリにSTAR_output_xxxという名前のディレクトリが生成されていると思います。その直下に、SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.bamという名前のBAMファイルがあると思います。

ここでは、そのBAMファイルから、PARalyzerと呼ばれるソフトウェアを用いてPeak callを行いたいと思います。

まず、PARalyzerではマッピングデータであるBAMファイルやGenome情報、各種パラメータの情報をまとめたiniファイルを用意する必要があります。

詳細については下記のURLを参照のこと。
https://ohlerlab.mdc-berlin.de/files/duke/PARalyzer/README.txt

iniファイルの準備

今回は、このiniファイルを生成するPythonスクリプトを用意したのでそれを用います。

$ python2 ./PARalyzer_ini_file_generator.py ./ SRR488740_PAR-CLIP_CstF-64 ./hg19.2bit ./PARalyzer_v1_5.ini
  • 1つ目の引数 BAMファイルの置いてあるディレクトリ。
  • 2つ目の引数: FASTQファイルのファイル名（.fastqを除いた部分）。
  • 3つ目の引数: 2bitのゲノム情報。
  • 4つ目の引数: 出力されるiniファイルの名前。

実際のiniファイルの中身。 自身でパラメータを設定したい場合は、下記の情報を書き換えたiniファイルを用意する必要があります。

BANDWIDTH=3
CONVERSION=T>C
MINIMUM_READ_COUNT_PER_GROUP=5
MINIMUM_READ_COUNT_PER_CLUSTER=5
MINIMUM_READ_COUNT_FOR_KDE=5
MINIMUM_CLUSTER_SIZE=8
MINIMUM_CONVERSION_LOCATIONS_FOR_CLUSTER=1
MINIMUM_CONVERSION_COUNT_FOR_CLUSTER=1
MINIMUM_READ_COUNT_FOR_CLUSTER_INCLUSION=5
MINIMUM_READ_LENGTH=13
MAXIMUM_NUMBER_OF_NON_CONVERSION_MISMATCHES=0
EXTEND_BY_READ
SAM_FILE=./STAR_output_SRR488740_PAR-CLIP_CstF-64_EndtoEnd/SRR488740_PAR-CLIP_CstF-64_4_STAR_result_Aligned.sortedByCoord.out.sam
GENOME_2BIT_FILE=./hg19.2bit
OUTPUT_DISTRIBUTIONS_FILE=SRR488740_PAR-CLIP_CstF-64_distribution.csv
OUTPUT_GROUPS_FILE=SRR488740_PAR-CLIP_CstF-64_groups.csv
OUTPUT_CLUSTERS_FILE=SRR488740_PAR-CLIP_CstF-64_clusters.csv
  • SAM_FILE: SAMファイルを指定。
  • GENOME_2BIT_FILE: 2bitのゲノム情報を指定。
  • OUTPUT_DISTRIBUTIONS_FILE: 出力ファイル。任意の名前を指定。
  • OUTPUT_GROUPS_FILE: 出力ファイル。任意の名前を指定。
  • OUTPUT_CLUSTERS_FILE: 出力ファイル。任意の名前を指定。

2bitのゲノム情報の取得

PARalyzerで2bitのゲノム情報が必要になるので、ダウンロードしておく。

$ wget http://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.2bit

PARalyzerの実行

PARalyzer v1.5から引数として「消費メモリ」と「iniファイル」を指定する必要がある。このとき、消費メモリはSGE (qsubによるジョブ)で指定・確保したメモリ量よりも小さい値を指定する必要がある。例えば、下記のコードでは、qsubによるジョブで「16G」のメモリを確保したケースで、PARalyzerの実行に「11G」のメモリを使用するように指定している。

$ PARalyzer 11G ./PARalyzer_v1_5.ini
  • 1つ目の引数: 先ほど作成したiniファイルを指定。

14. BED, FASTAファイルの作成

下記のURLに必要なスクリプト群が置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64 * PARalyzer2BEDandFASTA.py

Annotation付けやモチーフ検索のために、BEDとFASTAファイルを用意する。 入力ファイルは、PARalyzerから得られたXXX_clusters.csvを利用する。

$ python2 ./PARalyzer2BEDandFASTA.py SRR488740_PAR-CLIP_CstF-64_clusters.csv SRR488740_PAR-CLIP_CstF-64_clusters.bed SRR488740_PAR-CLIP_CstF-64_clusters.fa
  • 1つ目の引数: PARalyzerから出力されたCSVファイルを指定。
  • 2つ目の引数: 出力されるBEDファイルの名前を指定。
  • 3つ目の引数: 出力されるFASTAファイルの名前を指定。

15. MEMEによるモチーフ検索

PAR-CLIPのデータから得られた、CstF64タンパク質が結合しているRNA配列のFASTAファイルを用意します。

MEMEに投げる配列の数が多い場合、計算にかなりの時間を要します（一日以上）。そのため、数として10,000配列未満に抑えたほうがよいです。

PARalyzerから得られたピーク数が20万を超えるので、これを1万程度まで絞り込む必要があります。今回は、CstF64がRNAの3'endに結合していることがわかっているので、3'UTR上で検出されたピークのみを解析対象としたいと思います。

さらに、3'UTR上で複数のピークが確認された場合は、ピークのリード数が多いものを各遺伝子の代表のピークとして選択するようにします。

Representative isoformの用意

必要なスクリプト

下記のURLに必要なスクリプト群が置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

  • prep_gtf_for_PAR-CLIP_analysis.sh
  • A_prep_gtf_for_PAR-CLIP.py
  • gtf2bed.pl
  • B_prep_longest_trx_id_list.py
  • bed12to3UTRbed.pl

内容

各遺伝子に対してRepresentative isoformに対するBEDファイルとFASTAファイルをそれぞれ用意します。

まず、GTFファイルから代表的なIsoformを絞り込みます。Gencodeに登録されている各GeneのIsoformsはナンセンス変異を持つIsoformやTrancateなIsoformなどを含んでいます。

そこで、ここではbasicというタグの付いた完全長のIsoformだけを抽出しています。また、StatusがPUTATIVEとなっているIsoformについても除外しています。

# Prepare Protein-coding RNAs (Status: Basic)
$ python A_prep_gtf_for_PAR-CLIP.py gencode.v19.annotation.gtf \
gencode.v19.annotation_only_basic_mRNAs.gtf \
gencode.v19.annotation_only_basic_mRNAs_gene_trx_list.txt
  • 1つ目の引数: もとのGTFファイルを指定。
  • 2つ目の引数: 代表的なIsoformのみに絞ったGTFファイルを出力。
  • 3つ目の引数: 遺伝子リストを出力。

GTFファイルをBEDファイルに変換します。

# Convert gtf to bed
$ perl ./gtf2bed.pl gencode.v19.annotation_only_basic_mRNAs.gtf gencode.v19.annotation_only_basic_mRNAs.bed
  • 1つ目の引数: 先ほど作成したGTFファイルを指定。
  • 2つ目の引数: BED形式に変換したファイルを出力。

各Geneに対して１つの代表Isoformを決めるために、最長のIsoformのみを抽出します。

# Extract longest isoforms
$ python2 ./B_prep_longest_trx_id_list.py ./gencode.v19.annotation_only_basic_mRNAs.bed \
./gencode.v19.annotation_only_basic_mRNAs_gene_trx_list.txt \
./gencode.v19.annotation_only_Rep_basic_mRNAs.bed
  • 1つ目の引数: 先ほど作成したBEDファイルを指定。
  • 2つ目の引数: 先ほど作成した遺伝子リストを指定。
  • 3つ目の引数: 最終的に得られる代表IsoformのBEDファイルを出力。

先ほど作成したBEDファイルから、各遺伝子の3'UTRの領域についてのBEDファイルとFASTAファイルを用意します。

# Extract 3UTR bed file
$ perl ./bed12to3UTRbed.pl ./gencode.v19.annotation_only_Rep_basic_mRNAs.bed \
./gencode.v19.annotation_only_Rep_basic_mRNAs_3UTR.bed \
./gencode.v19.annotation_only_Rep_basic_mRNAs_non-3UTR.bed

# Get fasta file
$ bedtools getfasta -name -s -split -fi ./hg19.fa \
-bed ./gencode.v19.annotation_only_Rep_basic_mRNAs_3UTR.bed \
> ./gencode.v19.annotation_only_Rep_basic_mRNAs_3UTR.fa
bed12to3UTRbed.plのスクリプト。
  • 1つ目の引数: 先ほど作成した代表IsoformについてのBEDファイルを指定。
  • 2つ目の引数: 各遺伝子の3'UTRの領域についてのBEDファイルを出力。
  • 3つ目の引数: 3'UTRを持たないRNAのデータを出力（ncRNAはすべてこちらに出力される）。
bedtools getfastaコマンド。
  • -name: BEDの名前の列の情報を、FASTAファイルの各配列の名前として記述。
  • -s: Strandの向きを考慮する。
  • -bed: FASTAファイルに変換したいBEDファイルを指定。 標準出力で、FASTAファイルが出力される。

Representative peaksの選定

必要なスクリプト

下記のURLに必要なスクリプト群が置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

  • compare_annotation_3UTR_with_peaks.sh
  • D_prep_fasta_file_for_MEME_rep_version.py

PARalyzerにより検出されたピーク情報と3'UTRのアノテーション情報を比較。3'UTRの領域とOverlapするピークを抽出します（中でも、ピークのリード数が最も多い代表のピークを選抜する）。

# Compare 3'UTR annotation with PAR-CLIP peaks
$ bedtools intersect -a ./SRR488740_PAR-CLIP_CstF-64_clusters.bed -b ./gencode.v19.annotation_only_Rep_basic_mRNAs_3UTR.bed -wa -wb -s -split > ./SRR488740_PAR-CLIP_CstF-64_on_mRNAs_plus_anno.bed

# make a fasta file for MEME
$ python2 D_prep_fasta_file_for_MEME_rep_version.py ./SRR488740_PAR-CLIP_CstF-64_on_mRNAs_plus_anno.bed \
./SRR488740_PAR-CLIP_CstF-64_on_mRNAs_plus_anno_rep_version.fa
bedtools intersectコマンド
  • -a: 比較したいBEDファイルを指定。
  • -b: 比較したいBEDファイルを指定。
  • -wa: -aで指定したBEDファイルの情報を結果に出力。
  • -wb: -bで指定したBEDファイルの情報を結果に出力。
  • -s: Strandを考慮する。
  • -split: スプライシングの情報を含む12列のBEDファイル（BED12）に対応する。 標準出力で、2つのBEDファイルを比較したファイルが出力される。
D_prep_fasta_file_for_MEME_rep_version.pyスクリプト。
  • 1つ目の引数: 先ほど作成したBEDファイルを比較したファイルを指定。
  • 2つ目の引数: MEMEに用いるピークの配列情報（FASTAファイル）を出力。

Markov Backgroundモデルの作成

MEMEに用いるBackgroundモデルを作成します。

$ fasta-get-markov -m 3 ./gencode.v19.annotation_only_Rep_basic_mRNAs_3UTR.fa \ ./SRR488740_PAR-CLIP_CstF-64_markov_background_for_MEME.txt
  • 1つ目の引数: 先ほど作成した代表Isoformの3'UTR領域のFASTAファイルを指定。
  • 2つ目の引数: MEMEで使用するBackgroundファイルを出力。

詳しくは下記のサイトを参照のこと。
fasta-get-markov http://meme-suite.org/doc/fasta-get-markov.html

モチーフ検索

MEMEを用いて、モチーフ検索を行います。

$ meme ./SRR488740_PAR-CLIP_CstF-64_clusters_on_mRNAs_plus_anno_rep_version.fa \
-dna -bfile ./SRR488740_PAR-CLIP_CstF-64_markov_background_for_MEME.txt \
-oc MEME_result_SRR488740_PAR-CLIP_CstF-64_rep -minw 5 -maxw 8 -nmotifs 3 -maxsize 1000000000 -mod zoops
  • 1つ目の引数: 先ほど作成したピークの配列情報をまとめたFASTAファイルを指定。
  • -dna: DNA配列として処理（RNA、Proteinも選択できる）。
  • -oc: 出力先のディレクトリを指定。
  • -minw: 予測する最小モチーフ配列長を指定。
  • -maxw: 予測する最大モチーフ配列長を指定。
  • -maxsize: 最大文字数を指定（ファイルが大きいと最大文字数を超過してエラーが起こる。適当に大きな数字を入れておく）。
  • -mod zoops: １つの配列中に予測されたモチーフ配列が0もしくは1以上含まれる。

memeのオプションの詳細については、こちらを参照のこと。
http://meme-suite.org/doc/meme.html

モチーフ検索の結果

先行研究の結果と一致することを確認します。Cstf64タンパク質は、RNAのGU-richなモチーフに結合することが知られているので、モチーフ検索の結果、GU-richな配列が抽出されると期待されます。

以下がMEMEで予測されたモチーフ配列です。 Cstf64タンパク質が結合するRNA配列。

logo1.png (7.4 kB) p-value: 8.6e-1435

Cstf64はpolyA付加シグナル配列の近傍に結合することが知られているため、polyA付加シグナルのモチーフ配列も検出されました。

logo2.png (7.2 kB) p-value: 7.2e-828

今回解析したPAR-CLIPのデータ元の論文では、Top 500のピークについてモチーフ検索を行っているが、今回は約8,000のピークについて大まかに解析している。それでも、GU-richな配列が取れてきた。

参考

以下はRBPDBで検索した結果。

RBPDB - Google Chrome 2017-01-18 10.28.42.png (101.9 kB) （RBPDB [http://rbpdb.ccbr.utoronto.ca/]より抜粋）

Cstf64タンパク質の結合サイトに関する論文。

  • Cañadillas, J. M. P., & Varani, G. (2003). Recognition of GU‐rich polyadenylation regulatory elements by human CstF‐64 protein. The EMBO journal, 22(11), 2821-2830. [PubMed]

  • Martin, G., Gruber, A. R., Keller, W., & Zavolan, M. (2012). Genome-wide analysis of pre-mRNA 3′ end processing reveals a decisive role of human cleavage factor I in the regulation of 3′ UTR length. Cell reports, 1(6), 753-763. [PubMed]

（おまけ）16. アノテーション付け

内容

さきほど作成したRepresentative isoformの情報から、各遺伝子領域（5'UTR, CDS, 3'UTR, intron）上に見られるCstf64のPAR-CLIPのピークを確認します。

必要なスクリプト

下記のURLに必要なスクリプト群が置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/PAR-CLIP_Tutorial_CstF64

  • extract_gene_symbol_type_from_gtf.py
  • I_prep_mRNA_lncRNA_symbol_refid_list.py
  • F_split_into_each_section.py
  • G_add_intron_infor.py
  • H_annotate_peaks.py
  • prep_anno_list_for_user.sh
  • peak_anno.sh

遺伝子リストの準備

Gencode v19に収載されているmRNAとlncRNAの遺伝子をGTFファイルからリストアップする。

# Extract gene symbol & ref ID from GTF file
$ python2 ./extract_gene_symbol_type_from_gtf.py \
  ./gencode.v19.annotation.gtf \
  ./gencode.v19.annotation_symbol_type_list.txt

# Prepare mRNA & lncRNA data list
$ python2 ./I_prep_mRNA_lncRNA_symbol_refid_list.py \
  ./gencode.v19.annotation_symbol_type_list.txt \
  ./gencode.v19.annotation_symbol_type_mRNA_list.txt \
  ./gencode.v19.annotation_symbol_type_lncRNA_list.txt
extract_gene_symbol_type_from_gtf.pyスクリプト

GTFファイルから遺伝子リストを作成。 * 1つ目の引数: GTFファイルを指定。 * 2つ目の引数: GTFファイルから作成した遺伝子リストを出力。

I_prep_mRNA_lncRNA_symbol_refid_list.pyスクリプト

先ほど作成した遺伝子リストからmRNAとlncRNAをそれぞれ抽出。

  • 1つ目の引数: 先ほど作成した遺伝子リストを指定。
  • 2つ目の引数: mRNAのみの遺伝子リストを出力。
  • 3つ目の引数: lncRNAのみの遺伝子リストを出力。

PAR-CLIPのピークと遺伝子領域を比較

「15. MEMEによるモチーフ検索 / Representative isoformの用意」のところで作成した、各遺伝子の代表Isoformを抽出したgencode.v19.annotation_only_Rep_basic_mRNAs.bedのBEDファイルから、5' UTR、CDS、3' UTRごとに分割したBEDファイルを作成する。

# Get region infor for each gene (5UTR, ORF, 3UTR)
python2 ./F_split_into_each_section.py \
  ./gencode.v19.annotation_only_Rep_basic_mRNAs.bed \
  ./gencode.v19.annotation_only_Rep_basic_mRNAs_each_region.bed
  • 1つ目の引数: 各遺伝子の代表Isoformを抽出したBEDファイルを指定。
  • 2つ目の引数: 5' UTR、CDS、3' UTRごとに分割したBEDファイルを出力。

さらに、Intron領域についてもBEDファイルを作成する。

# Add intron region
$ python G_add_intron_infor.py  \
  ./gencode.v19.annotation_only_Rep_basic_mRNAs_each_region.bed \
  ./gencode.v19.annotation_only_Rep_basic_mRNAs_each_region_intron_plus.bed
  • 1つ目の引数: 先ほど作成した5' UTR、CDS、3' UTRごとに分割したBEDファイルを指定。
  • 2つ目の引数: 5' UTR、CDS、3' UTR、Intronごとに分割したBEDファイルを出力。

各遺伝子の5' UTR、CDS、3' UTR、Intronの領域とOverlapするPAR-CLIPのピークをbedtoolsを使って見つける。

# Compare mRNA region with CLIP peaks
$ bedtools intersect -a ./SRR488740_PAR-CLIP_CstF-64_clusters.bed.bed \
  -b ./gencode.v19.annotation_only_Rep_basic_mRNAs_each_region_intron_plus.bed \
  -wa -wb -loj > ./SRR488740_PAR-CLIP_CstF-64_clusters.bed_with_mRNA_region.tmp
  • -a: 先ほど作成した5' UTR、CDS、3' UTR、Intronごとに分割したBEDファイルを指定。
  • -b: 5' UTR、CDS、3' UTR、Intronごとに分割したBEDファイルを出力。
  • -wa: -aに指定したBEDファイルの情報を出力。
  • -wb: -bに指定したBEDファイルの情報を出力。
  • -loj: 通常、-aと-bで指定したBEDファイルを比較して、overlapした領域（各行）の情報しか出力しない。このオプションを付けると、-aに指定したBEDファイルの領域とoverlapしていなかった場合NULLを返すことで、-aで指定したBEDファイルは全行結果として出力されるようになる。

標準出力で比較結果が出力される。

bedtoolsの出力ファイルを整理して、各遺伝子についてPAR-CLIPのピークがどの遺伝子領域（5' UTR、CDS、3' UTR、Intron）に位置しているのかtab-delimited tableとしてまとめる（Excelにコピペして見れる状態にする）。

# Annotate peaks
$ python2 ./H_annotate_peaks.py SRR488740_PAR-CLIP_CstF-64_clusters.bed_with_mRNA_region.tmp \
  ./gencode.v19.annotation_symbol_type_mRNA_list.txt \
  ./SRR488740_PAR-CLIP_CstF-64_clusters_with_mRNA_region.txt
  • 1つ目の引数: 先ほど作成したbedtoolsの出力ファイルを指定。
  • 2つ目の引数: 遺伝子リストの準備の項で準備した、mRNAの遺伝子リストを指定。
  • 3つ目の引数: ピークの位置情報を加えた遺伝子リストを出力。

以上になります。