Palmsonntagmorgen

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

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は遺伝子配列長の計測にも利用可能です。