Zenn
📖

【ゲノム解析】BAMファイルから1座標の塩基割合を求める

2023/09/12に公開
R
Ubuntu
WSL 2
bioinformatics
tech

NGSリードを参照配列にマッピングして得られるBAMファイルをIGVで見ると、1座標単位でATGCの数と割合が確認できる。

BAMファイルをIGVで開いた画面

この集計を1ファイルずつ1座標ずつやってられないので、CUIで回したいところだがちょうどいいものを見つけるのが大変だったので記録しておく。最終的にはbam-readcountかRsamtoolsが使用可能であった。

本記事で使用したBAMファイルは、ヒト細胞のWES (whole exome sequence)リードを参照ゲノムhg19にbwa memでmappingし、Picard toolでPCR重複を除いたものである。

失敗した方法

bedtools intersect

bedtools intersect

[参考]
https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

bedtools intersectで2つのinputの共通部分を返す。探索したい塩基を記載したBEDファイルと、マッピング結果のBAMファイルを使用する。

BEDファイル

-a、-bに2つのファイルを指定するが、どちらの順番でも受け付ける。-aには参照元、-bにクエリファイルという順番で良いだろう。

bedtools intersectの出力は-aに指定したファイル形式が出力になってしまう。(-aにBAMファイルを指定すると出力もBAMファイルになる。)
-bedオプションを入れておきBEDファイル形式での出力を指定する。

bedtools intersect -a bamファイルパス bedファイルパス  -bed > out.bed


out.bed

出力は興味のある座標にまたがるNGSリードの情報（リードのIDやStrand、start, stop座標などなど）になっており、その座標の塩基はわからなかった。

bcftools mpileup

bcftools mpileup

[参考]
https://www.biostars.org/p/425139/
https://github.com/samtools/bcftools/issues/992
https://www.biostars.org/p/47945/#47951

bcftools mpileupに使用するBEDファイルで1塩基領域を示す場合は、2列目と3列目を同じ値にするとエラーが出る。
2444414の塩基が欲しければ、2列目は244413、3列目を2444414とする。

bcftools mpileup \
--regions-file bedファイルパス \
bamファイルパス \
--fasta-ref 参照ゲノム(fasta)ファイルパス \
--annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
> out.txt

100行以上のheaderに続いたのち、座標の参照配列、変異配列、、、のような情報が得られる。

out.txt

bcftoolsは遺伝子変異を扱うためのツールなので、単純に座標の塩基割合を知りたいだけだと余分な情報が多い。
また得られる情報は慣れていないと難解。

bcftools queryに渡せばまだ見やすくはなるが、それでも直感的ではない。

bcftools mpileup \
--regions-file bedファイルパス \
bamファイルパス \
--fasta-ref 参照ゲノム(fasta)ファイルパス \
--annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
| bcftools query --format '%CHROM\t%POS\t%REF\t%ALT[\t%AD]\n' \
> out.txt
igvtools count

igvtools count

IGVで見えているデータなのでIGVデータを操作するigvtoolsを使う方法もある。
igvtools countでは指定したウィンドウサイズでの塩基情報を集計できるので、ウィンドウサイズを1にすれば1塩基ずつの情報が得られる。
しかし、特定の領域を指定する方法が無かったので、bamファイル中の全塩基を集計してしまう。処理に時間がかかるし、出力は巨大なので多検体で行うのは現実的ではない。

[参考]
https://software.broadinstitute.org/software/igv/igvtools
https://ncrna.jp/blog/item/65-igv-toolkit

igvtools count \
-w 1 \
--bases \
bamファイルパス \
出力ファイルパス \
参照ゲノムのfastaファイルパス

【bam-readcount】

https://github.com/genome/bam-readcount

BAMファイルから1塩基単位での計測結果を取得するツール。

インストール

GitHubページをcloneしてソースコードからbuildする方法とdockerを使う方法が紹介されている。

ここではソースコードからbuildして使用する。うまくいかなければdockerを使う方法を検討してほしい。

GitHubページをclone
git clone https://github.com/genome/bam-readcount
ソースコードからbuild
# ディレクトリ移動
cd bam-readcount

# ディレクトリ作成して移動
mkdir build
cd build

# ソースファイルをコンパイル
cmake ..
make

実行ファイルは~/bam-readcount/build/bin/bam-readcountに入っていてこのままでは使いにくいので、ここにパスを通すか、パスが通ったディレクトリにシンボリックリンクを作成しておく。
私の場合、この手のツールを入れるtoolsフォルダを作ってパスを通しているのでここにシンボリックリンクを入れておいた。

ln -s ~/bam-readcount/build/bin/bam-readcount ~/tools/

whichコマンドでbam-readcountのパスが表示されていれば成功。

which bam-readcount

実行

参照ゲノムfastaファイル、BAMファイル、クエリとなるBEDファイルを使って特定座標の塩基の計測を実行する。

--site-list引数でBEDファイルを指定する。1座標だけ欲しい場合は2列目と3列目の値を同じにしておいて良い。

BEDファイル

bam-readcountの実行
bam-readcount \
--reference-fasta 参照ゲノムのfastaファイルパス \
BAMファイルパス \
--site-list BEDファイルパス \
--max-warnings 0  \ 
> out.tsv

出力はタブ区切りなので、txtかtsvあたりの拡張子のファイルにリダイレクトで書き出しておくとよい。
塩基の数以外の情報も多いが、1番染色体の2444414番目の塩基にまたがるリードが49本あり、Aが43、Gが6つであることが読み取れる。

out.tsv

結果の見方
https://kazumaxneo.hatenablog.com/entry/2018/07/05/230927

【Rsamtools】

Rsamtoolsを使えばRでBAMファイルを操作することができる。

https://bioconductor.org/packages/release/bioc/html/Rsamtools.html

# BiocManager::install("Rsamtools")
library(Rsamtools)

BAMファイルの登録

まずはBamFile()機能で参照するBAMファイルを登録する。同フォルダにbaiファイルもあれば同時に認識される。

bamfile <- BamFile(file = "BAMファイルパス")

scan parameterの作成

scanBam()機能を使ってこのBAMファイルからクエリとなる情報だけを取得していく。そのためのクエリをScanBamParam()機能を使って用意する。
ScanBamParam()ではwhich=引数に目的の座標、what=引数に取得するリードの情報の項目を指定する。

座標情報

which=引数に指定する座標情報をIRange()、GRanges()を使って作成する。染色体番号と塩基番号を使用するが、今回は1塩基だけが対象なので、IRanges()のstartとendは同じ値を使用している。

chr <- "1"
pos <- 2444414

which <- GRanges(seqnames = chr, ranges = IRanges(start = pos, end = pos))


whichの中身

計測項目

what=引数で指定可能な計測項目はscanBamWhat()機能で確認できる。


what=に指定できる項目名

ここでは対象座標を跨ぐリードのmapping位置の情報と、リードの塩基配列情報を取得する。

what <- c("pos", "seq")

scanBAM()

上記までのscan parameterを使って、座標を跨ぐリードとその情報を取得する。

reads <- scanBam(file = bamfile,
                 param = ScanBamParam(which = which, what = what)
                 )

1:2444414-2444414を跨ぐ49本のリードの情報がリストとして得られる。posにはリードの5'側座標、seqは座標を跨ぐリードの塩基配列である。

readsの中身

リードの保存順は、各リードの5'側座標の昇順である。

readsの保存順をIGVで見た場合

目的座標の塩基を取得

readsリストのseqから目的座標の塩基だけの情報にする。

# 各リードの何番目が目的座標に当たるかを計算
correct_pos <- pos - reads$`1:2444414-2444414`$pos + 1

# 配列情報のみを取り出す
sequences <- reads$`1:2444414-2444414`$seq

# 結果の保存先ベクトル
base <- c()
for( i in seq(length(sequences))){
  # 配列をベクトル化
  sequence <- as.vector(sequences[i])
  # 1文字ずつ異なる要素として保存
  sequence <- unlist(strsplit(as.vector(sequences[i]), ""))
  # 目的座標の要素を取得して、baseに追記
  base <- c(base, sequence[correct_pos[i]])
}

Discussion