macでインフォマティクス

macでインフォマティクス

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

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

 

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

 

bamCoverage

bamCoverage — deepTools 2.3.1 documentation 

 

deeptoolsのインストール 

sudo pip install deeptools

#Anaconda環境ならpipでなくcondaを使う。
conda install -c bioconda deeptools

インストール後にbamCoverageを実行すると、python2.7環境であったが

 Reason: Incompatible library version: pyBigWig.so requires version 9.0.0 or later, but libcurl.4.dylib provides version 7.0.0

 が出た。

python - error: Setup script exited with error: command 'x86_64-linux-gnu-gcc' failed with exit status 1 - Stack Overflow

 を参考に修復した。

 

 実行方法

カバレッジ情報をウィンドウサイズ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:plainmapping 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で導入できます。

 

 

引用

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.