Home
bash 応用, FASTX-Toolkit, シェルスクリプト
シェルスクリプト FASTAファイルの分割 fastx_toolkit, sed, コマンドの活用例

シェルスクリプト FASTAファイルの分割 fastx_toolkit, sed, コマンドの活用例

バイオインフォ道場、くまぞうです。

バイオインフォマティクスでは、大きなファイルを扱う機会が多いです。小さなファイルに分割すれば、サーバーで分散処理して大幅な時間短縮を行うことが出来るかもしれません。いくつかのツールを組合せてFASTAファイルを小さく分割してみます。

スポンサーリンク

FASTAファイルについて

FASTAファイルは、シーケンスデータの記述形式として良く使われるファイルフォーマットです。エントリー（ヘッダー行）と対応するシーケンス文字列の繰り返しで構成されています。ヘッダー行の先頭は>です。続いて、シーケンスを識別するための情報を記述します。これらの間に空白を入れてはいけません。また、塩基配列は、長すぎる場合は適当に改行されます。シーケンスデータの各行は８０文字未満が推奨されています。

FASTAファイルの例（Arabidopsisから一部抜粋）

>AT1G51370.2 | Symbols:  | F-box/...
ATGGTGGGTGGCAAGAAGAAAACCAAGATATGTGACAAAGTGTCACATGAGGAAGATAGG...
TTTGATATCTGAAATACTTTTTCATCTTTCTACCAAGGACTCTGTCAGAACAAGCGCTTT...
TTTGGCAATCGGTTCCTGGATTGGACTTAGACCCCTACGCATCCTCAAATACCAATACAA...
...

FASTAファイルの分割

FASTAファイルはシーケンス文字列が何行も続くので、単純に行数で分割するのは難しいです。分割後のファイルに、ヘッダー行と対応するシーケンス文字列がきちんとおさまるように分割したいと思います。まず、FASTAファイルを「ヘッダーとシーケンス文字列」の１行にします。次に、分割数で行数を計算し、Linuxコマンドのsplitで分割します。最後に、フォーマットを調整を行います。

FASTX-Toolkitの活用

FASTX-Toolkitは、配列情報に関する操作を提供している便利なツールです。今回はFASTAファイルを「ヘッダーとシーケンス文字列」の１行（TAB区切り）にまとめる作業と整形に用います。まずは、１行にまとめるために、fasta_formatterコマンドを使います。オプションは-tです。注意点は、この作業を通じて、ヘッダー行の>が消えることです。

ヘッダーとシーケンス文字列の１行化

$ cat my_fasta.fasta | fasta_formatter -t > my_fasta_line.fasta

$ head my_fasta_line.fasta
AT1G51370.2 | Symbols:  | F-box/...ATGGTGGGTGGCAAGAAGAA
...

Linux splitコマンドの活用

FASTAファイルが「ヘッダーとシーケンス文字列」の１行になったので、Linuxコマンドのsplitを使って行数で分けます。行数で分ける場合のオプションは-lです。デフォルトの設定では、行数で分けられたファイルはxaa, xab, xac….というファイル名で保存されます。オプションで変更することが可能です。

ファイルを行数で分割

$ wc -l < my_fasta_line.fasta
300
$ cat my_fasta_line.fasta | split -l 100  #100行ずつに分割

$ ls
xaa xab xac
$wc -l xa*
100 xaa
100 xab
100 xac

sed の活用

FASTAのヘッダー行の先頭に>を挿入し、ヘッダー行とシーケンス文字列を改行で分けます。今回は、「数字が含まれる行」をヘッダー行と仮定して処理しました。

フォーマットの手直し

cat xaa | sed -e '/[0-9]/s/.*/>&/' -e 's/\t/\n/' > my_fasta1.fasta

FASTX-Toolkitで仕上げる

最後に見た目を整えます。FASTX-Toolkitのfasta_formatterでシーケンス文字列の折り返し長さを指定して、きれいに整えます。オプションは-wです。折り返す文字数とともに指定します。

cat my_fasta1.fasta | fasta_formatter -w 60 > my_fasta1_60.fasta  # シーケンス文字列を60文字で折り返し

FASTAファイルの分割を自動化する

ここまでの処理を、シェルスクリプトで自動化してみました。FASTAファイルと分割数（デフォルト４）を引数にしてファイルを分けます。処理の都合上、オプションの追加・変更を行っています。

fasta_split.sh

#! /bin/bash

FILE=$1  # ファイル名（引数）
DIVN=$2  # 分割数（引数）

# 分割数のデフォルト値設定
DIVN=${DIVN:-4}

# ファイルチェック
if [ ! -e ./$1 ]; then
    echo "Not Found:$1"
    exit 1
fi

echo "$1 -> ${DIVN}分割"

# ヘッダーと配列情報を１行にまとめる
cat ${FILE} | fasta_formatter -t > ${FILE}_tab

# 分割数からファイル行数を計算
LINE=`wc -l < ${FILE}_tab`
NUM=`expr ${LINE} / ${DIVN}`
CAL=`expr ${NUM} \* ${DIVN}`

if [ $LINE -gt $CAL ]; then
    NUM=`expr $NUM + 1`
fi  

# 分割
cat ${FILE}_tab | split -l $NUM --numeric-suffixes=1
rm -f ${FILE}_tab

# 分割したファイルを整形（>の挿入・改行・幅調整）
for x in `ls x[0-9]*[0-9]`; do
    ID=${x:1}
    NAME=${1/\./$ID.}.pre
    cat $x | sed -e "/[0-9]/s/.*/>&/" -e 's/\t/\n/' > ${NAME}
    cat ${NAME} | fasta_formatter -w 60 > ${NAME:0:-4}
    echo ${NAME:0:-4}
    rm -f $x
    rm -f $NAME
done

$ chmod +x fasta_split.sh
$ ./fasta_split.sh myfasta.fa 4
myfasta.fa -> 4分割
myfasta01.fa
myfasta02.fa
myfasta03.fa
myfasta04.fa

$ ls
myfasta.fa 
myfasta01.fa
myfasta02.fa
myfasta03.fa
myfasta04.fa

スポンサーリンク