Palmsonntagmorgen

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

常染色体と性染色体のみのゲノム配列ファイル 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は遺伝子配列長の計測にも利用可能です。