bamCoverageを使いカバレッジトラックを作成する

2019 9/13 インストール追記

2024/.04/15更新

deeptoolsはRNA-seq解析やchip-seq解析に特化したアライメントのカウント分析ツール(webサーバ)である。ヒートマップ出力などの機能を持ち、ツールの中にあるbamCoverageを使うと、bamのカバレッジ情報をwig形式などで出力してカバレッジトラックをIGVなどに表示させることができる。

bamCoverage

bamCoverage — deepTools 2.3.1 documentation

インストール

Github

https://github.com/deeptools/deepTools

#conda(link)
mamba create -n deeptools python=3.8 -y
conda activate deeptools
mamba install -c bioconda deeptools

#pip(PyPI)
pip install deeptools

実行方法

カバレッジ情報をウィンドウサイズ10で出力する。マッピングクオリティが10以下リードはカウントしない。

bamCoverage --bam input.bam -o output.bw -of bigwig --binSize=10 --minMappingQuality 10 --numberOfProcessors=max

--minMappingQuality INT　If set, only reads that have a mapping quality score of at least this are considered. (default: None)
-o　outFileName
-of=bigwig　Output file type. Either “bigwig” or “bedgraph”. Possible choices: bigwig, bedgraph
--binSize=50　Size of the bins, in bases, for the output of the bigwig/bedgraph file.
--numberOfProcessors=max　Number of processors to use. Type “max/2” to use half the maximum number of processors or “max” to use all available processors.

他にも様々なオプションがある。詳細はこちらを参考にしてください。

出力されたbigwigをIGVに読み込ませる。

igv -g reference.fa input.bam,output.bw

f:id:kazumaxneo:20171103214125j:plain mapping qualityが30以下をカウント対象から外して計算したため、図の中央左の白いリード付近のカバレッジはほぼゼロになっている。これは複数箇所にアライメントされたリードはmapping quality scoreが非常に低くなる（ゼロになる）ため、上の計算条件ではカバレッジのカウント対象から外れるためである。

一方で、塩基置換が複数おきても曖昧なマッピングではないならばmapping quality scoreは影響を受けないので、MQ≥10で計算しても中央右のSNPs集中部位のカバレッジはほぼ落ち込んでいない。

参考　Understanding Qualities

bedで出力する。

bamCoverage --bam input.bam -o output.bw -of bedgraph --binSize=10 --minMappingQuality 10 --numberOfProcessors=max

user$ head -20 test.bw

chr 0 50 0

chr 50 70 2

chr 70 90 4

chr 90 100 6

chr 100 110 7

chr 110 120 9

chr 120 130 12

chr 130 140 15

chr 140 150 16

chr 150 170 17

chr 170 180 20

chr 180 190 22

chr 190 210 24

chr 210 220 25

chr 220 250 23

chr 250 280 25

chr 280 290 23

chr 290 300 21

chr 300 310 18

chr 310 320 19

igvtoolsでも可能。

igvtools count input.bam output.wig reference.fa --minMapQuality 30 -w 10

igvtools（リンク）はbrewやconda（conda install -c bioconda -y igvtools）で導入できます。

引用

deepTools: a flexible platform for exploring deep-sequencing data

Fidel Ramírez, Friederike Dündar, Sarah Dieh, Björn A. Grüning, and Thomas Manke

Nucleic Acids Res. 2014 Jul 1; 42(Web Server issue): W187–W191.

mapping quality scoresについて

Mapping short DNA sequencing reads and calling variants using mapping quality scores

Heng Li,1 Jue Ruan,2 and Richard Durbin1,3

Genome Res. 2008 Nov;18(11):1851-8. doi: 10.1101/gr.078212.108. Epub 2008 Aug 19.

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

bamCoverageを使いカバレッジトラックを作成する