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

スポンサーリンク

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

とにかく早く問題解決したい人はこちら>>直接、データ解析相談

FASTAファイルについて

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

FASTAファイルの例(Arabidopsisから一部抜粋)

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

FASTAファイルの分割

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

FASTX-Toolkitの活用

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

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

$ cat my_fasta.fasta | fasta_formatter -t > my_fasta_line.fasta
$ head my_fasta_line.fasta
AT1G51370.2 | Symbols:  | F-box/...<tab区切り>ATGGTGGGTGGCAAGAAGAA
...
</tab区切り>

Linux splitコマンドの活用

FASTAファイルが「ヘッダーとシーケンス文字列」の1行になったので、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-Toolkitfasta_formatterでシーケンス文字列の折り返し長さを指定して、きれいに整えます。オプションは-wです。折り返す文字数とともに指定します。

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

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

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

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}分割"

# ヘッダーと配列情報を1行にまとめる
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

論文・データ解析の支援

個別に直接相談できるのが良いと、たくさんの方にご利用頂いています。

とにかく早く問題解決したい人はこちら>>直接、データ解析相談
とにかく早く問題解決したい人はこちら>>直接、論文相談

関連記事

IGV 使い方 遺伝子位置や好きな位置への移動
IGVはリードのマッピング状況を視覚的に確認できるツールです。ゲノムデータ(fastaファイル)とアライメントデータ(bamファイル、及びbaiファイル)を準備すれば簡単に表示できますが、更に遺伝子の位置情報を読み込むともっと使いやすくなり...
IGV 使い方 インストール〜便利な使い方まで | リファレンス・マッピングデータ・アノテーションを読み込んで表示しよう
IGV(Integrative Genomics Viewer)は、マッピングデータを視覚的に確認することができる大変便利なツールです。とにかく早く問題解決したい人はこちら>>直接、データ解析相談IGVとはIGVは、米Broad Insti...
IGV 使い方 bowtieのマッピングデータをsamtoolsで変換してIGVに表示するまで
IGVを使って、bowtieのマッピング結果表示してみます。samtoolsを使ったbamへの変換やsortなど、IGVの表示に必要なステップを1つずつ確認しながら説明します。とにかく早く問題解決したい人はこちら>>直接、データ解析相談bo...
マッピング率が低いときの対応方法 - Trimmomatic 除去アダプターの登録方法
マッピング率が極端に低い場合は、使用しているシーケンサ特有のアダプター配列がリードに残っていないか確認するとうまくいくケースがあります。特に、ライブラリーサイズが小さくなった場合は、シーケンサ内部のソフトで特有の配列が除去しきれずに残るケー...
fastqc memory error - java実行時により大きなメモリを割り当てる方法
bamファイルに対してfastqcを実行したら、メモリが足りないというエラーが発生しました。実行環境は十分にメモリーを積んでいたので、fastqc実行時のjavaのメモリーオプションを強制的に書き換えて実行したら、正常終了しました。とにかく...
FastQC 使い方 Fastq リードの評価
バイオインフォ道場、くまぞうです。バイオインフォマティクスの解析を行うにあたって、次世代シーケンサから得られたリード情報の品質をチェックすることはとても大事なことです。代表的なツールFastQCの使い方をまとめます。とにかく早く問題解決した...
FastQC 使い方 コマンドラインでテキスト形式の結果を活用する
バイオインフォマティクスの解析を行う際、一番最初に行うのはリードの品質チェックです。代表的なツールはFastQCで結果を視覚的に確認できて非常に便利ですが、たくさんデータを使う場合はコマンドライン(command line)を活用して結果を...
FastQC - インストール
FastQCは、NGSのリードの品質をチェックするツールです。解析の直前で使用し、リードの品質に問題がなければ解析に進みます。もし、品質に問題があるようであれば、低品質の情報をトリミングツールで除去したり、対象となるデータを解析から外すなど...
blast 使い方 完全一致の連続領域を探索するには?
blastn で完全一致の連続領域を探索する方法blastn オプション word_sizeblastnのパラメータ-word_sizeは、完全一致の初期長さを設定します。長さを指定すると、「完全一致が連続する場所」で「指定長さを満たす」部...
blast 使い方 エラーと解決 tblastx
バイオインフォマティクスの解析では、様々なツールを使います。たまにはエラーに遭遇してしまい、解決するのに時間がかかったりします。BLASTは、大変使用頻度の高いツールの1つですが、見慣れないエラーに遭遇したので、エラーと対応方法をまとめます...
blast 使い方 makeblastdb
blastを実行するには、事前にデータベースを構築する必要があります。データベースの対象となるfastaファイルを準備し、makeblastdbコマンドを実行してデータベースを作ります。とにかく早く問題解決したい人はこちら>>直接、データ解...
blast+ 使い方 best hitの算出 awkとoutfmt7
blastはよく使うツールの1つです。色んな機能があるし、オプションもたくさんあるので使いこなすのが大変です。よく苦労するのが「best hitを1個だけ取り出す」ことですが、実は「outfmt7」を使うと簡単に得ることができます!使用した...
blast 使い方 インストール blast+ v2.3.0
今回はblast+のインストールに挑戦です。プリコンパイル版です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談インストールマニュアルに従い、プリコンパイル版をインストールします。場所# cd /usr/localダウンロード...
bowtie メモリー関連のエラー Warning: Exhausted best-first chunk memory for read
bowtieを--bestモードで実行したところ、メモリーに関するワーニング(警告)が大量に表示されました。エラーではないので処理は進んでいくのですが、「マッピング結果に影響を与えないか?」という点がとても気になりました。そこで、--chu...
bowtie 使い方 マッピングとオプション
bowtie は、代表的なマッピングソフトです。高速にマッピング処理を行うことが出来ます。ただし、ギャップは考慮しません。とにかく早く問題解決したい人はこちら>>直接、データ解析相談Bowtie - 使い方まずは、インデックスファイルを作る...
bowtie2 使い方 アライメント情報(alignment summary)をファイルに保存する
bowtie2のアライメント情報(alignment summary)をファイル保存する方法について考えます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談bowtie2 アライメント情報と標準エラー出力Linuxのコマンド結...
bowtie v1.1.2 - インストール
今回はbowtie v1.1.2のインストールに挑戦です。プリコンパイル版です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談インストール場所cd /usr/localダウンロードwget 展開unzip bowtie-1.1...
bowtie2 v2.2.6 使い方 インストール
今回はbowtie2 v2.2.6のインストールに挑戦です。プリコンパイル版です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談インストールマニュアルに従い、プリコンパイル版をインストールします。cd /usr/localwg...
gffread を使った transcripts fasta 転写物の配列取得
ゲノム配列から転写物の配列を取得転写物(transcript)の配列を取得したい場合、染色体上の位置(開始位置〜終了位置)をたよりに配列情報を抜き出します。位置情報を管理する方法としてよく使われるのがGFF3やGTF形式のファイルです。領域...
VCFtools 使い方 VCFファイルからの抽出(BEDファイル対応領域)
SNP Callerには、通常、領域指定オプションがあります。しかし、既に算出されたVCFデータが領域外のデータを含む場合、領域を参照しながら整理し直すのは面倒な作業です。そんなとき、vcftoolsを使えば、計算した後からでも簡単に領域の...
htseq-count 使い方 gene単位・transcript単位の数え方
RNA-Seq解析では、ゲノムに張りついたリードの数を数えます。gene単位・transcript単位の数え方は、オプションで切り替えることができます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談htseq-count 使い...
tophat2 使い方 エラーと対策
tophat2でマッピングしていたらエラー終了しました。一応解決したので報告します。とにかく早く問題解決したい人はこちら>>直接、データ解析相談Tophat2 マッピングで発生したエラーtophat2のマッピングをバックグラウンドで実行した...
tophat2 使い方 RNA-Seq マッピング
今回は、tophat2を使って、RNA-seqデータのマッピングに挑戦です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談tophat2とはtophat2は、代表的なマッピングツールの1つです。RNA-Seqのデータをスプライ...
tophat2 v2.1.0 - インストール
今回は、Tophat2 v2.1.0 のインストールに挑戦です。プリコンパイル版です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談準備Tophat2が内部で利用する解析ツールを事前にインストールします。bowtie2samt...
SRA_Toolkit 使い方 データ保存フォルダの変更
SRA Toolkiltのprefetchでデータをダウンロードする際の保存フォルダは変更可能です。vdb-configを使って設定します。とにかく早く問題解決したい人はこちら>>直接、データ解析相談SRA_toolkit デフォルトフォル...
SRA Toolkit 使い方 公開データのダウンロードとsra fastq変換
公開されているFastqデータを活用したい場合は、SRA_toolkitを使うと便利です。データのダウンロードやFastqへの展開がコマンド1つで簡単に実行できます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談SRA Too...
UCSCのゲノム配列ファイル .2bitファイルの展開
UCSCのゲノム配列ファイルは、.2bit形式で圧縮されているものがあります。「twoBitToFa」というツールを使って展開します。とにかく早く問題解決したい人はこちら>>直接、データ解析相談.2bit形式について塩基を2ビットで表現して...
fastx_toolkit 使い方 FastqからFastaへ変換
Fastq・Fastaファイルともに、バイオインフォマティクスで良く使うファイル形式です。FastqからFastaへ変換についてまとめます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談FastqからFasta形式への変換よく...
FASTX-Toolkit 使い方 インストール
FASTX-Toolkitは、次世代シーケンサのFASTA・FASTQの前処理に関連するツールを集めたものです。トリミングやクリッピング、クオリティ処理に関するツールが提供されています。インストール方法について手順を整理します。とにかく早く...
samtools 使い方 faidx Fastaから部分配列
塩基配列の変異に注目したSNP解析や配列パターンに注目した解析などでは、部分配列を抜き出して詳しく調べる場合があります。今回は、開始位置-終了位置を指定してFASTAファイルから注目する塩基配列の一部(部分配列)を取り出します。開始・終了の...
samtools 使い方 mpileup ( calling SNPs ) & annotation
samtoolsを使ったVariant callingについてまとめます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談解析作業の流れマッピングCalling SNPsフィルタリングアノテーション付与主なツールCalling ...
samtools 使い方 インストール
SAM tools のインストールについて、バージョン毎にまとめました。手順は同じ。指定バージョンを探すのが面倒なので、それぞれ整理しました。とにかく早く問題解決したい人はこちら>>直接、データ解析相談インストールインストール情報に従い、ダ...
hisat2 使い方 | index 作成
HISAT2のindexファイルを、hisat2-buildコマンドで自ら生成します。とにかく早く問題解決したい人はこちら>>直接、データ解析相談HISAT2は、高速かつ正確に次世代シーケンシング(NGS)データをリファレンスゲノムにマッピ...
HISAT2 使い方 インストールとマッピング
HISAT2 インストールHISAT2のサイトにインストールが簡単なバイナリー版が準備してあるので、ダウンロードして使います。使用したのは、「Linux x86_64 binary」版です。HISAT2HISAT2 インストール手順# ホー...
SnpEff 使い方 アノテーション
samtoolsなどで予測されたVariantに、snpEffを使ってアノテーションをつけてみます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談snpEffsnpEffとは、予測されたVariantにアノテーションやその影響...
Miniconda 使い方 | インストール と 環境作成・環境利用
Miniconda install公式サイトから、最新版 Miniconda installer(各OS毎)のシェルファイルをダウンロードします。Miniconda$ wget $ bash Miniconda3-latest-Linux-...
Trinity 使い方 RNA-seq de novo アセンブル
Trinityは、トランスクリプトーム解析で用いられる「de novo アセンブルツール」です。たくさんのメモリとマルチコアの環境を要求します。100万のペアエンドリードに対して1GB程度のメモリが推奨されています。とにかく早く問題解決した...
Trinity Docker 面倒な「install作業なし」で一発起動!
Trinityを、一番「簡単に動かす方法」だと思います。とにかく早く問題解決したい人はこちら>>直接、データ解析相談Trinity in Docker 使い方アセンブルツール Trinity(de novo transcriptome as...
Trinity v2.2.0, v2.4.0 インストール
Trinity v2.2.0、v2.4.0のインストールに挑戦です。ソースダウンロードして、make作業を行なってインストールします。とにかく早く問題解決したい人はこちら>>直接、データ解析相談インストールmake関連場所cd /usr/l...
Trinotate インストール | conda/biocondaを使って簡単インストール
trinotateは、trinityのアセンブル結果にアノテーションをつけるツールです。とにかく早く問題解決したい人はこちら>>直接、データ解析相談trinotate installTrinotate便利なツールなのですが、インストールが必...
fastp インストール・簡単な使い方・よく使うオプション
Fastqファイル用の高速前処理ツールです。クオリティの評価(トリミング前後)やトリミング処理を行います。とにかく早く問題解決したい人はこちら>>直接、データ解析相談fastp とは?fastp : fastp 使い方簡単な使い方シングルエ...
macs2 install & 使い方
MACS2、良く使われている ChIP-Seq 解析用の peak caller です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談macs2色々なサイトで、macs2 のインストール方法や使い方がまとめられていますが、バー...
freebayes インストール 簡単な使い方
freebayesは、代表的なSNP検出ツールの1つです。ベイズ的アプローチでIndelを検出します。インストールと簡単な使い方についてまとめます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談freebayes とはこちらで...
htseq-count 使い方 よく使うオプション
htseq-countは、RNA-Seq解析などでゲノムに張りついたリードの数を数えるツールです。よく使うオプションをまとめます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談htseq-count 使い方基本的な使い方hts...
HTSeq htseq-count 使い方 インストール v0.7.2
htseq-countは、RNA-Seq解析などでリードをゲノムにマッピングしたとき、どこにどのくらい張り付いたかをカウントするプログラムです。HTSeqのインストールには、いくつか事前にインストールしておくべきものがあります。とにかく早く...
STAR RNA-seq aligner 使い方 インストール・index・マッピング
STARは、RNA-Seq用のマッピングソフトです。非常に高速なマッピングが可能です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談STARとは?STAR インストールダウンロードした圧縮ファイルにプリコンパイル版が含まれてい...
塩基配列の表示幅を調整する - Fastaファイルの整形 -
Fastaファイルの整形についてまとめます。とにかく早く問題解決したい人はこちら>>直接、データ解析相談フォーマットを整えるFASTX_Toolkitは、Fastq・Fasta形式のファイルについて前処理をしたり、ファイル整形をするのに便利...
Fastq Fasta ファイル形式
Fastq・Fastaファイルともに、良く使うファイル形式です。とにかく早く問題解決したい人はこちら>>直接、データ解析相談ファイル形式Fastq形式1本の配列は4行構成で記述されます。1行目:「@」で始まる配列のID。以降に、説明を記述す...
タイトルとURLをコピーしました