タブ区切りゲノムポジションファイルにindexをつけて素早く問い合わせる tabix

2024/02/04追記、02/06追記

　コマンドラインやゲノムビューアで局所的なゲノム特徴を調べる場合、指定された領域に重なる特徴を検索するインターバルクエリを頻繁に実行する必要がある。インターバルクエリを数回しか行わない場合には、データファイル全体を読み込むこともできるし、インターバルクエリを頻繁に行う場合（ゲノムビューアなど）には、ファイルをメモリにプリロードしておくこともできる。しかし、ファイル全体を読むことは、巨大なデータセットを考えると、どちらの戦略も非効率的になる。効率性の問題に対する解決策は、データファイル用のデータベースを構築することであるが(Kent et al., 2002; Stein et al., 2002)、これも最適ではない。データベースの設定とスキーマの設計の技術的な複雑さもまた、平均的なエンドユーザーにとってのデータベースアプローチの採用を妨げている。

　このような状況下で、データ圧縮とリモートファイルアクセスをサポートしながら、巨大なデータセットへの効率的なランダムアクセスを実現するために、bigBed/bigWig (Kent et al., 2010)やBAM (Li et al., 2009)を含むいくつかの特殊なバイナリフォーマットが最近開発され、ルーチンのデータ処理とデータの可視化に大きく貢献している。現在のところ，これらの高度な索引付け技術は，BED，Wiggle，BAMにしか適用されていない．それにもかかわらず、ほとんどのTABで区切られた生物学的データフォーマット（PSL、GFF、SAM、VCF、UCSCデータベースダンプなど）には染色体の位置が含まれているため、これらすべてのフォーマットに対応したインデックスを作成する汎用的なツールが可能であることは想像に難くない。そのツールがTabixである。ここでは、Tabixのインデックス作成の技術的な進歩を説明し、アルゴリズムを説明し、生物学的データに対する性能を評価する。

manual

http://www.htslib.org/doc/tabix.html

current release

http://www.htslib.org/download/

インストール

tabixとbgzipは今ではHTSlib projectの一部になっているので、HTSlibをビルドした。condaでも導入可能。

Github

#bioconda (link) 
conda install -c bioconda tabix
conda install -c bioconda/label/cf201901 tabix

#from source (v1.11)
wget https://github.com/samtools/htslib/releases/download/1.11/htslib-1.11.tar.bz2
tar -jxvf htslib-1.11.tar.bz2
cd htslib-1.11/
make -j

> ./tabix

Version: 1.11

Usage: tabix [OPTIONS] [FILE] [REGION [...]]

Indexing Options:

-0, --zero-based coordinates are zero-based

-b, --begin INT column number for region start [4]

-c, --comment CHAR skip comment lines starting with CHAR [null]

-C, --csi generate CSI index for VCF (default is TBI)

-e, --end INT column number for region end (if no end, set INT to -b) [5]

-f, --force overwrite existing index without asking

-m, --min-shift INT set minimal interval size for CSI indices to 2^INT [14]

-p, --preset STR gff, bed, sam, vcf

-s, --sequence INT column number for sequence names (suppressed by -p) [1]

-S, --skip-lines INT skip first INT lines [0]

Querying and other options:

-h, --print-header print also the header lines

-H, --only-header print only the header lines

-l, --list-chroms list chromosome names

-r, --reheader FILE replace the header with the content of FILE

-R, --regions FILE restrict to regions listed in the file

-T, --targets FILE similar to -R but streams rather than index-jumps

-D do not download the index file

--cache INT set cache size to INT megabytes (0 disables) [10]

--separate-regions separate the output by corresponding regions

--verbosity INT set verbosity [3]

実行方法

１、indexing

#vcf
bgzip -@ 8 input.vcf
tabix -p vcf input.vcf.gz

#bed
bgzip -@ 8 input.bed
tabix -p bed input.bed.gz

#sam（他と同じくcoordinate sortされている必要がある。Pipingしてsamtools sortしておく）
bwa mem ref.fa pair1.fq.gz pair2.fq.gz |samtools sort -O SAM - > input.sam
bgzip -@ 8 input.sam
tabix -p sam input.sam.gz

-p, --preset <STR> gff, bed, sam, vcf
-C, --csi generate CSI index for VCF (default is TBI)

２、querying

tabix input.vcf.gz chr5:1000-9000 > output

#headerも含める
tabix -h input.vcf.gz chr1:1-100 > output

#含まれるchr名だけを取得
tabix -l input.vcf.gz

-h, --print-header print also the header lines
-H, --only-header print only the header lines
-l, --list-chroms list chromosome names

1000ゲノムのQ&Aに書かれているように、curlなどでダウンロードしなくてもFTPのファイルを扱うことが可能。

tabix -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz chr2:39967768-39967768 > output

#samtools viewでも可能。bamの特定領域だけ返す
samtools view -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00154/alignment/HG00154.mapped.ILLUMINA.bwa.GBR.low_coverage.20101123.bam 17:7512445-7513455 > output

補足

samtools faidxでゲノムのポジションを取得する

#chr1全長
samtools faidx genome.fa chr1 > chr1.fasta

#chr2の1000-10000
samtools faidx genome.fa chr2:1000-10000 > chr2_1000-10000.fasta

#抽出領域が多い時はリストファイルを指定（１行１形式）
samtools faidx genome.fa -r target_list > target.fasta