アメリエフのブログ

Amelieff Staff Blog

Pythonで全ゲノムFASTA=>染色体毎FASTAを作成してみた

こんにちは、iijm_lです🐰
12月になりました、今年ももうすぐ終わっちゃいますね🥺
来年もいい年になりますように...🍊(もう既に年末気分)

f:id:iijm-l:20201130165906p:plain


さて今回紹介するのは、

ゲノム配列ファイルから染色体毎の配列ファイルを作成する方法

です🌟


稀に、全ゲノム配列FASTAだけでなく、染色体ごとのFASTAが必要な解析ツールってありますよね〜。
先日気軽な感じで、
🐰「ツールの機能知りたいから、とりあえずテスト解析してみよ〜🥳」
と思って実行してみたら、
💻「染色体ごとのFASTAが無いよ、用意して😡
と怒られてしまいました。

全ゲノムのFASTAならある、、、そんな時はpythonですぐに解決できます💫

動作確認環境

  • python 3.7.9
  • Biopython 1.78

参考スクリプト

BiopythonのSeqIOを使って、FASTAをリードごとに分割し出力させます。

#!/bin/env python

import os
from Bio import SeqIO as SIO

genomefile = '/home/genome/hg38/genome.fa'
outdir = '/home/genome/hg38/byChr/'

# make outdir
try:
    os.makedirs(outdir)
except FileExistsError:
    pass

# separate by chrom number
for record in SIO.parse(genomefile, 'fasta'):
    if not record.id.startswith('chr'):
        continue
    outfile = os.path.join(outdir, record.id + '.fa')
    SIO.write(record, outfile, 'fasta')


変数genomefileoutdirに指定するファイルパスは環境に合わせて指定してくださいね😉
ついでに出力ディレクトリ作成機能も付けてみました。
書き出し権限があるディレクトリを指定しないと、エラーを吐くので注意が必要です。


いざ、実行!!!

$ python sepbychr.py 


1分ほどかかって正常終了しました⏱
果たしてファイルはできているのでしょうか...?

$ ls /home/genome/hg38/byChr/
chr1.fa   chr13.fa  chr17.fa  chr20.fa  chr4.fa  chr8.fa   chrY.fa
chr10.fa  chr14.fa  chr18.fa  chr21.fa  chr5.fa  chr9.fa
chr11.fa  chr15.fa  chr19.fa  chr22.fa  chr6.fa  chrMT.fa
chr12.fa  chr16.fa  chr2.fa   chr3.fa   chr7.fa  chrX.fa


成功しました👏🏻
※染色体の種類は、使用したリファレンスゲノムに依存しますので、上記と異なる場合がございます。

まとめ


BiopythonのSeqIOを使えば、FASTA形式のデータが簡単にparseできます✨
FASTQ形式やその他いろいろな形式がparse/outputできるので、非常に便利ですよ!