Tips : [Python] fastaファイルから配列を無作為復元抽出する

　大量の配列から無作為抽出したくなるときは時々あります。
　今回はそんなときに使うコードを紹介します。

　やり方は単純にファイルを開き、配列を乱数で選んでアクセスしていきます。

配列を乱数で選んでアクセス

# 配列取得
seqnum=random.randint(1,seqnumMAX)
name = linecache.getline(infn, seqnum*2-1) #配列名
seq = linecache.getline(infn, seqnum*2)　 #配列
# 書き込み
outfdl    =   open(outfn, 'w')
outfdl.writelines(name+seq)
outfdl.close()
linecache.clearcache() #キャッシュをクリア

　このとき読み込みにはlinecache.getline（入力ファイル名, 行数）を使います。
　これを使うと指定した行の読み込みを内部で最適化してくれます。
　Bioinformaticsでは大きな配列ファイルを扱いがちなので便利かもしれません。

　この処理を取得したい回数分だけループしていけば無作為復元抽出ができます。

　以上をちゃんとした処理としてまとめると以下のようになります。

randsampleFasta.py

# !/usr/local/bin/python3
# -*- coding: utf-8 -*-

"""
fastaファイルから配列を無作為復元抽出する。
"""
__author__  = "Kazuki Nakamae <[e-mail]>"
__version__ = "0.00"
__date__    = "29 May 2017"

import sys
import mmap
import linecache
import random

def randsampleFasta(infn,outfn,n):
    """
    fastaファイルから配列を無作為復元抽出する。
    @param  infn 読み込むfastaファイル名{string} 
    @param  outfn 出力するfastaファイル名{string} 
    @param  n 復元抽出する配列数{int} 
    """

    infdl     =   open(infn, 'r')
    #入力ファイルの配列数確認
    print('Checking how many sequences are in a FASTA format sequence file....(1/2)')
    buf = mmap.mmap(infdl.fileno(), 0, prot=mmap.PROT_READ)
    seqnumMAX= 0
    readline = buf.readline
    while readline():
        seqnumMAX += 1
    seqnumMAX = int(seqnumMAX / 2)
    infdl.close()
    print(str(seqnumMAX)+' sequences')
    #配列を無作為復元抽出
    print('random sampling from a sequence with replacement....(2/2)')
    outfdl    =   open(outfn, 'w')
    random.seed(a='hoge', version=2)    #シード設定
    seqi=1
    while seqi<=n:
        seqnum=random.randint(1,seqnumMAX)
        name = linecache.getline(infn, seqnum*2-1)
        seq = linecache.getline(infn, seqnum*2)
        outfdl.writelines(name+seq)
        seqi += 1
    outfdl.close()
    linecache.clearcache() 
    print('done.')

if __name__ == '__main__':
    argvs = sys.argv  # コマンドライン引数
    argc = len(argvs) # 引数の数

    if (argc != 4):   # 引数チェック
        print("USAGE : python3 randsampleFasta.py <INPUT.fa> <OUTPUT.fa> <NUMBER_OF_SEQUENCES>")
        quit()

    randsampleFasta(argvs[1],argvs[2],int(argvs[3]))
quit()

実施例

　入力するファイル

test.fa

>1
CCGTATTGGAAAGCTC
>2
AGGATTATCGGATACT
>3
ATCCGGACGGGGGGTT
>4
GACCTCGTTATCATCC
>5
AGTCAGGTTACCCGCA

　Bash上での入力

入力

python3 randsampleFasta.py test.fa out.fa 4

　Bash上での出力

標準出力

Checking how many sequences are in a FASTA format sequence file....(1/2)
5 sequences
random sampling from a sequence with replacement....(2/2)
done.

　出力ファイル

out.fa

>3
ATCCGGACGGGGGGTT
>2
AGGATTATCGGATACT
>4
GACCTCGTTATCATCC
>3
ATCCGGACGGGGGGTT

コードの利用

自由です。