Dockerコンテナ上でNGSのソフトウェアを実行してみる(Windows10・MacOSXローカル環境編)

きっかけ

近年では、BiocondaのようなBiology関連ソフトウェアに対するパッケージマネージャが登場し、NGS関連のソフトウェアのインストールがグッと簡単になりました(Biocondaが登場する前は、いちいちソースコードからコンパイルする必要があったり、各ソフトウェア固有のインストール方法を理解する必要があったり、とにかく手間がかかっていました…。NGS関連のソフトは依存関係が複雑で本当にやっかい)。

ただ、いろんな人のNGS関連のブログを見ていると、環境によってはBiocondaでうまくインストールできないという事例もあるようです(Biocondaレシピで想定されていなかったパッケージへの依存によるもの?)。確かに、Biocondaは非常に便利なんですが、他の環境でも同じ解析環境を構築しようとすると、再現よく環境を構築できないケースがあることがわかってきました。

そこで、ソフトウェアごとに仮想環境を構築し、各ソフトウェアの動作が保証された環境を色々な解析に使いまわすことで、上記の問題を解決しようと考えました。

そこで、Dockerの出番です。Dockerはコンテナ型仮想化技術の1つであり、ホストOS上で単一のプロセスとして起動します。ホストOSと同じカーネルを共有するため、VMWare, VirtualBoxなどのホスト型仮想化技術と比較して、軽量であるとされています。

「Docker」を全く知らない人のために「Docker」の魅力を伝えるための「Docker」入門
https://qiita.com/bremen/items/4604f530fe25786240db

前置きが長くなりましたが、今回は、Dockerを土台としたシステム基盤をつくることを考えるために、まずはローカル環境でDockerコンテナを動かしてみた (FastQCを実行して見る)という内容になります。

Dockerをベースにした解析パイプラインの活用事例

パッと目についた事例をあげてみました。他にもたくさんあると思います。

High-Throughput Genomics on AWS - LFS309 - re:Invent 2017

AWSのカンファレンス re:Invent 2017で紹介されたNGS解析のパイプライン
https://www.slideshare.net/AmazonWebServices/highthroughput-genomics-on-aws-lfs309-reinvent-2017

Broad Institute: Genomes-in-the-cloud

Broad Instituteで利用されているゲノム解析用のパイプライン
https://hub.docker.com/r/broadinstitute/genomes-in-the-cloud/

BioContainers

バイオインフォマティクス関連のDockerイメージファイルをまとめたもの
https://biocontainers.pro/

ベースとなるDockerイメージの選択

NGSのデータ解析をする上では、必要最小限のLinuxディストリビューションのDockerイメージを使用した方がいいと思います。1つはイメージサイズが小さいほど、Dockerイメージの起動や破棄、ダウンロードなどの取り回しの面で有利だから。もう1つは、軽量コンテナの方が不要なパッケージが少なく、セキュリティ上、潜在的な攻撃を受けにくいからです。

また、ベースイメージを選定する上で、OSにインストールマネージャーがついているかどうかや、ITベンダーなどにより定期的にメンテナンスされているかどうかもポイントになります。

以下では、主要なLinuxディストリビューションについて簡単に比較をして見たいと思います。

Dockerイメージのサイズ

以下が主要なLinuxディストリビューションのDockerイメージのサイズになります。

Linux OS Compressed size
debian8.10-slim 30 MB
debian 8.10 53 MB
ubuntu 18.04 35 MB
centos 7.4.1708 73 MB

上記のデータは、2018-03-13時点のデータです。少し前までは、Debianが軽量なDockerイメージだったので、Debian一択なのかと思っていましたが、Ubuntuも直近ではかなりスリム化を果たしていますね(ちなみに、Anacondaやその他NGS関連のDockerイメージはDebianをベースに作られたものが多い気がします)。

脆弱性レポート

各Dockerイメージのレポジトリのtagタブから、各バージョンの脆弱性レポートを確認することができます。

image.png (328.0 kB)

スキャンされた各イメージをクリックすると、詳細を見ることができます。以下では、各LinuxディストリビューションのDockerイメージの脆弱性レポートを見ていくことにします。

debian8.10-slim

image.png (281.9 kB)

debian 8.10

image.png (283.9 kB)

ubuntu

image.png (281.8 kB)

centos 7.4.1708

image.png (308.9 kB)

脆弱性レポートまとめ

上記で示した脆弱性レポートの結果をまとめて見ました。いくつのコンポーネント脆弱性を抱えているかカウントした表が以下の通りになります。

Linux OS Vulnerability (Components)
debian8.10-slim 6
debian 8.10 6
ubuntu 18.04 0
centos 7.4.1708 16

上記のデータは、2018-03-13時点のデータです。Ubuntu脆弱性が0なのが際立っている気がします(本当かなって結果になりました)。

現状では、Debian-slimかUbuntuを選択するのが良いのかなと思いました。

Dockerのインストール

自身のパソコンにDockerをインストールします。 以下のURLから、Docker (Community Edition)をダウンロードできます。
https://www.docker.com/community-edition#/download

MacOSXにDockerをインストールして使う場合は以上でインストール完了です。
Windows10環境の場合、既存のコマンドプロンプトPowershellがどうも使いにくい(個人的な感想ですが)ので、Windows Subsystem for Linux (WSL)を導入して、BashシェルからDockerを扱えるようにしたいと思います。

Windows Subsystem for Linux (WSL)経由でホストOS (Windows10)のDockerを動かす

1. WSLをインストールする

注意: Dockerを使用するにはHyper-Vを利用する必要があるため、OSがWindows10 Proであることが必須です。

Miscrosoft storeからWSLをダウンロード・インストールできるようになったので、ストアからインストールしましょう。今回はUbuntuを導入しました。

image.png (22.4 kB)

WSLを使用するためには、「プログラムと機能」から「Windowsの機能の有効化または無効化」を開き、Windows Subsystem for Linuxにチェックをつけます。

image.png (73.0 kB)

2. DockerをWSLにインストールする

以下のUbuntuへのDockerのインストール方法は、Dockerの公式ドキュメントを参考にしました。詳しくは、そちらを参照のこと。

Get Docker CE for Ubuntu
https://docs.docker.com/install/linux/docker-ce/ubuntu/#install-using-the-repository

WSLのBash shellを起動し、まずaptパッケージをアップグレードしておきます。

$ sudo apt-get update

次のDockerを動かすために必要な各種パッケージをインストールします。

$ sudo apt-get install \
    apt-transport-https \
    ca-certificates \
    curl \
    software-properties-common

Dockerの公式のGPGキーを加え、キー情報の確認を行っています。

$ curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo apt-key add -
$ sudo apt-key fingerprint 0EBFCD88

pub   4096R/0EBFCD88 2017-02-22
      Key fingerprint = 9DC8 5822 9FC7 DD38 854A  E2D8 8D81 803C 0EBF CD88
uid                  Docker Release (CE deb) <docker@docker.com>
sub   4096R/F273FCD8 2017-02-22

Dockerをインストールします。

$ sudo add-apt-repository \
   "deb [arch=amd64] https://download.docker.com/linux/ubuntu \
   $(lsb_release -cs) \
   stable"
$ sudo apt-get update
$ sudo apt-get install docker-ce

3. WSL側とホストOS (Windows10)側のDockerを連携させる

Docker for Windowsのデーモンに接続するための設定を行います。今の状態では、あくまでホストOS側(Windows10)とWSL側で別々にDockerをインストールしただけなので、WSL側から、ホストOS側のDockerのデーモンプロセスを操作できるのように連携させる必要があります。

以下のように、WSL側のDockerクライアントからホストOS側のデーモンを見れるように、環境変数DOCKER_HOSTを設定します。

$ echo "export DOCKER_HOST='tcp://0.0.0.0:2375'" >> ~/.bashrc
$ source ~/.bashrc

最後に、Windows10側のDockerの設定画面を開き、Expose daemon on tcp://localhost:2375 without TLSにチェックを付けます。これで、WSL側からホストOS側のDockerのデーモンを操作できるようになりました。

image.png (142.5 kB)

参照

Docker for Windowsで快適な環境を得るまでの そこそこ長い闘い
https://qiita.com/YukiMiyatake/items/73c7d6c4f2c9739ebe60
WSL(Bash on Windows)でDockerを使用する
https://qiita.com/yoichiwo7/items/0b2aaa3a8c26ce8e87fe
Docker for WindowsをWSLから使う時のVolumeの扱い方
https://qiita.com/gentaro/items/7dec88e663f59b472de6

VSCodeの総合ターミナルでWSLを使う

VSCodeでは、エディタ内でターミナルを開くことができるので、こちらの標準のターミナルをWSLに変更しておくと便利です。

まず。設定画面から、terminal.integrated.shell.windowsで検索します。

対象の設定に対して "terminal.integrated.shell.windows": "C:\\Windows\\System32\\bash.exe"のようにWSLのBash.exeのパスを指定します。

image.png (111.0 kB)

総合ターミナルを開いて確認すると、めでたくBashシェルが使えるようになっているはずです!

image.png (30.8 kB)

Dockerファイルの作成

まず、Linux OSのDockerイメージをベースイメージにして、その上にNGS解析のソフトウェアをインストールしたDockerイメージを作成します。

1. Debian9-slimにBiocondaをインストールしたDockerイメージを作成する

Linux OSのDebina9-slimをベースイメージにして、BiocondaをインストールしたDockerイメージを作成します。

以下のDockerfileは、condaの公式Dockerfileを参考しました。
https://hub.docker.com/r/conda/miniconda2/~/dockerfile/

&&で処理を連続して行っているのは、レイヤーを減らすためです。また、中間ファイルを削除することでDockerのイメージファイルのサイズを小さくしています。

Dockerfile

# Base image
FROM debian:9-slim

# Install Miniconda2 & Bioconda
RUN apt-get -qq update && apt-get -qq -y install curl bzip2 \
  && curl -sSL https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -o /tmp/miniconda.sh \
  && bash /tmp/miniconda.sh -bfp /usr/local \
  && rm -rf /tmp/miniconda.sh \
  && conda install -y python=2 \
  && conda update conda \
  && conda config --add channels r \
  && conda config --add channels defaults \
  && conda config --add channels conda-forge \
  && conda config --add channels bioconda \
  && apt-get -qq -y remove curl bzip2 \
  && apt-get -qq -y autoremove \
  && apt-get autoclean \
  && rm -rf /var/lib/apt/lists/* /var/log/dpkg.log \
  && conda clean --all --yes

ENV PATH /opt/conda/bin:$PATH

次に、DockerfileからDockerイメージを作成します。

$ docker build -t debian9-slim-bioconda2:latest -f Dockerfile .

2. FastQCをインストールしたDockerイメージを作成する

先程作ったDockerイメージをベースにして、FastQCのインストールしたイメージファイルを作成します。

Dockerfile

# Base image
FROM debian9-slim-bioconda2:latest

# Application entry point
# Dejavu fonts are not included with default::openjdk(fastqc recipe).
# They are already included with conda-forge::openjdk.
# https://github.com/conda-forge/staged-recipes/issues/3164 
RUN conda update conda \
  && conda install fastqc=0.11.6 \
  && conda install -c conda-forge openjdk \
  && conda clean --all --yes
COPY ./fastqc/run_fastqc_local.py /run_fastqc_local.py

ENTRYPOINT ["python", "/run_fastqc_local.py"]

Pythonのコードをエントリポイントに設定します。以下が、そのソースコードになります。

from __future__ import print_function

import os
import shlex
import subprocess
from argparse import ArgumentParser


def run_fastqc(fastq_path, cmd_args, working_dir):
    """
    Runs fastqc  
   :param fastq_path: Local path to fastq file  
   :param cmd_args: Additional command-line arguments to pass in  
   :param working_dir: Working directory  
   :return: Path to the result folder
    """
    # Prepare fastqc result path
    fastqc_result_dir = os.path.join(working_dir, '')

    # Prepare fastqc log path
    fastqc_log_path = os.path.join(fastqc_result_dir, "log_fastqc.txt")

    # Prepare fastqc command
    cmd = "fastqc -o {0} {1} {2}".format(
        fastqc_result_dir, cmd_args, fastq_path)
    print(cmd)

    # Execute fastqc
    with open(fastqc_log_path, 'w') as log_file:
        subprocess.check_call(shlex.split(cmd), stdout=log_file)

    return fastqc_result_dir


def main():
    argparser = ArgumentParser()

    argparser.add_argument('--fastq_path', type=str,
                           help="fastq sequence file", required=True)
    argparser.add_argument('--cmd_args', type=str,
                           help="arguments/options for fastqc", default="")
    argparser.add_argument('--working_dir', type=str, default="data/")

    args = argparser.parse_args()

    print("Creating working directory...")
    if not os.path.exists(args.working_dir):
        os.mkdir(args.working_dir)

    print("Running fastqc...")
    local_result_dir = run_fastqc(
        args.fastq_path, args.cmd_args, args.working_dir)

    print('Save the result in {0}.'.format(local_result_dir))

    print("successfully Completed !!")


if __name__ == '__main__':
    main()

次に、DockerfileからDockerイメージを作成します。Pythonのコードも同じディレクトリに置いておきます。

$ docker build -t debian9-slim-bioconda2-fastqc:latest -f Dockerfile .

以下のコマンドで作成したDockerイメージをチェックできます。

$ docker images

image.png (53.4 kB)

Dockerイメージの実行

1. FASTQファイルの準備

適当なFASTQファイルを持っていない場合、以下の記事を参考にしてテストに使うFASTQファイルを用意してください。

imamachi-n.hatenablog.com

2. Dockerコンテナを起動する

FASTQ_FILEにFASTQファイル、FASTQファイルが置いてあるLOCAL_WORKING_DIRに作業ディレクトリを指定します。また、DOCKER_IMAGEに先ほど作成したDockerイメージを、DOCKER_WORKING_DIRにDockerコンテナ内の適当な作業ディレクトリを指定します。

#!/bin/bash
#LOCAL_WORKING_DIR="/Users/imamachinaoto/Desktop/NGS/data" # MacOS
LOCAL_WORKING_DIR="/c/Users/imama/Desktop/NGS/fastqc/data"    # Windows10(WSL)
DOCKER_WORKING_DIR="/data"
DOCKER_IMAGE="debian9-slim-bioconda2-fastqc:latest"
FASTQ_FILE="RefSeq_hg38_2015-08-19_NM_ultraSlim.fq"

docker run --rm -it \
-v ${LOCAL_WORKING_DIR}:${DOCKER_WORKING_DIR} \
${DOCKER_IMAGE} \
--fastq_path ${DOCKER_WORKING_DIR}/${FASTQ_FILE} \
--working_dir ${DOCKER_WORKING_DIR}

ちなみに、DockerでWSL環境からホストOS(Windows10)のフォルダを見るとき、Cドライブは/mnt/cではなく、/cです。Bashシェルからアクセスするときは、/mnt/c(WSLにマウントされているCドライブを見る形)なのでパスの違いに注意が必要です。Windowsだと、これで詰まりました…。

コンテナでデータを管理する
http://docs.docker.jp/engine/userguide/dockervolumes.html

上記シェルスクリプト(もしくは以下で説明するmakeコマンド)を実行することで、FastQCが実行されます。

image.png (191.3 kB)

ホストOSのディレクトリを、Docker上のディレクトリに紐づけておくこと(ボリュームのマウント)で、以下のようにDockerコンテナが破棄されたあとでも、データを永続化させることができます。

image.png (36.4 kB)

Dockerコマンドをまとめたmakefileの作成

Dockerコマンドはイメージファイルの作成やコンテナ起動時など、繰り返し使用します。なので、今回はmakefileにコマンド類をまとめて、make buildmake startなどを用意することで、コマンド操作が簡単になります。

Windows10やWSLには、makeコマンドがインストールされていないため、あらかじめインストールしておきます。

Windows10でmakeファイルを実行する

WSLを使わない場合のmakeコマンドのインストール方法を示します。Windows10にはコマンドプロンプトからmakeコマンドを使えないので、まずmakeコマンドをインストールします。

1. makeコマンドをコマンドプロンプトに導入する(参考)

以下のURLから、makeをダウンロードします。
http://gnuwin32.sourceforge.net/packages/make.htm

image.png (225.1 kB)

C:\Program Files (x86)\GnuWin32\binにmake.exeが保存されているので、パスを通します。

image.png (56.5 kB)

2. WSL環境にmakeコマンドを導入する

WSLのコンソールから、以下のコマンドを実行するだけです。

$ sudo apt-get install make

Windowsでmakeコマンドを使う
https://qiita.com/tokikaze0604/items/e13c04192762f8d4ec85

makefileの例

NAME=debian9-slim-bioconda2-fastqc
TAG=0.11.6.local
UNDERSCORE_TAG=0_11_6_local

LOCAL_WORKING_DIR=/c/dev/NGS/fastqc/data
DOCKER_WORKING_DIR=/data
DOCKER_IMAGE="fastqc-aws-debian9-slim-bioconda2:0.11.6.local"
FASTQ_FILE=RefSeq_hg38_2015-08-19_NM_ultraSlim.fq

all: build push

build:
      docker build -t $(NAME):$(TAG) -t $(NAME):latest -f Dockerfile .
build-no-cache:
      docker build --no-cache -t $(NAME):$(TAG) -t $(NAME):latest -f Dockerfile .

push:
      docker push $(NAME):$(TAG)
      docker push $(NAME):latest

start:
      docker run --rm -it -v ${LOCAL_WORKING_DIR}:${DOCKER_WORKING_DIR} ${DOCKER_IMAGE} --fastq_path ${DOCKER_WORKING_DIR}/${FASTQ_FILE} --working_dir ${DOCKER_WORKING_DIR}

rm:
      docker rm `docker ps -aq`

rmi:
      docker rmi $(NAME):$(TAG) $(NAME):latest

つまづいたポイント

makefileの作成

インデントはスペースではなく、タブにしないとエラーになる。 Visual Studio Codeの場合、makefileと認識されたファイルでは、インデントをタブにするように設定されている。しかし、ファイル名がMakefileだと下記の設定が適応されない。makefileと小文字のmにする必要がある。

image.png (76.5 kB)

Docker関連の参考

docker コマンド チートシート
https://qiita.com/voluntas/items/68c1fd04dd3d507d4083
dockerでなrepositoryを消す
https://qiita.com/lirispp/items/06fc74c5bbc64fddf9ab

ART: シミュレーション用のRNA-seqデータを作成する

テスト用のFASTQファイルの作成

ART (a next-generation sequencing read simulator)を用いて、テスト用のFASTQファイルを作成します。

ARTのインストール

※ Biocondaをインストールしていることを前提にしています。
ARTをインストールします。

$ conda intall art

ARTが正しくインストールされたか確認します。

$ art_illumina --help

RefseqのGTFファイルをダウンロード

ARTを使ってテスト用のFASTQファイルを作るのですが、今回はRNA-seqのデータを作ろうと思います。

そこで、元となる配列を用意したいと思います。 TranscriptomeのデータをGTFファイルと取得し、それをFASTQのファイルに変換したいと思います。

GencodeのGTFファイルでもいいですが、ここでは必要なTranscriptomeデータが含まれておりコンパクトなRefseqのデータを使います。

illuminaのiGenomeサイトから、hg38のデータを丸ごとダウンロードします。この中に、RefSeqのGTFファイルも含まれています。

image.png (182.7 kB)

ダウンロードしたら、解凍します。

$ tar zxvf Home_sapiens_UCSC_hg38.tar.gz

解凍が終わったら、Homo_sapiensというディレクトリができるはずです。 Homo_sapiens/UCSC/hg38/Annotation/Archives/archive-2015-08-14-08-18-15/Genes/にあるgenes.gtfというファイルがRefSeqのGTFファイルになります(ちょっと古いデータですが、実際の解析に用いるわけではないので、良しとしましょう)。

これだけ適当なディレクトリに移動もしくはコピーしておき、それ以外のファイルは削除します。

このままだと、何のファイルかわからないのでリネームしておきます。

$ mv genes.gtf RefSeq_hg38_2015-08-19.gtf

GTFファイルからBEDファイルを作成

以下のスクリプトを使って、GTFファイルをBEDファイルに変換する。 gtf2bed.pl (3.8 kB)

$ perl gtf2bed.pl RefSeq_hg38_2015-08-19.gtf > RefSeq_hg38_2015-08-19.bed

このままだと、データが大きすぎるのでmRNAに限定し、かつRepresentative transcriptsをできるだけ抽出するようにします。
ExtractMrna.py (360 B)

$ python ExtractMrna.py RefSeq_hg38_2015-08-19.bed RefSeq_hg38_2015-08-19_NM_slim.bed

ヒトゲノム配列をダウンロード

image.png (1.5 MB)

image.png (1.5 MB)

image.png (300.8 kB)

image.png (71.1 kB)

image.png (273.5 kB)

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
$ gunzip hg38.fa.gz

BEDファイルからFASTAファイルを作成

BEDファイルからFASTA ファイルへ変換するために、Bedtoolsを使います。 まず、Bioconda経由でBedtoolsをインストールします。

$ conda install bedtools

BEDファイルをFASTAファイルに変換します。 ここで、先ほどダウンロードしたhg38.faファイルを使います。

$ bedtools getfasta -name -s -split -fi hg38.fa -bed RefSeq_hg38_2015-08-19_NM_slim.bed > RefSeq_hg38_2015-08-19_NM_slim.fa

-name: BEDファイルのname列の値を各配列の名称とする。
-s: ゲノム上での向き (Strandness)を考慮する(Transcriptomeなのでもちろん考慮します)。
-split: 12BED (12列のデータがあるBEDふぁいる)がInputである場合に指定 。
-fi: Referenceになる配列(ゲノム配列)のFASTAファイルを指定。
-bed: TranscriptomeのBEDファイルを指定。

このコマンドの実行には、FASTAファイルのインデックスファイルが必要になりますが、ファイルがない場合、自動的に生成してくれるので問題ないです。

> index file hg38.fa.fai not found, generating...

生成されたTranscriptomeのFATSAファイルは、標準出力で出力されます。

TranscriptomeのFASTAファイルからFASTQファイルを生成

今回は、illuminaのFASTQシーケンスデータ (RNA-seq, single-end)を生成したいと思います。ここでは、HiSeq2000のデータを作ります。

$ art_illumina -ss HS20 -i RefSeq_hg38_2015-08-19_NM_slim.fa -o RefSeq_hg38_2015-08-19_NM_slim -l 36 -f 2 --noALN
$ art_illumina -ss HS20 -i RefSeq_hg38_2015-08-19_NM_slim.fa -o RefSeq_hg38_2015-08-19_NM_slim_paired -l 36 -f 1 -p -m 300 -s 10 --noALN

-i --in: fasta形式の入力ファイルの指定。
-o --out: fastq形式の出力ファイルの指定 (拡張子はつけなくて良い)。
-ss --seqSys: シーケンサーの種類を選択。
指定する引数(System ID) - シーケンサー名 (リード長)

  • GA1 - GenomeAnalyzer I (36bp,44bp),
  • GA2 - GenomeAnalyzer II (50bp, 75bp)
  • HS10 - HiSeq 1000 (100bp)
  • HS20 - HiSeq 2000 (100bp)
  • HS25 - HiSeq 2500 (125bp, 150bp)
  • HS10 - HiSeq 1000 (100bp)
  • HS20 - HiSeq 2000 (100bp)
  • HS25 - HiSeq 2500 (125bp, 150bp)
  • HSXn - HiSeqX PCR free (150bp)
  • HSXt - HiSeqX TruSeq (150bp)
  • MinS - MiniSeq TruSeq (50bp)
  • MSv1 - MiSeq v1 (250bp)
  • MSv3 - MiSeq v3 (250bp)
  • NS50 - NextSeq500 v2 (75bp)

-l --len: リード長の指定。
-f --fcov: リードカバレージの指定。
-p --paired: Paired-endのデータとして出力(指定しない場合、Single-endのデータとして出力)。
m --mflen: フラグメントの平均サイズの指定(Paired-endのときのみ)。
-s --sdev: フラグメントのサイズの標準偏差の指定(Paired-endのときのみ)。
-sam --samout: SAMファイルを出力。
-na --noALN: アライメントファイル (ALN)を出力しない。
別にリード長は、指定された長さよりも短くても大丈夫です。

もっとファイルサイズを小さくしたい場合は、以下のように特定の行を抽出したファイルを作成します。

$ sed -n 1,200000p RefSeq_hg38_2015-08-19_NM_slim.fq > RefSeq_hg38_2015-08-19_NM_ultraSlim.fq

Seleniumコトハジメ - PubMed検索とUCSC Genome browserへのCustom Trackのアップロードを自動化してみた

はじめに

Seleniumとは、Webアプリケーションの画面操作を自動化するためのツールです。主に、Webアプリケーション開発のUIテスト時に利用するものです。また、キー入力などの操作を自動化できることから、Webスクレイピングのツールとして活用される例もあるようです。

Seleniumは、Apatch 2.0ライセンスで公開されているオープンソースのツールです。
https://www.seleniumhq.org/about/license.jsp

Seleniumの内部で使われているWebDriverは、Webブラウザを外部から遠隔操作するためのインターフェースであり、W3Cによって標準化されています(厳密には、標準化された内容と乖離する部分あり?)。
https://www.w3.org/TR/webdriver/

Selenium2では、Selenium RCと呼ばれる古い実装が残っていましたが、2016年に公開されたSelenium3では、完全にWebDriverの実装に移行しました。古い記事では、Selenium RCでの操作方法が記載されており、これは現在では動作しない方法になるので注意が必要です。

本来はUIテスト用のツールですが、今回は研究室時代にめんどくさいと思っていた作業を自動化してみました。1つはPubMedによるAdvanced Search、もう1つはUCSC Genome BrowserにCustom trackをアップロードする作業です。繰り返し同じことを繰り返す作業に対して、こういった自動化の方法は役に立ちますね。

Selenium実行サンプル

PubMed検索の自動化

selenium-test.gif (753.1 kB)

UCSC Genome BrowserへのCustom Trackのアップロード

selenium-test2.gif (834.3 kB)

実行環境・ツール

Mac OSX Sierra (v10.12.6)

  • Selenium Client for Java (v3.10.0)
  • ChromeDriver (v2.35)
  • GeckoDriver (v0.19.0)
  • SafariDriver + Safari (v11.0.3)

Windows10 (Build 16299)

Seleniumの導入

Seleniumを動かすのに必要なものは以下の通りです。

  • Selenium Client & WebDriver Language Bindings
  • 各種ブラウザのWebDriver

今回はJavaでの例を示したいと思いますが、javaであれば、MavenやGradleを使って導入する方法もあります。ただ、今回はJARファイルをダウンロードして1からプロジェクトを作っていきます。

Selenium Client & WebDriver Language Bindings

Seleniumは、JavaC#RubyPythonJavascript (Node)を公式でサポートしており、それぞれの言語についてSeleniumAPIが使えるライブラリ(Selenium Client)を用意しています。

まず、以下のサイトからJavaSelenium Clientのファイルをダウンロードします。 selenium-java-3.xx.xx.zipというファイルをダウンロードできるはずです。
https://www.seleniumhq.org/download/

image.png (195.3 kB)

各種ブラウザのWebDriver

続いて、各種ブラウザのWebDriverをダウンロードします。私の使っているPCはMac OSXなので、とりあえずFireFoxGoogle ChromeSafariのWebDriverをそれぞれダウンロードすることにします(自身のパソコンにそれぞれのブラウザがすでにインストール済みでであることを前提とします)。

image.png (310.2 kB)

余談ですが、WebDriverはSeleniumを開発しているグループが作成しているわけではなく、Webブラウザの開発元であるGoogleMozillaAppleMicrosoftなどがそれぞれ開発をしています。そのため、基本的に最新バージョンのWebブラウザに対応するWebDriverを、これらのITベンダーが提供してくれる状態にあります。こういった点も、Seleniumを使うメリットですね。

Mozilla GeckoDriver (FireFox) (Selenium3.5以降でのみ動作)

以下からダウンロードします。
https://github.com/mozilla/geckodriver/releases

自身のパソコンのOSに合ったWebDriverをダウンロードしてきます。

image.png (251.7 kB)

以下のコマンドでファイルを解凍します。

$ tar zxvf geckodriver-v0.19.1-macos.tar.gz

geckodriverという名前のWebDriverの実行ファイルが解凍されます。

Google Chrome Driver (Google Chrome)

以下からダウンロードします。
https://sites.google.com/a/chromium.org/chromedriver/

image.png (410.3 kB)

WebDriverをインストールする際は、Google Chromeのバージョンに気をつけてください。最新のものだと、対応しているバージョンはv62-64に限定されます(もちろん、以前のバージョンのWebDriverについてもダウンロード可能です)。自分のパソコンにインストールされているGoogle Chromeのバージョンを確認し、古いバージョンであれば、そのバージョンをサポートしているWebDriverをダウンロードしてください。

image.png (334.8 kB)

自身のパソコンのOSに合ったWebDriverをダウンロードしてきます。 image.png (175.0 kB)

以下のコマンドでファイルを解凍します。

$ unzip chromedriver_mac64.zip

chromedriverという名前のWebDriverの実行ファイルが解凍されます。

SafariDriver (Safari)

Mac OSX内にすでにインストール済みです。
/usr/bin/safaridriverの場所に存在します。

Safariの設定を変更する必要があります。 まず、Safariの環境設定を開きます。 image.png (67.6 kB)

次に、「メニューバーに'開発'メニューを表示」にチェックをつけます。 image.png (362.8 kB)

すると、メニューバーに「開発」が追加されるので、そこから「リモートオートメーションを許可」をクリックします。これで、ブラウザ側の設定は完了です。 image.png (1.1 MB)

リモートオートメーションを許可していない場合、Javaのテストコードを実行すると、以下のようなエラーが出力されます。このエラーが出た場合は、上記の設定ができてないので、設定を変更してください。

org.openqa.selenium.SessionNotCreatedException: Could not create a session: You must enable the 'Allow Remote Automation' option in Safari's Develop menu to control Safari via WebDriver. (WARNING: The server did not provide any stacktrace information)

Microsoft Edge

以下からダウンロードします。
https://developer.microsoft.com/en-us/microsoft-edge/tools/webdriver/

image.png (389.2 kB)

自身のパソコンのOSのビルドに合ったWebDriverをダウンロードしてきます。 image.png (120.6 kB)

Internet Explorer 11

設定が複雑なので、必要でない限り使わないほうが良いと思います。

こちらに詳しい使用方法が記載されているので参照してください。 保護モードの設定が、「インターネット」「イントラネット」「信頼済みサイト」「制限付きサイト」の4つで統一されていることがかなり重要です(この設定が揃っていないとWebDriverが動作しません)。

seleniumでinternet-explorer11を動かす方法

http://bitwave.showcase-tv.com/seleniumでinternet-explorer11を動かす方法/

トラブルシューティング(一部)

画面は拡大せずに、100%でないと動作しません。以下のエラーが出た場合は、拡大率を100%にしましょう。

org.openqa.selenium.SessionNotCreatedException: Unexpected error launching Internet Explorer. Browser zoom level was set to 110%. It should be set to 100%

image.png (26.9 kB)

Seleniumを動かしてみる

以下では、Eclipseがインストール済みであることを前提として説明を行います。 素のEclipseをインストールするのではなく、以下のAll-in-Oneのパッケージ(Pleiades All in One)をインストールすることをお勧めします。
http://mergedoc.osdn.jp/

1. プロジェクトファイルの作成

Seleniumを動かすためのプロジェクトファイルを作成します。 [新規] -> [Javaプロジェクト]をクリックします。 image.png (70.3 kB)

適当なプロジェクト名をつけて、「次へ」をクリックします。 image.png (51.9 kB)

先ほどダウンロードしたSelenium Clientを解凍し、libsフォルダのjarとclient-combined-3.10.0.jarclient-combined-3.10.0-sources.jarをそれぞれプロジェクトファイルに取り込みます。

image.png (86.4 kB)

まず、「ライブラリ」タブを選択し、「外部JARの追加」をクリックします。 image.png (41.3 kB)

上述のSelenium ClientのJARファイル一式を追加します。 最後に、「完了」ボタンをクリックし、プロジェクトを生成します。 image.png (75.1 kB)

すると、参照ライブラリに先ほど取り込んだSelenium ClientのJARファイルが参照されて取り込まれていることが確認できると思います。

image.png (51.5 kB)

次に、WebDriverを保存するためのフォルダを用意します。 [新規] -> [フォルダー]をクリックします。

image.png (117.2 kB)

フォルダー名を「exe」として、「完了」ボタンをクリックします。 image.png (25.0 kB)

作成した「exe」フォルダ内に、先ほどダウンロードしたWebDriverをドラック&ドロップでインポートします。(ここでは、chromedriver.exeをインポートした例を示した。)

image.png (13.3 kB)

Mac OSXの場合、実行権限を付与しておきます。 image.png (58.4 kB)

$ chmod 755 geckodriver

続いて、Seleniumを実行するためのテストクラスを作成していきます。 プロジェクトの「src」フォルダを右クリックし、[新規] -> [Junitテスト・ケース]をクリックします。

image.png (80.9 kB)

適当なパッケージ名とクラス名を記入し、setUp()tearDown()にチェックを付けます。

image.png (47.3 kB)

JUnit4がビルドパスに入っていない場合、追加するかどうか聞かれるので、「OK」をクリックしてビルドパスに加えます。

image.png (20.0 kB)

2. PubMedUCSC Genome Browserを操作するコードを書いてみる

ページオブジェクトパターン

下図にページオブジェクトパターンと呼ばれる、一般的なクラスの設計方法を示しました。この方法を取り入れることで、メニューやボタン、検索入力欄など、Webページへの操作を共通クラスにまとめることでメンテナンス性の高いコードを書くことができます。

image.png (538.4 kB)

ロケータの調べ方

selenium用のコードを書くために、まずロケータ(ボタンやメニューなどを操作するための目印)を調べる必要があります。例えば、FireFoxを使うと、要素(下の例だと"Advanced")を右クリックして「要素を調査」をクリックすることでロケータを調べることができます。

image.png (353.4 kB)

どのタグが使われているかなどの情報を調べることができます(以下の図だと、INPUTタグのID要素が"Submit"であることがわかります。これを目印にしてSeleniumでWebブラウザを操作していきます)。

image.png (682.9 kB)

コーディング

以下のような構成で、プログラムを用意しました。 pageパッケージにページオブジェクトクラスを、testパッケージにJUnitのテストケースクラスを配置するようにしました。

image.png (36.0 kB)

WebDriverの呼び出し

System.setProperty()の2番目の引数は、ダウンロードしたWebDriverを指定してください。OSによって、拡張子が異なるので注意する必要があります。以下でいくつか例を示しました。

Google Chrome (Mac OSX)
System.setProperty("webdriver.chrome.driver", "exe/chromedriver");
WebDriver driver = new ChromeDriver();
Firefox (Mac OSX)
System.setProperty("webdriver.gecko.driver", "exe/geckodriver");
WebDriver driver = new FirefoxDriver();
Safari (Mac OSX)
WebDriver driver = new SafariDriver();
Microsoft Edge (Windows10)
System.setProperty("webdriver.edge.driver", "exe/MicrosoftWebDriver.exe");
WebDriver driver = new EdgeDriver();
Internet Explorer (Windows10)
System.setProperty("webdriver.ie.driver", "exe/IEDriverServer.exe");
WebDriver driver = new InternetExplorerDriver();

ソースコード

ソースコードの内容について言及しません。 主に、「Selenium実践入門」と言う本と下記のWebサイトを参考にしました。

PubMedSearchTest.java
package test;

import java.util.concurrent.TimeUnit;

import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.openqa.selenium.WebDriver;
import org.openqa.selenium.chrome.ChromeDriver;

import page.PubMedAdvancedSearchPage;

public class PubMedSearchTest {

    private WebDriver driver;

    @Before
    public void setUp() throws Exception {
        System.setProperty("webdriver.chrome.driver", "exe/chromedriver.exe");
        driver = new ChromeDriver();
        //System.setProperty("webdriver.edge.driver", "exe/MicrosoftWebDriver.exe");
        //driver = new EdgeDriver();
        //System.setProperty("webdriver.ie.driver", "exe/IEDriverServer.exe");
        //driver = new InternetExplorerDriver();

        driver.manage().timeouts().implicitlyWait(5, TimeUnit.SECONDS);
    }

    @After
    public void tearDown() throws Exception {
        //driver.quit();
    }

    @Test
    public void pubmedAdvancedSearchTest() {
        // PubMed Top
        PubMedAdvancedSearchPage pmPage = new PubMedAdvancedSearchPage(driver);

        // Go to Advanced Search page
        pmPage.goToAdvancedSearch();

        // Set keywords
        pmPage.setJournal("Nature communications");
        pmPage.setJournal("Nature biotechnology");
        pmPage.setSearchSelect(1, "OR");
        pmPage.setJournal("Nature methods");
        pmPage.setSearchSelect(2, "OR");
        pmPage.setJournal("Nature genetics");
        pmPage.setSearchSelect(3, "OR");
        pmPage.setJournal("Molecular cell");
        pmPage.setSearchSelect(4, "OR");
        pmPage.setJournal("eLife");
        pmPage.setSearchSelect(5, "OR");
        pmPage.setJournal( "PLoS biology");
        pmPage.setSearchSelect(6, "OR");
        pmPage.setJournal("Genome research");
        pmPage.setSearchSelect(7, "OR");
        pmPage.setJournal("Genes & development");
        pmPage.setSearchSelect(8, "OR");
        pmPage.setJournal("Nature cell biology");
        pmPage.setSearchSelect(9, "OR");

        // Click Search button
        pmPage.clickSearchButton();
    }

}
UCSCGenomeBrowserTest.java
package test;

import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.openqa.selenium.WebDriver;
import org.openqa.selenium.chrome.ChromeDriver;

import page.UCSCGenomeBrowserUploadPage;

public class UCSCGenomeBrowserTest {

    private WebDriver driver;

    @Before
    public void setUp() throws Exception {
        System.setProperty("webdriver.chrome.driver", "exe/chromedriver.exe");
        driver = new ChromeDriver();
        //System.setProperty("webdriver.edge.driver", "exe/MicrosoftWebDriver.exe");
        //driver = new EdgeDriver();
        //System.setProperty("webdriver.ie.driver", "exe/IEDriverServer.exe");
        //driver = new InternetExplorerDriver();
    }

    @After
    public void tearDown() throws Exception {
    }

    @Test
    public void test() {
        UCSCGenomeBrowserUploadPage gmPage = new UCSCGenomeBrowserUploadPage(driver);
        gmPage.goToCustomTrack();

        gmPage.setGenomeRefSelect("hg19");

        // Upload file
        gmPage.UploadFile("C://Users/imama/Desktop/test1.bed");
        gmPage.submit();

        // Back to custom track page
        gmPage.goBackToCustomTrack();

        // Upload file
        gmPage.UploadFile("C://Users/imama/Desktop/test2.bed");
        gmPage.submit();

        // Back to custom track page
        gmPage.goBackToCustomTrack();

        // Upload file
        gmPage.UploadFile("C://Users/imama/Desktop/test3.bed");
        gmPage.submit();

    }

}
PubMedAdvancedSearchPage.java
package page;

import java.util.concurrent.TimeUnit;

import org.openqa.selenium.By;
import org.openqa.selenium.WebDriver;
import org.openqa.selenium.WebElement;
import org.openqa.selenium.support.ui.ExpectedCondition;
import org.openqa.selenium.support.ui.ExpectedConditions;
import org.openqa.selenium.support.ui.Select;
import org.openqa.selenium.support.ui.WebDriverWait;

public class PubMedAdvancedSearchPage {

    /* WebDriverクラスのインスタンス */
    private WebDriver driver;

    // 明示的な待機
    WebDriverWait wait;

    // キーワードカウンタ
    private Integer keywordCounter;

    // コンストラクタ
    public PubMedAdvancedSearchPage(WebDriver driver){
        this.driver = driver;
        this.keywordCounter = 0;
        this.wait = new WebDriverWait(driver, 10);
    }

    // Advanced saerch画面へ遷移
    public void goToAdvancedSearch() {
        driver.get("https://www.ncbi.nlm.nih.gov/pubmed/");
        driver.findElement(By.linkText("Advanced")).click();

        // Goto PubMed Advanced search
        String previousURL = driver.getCurrentUrl();
        System.out.println(previousURL);
        driver.manage().timeouts().implicitlyWait(30, TimeUnit.SECONDS);
        ExpectedCondition<Boolean> e = new ExpectedCondition<Boolean>() {
              public Boolean apply(WebDriver d) {
                return (d.getCurrentUrl() != previousURL);
              }
        };
        wait.until(e);
        String currentURL = driver.getCurrentUrl();
        System.out.println(currentURL);

        // ページタイトルが変更されるまで待ち
        wait.until(ExpectedConditions.titleContains("Advanced search"));
    }

    // Set journal
    public void setJournal(String journal) {
        // FieldにJournalを選択
        Select select = new Select(driver.findElement(By.id("ff_" + keywordCounter)));
        select.selectByValue("Journal");

        // Journal名を入力
        WebElement searchBox = driver.findElement(By.id("fv_" + keywordCounter));
        searchBox.sendKeys(journal);

        // カウントアップ
        keywordCounter += 1;
    }

    // Set search select
    public void setSearchSelect(Integer id, String value) {
        Select select = new Select(driver.findElement(By.id("fop_" + id)));
        select.selectByValue(value);
    }

    // Click Search button
    public void clickSearchButton() {
        driver.findElement(By.id("search")).click();
    }
}
UCSCGenomeUploadPage.java
package page;

import java.util.concurrent.TimeUnit;

import org.openqa.selenium.By;
import org.openqa.selenium.WebDriver;
import org.openqa.selenium.WebElement;
import org.openqa.selenium.support.ui.ExpectedConditions;
import org.openqa.selenium.support.ui.Select;
import org.openqa.selenium.support.ui.WebDriverWait;

public class UCSCGenomeBrowserUploadPage {

    /* WebDriverクラスのインスタンス */
    private WebDriver driver;

    // 明示的な待機
    WebDriverWait wait;

    // コンストラクタ
    public UCSCGenomeBrowserUploadPage(WebDriver driver){
        this.driver = driver;
        this.wait = new WebDriverWait(driver, 10);

        driver.manage().timeouts().implicitlyWait(5, TimeUnit.SECONDS);
    }

    // Go to Custom track page
    public void goToCustomTrack() {
        driver.get("https://genome.ucsc.edu/");
        driver.findElement(By.id("myData")).click();

        wait.until(ExpectedConditions.visibilityOfElementLocated(By.id("customTracksMenuLink")));
        driver.findElement(By.id("customTracksMenuLink")).click();
    }

    // Upload file
    public void UploadFile(String filePath) {
        WebElement upload = driver.findElement(By.name("hgt.customFile"));
        upload.sendKeys(filePath);
    }

    // Submit
    public void submit() {
        driver.findElement(By.id("Submit")).click();
    }

    // Return Custom track page
    public void goBackToCustomTrack() {
        WebDriverWait wait2 = new WebDriverWait(driver, 10000);
        wait2.until(ExpectedConditions.visibilityOfElementLocated(By.id("addTracksButton")));
        driver.findElement(By.id("addTracksButton")).click();
    }

    public void setGenomeRefSelect(String genomeRef) {
        Select select = new Select(driver.findElement(By.id("db")));
        select.selectByValue(genomeRef);
    }
}

スクリーンショットの取り方について(トラブルシューティング

余談ですが、Java8の場合、FileUtils(Java7以前)ではなく、FileHandler(Java8以降)なので間違いがないように。 https://github.com/SeleniumHQ/selenium/issues/4919

File tmpFile = ((TakesScreenshot)driver).getScreenshotAs(OutputType.FILE);
try {
    FileHandler.copy(tmpFile, new File("/Users/imamachinaoto/Desktop/test.png"));
} catch (IOException e) {
    // TODO 自動生成された catch ブロック
    e.printStackTrace();
}

参考

[初心者向け] JavaSeleniumを動かす
https://qiita.com/tsukakei/items/41bc7f3827407f8f37e8

selenium get current url after loading a page
https://stackoverflow.com/questions/16242340/selenium-get-current-url-after-loading-a-page

Selenium WebDriverのwaitを活用しよう
http://softwaretest.jp/labo/tech/labo-294/

Best Practices & Tips: Selenium File Upload
https://saucelabs.com/resources/articles/best-practices-tips-selenium-file-upload

bioRxivでバズってる論文を抽出してSNSにメッセージを流してみた

きっかけ

近年、NatureやCell, Scienceに論文が投稿される前に、preprintサーバであるBioRxivに投稿する流れが増えてきています(2015年以降、急激に増加。ちなみに、Natureなどの多くの主要ジャーナルが、preprintサーバへ事前に投稿することを許可しているので、BioRxivへの投稿はいわゆる二重投稿に引っかかりません)。

また、BioRxivへの投稿が増えている背景には、オープンイノベーションを推進するという目的だけでなく、論文のPriorityを確保することや研究の宣伝を行う目的なども相まってのことだと考えられます。

とにかく、これからのサイエンスでは最新の研究動向を理解する上で、BioRxivの論文をフォローしていくことはますます重要になってくるでしょう。しかしながら、その論文の本数が多すぎます…。

直近で、preprintサーバに1ヶ月あたり1673本の論文が登録されています(1日あたり約50〜60本)。今後、さらに数が伸びることが予想されます。これを人力で全て確認するのは不可能です…。 image.png (218.9 kB) 元データはこちら。
http://www.prepubmed.org/monthly_stats/

BioRxivに登録されている論文はとにかく玉石混交なので、注目されている論文だけを抽出することはできないか考えました。そんな中で見つけたのが、こちらの記事です。

A Twitter bot to find the most interesting bioRxiv preprints
https://gigabaseorgigabyte.wordpress.com/2017/08/08/a-twitter-bot-to-find-the-most-interesting-biorxiv-preprints/

こちらの記事では、 Altmetric scoreという指標(論文の注目度を数値化したもの)を用いて、BioRxivに投稿された論文の中で、特に注目度の高い記事を抽出し、Tweetするという仕組みを作っていました。

ちなみに、Altmetric scoreとは、TwitterFacebookなどのSNSやNews、Mendeleyへの登録数などの様々な指標を元に算出した論文の注目度を示すスコアです(詳しく解析したわけではないですが、個人的な印象として、BioRxivに登録されている論文のAltmetric scoreは、Tweet数に引っ張られている印象を受けましたが…)。

上記の記事で紹介されているスクリプトをそのまま使うのも手でしたが、後々PubMedに登録されている論文に対しても応用したかったのと、データのストアをSQLデータベースではなく、Pickleで管理しているなど、実装部分で個人的に作り変えたい部分があったので、1から自分で作ってみることにしました。

そんなわけで、上記の記事を参考に、Raspberry Pi上で動作するBioRxivのキュレーションシステムを作って見たというのが今回の内容になります。

サンプル

SNSへのメッセージ出力として、SlackとTwitterに対応させました(現状、個人的に興味のある「Genomics」と「Bioinformatics」の分野に関する論文しかキュレーションの対象にしていません)。

サンプルのTwitter botはこちらになります。
https://twitter.com/BioRxivCurator

image.png (546.9 kB)

ソースコード

BioRxivCurator
https://github.com/Imamachi-n/BioRxivCurator

内容

全体像

システム全体の流れは以下の通りです(内部的なデータ処理の流れとはちょっと違いますが、概略図ということで解釈してもらえば良いです)。 image.png (218.7 kB)

RSS feedで取得した論文データについて、Altmetric scoreを取得し、データをSQLiteデータベースに格納しておきます。Altmetric scoreが一定以上の論文(注目度の高い論文)について、定期的にTwitterやSlackにメッセージを配信します。

また、一定のスコアに達しない論文については、情報をデータベースに保持し続けておき、定期的にAltmetric scoreを確認します。データの保持期限は1ヶ月とします(多くの場合、1ヶ月以降はあまりAltmeric scoreが変動しない傾向にあるため)。それを過ぎた論文データは破棄します(厳密には、フラグ管理により論理削除する形をとります)。

Pythonで実装し、プログラム自体はRaspberry Pi3 ModelB上で動かします。手作業でプログラムを動かすのは面倒なので、ジョブスケジューラであるcronを使って毎日論文を自動でチェックする形にしました。 image.png (171.0 kB) https://commons.wikimedia.org/w/index.php?curid=47497384

Pythonによる実装

まず、Pythonによる実装例を見ていきます。かいつまんで重要な箇所を説明しているだけなので、詳しくは上述のソースコードをご覧ください。

RSS feedの取得

BioRxivのRSS feedはこちらから取得できます。直近の30 articlesを取得できます。
https://www.biorxiv.org/alertsrss

論文の絞り込みとしては、
http://connect.biorxiv.org/biorxiv_xml.php?subject=all
で全体から取得するか、
http://connect.biorxiv.org/biorxiv_xml.php?subject=キーワード
で特定の分野の論文のみを取得する方法があります。

後者の場合、+でつなげて
http://connect.biorxiv.org/biorxiv_xml.php?subject=genomics+bioinformatics
のように複数の分野からRSSを取得することも可能です。

PythonRSS feedから得たデータをパースするために、今回はfeedparserライブラリを使いました。

import feedparser

feed = feedparser.parse(
        "http://connect.biorxiv.org/biorxiv_xml.php?subject={0}".format("+".join(subjects)))

 for pub in feed["items"]:
    <処理を書く>
     # pub["dc_identifier"] # DOI
     # pub["title"]          # 論文のタイトル
     # pub["link"].split('?')[0]  # 論文のURL
     # pub["updated"]        # 更新日

以上のように、feedparser.parse()に先ほど取得したRSS feedのURLを渡すことで、ディクショナリ型にパースされたRSSデータを取得できます。

各論文のデータは、itemsごとに格納されているので、Forループで各論文の情報を取得します。例えば、dc_identifierにはDOI、titleには論文のタイトル、linkには論文のリンク、updatedには更新日が格納されています。DOIはAltmetric scoreの取得の際に必要となるので、必ず取得しておきます。

実際にどんなデータが入っているか確認したい場合は、RSS feed URLに直接アクセスして、中身のXMLファイルを覗いてみるといいです。 ブラウザによって見え方が異なります。以下の例では、Google Chromeを使って表示しています。 image.png (389.8 kB)

Altmetric scoreの取得

RSS feedから取得したDOIを使って、Altmetric scoreを取得します。 Altmetric scoreの取得には、AltmetricのAPIを使用します。
https://api.altmetric.com/

ソースコード中では、Altmetric APIの薄いラッパークラスを作っていますが、要はAltmetric scoreの取得先URLに対してGETしてるだけです。

import requests
import json

# GET request
doi = "10.1038/480426a"
request = requests.get("https://api.altmetric.com/v1/doi/{0}".format(doi))
response = json.loads(request.text)
# response["context"]['journal']['pct']    # Altmetric scoreを「0〜100」に正規化した値。
# response["score"]            # Altmetric score

Altmetric APIは、
リクエストURL: https://api.altmetric.com/v1/doi/DOI
に対してGETすると、JSON形式でデータを返してくれます。

そこで、JSON形式のデータをjson.load()を使ってディクショナリ型に変換します。 あとは、キーを指定して欲しいデータを取り出すだけです。

こちらも、実際にどんなデータが入っているか確認したい場合は、リクエストURLに直接アクセスして、中身のJSONファイルを覗いてみるといいです。
ブラウザによって見え方が異なります。以下の例では、Firefoxを使って表示しています。Firefoxではこのように、 JSON形式のデータを見やすく表示してくれます。 image.png (298.8 kB) 今回は、赤枠で囲った値を取得しています。

SQLiteデータベースへの格納

順番は前後しますが、ここでまとめてSQLiteデータベースの処理を説明します。

1. テーブルの作成(CREATE TABLE

まず、テーブルを作成します。 ここでは、CREATE TABLE IF NOT EXISTSとすることで、システムの初回起動時にのみテーブルが作成される仕組みになっています。 また、論文のDOIをprimary keyに指定します(各論文に対してユニークな文字列が割り振られるため)。

import sqlite3

sqlite3_file = "./storeAltmetrics.sqlite3"
with sqlite3.connect(sqlite3_file) as conn:
        c = conn.cursor()

        # Create biorxiv_altmetrics_log table
        sql = """CREATE TABLE IF NOT EXISTS biorxiv_altmetrics_log
                (doi TEXT,
                 title TEXT,
                 link TEXT,
                 update_date TEXT,
                 altmetric_score INTEGER,
                 altmetric_pct INTEGER,
                 altmetric_flg INTEGER,
                 PRIMARY KEY(doi)
                )"""
        c.execute(sql)
        conn.commit()

基本的に、Javaなどの言語と書き方は似てます。 sqlite3.connect()でデータベースに接続します。ここでは、sqlite3_fileが格納先のデータベースのファイルだとします。 続いて、cursorオブジェクトを作成し、execute()メソッドでSQL文を実行します。最後に、変更点を保存するために、Connectionオブジェクトのcommit()メソッドを実行します(コミット)。

テーブルのフィールドについてです。 RSS feedから取得したデータをdoi, title, link, update_dateにそれぞれ格納します。 また、Altmetric APIから取得したスコアを、altmetric_score, altmetric_pctにそれぞれ格納します。

フラグ 意味
0 PCTが90未満
1 PCTが90以上
-1 PCTが90未満で1ヶ月が経過(解析対象から除外)
2. RSS feedから得た論文データの格納(INSERT INTO

RSS feedから取得した論文のデータ(論文のタイトル、論文のURL、DOI、更新日)をデータベースのテーブルに格納します。 ここでは、INSERT OR IGNORE INTOとすることで、テーブルに格納されていない論文についてのみ登録する形になっています。

with sqlite3.connect(sqlite3_file) as conn:
            c = conn.cursor()

            # Insert article info into biorxiv_altmetrics_log if not already exists
            sql = """INSERT OR IGNORE INTO biorxiv_altmetrics_log
                     VALUES(?,?,?,?,?,?,?)"""
            doi_info = [tuple([p.doi, p.title, p.url, p.date, 0, 0, 0])
                        for p in RSS_data_list]
            c.executemany(sql, doi_info)
            conn.commit()

ここでは、RSS_dataクラスのインスタンスRSS_data_listに論文のデータが格納されていることを想定しています。 データを渡す際には、タプル型である必要があります。

("value1,")
("value1", "value2", "value3")

また、タプル型のデータのリストを作成し、cursorオブジェクトのexecutemany()メソッドを使うことで、一気に複数のレコードをテーブルに格納することができます。最後は、connectionオブジェクトのcommit()メソッドを使って、変更点をコミットします。

ちなみに、突然出てきたRSS_data_listですが、下記の通りに定義しています。

class RSS_data(object):
    def __init__(self, doi, title, url, date):
        self.doi = doi
        self.title = title
        self.url = url
        self.date = date

RSS_data(doi=pub["dc_identifier"],
         title=pub["title"],
         url=pub["link"].split('?')[0],
         date=pub["updated"]))

またここでは、pubは上述したfeedparser.parse()メソッドで取得したRSS feedから得たデータをディクショナリ型で格納したオブジェクトだとします。

3. Altmeric scoreを取得する論文データの抽出(SELECT

続いて、過去に取得した論文データを含めて、一定のAltmetric scoreを超えていない論文のリストを取得します。
先ほど説明した通り、altmetric_flgが「0」である論文がAltmetric scoreを検索する対象となります。

with sqlite3.connect(sqlite3_file) as conn:
            c = conn.cursor()

            # Select target doi from biorxiv_altmetrics_log
            sql = """SELECT doi, title, link, update_date from biorxiv_altmetrics_log
                     WHERE altmetric_flg = 0"""
            c.execute(sql)

            # Store doi data as target_doi_data object
            target_doi_list = []
            for doi_info in c.fetchall():
                target_doi_list.append(target_doi_data(
                    doi=doi_info[0], title=doi_info[1], url=doi_info[2], date=doi_info[3]))

cursorオブジェクトのexecute()メソッドでSQL文を実行した後、検索した結果はcursorオブジェクトに格納されます。fetchall()メソッドを使うことで、条件に一致した全てのレコードをリストとして取得できます。

class target_doi_data(object):
    def __init__(self, doi, title, url, date):
        self.doi = doi
        self.title = title
        self.url = url
        self.date = date

取得したレコードの情報は、上記target_doi_dataオブジェクトのリストとして格納しておきます。

4. Altmetric scoreの値の格納(UPDATE

最後に、Altmetric scoreの更新を行います。APIから取得したスコアをテーブルに書き込みます。最後に忘れずに、変更をコミットします。

with sqlite3.connect(sqlite3_file) as conn:
            c = conn.cursor()

            # insert altmetric score into biorxiv_altmetrics_log
            sql = """UPDATE biorxiv_altmetrics_log
                     SET altmetric_score = ?,
                          altmetric_pct = ?,
                          altmetric_flg = ?
                     WHERE doi = ?"""
            c.execute(sql, tuple([altmetrics_data.altmetric_score,
                                  altmetrics_data.pct, altmetrics_data.flg,
                                  doi]))
            conn.commit()

Slackのアクセストークンを取得する

SlackのAPIを介してデータのやり取りをします。そのため、Slackのアクセストークンを取得する必要があります。以下ではその取得方法について説明します。

基本的に以下の記事を参考にしました。
https://qiita.com/ykhirao/items/3b19ee6a1458cfb4ba21

1. アプリの登録

まず、以下のURLにアクセスします。Slackにログインしていない場合、ログインしてください。
https://api.slack.com/apps

Create an Appをクリックします。 image.png (250.5 kB)

Create a Slack Appという画面が出てくるので、
App Name: 好きな名前を入力する。
Development Slack Workspace: Slackのメッセージを流したいワークスペースを選択。
を入力して、Create Appをクリックします。 image.png (194.2 kB)

2. スコープの設定

登録したアプリの権限の範囲を設定します。 Install your app to your workspaceを開くと、権限のスコープを設定してくださいとメッセージが出ているので、permission scopeをクリックします。 image.png (280.0 kB)

今回は、Slackへメッセージを送るだけなので、ChatのSend messages as xxxxxxを選択します。 image.png (233.3 kB)

すると、選択したPermission Scopeが表示されます。
状態を保存するために、Save Changesを必ずクリックしてください。 image.png (181.5 kB)

3. アプリをSlackにインストールする

最後に、アプリをSlackにインストールします。 image.png (168.1 kB)

先ほど設定したメッセージを送る権限を付与していいか聞かれるので、Authrizeをクリックします。 image.png (81.0 kB)

インストールが完了すると、OAuth Access Tokenを取得できます。 この文字列を使うので、コピーして保存しておきます。 image.png (119.7 kB)

Slackへのメッセージ配信

Slackへのメッセージ送信は、公式のPythonライブラリであるSlackClientを使います。
https://slackapi.github.io/python-slackclient/basic_usage.html#sending-a-message

使い方は簡単で、先ほど取得したSlackのアクセストークslack_tokenを指定し、spi_call()メソッドに、送り先のチャンネルchannelと送信するテキストtextを指定するだけです。エラーハンドリングについては、上記のURLの記事を参照してください。

slack_token = "xxxxxxxx"
channel = "@imamachi"
text = "Hello World !!"

sc = SlackClient(slack_token)
    response = sc.api_call(
        "chat.postMessage",
        channel=channel,
        text=message
    )

Twitterのアクセストークンを取得する

Slack同様、TwitterAPIを介してデータのやり取りをします。そのため、Twitterについてもアクセストークンを取得する必要がありmす。以下ではその取得方法について説明します。

1. Twitterアプリの作成

まず、Twitter Application Managementにアクセスします。アカウントにログインしていない場合は、ツイート対象となるアカウントにログインしてください。
https://apps.twitter.com/

Create New Appをクリックします。 image.png (84.9 kB)

Name: アプリ名を記載する。
Description: アプリの説明を記載する。
Website: ソースコードやアプリが入手できるURLを入力。
Create your Twitter applicationをクリックし、アプリを作ります。 image.png (343.5 kB)

image.png (55.1 kB)

すると、次のようなページに遷移します。 image.png (223.4 kB)

2. アクセストークンの取得

Keys and Access Tokensタブをクリックし、 Consumer Key (API Key)Consumer Secret (API Secret)を確認します。 これらはAPIでのやり取りに必要となるので、コピーして保存しておきます。 image.png (318.7 kB)

続いて、Your Access TokenでCreate my access tokenをクリックします。 image.png (65.5 kB)

Access TokenAccess Token Secretを取得します。 これらも同様に、APIでのやり取りに必要となるので、コピーして保存しておきます。 image.png (85.4 kB)

Twitterへのメッセージ配信

Twitterへのツイートは、tweepyというライブラリを使います。 こちらも使い方は非常に簡単で、先ほど取得した各種の値を設定し、メッセージをupdate_status()メソッドに指定するだけです。

こちらも参照。
https://review-of-my-life.blogspot.jp/2017/07/python-cloud9-tweepy.html

consumer_key = "xxxxxxxx"
consumer_secret = "xxxxxxxx"
access_token = "xxxxxxxx"
access_token_secret = "xxxxxxxx"
message = "Hello World !!"

auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_token, access_secret)
api = tweepy.API(auth)
api.update_status(message)

追記: 余談ですが、TweepyはメインのContributorがいなくなり、他のmaintainerに引き継がれています。アップデートが最近されていないので、TwitterAPIの仕様変更で動かなくなるリスクがあります。
https://qiita.com/utgwkk/items/4beed333e8262c675028

Raspberry Pi3の用意

以下のサイトを参考にRaspberry Pi3のOSをインストールしたMicro SDカードを用意します。
https://getpocket.com/redirect?url=http%3A%2F%2Fhirazakura.hatenablog.com%2Fentry%2Fraspberrypi%2Fsetup%2Ffirst&formCheck=3f0b1d005779f9a1ac22bfac92159bea

続いて、アップデートをかけて、リブートします(再起動)。

$ sudo apt-get update
$ sudo apt-get upgrade
$ sudo rpi-update
$ sudo reboot

テキストエディタが入っていないので、VSCodeをインストールします。
https://code.headmelted.com/ 管理者権限でないとインストールできないので注意してください。

$ wget -O - https://code.headmelted.com/installers/apt.sh 
$ chmod 755 apt.sh
$ sudo apt.sh

以下でVisual Studio Codeを実行できます。

$ code-oss

SQLiteのクライアントアプリをインストール

SQLiteデータベースをGUIで見たい場合、以下のクライアントアプリがオススメです。
http://sqlitebrowser.org/

Raspberry Pi上では、以下のコマンドでインストール可能です。

$ sudo apt-get install sqlitebrowser

cronの設定

Raspberry Piで作成したPythonスクリプトを定期的に自動で実行するために、cronを使います。 cronはUnix OS系に標準で搭載されているジョブスケジューリングのプログラムです。

cronについての記事はこちらを参照のこと。
https://qiita.com/hikouki/items/e744b3a4d356d2af12cf

ここでは、簡単な設定例を示したいと思います。 まず、cronのconfigureファイルcron.confを作成します。 スペース区切りで日時を指定、最後に実行コマンドを指定します。 行の最後に改行が必要なので、入れ忘れないように。

分 時 日 月 曜日 <実行コマンド> [改行]

あと、念のため、エラーになった時にエラーログを確認できるように、標準出力でエラーログを取っておきます。

0 11 * * * bash /home/pi/Desktop/BioRxivCurator/src/startup4RaspberryPi.sh > /home/pi/Desktop/error.txt 2>&1

cronのconfigureファイルを作成したら、以下のコマンドを実行します。

$ crontab ./cron.conf

設定内容を確認したい場合は、

$ crontab -e

で設定内容を確認・編集できます。

間違って登録してしまった場合、

$ crontab -r

で全ての設定を削除することができます(このコマンドの実行には注意が必要です)。

cron使用時の注意点

だいたいのエラーは、上記pythonへのパスに起因した問題です。

相対パスに注意

cronが実行するコマンドは、実行ユーザーのホームディレクトリで実行されます。

特に今回の場合、pythonのコードの中に相対パスを設定している箇所があるので、実行コマンド中で、cdでカレントディレクトリをソースコードの置いてあるディレクトリに移動しておきます。

$ cd /home/pi/Desktop/BioRxivCurator/src

Minicondaを使用している場合

python絶対パスで指定する必要があります。そうしないと、cronは標準でインストールされているpythonを使おうとします。

/home/pi/miniconda3/bin/python /.main.py --yaml_setting_file ./production.yaml

原因がわからないエラーの場合

cronで実行するコマンドの最後に、エラー出力先を指定しておくと良いです。

0 11 * * * bash /home/pi/Desktop/BioRxivCurator/src/startup.sh > /home/pi/Desktop/error.txt 2>&1

こうすることで、指定したファイルerror.txtにエラーログが出力される。

また、/var/log/syslog内にもcron実行時のログが残るので、そちらも参照すると原因の特定につながります。

cron.confが登録できない

コマンドの打ち間違え(こんなミスをするのは私だけかもしれませんが…)。

$ cron ./cron.conf
cron: can't lock /var/run/crond.pid, otherpid may be 3505: Resource temporarily unavailable

正しくは、cronではなくcrontabです。

$ crontab ./cron.conf

どうしようもなくなったら、とりあえず、リブート(再起動)してみましょう。

$ reboot

Windows Mixed RealityのImmersive Headsetの初期設定

概要

Acerから出たAcer Windows Mixed Reality Development Editionが届いたので、セットアップ方法をまとめました。

IMG_20170826_145612.jpg (2.6 MB)

設定方法

(1)HDMIケーブルとUSB3.0ケーブルをパソコンに差し込みます。

(2)しばらく待っていると、次のような画面が出てきます。 image.png (562.6 kB)

(3)次へを押していくと、接続しているパソコンのスペックが動作要件を満たしているかどうかチェックが行われます。問題なければ、次へを押します。 image.png (188.2 kB)

(4)HMDの説明が入り、ヘッドフォンジャックの位置などが示されます。 image.png (158.5 kB)

(5)身長を聞かれます(メートル法じゃない…)。 image.png (176.0 kB)

フィートとかインチでは身長がわからないので、メートル法を指定して自分の身長を入力します。 image.png (170.8 kB)

(6)フロアの設定が始まります。 image.png (477.6 kB)

HMDかぶって部屋の中央に立ち、開始ボタンを押します。 image.png (192.3 kB)

10秒間待つと、フロアが検出されましたと表示されます。 image.png (182.4 kB)

(7)続いて、フロアの境界を検知します。 image.png (664.1 kB)

できるだけ障害物を除いて、 image.png (669.3 kB)

腰のあたりにHMDを持ち、フロアの境界(壁)をスキャニングしていきます。 image.png (238.4 kB)

一周りすると、いかのように境界が設定されます(あんまりうまくない例です)。 image.png (262.3 kB)

(8)以上で基本的な設定は終わりで、 image.png (176.0 kB)

複合現実ポータルに入れます。 image.png (779.4 kB)

右クリックで、場所の移動。 キーボードのWindowsキーを押すと、メニューが出現します。 全体的な操作感は、HoloLensと似ています。 image.png (902.6 kB)

簡単ではありますが、以上になります。

Ribo-seqデータ解析 (RiboTaper, short ORF prediction)

概要

Ribo-seq (Ribosome Profiling)とは、リボソームが結合しているRNA領域をシークエンスする手法です。この方法を利用することで、mRNAのどの領域が翻訳されているか、つまりはORFを予測することができます。また、特定の刺激などによる翻訳活性の変化を捉えることができます。

ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、Ribo-seqのデータ解析の詳細をStep by Stepで説明します。ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。

使用するデータ

SRA160745 (Gao X et al. Nat. Methods 2015)

https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=SRA160745

使用するソフトウェア

解析ワークフロー

今回の解析の流れを以下にまとめました。

PowerPoint スライド ショー - [2016-08-22_NGS_workflow.pptx] 2017-03-30 10.21.48.png (140.1 kB)

サンプルコード

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

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

権限の問題

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

$ chmod 755 <hoge.sh>

改行コードの問題

WindowsLinuxで改行コードが異なるため、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

まず、Biocondaで仮想環境をつくります。ここでは、python2.7をベースに仮想環境を構築します。

conda create -n ribotaper python=2.7
  • -nで仮想環境の名前を指定します。
  • pythonに2.7を指定することで、python2.7を仮想環境にインストール

source activateに続けて、先ほど作成した仮想環境の名称を指定することで、仮想環境に入ることができます。

source activate ribotaper

conda info --envsで現在どの環境にいる確認できるので、前後でどの環境にいるか確認してみましょう。

$ conda info --envs
  # conda environments:
  #
  ribotaper                /home/akimitsu/miniconda2/envs/ribotaper
  root                  *  /home/akimitsu/miniconda2

$ source activate miso
$ conda info --envs
  # conda environments:
  #
  ribotaper             *  /home/akimitsu/miniconda2/envs/ribotaper
  root                     /home/akimitsu/miniconda2

最初はrootという環境にいたのが、source activate ribotaperの後、ribotaperという仮想環境に移動していることがわかると思います。

RiboTaperの実行は、このribotaperと名付けた仮想環境上で行うこととします。また、その他の必要なソフトウェアに関しても同様の環境上にインストールしていきます。

それでは、RiboTaperに必要なパッケージやソフトウェアをBiocondaを利用してあらかじめインストールしておきましょう。

注意!!

samtoolsとbedtoolsはバージョンに依存しているので、最新のものをインストールしてはならない。

$ conda install fastqc
$ conda install cutadapt
$ conda install fastx_toolkit
$ conda install star
$ conda install samtools
$ conda install rseqc
$ conda install bedtools
$ conda install subread

# For RiboTaper
$ conda install samtools=0.1.19
$ conda install bedtools=2.17.0
$ conda install r-seqinr
# conda install r-ade4
$ conda install r-multitaper
$ conda install r-xnomial

# DOWNGRADED due to dependency conflicts (R(v3.2.2-0))
$ conda install r-domc
# conda install r-iterators
# conda install r-foreach

RiboTaperのソースコードをダウンロードして、ファイルを解凍します。

$ wget https://ohlerlab.mdc-berlin.de/files/RiboTaper/RiboTaper_v1.3.tar.gz
$ tar zxvf RiboTaper_v1.3.tar.gz

RiboTaper_v1.3ディレクトリにパスを通します。下記のように、ホームディレクトリにある.bashrcにパスを書き加えます。

export PATH=/home/akimitsu/software/RiboTaper_v1.3/scripts:${PATH}

1. SRAダウンロード

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

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

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

まず、下記のURLにアクセスすると、
https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=SRA160745

DDBJ DRAのサイトはこんな感じ。NCBI GEOと比較してIDが羅列してあるだけなので、どのIDがどのサンプルに対応しているかぱっと見わかりません…。

SRA160745 - DRA Search と 2 ページ - Microsoft Edge 2017-03-30 18.23.54.png (155.3 kB)

下記のサイトでDRAに登録されているデータの詳細をリスト化してみることができます。

DBCLS SRA: Survey of Read Archives
http://sra.dbcls.jp/v1/

早速、SRA#の欄にSRA160745と入力して、検索してみましょう。

DBCLS SRA と 7 ページ - Microsoft Edge 2017-03-30 19.15.36.png (109.3 kB)

Expsの項目をクリックします。

STUDIES と 10 ページ - Microsoft Edge 2017-03-30 19.27.42.png (36.5 kB)

すると、各IDに対してタイトルなどの情報を見ることができます。

今回使用するのは、No16とNo19のデータになります。Exp#の項目のSRX740748とSRX740751をそれぞれクリックします。

EXPERIMENTS と 10 ページ - Microsoft Edge 2017-03-30 19.28.11.png (137.1 kB)

DDBJ DRAのページに移動します。NavigationのRunのSRAボタンをクリックします。

SRX740748 - DRA Search と 10 ページ - Microsoft Edge 2017-03-30 19.32.40.png (42.9 kB)

SRR1630831のURLをコピーします、

_ddbj_database_dra_sralite_ByExp_litesra_SRX_SRX740_SRX740748_SRR1630831 のインデックス - Google Chrome 2017-03-30 19.33.17.png (29.7 kB)

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

$ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/SRX/SRX740/SRX740748/SRR1630831/SRR1630831.sra
$ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/sralite/ByExp/litesra/SRX/SRX740/SRX740751/SRR1630838/SRR1630838.sra

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

まず、SRA IDだけだったファイル名にサンプルの情報を付け加えます(SRRのIDだけだと何のファイルかわからなくなるので)。mvコマンドを使ってリネームを行います。

$ mv SRR1630831.sra Ribo-Seq_HEK293_Cell_Control_SRR1630831.sra
$ mv SRR1630838.sra RNA-Seq_HEK293_Cell_SRR1630838.sra 

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

$ fastq-dump Ribo-Seq_HEK293_Cell_Control_SRR1630831.sra
$ fastq-dump RNA-Seq_HEK293_Cell_SRR1630838.sra 

カレントディレクトリにFASTQファイルが展開されているはずです。

3. Quality check

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

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

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

実行例

$ mkdir fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831
$ fastqc -t 8 -o ./fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831 ./Ribo-Seq_HEK293_Cell_Control_SRR1630831.fastq -f fastq

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

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

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

FastQCから出力されるデータが、fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831ディレクトリの中のRibo-Seq_HEK293_Cell_Control_SRR1630831_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を参照して得ます。

今回は、論文から情報を得ました。Ribo-seqで用いられるアダプター配列はほとんど同じで、下記の論文に記載されている、TCGTATGCCGTCTTCTGCTTGという配列です。

Hafner, et al. Methods, 44:3-12 (2008);
https://www.ncbi.nlm.nih.gov/pubmed/18158127

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

$ less Ribo-Seq_HEK293_Cell_Control_SRR1630831.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)

上図では、今回のサンプルとは異なるサンプルを見ています。

実行例

$ cutadapt  -m 10 -a TCGTATGCCGTCTTCTGCTTG \
  Ribo-Seq_HEK293_Cell_Control_SRR1630831.fastq > Ribo-Seq_HEK293_Cell_Control_SRR1630831_1_trimmed_adapter.fastq 2>> ./log_Ribo-Seq_HEK293_Cell_Control_SRR1630831.txt
  • -f: フォーマットの指定。
  • -m: トリミング後のリードの最低配列長。例えば、10と指定すると、トリミング後にリード長が10未満のリードは排除する。
  • -a: 3'側にライゲーションされている、指定したアダプター配列を除去する。

5. Quality filtering

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

コード例

$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_1_trimmed_adapter.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o Ribo-Seq_HEK293_Cell_Control_SRR1630831_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由来のリードを除く

リボソームをIPしている関係で、Ribo-seqのFASTQファイルにはある程度リボソームRNA由来のコンタミしています。そのため、ゲノムにマッピングする前に、リボソームRNA由来のリードを除去します。

まず、リボソームRNAの配列を用意します。 NCBIの登録されているX12811.1| Human 5S DNAU13369.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としてまとめます。

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

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

実行例

$ mkdir ./STAR_Index_rRNA_contam
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./STAR_Index_rRNA_contam \
     --genomeFastaFiles ./contam_Ribosomal_RNA.fa

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

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

rRNAを除く

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

実行例

$ STAR --runMode alignReads --runThreadN 8 --genomeDir ./STAR_Index_rRNA_contam \
--readFilesIn ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_2_filtered.fastq \
--outSAMunmapped Within --outFilterMultimapNmax 30 --outFilterMultimapScoreRange 1 \
--outFileNamePrefix ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq \
--outSAMattributes All --outStd BAM_Unsorted --outSAMtype BAM Unsorted \
--outFilterType BySJout --outReadsUnmapped Fastx --outFilterScoreMin 10 \
--alignEndsType EndToEnd > ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_repbase_rrna_comtam.bam
$ mv ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastqUnmapped.out.mate1 ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_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に由来するリードを除いた後、フィルタリング後のFASTQファイルをFastQCを用いて再確認します。

実行例

$ mkdir fastqc_${file}_filtered
$ fastqc -o ./fastqc_Ribo-Seq_HEK293_Cell_Control_SRR1630831_filtered ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq -f fastq

8. Mapping

Indexファイルの作成

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

まず、Human Genome(hg38)のFASTAファイルのダウンロードしてきます。 各染色体ごとにファイルが用意されているので、catコマンドでダウンロードしたすべてのファイルをマージして、hg38.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/hg38/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/24.html

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

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

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

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

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

ここまでで、ゲノムとトランスクリプトームのIndexファイル用意できたので、いよいよSTARを用いて各シークエンスリードマッピングします。

実行例

$ mkdir hg38_Genome_Gencode_v24
$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir hg38_Genome_Gencode_v24 \
     --genomeFastaFiles ./hg38.fa --sjdbGTFfile ./gencode.v24.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ファイルが用意できたので、ゲノム(hg38)とトランスクリプトーム(Gencode v24)にリードをマッピングします。

実行例

$ mkdir STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831
$ STAR --runMode alignReads --runThreadN 8 --genomeDir ${indexFile} \
  --readFilesIn ./Ribo-Seq_HEK293_Cell_Control_SRR1630831_3_rm_repbase_rrna.fastq \
  --outFilterMultimapNmax 8 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --sjdbScore 1 \
  --outFilterMismatchNmax 4 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 \
  --outFileNamePrefix ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_ \
  --outSAMattributes All --outSAMtype BAM SortedByCoordinate --limitBAMsortRAM ${maxRAM} \
  --outFilterType BySJout --outReadsUnmapped Fastx
STARのパラメータ。
  • --outSAMtype BAM SortedByCoordinate: SAMファイルでなく、BAMファイルとしてマッピングされた情報を直接出力する。出力ファイルは染色体の座標でソートする。
  • --limitBAMsortRAM 32000000000: BAMファイルをソートするときに割り当てられるメモリの上限を設定。

9. RiboTaperによるORF予測

80Sリボソームと結合するRNA領域の長さは28-30nt程度であることが知られており、実際にRibo-seqのリード長も同じ長さを取っています。

また、コドンは3塩基ごとに読まれていくので、マッピングされたシークエンスリードの分布には、ORFの読み枠に合わせてバイアスがかかっていることが知られています。

例えば、RibosomeのP-siteにmRNAのAUG配列がロードされている場合、シークエンスリードの5'末から12nt下流の位置を1とすると、AUGのAがそれに該当します(13-15nt目がリボソームのP-siteの位置になります)。

------------AUG---------------
<-- 12nt -->1<----15/16nt---->

Ribo-seqのデータでは、先ほど定義したP-siteの相対的な位置(1: 5'末から13nt目)に、コドンの最初の塩基が来るケースが多いです(これがマッピングされたリードの偏り)。 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4193932/

RiboTaperでは、上記のようなRibo-seqにおけるシークエンスリードの分布の偏りから、リードがマッピングされた領域がORFとしてリボソームに読まれて翻訳されている領域かどうか判断しています。

また、RiboTaperでは、このシークエンスリードマッピングバイアスを事前にチェックするプログラムが用意されているので、それを使ってデータのクオリティに問題がないか最初に確認しましょう。

# RiboTaper用のアノテーション情報の準備

まず、ORFなどのアノテーション情報をRiboTaper用に用意する必要があります。RiboTaper側で用意されているcreate_annotations_files.bashシェルスクリプトを用いて、GencodeなどのGTFファイルからRiboTaper用のアノテーション情報を作成します。

$ create_annotations_files.bash  \
  ./gencode.v24.basic.annotation.gtf \
  ./hg38.fa \
  true \
  true \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \
  /home/akimitsu/software/RiboTaper_v1.3/scripts
  • 1つ目の引数: アノテーション情報のGTFファイルを指定。
  • 2つ目の引数: ヒトゲノム(hg38)のFASTAファイルを指定。
  • 3つ目の引数: CCDS annotationを利用するかどうか(Gencodeのデータを推奨)。
  • 4つ目の引数: Appris annotationを利用するかどうか(Gencodeのデータを推奨)。
  • 5つ目の引数: 保存先のディレクトリを指定。
  • 6つ目の引数: 使用するBedtoolsのパスを指定。
  • 7つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。

CCDSとApprisについては下記のサイトを参照のこと。 https://www.gencodegenes.org/gencode_tags.html

# Ribo-seqデータのチェック

冒頭で説明したとおり、まず、マッピングしたRibo-seqのBAMファイルをチェックします。

実行例

$ create_metaplots.bash \
  ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24/start_stops_FAR.bed \
  metaplots_result \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir/ \
  /home/akimitsu/software/RiboTaper_v1.3/scripts/
  • 1つ目の引数: マッピングしたRibo-seqのBAMファイルを指定。
  • 2つ目の引数: 先ほど作成したRiboTaper用のアノテーション情報のstart_stops_FAR.bedというファイルを指定。
  • 3つ目の引数: 保存するファイルの名称を指定。
  • 4つ目の引数: 使用するBedtoolsのパスを指定。
  • 5つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。

実行すると、リード長ごとに、開始コドンと終止コドンの周辺のリードの分布を図にしてくれます。

緑色のバーがORFの読み枠(例えば、AUGのA)と一致しているリードの数になります。結果の図を見たときに、良い例のように緑色のバーにリードの分布が偏っているのが使えるシークエンスデータになります。

良い例

metaplots_result_28.pdf - Adobe Acrobat Reader DC 2016-09-01 15.47.32.png (103.9 kB)

一方で、悪い例のようにリードの分布にバイアスがイマイチかかっていないサンプルは、ノイズの大きいデータになるので、RiboTaperによる解析には使えません。

悪い例

metaplots_result_28.pdf - Adobe Acrobat Reader DC 2016-09-01 15.48.43.png (107.0 kB)

データのクオリティに問題がないと判断できたら、RiboTaperのシェルスクリプトを実行します。

RiboTaperの実行

RiboTaperには計算上のボトルネックが存在したり、プログラムが不安定で途中で処理が止まってしまったりします。

今回は、スクリプトを修正したり、処理を細分化して実行時にエラーが起こらないようにしたファイルを用意したので、そちらを使用した場合の実行例を示したいと思います。

とりあえず、下記のサイトにおいてあるscriptsディレクトリを、自分の環境のRiboTaper_v1.3の直下に上書きすればOKです。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Ribo-seq_RiboTaper
(注: RiboTaperのソースコードGNU General Public Licenseのもとに配布されています。)

参考(スクリプトの修正箇所)

念のため、スクリプトの修正箇所について言及しておきます。 修正を加えている点はすべて、Rスクリプトを実行している箇所です。

もとのRスクリプトでは、1行目が、

#!/usr/bin/Rscript

となっており、スクリプトの実行はRscriptを利用したコマンドライン上での実行を行っていません。 例えばこんな感じ。

./hoge.R

そのため、Rscript/usr/bin/Rscriptに存在する環境下では機能しますが、Biocondaなどでローカル環境にRをインストールした場合には機能しません。

そこで、以下のように書き換えます。

Rscript ./hoge.R

こうすることで、/usr/bin/RscriptにRscriptがなくても、Rスクリプトが機能します。

修正箇所は、以下のとおりです。

CCDS_orf_finder.R (821行目)
# syst_scr<-paste(scr,"bed_tocheck_ccds",sep = " ")
scr2<-paste("Rscript",scr,sep=" ")
syst_scr<-paste(scr2,"bed_tocheck_ccds",sep = " ")
NONCCDS_orf_finder.R (738行目)
# syst_scr<-paste(scr,"bed_tocheck_nonccds",sep = " ")
scr2<-paste("Rscript",scr,sep=" ")
syst_scr<-paste(scr2,"bed_tocheck_nonccds",sep = " ")
Ribotaper_ORF_find.R (105-121行目)
# $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores
Rscript $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores

echo "NONCCDS ORF finding..."

# $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores
Rscript $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores

# Groups ORFs and creates BED files + protein fasta database

echo "Grouping ORFs and creating protein fasta database..."

# $scripts_dir"/create_protein_db.R"
Rscript $scripts_dir"/create_protein_db.R"

# makes summary plot for the found ORFs

echo "Summarizing ORF finding results"

# $scripts_dir"/ORF_final_results.R"
Rscript $scripts_dir"/ORF_final_results.R"
create_annotations_files.bash (220行目, 225行目, 235行目)
# $scripts_dir_full"/write_startstops.R"
Rscript $scripts_dir_full"/write_startstops.R"
# $scripts_dir_full"/gtf_to_start_stop_tr.R" 
Rscript $scripts_dir_full"/gtf_to_start_stop_tr.R" 
# $scripts_dir_full"/genes_coor.R"
Rscript $scripts_dir_full"/genes_coor.R"
create_metaplots.bash (68行目)
# $scripts_dir"metag.R" $3
Rscript $scripts_dir"metag.R" $3
Ribotaper.sh (135-153行目, 167-177行目)

今回は、このスクリプトを分割して実行していますが、シェルスクリプトの中身は以下のように修正しています。

# $scripts_dir"/tracks_analysis.R" ccds $scripts_dir $n_of_cores 
Rscript $scripts_dir"/tracks_analysis.R" ccds $scripts_dir $n_of_cores 

echo "Running calculations exons_ccds..."

# $scripts_dir"/tracks_analysis.R" exonsccds $scripts_dir $n_of_cores 
Rscript $scripts_dir"/tracks_analysis.R" exonsccds $scripts_dir $n_of_cores 

echo "Running calculations nonccds..."

# $scripts_dir"/tracks_analysis.R" nonccds $scripts_dir $n_of_cores 
Rscript $scripts_dir"/tracks_analysis.R" nonccds $scripts_dir $n_of_cores 

# annotates the exons relative to ccds regions TO BE ADAPTED, CHECK WHICH FILES THEY NEED.

echo "Annotate exons..."

# $scripts_dir"/annotate_exons.R" $annot_dir $scripts_dir $n_of_cores 
Rscript $scripts_dir"/annotate_exons.R" $annot_dir $scripts_dir $n_of_cores 

echo "Making quality plots..."

# $scripts_dir"/quality_check.R" $annot_dir
Rscript $scripts_dir"/quality_check.R" $annot_dir
# $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores
Rscript $scripts_dir"/CCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores

echo "NONCCDS ORF finding..."

# $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores
Rscript $scripts_dir"/NONCCDS_orf_finder.R" $annot_dir $scripts_dir $bedtools_dir $n_of_cores

# Groups ORFs and creates BED files + protein fasta database

echo "Grouping ORFs and creating protein fasta database..."

# $scripts_dir"/create_protein_db.R"
Rscript $scripts_dir"/create_protein_db.R"

今回実行に使うRiboTaperのシェルスクリプトは以下の通りに分割しておきます。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/Ribo-seq_RiboTaper/scripts

  • Ribo-seq-RiboTaper-3-define_ORF_STEP1.sh
  • Ribo-seq-RiboTaper-3-define_ORF_STEP2-1.sh
  • Ribo-seq-RiboTaper-3-define_ORF_STEP2-2.sh
  • Ribo-seq-RiboTaper-3-define_ORF_STEP2-3.sh
  • Ribo-seq-RiboTaper-3-define_ORF_STEP3.sh
  • Ribo-seq-RiboTaper-3-define_ORF_STEP4.sh

上記が、実行するためのシェルスクリプトになります。これらを、ribotaperscriptsディレクトリにあらかじめコピーしておきます。

実行時に注意すべき点は、RiboTaperはマルチコアで処理させる場合、使用するCPUのコア数によりますが、実行時に50-60GB程度のメモリを確保したほうがよいです(メモリが不足すると途中で処理が止まってしまいます)。

以下では、コマンドライン上で1つずつ実行したときの例を示します。 まず、最初に作成したribotaperの仮想環境に移行します。

source activate ribotaper
STEP1(所要時間: 2.5hr)
Ribotaper_step1.sh \
  ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \
  ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \
  27,28,29 \
  12,12,12 \
  /home/akimitsu/software/RiboTaper_v1.3/scripts \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \
  8
  • 1つ目の引数: Ribo-seqデータのBAMファイルを指定。
  • 2つ目の引数: RNA-seqデータのBAMファイルを指定。
  • 3つ目の引数: 先ほど作成したRiboTaper用のアノテーション情報のディレクトリを指定。
  • 4つ目の引数: 使用するリード長(ノイズでなく、Ribosomeが結合していたリードのみを選抜)を指定(カンマ区切り)。先ほどのRibo-seqデータのチェック時に、緑色のバーに偏ってリードがマッピングされたリード長のみを選抜してくる。
  • 5つ目の引数: P-sitesの計算に使用するOffset値を指定(カンマ区切り)。
  • 6つ目の引数: 使用するRiboTaper用のスクリプトのパスを指定。
  • 7つ目の引数: 使用するBedtoolsのパスを指定。
  • 8つ目の引数: 使用するCPUのコア数。

実行後、以下のファイル群が得られるはずです。

P_sites_all
Centered_RNA
RIBO_unique_counts_ccds
RIBO_best_counts_ccds
RNA_unique_counts_ccds
RNA_best_counts_ccds
RIBO_unique_counts_exonsccds
RIBO_best_counts_exonsccds
RNA_unique_counts_exonsccds
RNA_best_counts_exonsccds
RIBO_unique_counts_nonccds
RIBO_best_counts_nonccds
RNA_unique_counts_nonccds
RNA_best_counts_nonccds

# data_tracks/
Centered_RNA_tracks_ccds
Centered_RNA_tracks_exonsccds
Centered_RNA_tracks_nonccds
P_sites_all_tracks_ccds
P_sites_all_tracks_exonsccds
P_sites_all_tracks_nonccds
Psit_Ribo_Rna_Cent_tracks_ccds
Psit_Ribo_Rna_Cent_tracks_exonsccds
Psit_Ribo_Rna_Cent_tracks_nonccds
RIBO_tracks_ccds
RIBO_tracks_exonsccds
RIBO_tracks_nonccds
RNA_tracks_ccds
RNA_tracks_exonsccds
RNA_tracks_nonccds
index_tracks_ccds
index_tracks_exonsccds
index_tracks_nonccds
STEP2-1, 2-2, 2-3(所要時間: 2.5hr)

STEP2以降も同様に実行するだけです。STEP2-1, 2-2, 2-3は、並行して実行して構いません。Sun Grid Engineなどの並列計算システムを利用している場合、計算時間を短縮するために、並行してジョブを実行するようにしましょう。

以下には、STEP2-1の実行例を示しました。

Ribotaper_step2-1.sh \
  ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \
  ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \
  27,28,29 \
  12,12,12 \
  /home/akimitsu/software/RiboTaper_v1.3/scripts \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \
  8

実行後、以下のファイル群が得られるはずです。

results_ccds
results_exonsccds
results_nonccds
STEP3(所要時間: 40min)
Ribotaper_step3.sh \
  ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \
  ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \
  27,28,29 \
  12,12,12 \
  /home/akimitsu/software/RiboTaper_v1.3/scripts \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \
  8

実行後、以下のファイル群が得られるはずです。

results_nonccds_annot
all_calculations_ccdsgenes_annot_new
quality_check_plots.pdf
STEP4(所要時間: 8-9hr)
Ribotaper_step4.sh \
  ./STAR_output_Ribo-Seq_HEK293_Cell_Control_SRR1630831/Ribo-Seq_HEK293_Cell_Control_SRR1630831_4_STAR_result_Aligned.sortedByCoord.out.bam \
  ./STAR_output_RNA-Seq_HEK293_Cell_SRR1630838/RNA-Seq_HEK293_Cell_SRR1630838_4_STAR_result_Aligned.sortedByCoord.out.bam \
  /home/akimitsu/software/RiboTaper_v1.3/annotation/hg38_Gencode_v24 \
  27,28,29 \
  12,12,12 \
  /home/akimitsu/software/RiboTaper_v1.3/scripts \
  /home/akimitsu/software/RiboTaper_v1.3/bedtools_dir \
  8

実行後、以下のファイル群が得られるはずです。

# ORFs_CCDS/
# ORFs_NONCCDS/
  • ORF_max
  • ORF_max_filt: filtered for multimapping and non-coding nonccds_coding_ORFs (see below) (recommended for further analysis)
  • protein_db_max.fasta: (for proteomics search, not filtered for multimapping)
  • translated_ORFs_sorted.bed: corresponding to ORF_max
  • translated_ORFs_filtered_sorted.bed: corresponding to ORF_max_filt filtered for multimapping and non-coding nonccds_coding_ORFs (see below) (recommended for further analysis)
  • ORFs_genes_found: tab-separated values for number of ORFs found and their corresponding genes
  • Final_ORF_results.pdf: pdf file with length/coverage distributions for different ORF categories.

太文字で書かれたファイルが予測されたORFになります。ファイルの中身を確認すると以下の通り。

f:id:biodata:20170331214311p:plain

詳細については言及しませんが、Categoryを確認すると、主に、

  • ORFs_ccds
  • uORF
  • ncORFS

があります。

ORFs_ccdsはmRNA由来のCanonical ORFを、uORFはmRNAの上流5'UTR中に存在するUpstream ORFを、ncORFsはNoncoding RNA中に予測されるORFを指しています。

PAR-CLIP解析データ(Peak calling, Motif discovery, Peak annotation)

概要

CLIP-seqとは、RNA結合タンパク質(RBP)の結合サイトを網羅的に同定する手法です。CLIP-seqには主に、HITS-CLIP、iCLIP、PAR-CLIP、eCLIPの4種類が存在します。今回は、中でも幅広く利用されている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>

改行コードの問題

WindowsLinuxで改行コードが異なるため、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.refhumrep.refhumsub.refをマージします。ファイル名はRepBase_human_v22_01.faとしておきます。

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

次に、リボソームRNAの配列を用意します。 NCBIの登録されているX12811.1| Human 5S DNAU13369.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 stringMD 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ファイルの1行目の情報を見ています。

  • 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に対して1つの代表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: 1つの配列中に予測されたモチーフ配列が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つ目の引数: ピークの位置情報を加えた遺伝子リストを出力。

以上になります。