シークエンスデータのファイル形式とBiopythonオブジェクト【Python】

2020-05-252020-10-25

ここでは、シークエンスデータを表す各種のファイル形式について説明します。BiopythonではSeqRecordオブジェクトとしてシークエンスデータを扱いますが、それぞれのファイル形式がSeqRecordオブジェクトにどのように変換されるのかを示します。

バイオインフォマティクス環境

以下の方法で構築した環境であることを前提として説明します。

Windows10 + WSL (WindowsユーザーのためのPythonを用いたゲノム解析環境)

※ 検証はしていませんが、プログラム自体は基本的にmacOS、Linuxでも動作すると思われます。

[toc]

フラットファイル形式

DDBJ/ENA/GenBankで使われている形式で、DDBJ形式やGenBank形式などとも呼ばれます。配列の情報と、その配列にまつわる様々な付加情報が含まれています。

DDBJに登録されているフラットファイル形式のフォーマットについては以下を参照してください。GenBankやENAに登録されているものとは記載の順序など若干の違いがありますが、書かれている情報としては同一です。

https://www.ddbj.nig.ac.jp/ddbj/flat-file.html

フラットファイル形式の構成

フラットファイル形式の主な項目の意味は次のようになっています。

LOCUS	Locus名、配列長、塩基配列の分子タイプ、塩基配列の分子形態、Division、データの最終公開日 ※ Locus 名は、データベース中でそのエントリのみが持つユニークな名前が付けられています。 ※ DDBJ では、1996年7月以降公開されたデータに、アクセッション番号と同じ名前を付けています。
DEFINITION	データの定義や遺伝子などに関する簡単な情報（一定の記述ルールに従ってデータバンク側で作成しています）
ACCESSION	INSD が発行する登録番号（アクセッション番号）
VERSION	アクセッション番号とバージョン番号の組み合わせ
DBLINK	BioProject, BioSample の accession number, Sequence Read Archive の Run accession number 他、特定データベースへのリンク
KEYWORDS	データの詳細種別 (EST, TSA, HTC, HTG, GSS, WGS, TPA など)、配列の特性、実験手法、ゲノム配列の完成度などに関する情報（必要に応じて）
SOURCE	配列データが由来する生物の学名 (と存在する場合は一般名)
ORGANISM	由来生物の生物名と系統関係
REFERENCE	そのエントリに関連する論文情報
COMMENT	FEATURES/Location/Qualifiers で記述できないその他の情報やコメントなど
FEATURES	登録配列の生物学的な特徴 ※ feature key （特徴を表す項目）、location （配列上の位置情報）および Qualifier （特徴をさらに特定する項目）を用いて記述しています。
ORIGIN	そのエントリの表す配列情報

フラットファイル形式のBiopythonでの扱い

SeqIO.read関数を用いて、フラットファイル形式からSeqRecordオブジェクトを取得することができます。例として、フラットファイル形式をTogoWSを用いて取得して、SeqRecordオブジェクトとするコードを示します。

from Bio import TogoWS, SeqIO
with TogoWS.entry('nucleotide', 'NC_045512.2') as handle:
    record = SeqIO.read(handle, "genbank")

SeqRecordオブジェクトの属性とフラットファイル形式の項目は次のように対応します。

SeqRecordオブジェクトの属性		フラットファイル形式の項目
annotations	←→	下記以外の項目
description	←→	DEFINITION
dbxrefs	←→	DBLINK
features	←→	FEATURES
id	←→	VERSION
letter_annotations	←→	なし
name	←→	Locus名
seq	←→	ORIGIN

例えば上記の例で、フラットファイル形式におけるDEFINITIONを取得するには次のようにします。

print(record.description)

FASTA形式

もともとは配列類似性検索ツールFASTAにおける入力配列形式として開発されたもので、現在は塩基配列やアミノ酸配列の情報を保持するために幅広い分野で用いられているファイル形式です。

FASTA形式の構成

FASTA形式では2つの行で1つのエントリーを表現します。

1行目：ヘッダー行

>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome

>エントリーのID エントリーの定義・説明

‘>’に引き続いて、エントリーのID(アクセッション番号など)が記述されています。
1番最初のスペースの後に、そのエントリーの定義・説明などが記述されています。

2行目以降：シークエンスデータ領域

ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCT
GTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACT
CACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATC
    ....

配列を表す文字列(核酸・アミノ酸を表すコード)が記述されています。
次の’>’が出てくるまでがシークエンスデータ領域です。（次の’>’から新しいエントリーになります）

FASTA形式のBiopythonでの扱い

FASTAファイルに含まれるエントリーが1つのみの場合

SeqIO.read関数を用いて、FASTA形式からSeqRecordオブジェクトを取得することができます。例として、TogoWSを用いてFASTA形式を取得して、SeqRecordオブジェクトとするコードを示します。

from Bio import TogoWS, SeqIO
with TogoWS.entry('nucleotide', 'NC_045512.2', 'fasta') as handle:
    record = SeqIO.read(handle, "fasta")

SeqRecordオブジェクトの属性とFASTA形式の項目は次のように対応します。

SeqRecordオブジェクトの属性		FASTA形式の項目
annotations	←→	なし
description	←→	そのエントリーの定義・説明など (フラットファイルのDEFINITION)
dbxrefs	←→	なし
features	←→	なし
id	←→	エントリーのID (フラットファイルのVERSION)
letter_annotations	←→	なし
name	←→	エントリーのID (フラットファイルのVERSION)
seq	←→	シークエンス (フラットファイルのORIGIN)

FASTAファイルに複数のエントリーが含まれる場合

SeqIO.read関数はFASTAファイルに複数のエントリーが含まれる場合は’more than one record found in handle’というエラーを出して停止してしまいます。その場合は、次のようにSeqIO.parse関数を用いて1つ1つのエントリーごとにSeqIO.read関数で読み込みましょう。

例えば、multi_seq.fastaというファイルに複数のエントリーが含まれる場合は、次のように読み込みます。

from Bio import SeqIO
records = SeqIO.parse('multi_seq.fasta', 'fasta')
for record in records:
    print('------------------------------')
    print(record)

SeqIO.parse関数によってSeqRecordオブジェクトをまとめたオブジェクトを取得できるので、forループもしくは組み込み関数のnext関数を用いて1つ1つのSeqRecordオブジェクトを取り出します。

FASTQ形式

FASTA形式にSequenc Quality値を加えたデータ形式です。NGSの出力するデータは配列情報に加えて、その精度を表すSequence Quality値がセットになっているので、FASTQ形式を用いることでNGSの出力データを表現することができ、NGS配列データ形式のデファクトスタンダードとなっています。

FASTQ形式の構成

FASTQ形式は次の4行で1セットとなり、1つのリード(シークエンス)を表します。

1行目：’@’と配列のID（スペース区切りで配列の説明を加えてもよい）
2行目：塩基配列
3行目：’+’（スペース区切りで配列のID・説明を加えてもよい）
4行目：2行目の塩基配列に対応するSequence Quality値

最小のFASTQ形式は以下のようになります。(Wikipediaより)

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

1行目はFASTA形式のヘッダー行に相当します。また、1行目の@以降と3行目の+以降は同じ文字列を入れることが多いようです。

4行目のSequence Quality値の表現は記号の羅列で暗号のようになっていますが、これはその記号のASCIIコード – 33がそのSequence Quality値を表しています。Sequence Quality値の数値そのままでは2桁の値になったときに塩基配列部分と文字数が一致しなくなってしまうので、1つのSequence Quality値を1文字で表現するためにこのような形になっています。なお、NGSのベースコールはPhredが広く普及しており、基本的にこのSequence Quality値はPhredクオリティスコアになります。

実際のFASTQ形式ファイルは上記の1～4行が一まとまりで1つのリードを表し、改行によって新たなリードとして認識されます。FASTQ形式では複数のリードが連続していき以下のようになります。

@DRR000001.265 3060N:7:1:502:2032 length=36
GTTTTTCCCCATTATTTATACCTCTGATAAAAGTAA
+DRR000001.265 3060N:7:1:502:2032 length=36
IIIIIIIIII<II@IGIHI3B3IA?1322+)--/:%
@DRR000001.572 3060N:7:1:620:2034 length=36
GGTGACAGCAGGATTACGGAAGACANNNNTNNGNNT
+DRR000001.572 3060N:7:1:620:2034 length=36
IIIIIIIIIIIIHAC=869-3852*!!!!+!!#!!0
@DRR000001.904 3060N:7:1:873:2032 length=36
GGCGGTTGTCAAAATAGGGATTCGATTTGCCGTTAA
+DRR000001.904 3060N:7:1:873:2032 length=36
IIIII*I>I6+9AI+F.:I138(.(,1<&&(%)*(&

FASTQ形式のBiopythonでの扱い

SeqIO.read関数を用いて、FASTQ形式からSeqRecordオブジェクトを取得することができます。

SeqRecordオブジェクトの属性		FASTQ形式の項目
annotations	←→	なし
description	←→	その配列のID・説明
dbxrefs	←→	なし
features	←→	なし
id	←→	配列のID
letter_annotations	←→	Sequence Quality値
name	←→	配列のID
seq	←→	配列

ただし、複数のリードが含まれたデータの場合は、そのままSeqIO.read関数で読み込もうとすると’more than one record found in handle’というエラーで止まってしまいます。複数のリードを含む場合は、SeqIO.parse関数を用いて1つ1つの分割して読み込みましょう。

NGSから出力されるデータや公共データベースからの入手されるデータは大量のリードから構成されているので、通常はSeqIO.parse関数を用いることになると思います。ここでは例としてまず以下のファイルを用意します。

@SEQ_ID test_seq
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+SEQ_ID test_seq
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
@SEQ_ID2 test_seq
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+SEQ_ID2 test_seq
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

ファイル名 seq.fastq として保存し、Biopythonで読み込みます。SeqIO.parse関数により複数のSeqRecordオブジェクトをまとめたオブジェクトとして取得できるので、forループにより一つ一つのSeqRecordオブジェクトを取り出します。

from Bio import SeqIO
records = SeqIO.parse('seq.fastq', 'fastq')
for record in records:
    print('------------------------------')
    print(record)

------------------------------
ID: SEQ_ID
Name: SEQ_ID
Description: SEQ_ID test_seq
Number of features: 0
Per letter annotation for: phred_quality
Seq('GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT', SingleLetterAlphabet())
------------------------------
ID: SEQ_ID2
Name: SEQ_ID2
Description: SEQ_ID2 test_seq
Number of features: 0
Per letter annotation for: phred_quality
Seq('GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT', SingleLetterAlphabet())

SRA形式

FASTQ形式を圧縮した形式で、Sequence Read Archive (SRA)で使われています。テキストファイルであるFASTQ形式をバイナリ形式で圧縮することで、ファイルサイズは10分の1程度にまで圧縮できますが、SRAデータのためだけの専用の形式であり、汎用の圧縮形式ではありません。

NGSの出力データは膨大な量のファイル容量が必要であり、データベースの容量を押さえるために、NCBIなどのデータベースではFASTQ形式では直接ダウンロードできずに、SRA形式でしかダウンロードできなくなっています。

SRA形式はNCBIが作成し配布しているSRA Toolkitのfasterq-dumpでFASTQ形式に変換できます。(fasterq-dumpで直接ダウンロードすることも可能です)

SeqRecordオブジェクトからファイルへの保存

SeqIO.write関数を用いて、SeqRecordオブジェクトを指定したフォーマットのファイルに保存することが可能です。

フラットファイル形式でダウンロードしSeqRecordオブジェクトとして格納したものを、FASTA形式でファイルに保存する例です。

from Bio import TogoWS, SeqIO
with TogoWS.entry('nucleotide', 'NC_045512.2') as handle:
    record = SeqIO.read(handle, "genbank")
SeqIO.write(record, 'seq.fasta', 'fasta')

このようにSeqRecordオブジェクトはフラットファイル形式やFASTA形式、FASTQ形式などの区別なく情報を格納できるので、保存するファイル形式も自由に変更できます。