Bioinformaticsで使うファイルフォーマットまとめ

TL;DR

バイオインフォマティクスをしていて、障壁になることの1つにファイルフォーマットが多すぎる、という問題があります。ツールを動かそうとすると出現するフォーマットが多く、どうやってこの形式のファイルを作ればいいんだということはよくあります。備忘録を兼ねて、よく使うフォーマットと関連するツールについてまとめておきます。基本的なフォーマットは網羅してるはずですが、また新しいフォーマットとかに出会えば追記していきます。

よく使うファイルフォーマット一覧

• fasta
• fastq
• sam/bam
• bed
• gtf/gff
• wig/bigwig
• vcf/gvcf

fasta

いろんな場面で使いますが多分一番最初に見ること人が多いファイルフォーマットです。 >で始まるID行と、配列データそのものを保存する行に分かれています。配列行では改行が許されています。配列行ではIUB/IUPACで規定されている塩基配列コードとアミノ酸コードを使用できます。詳しくはWikipediaとかを参照してください。

fasta sample

例としては、以下のようなフォーマットになります。

>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

使い道としては、

• blastdbの作成
• マルチプルアラインメントの作成
• 系統解析
• bowtie2やSTARなどのmapping toolのインデックスの作成
• bedファイルなどから配列データへのアクセス

などが主な使い道でしょうか。

fastaフォーマットを扱うツール

fastq

NGS解析で一番最初に作成されるファイルフォーマット。厳密には画像データが一次データですが、多分シーケンサーを持っていてそこからデータを直接扱っていない限りはこれ以前のファイルは見ないんじゃないんでしょうか。

fastqファイルには、NGSで読まれたリードの名前を示す@から始まるヘッダ行、配列、配列のクオリティが記載されています。また、配列と配列クオリティを分けるために+から始まるヘッダ行が配列と配列クオリティの間に置かれています。fastaフォーマットとは違い、配列、配列クオリティ行内では改行が許されていません。

fastq sample

例えばNCBIのSRAに存在するfastqは以下のようなフォーマットになります。

@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC
+SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC

配列には、AGCTNのみが許されており、配列クオリティには、Phredクオリティスコア（下の式）というものが使われています。基本的に高いほどシーケンサーのエラーである可能性が低いです。最近のバージョンではサンガーの式が使われていますが、Wikipediaによるとオッズ比などが使われていることもあるそうです。実際には数字ではなくASCIIコードで33から126の文字としてエンコーディングされます。このエンコーディングはSAM/BAMフォーマットでも共通のものです。

$Q = -10log_{10}p$

このファイルフォーマットはクオリティコントロール程度にしか使われず、基本的にはSAM/BAMに変換してから扱うことが多い印象です。最近ではRNA-seqなどにはSAM/BAMを介さずそのまま発現量測定などをすることもあります。

クオリティコントロールツール

クオリティコントロールには以下のようなツールがよくつかわれている気がします。他にもいろいろあります。

• fastQC
• fastp
• Trimmomatic

マッピング・定量ツール (fastq -> SAM/BAM)

SAM/BAMに変換する際には、以下のようなMapping Toolが使われていることが多いように思えます。RNA-seqの際にはイントロン等を考慮する必要があるので、DNAを読むときとは別に処理が必要になり、専用のMapping Toolを使う必要があります。Bisulfite SequenceなどはDNAですが、処理が特殊なので専用のMapping Toolが必要です。

bulk NGS sequence

その他にリードの分割やダウンサンプリングなどを行いたい場合にはfastaで紹介したようなseqkitなどが有用です。マージはcatとかでいいです。小ネタとしてgzip形式のものでもcatでマージできます。

bulkのRNA-seqでは以下のようなツールでSAM/BAMを介さずそのまま発現量テーブルを作成でます。また、これらのツールのほうが精度は高いらしいです。

この辺はCSVとかわかりやすい形式ではなく、よくわからない形式で出力されるのでR packageのtximportなどを使ってテーブル形式に変換します。変換の仕方などはこちらが参考になります。

scRNA-seq

scRNA-seqを扱う場合には、それ専用のツールがまたいろいろありますが、代表的なものとしては以下のようなものがあります。

miRNA-seq

miRNAの定量はイントロンとかないのでDNAと同じ感じでもいいのですが、isomirみたいな概念もあり、なんか色々ツールがあったりします。BAMとかじゃなくて独自形式に変換されていくものが多いです。独自形式を統一するための概念としてmirGFFというものが提案されていますが、かなり未成熟な印象があります。詳しくはmirGFF formatのところで書きます。

TCGAとかで扱われているmiRNA-seq解析はまた別のパイプラインが使われていたりします。

SAM/BAM/CRAM

マッピングを行ったあと扱うようになるファイルフォーマットです。BAMはSAMをバイナリ化したものでフォーマットとしては同一です。CRAMはfasta情報を使って更に圧縮率を上げることができるフォーマットです。あまりSAMのまま扱うことはなく、BAM/CRAMに変換されることが多いです。リードのヘッダ、配列、クオリティ、マッピング位置などほぼすべての情報が格納されている。情報が膨大なので、フォーマットの詳細はマニュアルを読んでほしいです。マニュアル以外の有用そうなリンクをまとめておきます。

SAM sample

例としてはこんな感じです。@から始まるヘッダ行とリードの情報が格納されているボディ部分に分かれています。

@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r002 0 ref 9 30 3S6M1P1I4M * 0 0 AAAAGATAAGGATA *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1

マニュアル以外の有用そうなリンク

• SAM Formatのcigar列の読み方(samtoolsとか)
• SamファイルのCIGAR string, MD tagから変異(mismatch&indels)を取得する
• Explain SAM Flags - GitHub Pages
• Strand-specific アラインメントの分割

基本操作

基本的にはsamtoolsを使えばたいていのことはできます。picardなども有用です。

可視化

どんなふうにリードが貼りついているのか、などを確認するのはクオリティコントロールなどの観点から重要です。IGVを使えば簡単に可視化できます。IGVは後述するGFF/GFTやbed、wig/bigwig、bedgraphなども可視化できるのでとりあえずインストールしておくべきツールです。

クオリティコントロール

duplicate readの除去や、マッピングクオリティによるフィルターなどを行うことがあります。基本的には先ほど上げたツールを使えば問題ないですが、少し複雑なフィルターがしたい時などには、deeptoolsのalignmentSieveコマンドが便利です。フラグメントサイズによるフィルターやStrand Specificなリードの抽出などを行えます。あとはbamUtilsを使えばリードのトリミングとかができます。

発現量の定量 (bam -> csv etc.,)

RNA-seqを行った後に行う代表的な解析は、発現量の定量です。ツールとしては色々ありますが、代表的そうなものを紹介します。cuffdiffなんかは有名ですが使用は推奨されていないようです。

遺伝子アノテーション (bam -> gff/gtf)

変異の検出 (bam -> vcf)

SAM/BAMフォーマットからVariant Call Format(VCF)に変換するステップと考えてもいいです。基本的にはSNVの検出を想定しており、SVなどは考慮していません。よく使われていそうなツールは以下のようなものがあります。

GATKに関する処理について、公式以外で役に立ちそうなリンクを貼っておきます。

※ samtools mpileupは推奨されていないようです。

ピークのコール (bam -> bed etc.,)

ChIP-seqなどではリードが集中した領域をピークとして扱うことが多いです。この時もBAMから何らかのフォーマットへの変換が行われます。たいていはbedに準ずる形式へと変換されます。代表的なツールは以下のようなものがあります。また、転写因子に関するChIPでは狭いピークが、ヒストン修飾などのChIPでは広いピークが見られます。これらは、検出方法が異なるので、ツールやオプションを使い分ける必要があります。たいていbedかそれに準ずる形式のファイルに変換されます。

htslib

samtoolsの本体です。SAM/BAMフォーマットを扱う際のAPIを提供しています。凝ったことをしたくなると使います。もともとはCで書かれていて、いろんな言語でWrapperが作成されています。個人的に知っているのは以下です。GitHubへのリンクを貼ります。

• pysam
• rust-htslib
• hts-nim

各種変換

bam -> fasta

変換する用途があまり思い浮かびませんが、samtoolsを使えばできます。

samtools fasta input.bam > output.fasta

bam -> fastq

samtoolsとかbedtoolsを使えば変換できます。あまり使うことはない気がします。k-merとか使って機械学習したいときとかにマップされたリードだけ使う、などの用途が考えられます。

# samtools
samtools fastq input.bam > output_single.fastq
# bedtools
bedtools bamtofastq -i input.bam -fq output_single.fastq

bam -> bed

samtoolsとawkでできる気はしますが、bedtoolsを使うと簡単です。

bedtools bamtobed -i input.bam > output.bed

bam -> bedgraph

bedtoolsのgenomecovかdeeptoolsのbamCoverageで変換できます。

# bedtools genomecov
bedtools genomecov -i input.bam -bg > output.bedgraph

# deeptools bamCoverage
bamCoverage -b input.bam -o output.bedgraph -of bedgraph

bam -> bigwig

deeptoolsのbamCoverageを使います。オプションとかは色々あって、ノーマライズなどもしてくれたりします。ChIPなどの解析の際にも、RPGC normalizeに対応しているので使えます。RNA-seqでもRPKMやCPMに対応しています。TPMにも対応してほしいです。

bamCoverage -b input.bam -o output.bw

bed

bedtoolsなどで扱います。Pythonなどではpybedtoolsのようなライブラリが提供されています。最初の三行(chrom, chromStart, chromEnd)が必須で、その他が自由なフォーマットです。一応ある程度は決まっていて、UCSCのFAQでは、

1. chrom: 染色体名
2. chromStart: スタート位置(0-index)
3. chromEnd: 終了位置
4. name: 遺伝子名など
5. score: 任意のスコア(track上での色の濃淡とかに反映される)
6. strand: strand (+, -)
7. thickStart: CDSの開始位置
8. thickEnd: CDSの終了位置
9. itemRgb: track上でのRGBカラー
10. blockCount: exonのブロック数
11. blockSizes: ブロックサイズ
12. blockStarts: exonの転写開始位置から見たスタート位置

という風に決まっているそうです。7行目以降は可視化する際に使われるパラメーターです。最初の三行のみのBEDをBED3、6行目までのBEDをBED6、12行目までのBEDをBED12と呼んだりします。三行目までのデータがあればbedtoolsで扱うことができます。またフォーマットは微妙に異なるのですが、GFFとかVCFもbedtoolsで扱えます。そういう意味では非常に基本的なフォーマットです。

bed -> fasta

bedtoolsのgetfastaを使うことで変換できます。興味のある領域のbedを作成した後、getfastaで配列を取得してMotif Enrichmentを行うなどの使用用途があります。

bedtools getfasta -fi genome.fasta -bed input.bed > output.fasta

bed -> bam

bedtoolsのbedToBamで変換できます。ただmutation情報などは失われます。

bedToBam -i input.bam -g genome > output.bam

bedgraph

Bedの亜種っぽい感じです。ProbabilityやTranscriptomeなど連続性のあるデータを表示させるために使われるフォーマットらしいです[参考]。あまり使ったことがありませんが、Bisulfite Sequenceの解析の際にMethylDackelというツールを使うと出てきました。MACS2のinput/outputでも使われてます。

gff/gtf

遺伝子のアノテーションなどは基本的にこのフォーマットでまとまっていることが多いです。GFFにはversion2とversion3があり、微妙にフォーマットが違います。また、GFF/GTFを扱うツールとしては以下のようなものがあります。

GFF format

1. chrom: 染色体番号や、Scaffold番号など
2. source: 何をもとに作られたか、どこのデータかなど
3. feature: CDS, exon, gene, five_prime_utrなど
4. start: featureの開始位置
5. end: featureの終了位置
6. score: なにかのスコア
7. strand: (+, -, .)。.は方向が不明な際に使われる。
8. frame: coding exonの場合はどのフレームなのかが書かれている。
9. attribute: 他のデータがセミコロン区切りで入力されている。gene_idや、parent_idなど。

GFF3 sample

X	Ensembl	Repeat	2419108	2419128	42	.	.	hid=trf; hstart=1; hend=21
X	Ensembl	Repeat	2419108	2419410	2502	-	.	hid=AluSx; hstart=1; hend=303
X	Ensembl	Repeat	2419108	2419128	0	.	.	hid=dust; hstart=2419108; hend=2419128
X	Ensembl	Pred.trans.	2416676	2418760	450.19	-	2	genscan=GENSCAN00000019335
X	Ensembl	Variation	2413425	2413425	.	+	.	
X	Ensembl	Variation	2413805	2413805	.	+	.

GTF format

基本的にGFFと同じですが、9行目が厳格に決められており、gene_idとtranscript_idを持たなければならない。また、#によるコメントが許されていなかったりする。

GTF sample

chr1    hg19_rmsk       exon    16777161        16777470        2147.000000     +       .       gene_id "AluSp"; transcript_id "AluSp";

chr1    hg19_rmsk       exon    25165801        25166089        2626.000000     -       .       gene_id "AluY"; transcript_id "AluY";

chr1    hg19_rmsk       exon    33553607        33554646        626.000000      +       .       gene_id "L2b"; transcript_id "L2b";

gff/grf -> fasta

gffread input.gtf -g genome.fa -w output.fa

gff/gtf -> bed

gffread input.gtf -g genome.fa -E --bed -o output.bed

wig/bigwig

bigwigはwigをバイナリ化したものです。UCSC genome browserで可視化するときに使用されている形式です。bigwigなどはbamとかと比べると本当に軽いので、可視化などが目的のときは最もおすすめできるフォーマットです。また、deeptoolsを使うことで、ヒートマップやPCA、相関解析などを行えます。deeptoolsは高機能なのでこれはこれで記事が書きたいです。あとオープンソースなのでGithubのコードを読むと勉強になります。Pythonで書かれています。

vcf/gvcf

変異情報が格納されているフォーマットです。samtoolsと同じところがフォーマットを決定しており、現状はver4.2です。非常に情報量が多いフォーマットなので、詳細はマニュアルを参照してください。gvcfは変異がコールされていない、という情報を加えて含んでいます。gatkのワークフローで登場しますが、あまり解析には使わないイメージです。詳しくは公式ページなどを参照してください。また、vcfを扱うツールとしては以下のようなものが有名です。ただvcfに関してはプログラム組んで動かした方が早い気もします。

プログラミング言語として扱えるパッケージはいろいろありますが、htslibのWrapper系列は大体対応しています。タブ区切りのファイルなので、Pythonならpandas等でも扱えます。他には、

• pyvcf
• vcfR

などが候補です。

vcf sample

##から始まるヘッダ行とそれ以外のボディ部分に分かれています。

##fileformat=VCFv4.2
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4
20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
20 1234567 microsat1 GTC G,GTCT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3

dbSNPs

よく知られているSNPsなどはデータベースとしてまとまっていて、これらは基本的にvcfフォーマットで配布されています。broadinstituteのgoogle cloud platformとかで配布されています。

たまに使うファイルフォーマット

• twobit
• mirGFF

twobit

deeptoolsを使ってGCBiasを補正するときに使いました。他に使ったことはないです。たぶんfastaをビット形式で扱っているので、効率がいいです。UCSCのツール群にあるfaToTwoBitを使えば作成できます。もう少し詳しい使い方などはこちらのサイトが詳しいです。

mirGFF

GFFに準拠したようなフォーマットでmiRNA系のNGSデータを統一的に扱うために策定されたフォーマットです。昔見たときはbioaxivだったんですが、最近論文になっているようです(Desvignes et al., 2020 Bioinfomatics)。様々な形式から相互変換ができるフォーマットでmirtopというパッケージを使って作成できます。miRNA-seqで使われているツールとしては以下のようなものがあります。なんか昔はエンコーディングが対応してなくて変換できない、とかだったんですが修正されているのでしょうか。

• seqbuster
• miRge2.0
• isomiR-SEA
• sRNAbench
• Prost!

これらのツールは独自形式のものを出力することが多いのですが、mirTopを通すことで、以下のようなフォーマットに変換できます。

• mirGFF3
• isomiRs
• VCF
• fasta
• count matrix

また、isomiRsというのはisomiRを考慮した解析するisomiRsというR packageで用いられている形式になります。

最後に

思ったよりすごい分量になってしまいました。間違いなどがあればご指摘いただけると幸いです。