BAMファイルのヘッダーを変更する

BAMファイルのヘッダーを変更したい場合には、samtools reheaderコマンドが利用できる。

samtools reheader

https://www.htslib.org/doc/samtools-reheader.html

samtools reheader – replaces the header in the input file

使用しているsamtoolsのバージョン

samtools version | head -2

samtools 1.21
Using htslib 1.21

BAMファイルのヘッダーを変更したい

インターネット上のどこかから誰か他の人がマッピングしたBAMファイルをダウンロードしてきた。 リードのカバレッジを計算してBIGWIGファイルに出力し、ゲノムブラウザで見ようとしたがうまく見られなかった。 ダウンロードしたBAMファイルのヘッダーを、samtools headで確認すると、 どうやらマッピングに使用した参照ゲノム配列の配列名が自分のゲノムブラウザのものと異なるようだった。

samtools head temp.bam

@HD VN:1.6  SO:coordinate
@SQ SN:1    LN:30427671
@SQ SN:2    LN:19698289
@SQ SN:3    LN:23459830
@SQ SN:4    LN:18585056
@SQ SN:5    LN:26975502
@SQ SN:M    LN:366924
@SQ SN:C    LN:154478

例えば、@SQ SN:1 LN:30427671とあるように第一染色体の配列名は1となっている。 これを@SQ SN:Chr1 LN:30427671のように変更して、配列名をChr1に変更したい。

samtools reheaderコマンド

samtools reheaderコマンドを使うと、SAM/BAMファイルのヘッダーを変更することができる。

コマンドのヘルプメッセージは以下。

samtools reheader

Usage: samtools reheader [-P] in.header.sam in.bam > out.bam
   or  samtools reheader [-P] -i in.header.sam file.cram
   or  samtools reheader -c CMD in.bam
   or  samtools reheader -c CMD in.cram

Options:
    -P, --no-PG         Do not generate a @PG header line.
    -i, --in-place      Modify the CRAM file directly, if possible.
                        (Defaults to outputting to stdout.)
    -c, --command CMD   Pass the header in SAM format to external program CMD.

コマンド（-cオプション）で変更する

先ほどのBAMファイルのヘッダーを変更してみよう。やり方としては以下の2通りある。

1. 正しいヘッダーを持つ別のBAMファイルを用意してそのヘッダーをコピーする。
2. ヘッダーを書き換えるコマンドを指定して変更する。

ここでは2つ目のコマンドで書き換える方法で変更してみる。 以下のように-cオプションにsedでヘッダーを書き換えるコマンドを指定した。

samtools reheader -c 'sed -e "s/@SQ\tSN:/@SQ\tSN:Chr/"' temp.bam | \
   samtools head

@HD VN:1.6  SO:coordinate
@SQ SN:Chr1 LN:30427671
@SQ SN:Chr2 LN:19698289
@SQ SN:Chr3 LN:23459830
@SQ SN:Chr4 LN:18585056
@SQ SN:Chr5 LN:26975502
@SQ SN:ChrM LN:366924
@SQ SN:ChrC LN:154478
@PG ID:samtools PN:samtools VN:1.21 CL:samtools reheader -c sed -e "s/@SQ\tSN:/@SQ\tSN:Chr/" temp.bam

うまくヘッダーを書き換えることができた。