Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

pyenvでPython環境を構築する

ざっくりですが。

Pythonツールのインストール

最近はPythonで書かれたツール類が増えてきています。 pip や conda でツールをインストールできるのはエラーも少なく、大変便利ですね。
ただ、pythonは2系と3系があり、python2.xでしか動かないツール、python3.xでしか動かないツールがあります。 違うバージョンで動かそうとするとエラーになったりして、ローカルPCなら適宜使い分ければいいのですが、 共有サーバ上で使い分けたりする場合はなかなかしんどいです。 (今そういう状態になっている)

ここでは、python2系と3系を使い分ける方法を簡単に記載します。
ここに書かれている内容は以下のエントリをベースにしていますので、 anaconda, pyenvとは何ぞや?という人はまず以下のエントリを参照ください。

データサイエンティストを目指す人のpython環境構築 2016

Win, MACのローカルPCの場合

Anacondaを直接インストールするのが一番簡単だと思います。
Anacondaも2系と3系がありますが、3系をインストールし、必要があれば中で2系の環境を構築するのがよいでしょう。
(構築方法は上記エントリを参照)

pyenvで環境構築

Linuxのターミナル上で動かす場合は、pyenvが一番良いのではないかと今のところ思っています。
システムにインストールされているpythonを色々いじるのはリスクが大きいことと、 pyenv内であれば、よくわからない状態になってしまった場合に新規入れ直しが簡単なためです。
anaconda3の中でpy2環境を構築するのでも良いのですが、一部のいけてないツールはpy2環境のbinを読みに行ってくれなかったりするので、pyenvで切り替える方がエラーが少なく済むかもしれません。
ただ、pyenvの環境設定は基本的にユーザ単位のようなので、共有サーバ上で使う場合、各ユーザが個別にpyenvをインストールし、設定作業をする必要があるようです。

pyenvのインストール

git がインストールされている前提です。

 $ git clone https://github.com/yyuu/pyenv.git ~/.pyenv
 $ echo 'export PYENV_ROOT="$HOME/.pyenv"' >> ~/.bashrc
 $ echo 'export PATH="$PYENV_ROOT/bin:$PATH"' >> ~/.bashrc
 $ echo 'eval "$(pyenv init -)"' >> ~/.bashrc
 $ source ~/.bashrc

pyenvでanacondaのインストール

ここではpyenv localを使う前提でanacondaの2系と3系をインストールしています。内部でpy2 環境を構築する場合はanaconda 3のみインストールします。

 $ pyenv install anaconda2-4.4.0
 $ pyenv install anaconda3-4.4.0
 $ pyenv rehash
 $ pyenv global anaconda3-4.4.0

ツール類のインストール

何しろnumpyとscipyは必ず必要になりますので最初に入れておきましょう。

 $ conda install -y numpy scipy

例えばDeepToolsをインストールしたい場合は以下のようにします。

 $ conda install -c bioconda deeptools

biocondaには色々とツールが入っているのですが、バージョンが古かったりエラーがあったりしますので、注意してください。

Python 2.x でしか動かないもの

  • MACS2 及び 内部でMACS2を動かすツール(Mangoなど)
  • Segway 2.0
  • Hi-C pro
  • RSeQC
  • qiime

など

Python 3.x で動くもの

  • DeepTools
  • HiClib
  • mageck-vispr
  • Snakemake

など

ENA,DDBJからfastqを取得する(2019/9/6追記)

前回のエントリでは、SRAからfastqを取得する方法を紹介しました。

SRAからfastqを取得する (2019/07/16 追記あり) - Palmsonntagmorgen

一方、ENA (European Nucleotide Archive) やDDBJから直接fastqファイルをダウンロードすることも可能です。

ENA Browser
http://trace.ddbj.nig.ac.jp/dra/index_e.html

SRA, DDBJ, ENAのデータベースは互いに同期されており、基本的にどのサイトからでも同一のファイルを取得できます。
2017/7/13の時点では、ENAからは fastq.gz 形式、DDBJからは fastq.bz2 形式でダウンロード可能です。
SRAを使うと、sraのダウンロード -> fastq への変換 -> fastq.gz への圧縮と三段階必要であるのに比べ、 上記サイトを使えば fastq.gz のダウンロードのみになるので、簡便かつ高速です。 なお私的には.bz2形式よりも.gz形式の方が使いやすいので、ここではENAからのダウンロードを紹介します。
DDBJからのダウンロードについては下記のサイトに詳述されていますので、そちらを参考にしてください(適当)。

https://bi.biopapyrus.jp/rnaseq/fastq-data/download-fastq-ddbj.html

(2019/9/6追記:DDBJは2017年4月7日からSRAのミラーリングを停止しているそうです。SRAのデータのダウンロードにはSRA本体かENAを使ってください。
 参考:https://www.ddbj.nig.ac.jp/news/ja/180918_2.html
(2021/5/7追記:2020年からDDBJのミラーリングは再開されたそうです。)

ENAからデータの取得

ここではサンプルデータとして SRR390728 を使います。
ENAのサイト上でSRR390728を検索すると、下記のような画面になります。

f:id:rnakato:20170713151016p:plain

"FASTQ files (FTP)" のリンクをクリックするとダウンロードが始まります。 SRR390728 は paired-end なので、2ファイルに分かれています。 singleとpairedを間違うこともないので、簡単ですね。
Linux ターミナル上ではwgetを使ってダウンロードをすることも可能です。 上記ページからFTPのリンクをコピーし、以下のようにコマンドを打ちます。

 $ wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR390/SRR390728/SRR390728_1.fastq.gz
 $ wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR390/SRR390728/SRR390728_2.fastq.gz

wget を使ったファイルの一括取得

サンプル数が多い場合にひとつひとつ手作業でダウンロードするのは面倒です。 そのような時はFTPリンクをファイルに保存し、wget -i を使うと簡単です。

ここではサンプルデータとして SRP000712 を使います。
ENA上で SRP000712 を検索すると、以下のように25サンプル含まれていることがわかります。

f:id:rnakato:20170713152534p:plain

画面上部の "TEXT" をクリックすると、各サンプルの詳細をテキスト形式でダウンロードできます。 この例では "PRJNA106023.txt" というファイル名で保存されます。
FTPのURLは10列目に記載されているので、10列目だけ抽出したファイルを作成しましょう。

 $ cut -f10 PRJNA106023.txt > ftp.txt

あとはこのファイルにwget -i を実行すれば順番にダウンロードされます。

 $ wget -i ftp.txt

サンプルによっては1サンプルにつき複数のURLが ; 区切りで記載されていたりしますが、 ; を改行に変換するなどして対応してください。

SRAからfastqを取得する (2019/07/16 追記あり)

更新の間が随分空いてしまいました。 その間に2つの学会に参加してきたのですが、海外の解析手法の進化具合にずいぶん衝撃を受けました。
が、ここは予定通り初歩的な作業から説明していきたいと思います。

今日はSRA(Sequence Read Archive) からfastqファイルを取得する方法です。

SRA Toolkit

まずはNCBIのサイトから SRA Toolkit をダウンロードします。

SRA Toolkit download

MAC用、Win用、Linux用と分かれているので、適切なものを選んでください。
例えばubuntu64なら、ダウンロードした圧縮ファイル(sratoolkit.2.8.2-1-ubuntu64.tar.gz)を解凍後、sratoolkit.2.8.2-1-ubuntu64/bin ディレクトリにパスを通せばOKです。
SRA Toolkitはかなり頻繁にアップデートされるので、ちょくちょく最新版をチェックしておくことをおすすめします。
2017/7/5時点での最新版は 2.8.2-1 です。

fastq-dumpでfastqをダウンロード (single-end)

ここではサンプルデータとして SRR390728 を使います。

 $ fastq-dump SRR390728

single-end データの場合、基本的にはこれだけでOKです。 Hiseqデータならば.fastq形式で、SOLiDでシーケンスされたcolor spaceデータならばcsfastq形式で出力されます。
圧縮形式(.fastq.gz)で出力したければ --gzip オプションを使います。

 $ fastq-dump SRR390728 --gzip

その他のオプションが知りたい場合は--help オプションを使用してください。

fastq-dumpでfastqをダウンロード (paired-end)

paired-endの場合、デフォルトで変換するとF3, F5リードがつながったSRR390728.fastq ファイルが出力されてしまいます。
F3, F5リードを分割して出力するには --split-files オプションを追加します。

 $ fastq-dump SRR390728 --split-files

これで、SRR390728_1.fastq, SRR390728_2.fastq と分割されたfastqファイルが生成されます。
なお、single-end データに --split-files オプションを使った場合、XXX_1.fastq ファイルのみ生成され、内容はオプションなしの場合と同一です。 ので、single-end とpaired-endデータが混在していて面倒という場合には、全てのファイルに --split-files オプションを付加しても問題ありません。

fastq-dumpで.sraから.fastqに変換

上記コマンドを実行した場合、fastq-dumpコマンドは内部で

(1) .sraファイルのダウンロード(~/ncbi/public/sraに保存される)
(2) .sraファイルの.fastqへの変換

を行っています。

fastq-dumpでダウンロードされた.sraファイルは.fastq変換完了後も削除されませんので、 何らかの理由で.fastqファイルへの変換に失敗した(paired-endなのに--split-files オプションをつけ忘れたなど)場合には、 ダウンロード済の.sraから.fastqの変換を行うと二度手間になりません。

 $ fastq-dump ~/ncbi/public/sra/SRR390728.sra

こうすると、 SRR390728.sra から SRR390728.fastq への変換のみが行われます。
逆に、.sraファイルをいつまでも保存していると容量を圧迫しますので、必要なくなった.sraファイルはまとめて削除しましょう。

(2019/07/16 追記) version 2.9.1 から、fasterq-dumpというコマンドが登場しています。fastq-dumpの後継で、マルチスレッドにすることで高速化しているようです。fastq-dumpは将来的に廃止されるそうですので、今後はこちらを使うとよいでしょう。なお--gzip オプションは消えているようですので、生成したfastqを pigzなどで圧縮するようにしましょう。

prefetchでsraをダウンロード

fastqファイルの取得・変換には時間がかかります。
大量のデータを一度に取得する場合には、 安全を期してまず.sraファイルを取得したいという場合もあるでしょう。
prefetchコマンドを使うと、sraファイルのダウンロードのみを行うことができます。

 $ prefetch SRR390728

このコマンドを実行すると、SRR390728.sra ファイルが ~/ncbi/public/sraに保存されます(出力先ディレクトリを何故指定できないのかは謎です…)。

--option-file コマンドを使うと、一度に複数の.sraファイルをダウンロードすることもできます。

 $ prefetch --option-file SRR_Acc_List.txt 

こうすると、SRR_Acc_List.txt に記載されているサンプル全てをダウンロードします。
これらのオプションはツールのバージョンアップと共に変わっていく可能性がありますので、新しいバージョンを使う場合には動作確認を忘れずに。

2bit genome を作成する

2bit genome はゲノム配列ファイルを2bit (バイナリ)形式で格納したものです。
2bit 形式はテキストエディタで開くことはできませんが、multifasta 形式よりも非常に高速にプログラムに読みこむことができるため、
ゲノム解析ツールを使う際にまれに 2bit genomeを要求されることがあります。
初めて見ると戸惑うかもしれませんが、UCSC genome browserで提供されているので、それをダウンロードして使えばOKです。
例えばhg38であれば http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ の hg38.2bit がそれに該当します。

2bit genome の作成

この2bit genomeファイルはUCSCから提供されているプログラムで作成することも可能です。
UCSC genome browserからダウンロードした ゲノム配列データにはコンティグ配列も含まれていますので(下記エントリ参照)、
それらを取り除いた2bitファイルが欲しいというような場合は自分で作成するとよいでしょう。

rnakato.hatenablog.jp

twoBitToFa と faToTwoBit を使った相互変換

2bit 形式 と multifasta 形式の相互変換には twoBitToFa と faToTwoBit のふたつのプログラムを使います。
まずはこの2つのプログラムをダウンロードし、実行権限を付与します。

 $ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
 $ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
 $ chmod +x faToTwoBit twoBitToFa  # 実行権限の付与

UCSC genome browserからhg38 の 2bit genomeをダウンロードしてみましょう。

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit
 $ less hg38.2bit   # lessコマンドでオープンを試みると
"hg38.2bit" may be a binary file.  See it anyway? 

hg38.2bitはバイナリファイルなので開けませんよ、と警告が出ます。

2bit genome → multifasta の変換

2bit genome → multifasta の変換には twoBitToFa を使います。

 $ ./twoBitToFa hg38.2bit hg38.fa
 $ less hg38.fa
>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
...

multifastaはテキストファイルなので、lessでオープンすることができます。
なお、ここでNが並んでいるのは、chr1の末端はテロメア領域であり配列が解読されていないためです。 このような未解読の領域は「ギャップ領域」と呼ばれます。
ヒトゲノム解読が完了したとは言っても、このようなギャップ領域はまだたくさん残っているということですね。

multifasta → 2bit genome の変換

こちらは faToTwoBit を使います。

 $ ./faToTwoBit hg38.fa hg38.2bit

hg38.faとhg38.2bitはファイル形式が異なるだけで内容は同一ですので、
このように、何度でも相互変換をすることができます。

なお、http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ に格納されている UCSCのコマンドツール群には他にも便利なものが色々ありますので、
全部ダウンロードして利用可能にしておくのがオススメです。

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する

UCSC genome browserからダウンロードした ゲノム配列データにはコンティグ配列なども含まれていますが、これらは通常ゲノムの解析には用いません。
そこでこれらを除去し、常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成します。
ここではhg38を用います。異なるbuildの場合は適宜読み替えてください。

hg38.fa.gz のダウンロード

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

生成された hg38.fa に含まれる配列を確認します。

 $ grep \> hg38.fa 
>chr1
>chr10
>chr11
>chr11_KI270721v1_random
>chr12
>chr13
>chr14
>chr14_GL000009v2_random
>chr14_GL000225v1_random
>chr14_KI270722v1_random
>chr14_GL000194v1_random
...
>chrUn_GL000218v1
>chrX
>chrY
>chrY_KI270740v1_random

コンティグ配列が含まれています。

染色体配列に分割

multifastaファイルをシングルfastaファイルに分割して出力する splitmultifasta.pl を使って、染色体配列ファイルを生成します。
github:DROMPA3のscripts フォルダからsplitmultifasta.pl をダウンロードしてください。
下記の git clone コマンドを実行してDROMPA3ごとダウンロードすることもできます:

 $ git clone https://github.com/rnakato/DROMPA3.git

DROMPA3/script/ ディレクトリの中に splitmultifasta.pl が入っています。
hg38.fa から染色体配列ファイルを生成します。

 $ splitmultifasta.pl hg38.fa # 分割
 $ ls *.fa                    # 確認
chr1.fa                    
chr10.fa                  
chr10_GL383545v1_alt.fa    
chr10_GL383546v1_alt.fa
chr10_KI270824v1_alt.fa  
...

染色体ごとのfastaファイルが生成されました。
生成されるファイル名は元ファイルのヘッダ名になります。

染色体ファイルの結合

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成します。

 $ s=""
 $ for i in $(seq 1 22) X Y  # 1~22番染色体とXとYを指定。ミトコンドリア(M)を追加してもよい
 $ do
 $   s="$s chr${i}.fa"
 $ done
 $ cat $s > genome.fa

生成された genome.fa を確認します。

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

無事目的のファイルを作成できました。

genome tableを作成する

genome table はゲノム中に存在する各染色体の名前とその長さをタブ区切りで記述したファイルで、DROMPA や bedtools などの解析ツールを使う時に必要になります。
UCSC genome browserの *.chrom.sizes ファイルをダウンロードしてもいいのですが、自分で自由に作成できる方が便利なので、makegenometable.plスクリプトを作成しました。

Download

github:DROMPA3のscripts フォルダからmakegenometable.pl をダウンロードしてください。
下記の git clone コマンドを実行してDROMPA3ごとダウンロードすることもできます:

 $ git clone https://github.com/rnakato/DROMPA3.git

DROMPA3/script/ ディレクトリの中に makegenometable.pl が入っています。

実行

makegenometable.pl の入力となるのは全染色体を含んだゲノム配列データ(multifasta形式)です。
例としてhuman genome build hg38のgenome tableを作成してみます。
まずUCSC genome browserのhg38.fa.gz をダウンロード・解凍します。

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

生成された hg38.fa に含まれる配列を確認します。

 $ grep \> hg38.fa 
>chr1
>chr10
>chr11
>chr11_KI270721v1_random
>chr12
>chr13
>chr14
>chr14_GL000009v2_random
>chr14_GL000225v1_random
>chr14_KI270722v1_random
>chr14_GL000194v1_random
...
>chrUn_GL000218v1
>chrX
>chrY
>chrY_KI270740v1_random

細かいコンティグなども含まれていますが、今回はこれをそのまま利用します。
作成されたhg38.fa に対してmakegenometable.pl を実行し、出力先をgenometable.txtに指定します。

 $ makegenometable.pl hg38.fa > genometable.txt

catでgenometable.txtを表示して以下のように表示されれば成功です。

 $ cat genometable.txt
chr1    248956422
chr10   133797422
chr11   135086622
chr11_KI270721v1_random 100316
chr12   133275309
chr13   114364328
chr14   107043718
chr14_GL000009v2_random 201709
chr14_GL000225v1_random 211173
chr14_KI270722v1_random 194050
chr14_GL000194v1_random 191469
...
chrUn_GL000218v1    161147
chrX    156040895
chrY    57227415
chrY_KI270740v1_random  37240

makegenometable.plは遺伝子配列長の計測にも利用可能です。