macでインフォマティクス

macでインフォマティクス

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

構造変化が起きた部位のマッピング状況を出力する samplot

 

samplotはbamやcramを入力として、SVの起こった領域の図を出力してくれるツール。vcfからの一括描画にも対応しているため、variant call format(VCF)を出力したら、そのままsamplotに送るようなスクリプトを書くことで、推定SV全てを目視で簡単に確認できるようになる(IGVなどのビューアで見ていく作業が必要なくなる)。

 

インストール

依存

  • numpy
  • matplotlib
  • pylab
  • pysam
  • statistics
  • bcftools

依存ライブラリは全てpipで導入できる。

本体  Github

https://github.com/ryanlayer/samplot

git clone https://github.com/ryanlayer/samplot.git 

python samplot/src/samplot.py -h

 $ python samplot.py -h

/Users/user/.pyenv/versions/anaconda2-4.2.0/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.

  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Usage: samplot.py [options]

 

Options:

  -h, --help            show this help message and exit

  -n TITLES             Plot title (CSV)

  -r REFERENCE          Reference file for CRAM

  -z Z                  Number of stdevs from the mean (default 4)

  -b BAMS               Bam file names (CSV)

  -o OUTPUT_FILE        Output file name

  -s START              Start range

  -e END                End range

  -c CHROM              Chromosome range

  -w WINDOW             Window size (count of bases to include), default(0.5 *

                        len)

  -d MAX_DEPTH          Max number of normal pairs to plot

  -t SV_TYPE            SV type

  -T TRANSCRIPT_FILE    GFF of transcripts

  -A ANNOTATION_FILE    bed.gz tabixed file of transcripts

  -a                    Print commandline arguments

  -H PLOT_HEIGHT        Plot height

  -W PLOT_WIDTH         Plot width

  --common_insert_size  Set common insert size for all plots

——

python samplot/src/samplot_vcf.sh -h

$ python samplot/src/samplot_vcf.sh -h

    usage: samplot_vcf.sh OPTIONS -v svs.vcf sample_1.bam sample_2.bam ...

 

    General options:

    -h      Show this message

    -T      Sorted and indexed transcript GFF file

    -A      CSV of sorted and indexed genome annotation BED files

    -o      Output directory (/Users/user/local)

    -O      Output type (default png)

 

    Path options:

    -B      BCFTOOLS path (/Users/user/.pyenv/shims/bcftools)

    -S      SAMPLOT path ()

——

 

 

ラン

まずサンプルデータのSVを描画してみる。data/に3つのbamと3つのCRAMが入っている(NA12878、NA12889、NA12890の一部の領域(家系はこちらを参照))。このbamのchr4:115928726-115931880のカバレッジを描画する。

samplot/src/samplot.py \ -n NA12878,NA12889,NA12890 \
-b samplot/test/data/NA12878_restricted.bam,samplot/test/data/NA12889_restricted.bam,samplot/test/data/NA12890_restricted.bam \
-o 4_115928726_115931880.png \
-c chr4 \
-s 115928726 \
-e 115931880 \
-t DEL
  • -b     Bam file names (CSV)
  • -o     Output file name
  • -c     Chromosome range
  • -s     Start range  
  • -e     End range
  • -t      SV type

3つのbamを-bフラグの後に","でつないで指定している。

出力の4_115928726_115931880.pngは下のような画像になる。

f:id:kazumaxneo:20180428170437p:plain 

カバレッジの分布がプロットされその上に、SVの信号を拾ったペアエンドのリードが線で表示されている。縦軸はインサートサイズも表しており、真ん中のbamはインサートサイズが非常に小さいペアエンド(=> deletionを示唆)がたくさんあるのがわかる。フォントのレンダリングは問題ないが、肝心のグラフのプロットの輪郭が眠いように見える。公式でもぼやけているように見えるので、このような仕様なのかもしれない。"-d 100"をつけると、100だけリードをランダムサンプリングして描画する。全リードを使うより早く実行できる。

 

 

vcfファイルの全コールを描画するならsamplot_vcf.shを使う。使うには、bcgtoolsとsamplot.pyのパスを指定する必要がある。

samplot/src/samplot_vcf.sh -B /user/local/bin/bcftools -S samplot/src/samplot.py -v input.vcf input.bam 

INFO/SVTYPEが定義されているバリアントのみ描画対象となる。順番にJSONファイルと絵(.png)が出力されてくる。

 

 

 公式に、遺伝子領域とリピートマスクされた領域をバーで表示する流れが記載されている。ここでもまとめておく。

まず遺伝子のアノテーションをダウンロードしてbedtoolsでポジションソートする。samtoolsのbgzipとbedtools sort(リンク)を使っている。

wget ftp://ftp.ensembl.org/pub/grch37/release-84/gff3/homo_sapiens/Homo_sapiens.GRCh37.82.gff3.gz


bedtools sort -i Homo_sapiens.GRCh37.82.gff3.gz | bgzip -c > Homo_sapiens.GRCh37.82.sort.gff3.gz
tabix Homo_sapiens.GRCh37.82.sort.gff3.gz

 

ゲノムのアノテーションとrepeatmaskerのリピート領域テキストファイルをダウンロードして、それぞれbed.gzに変換する。

wget http://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDukeMapabilityUniqueness35bp.bigWig
bigWigToBedGraph wgEncodeDukeMapabilityUniqueness35bp.bigWig wgEncodeDukeMapabilityUniqueness35bp.bed
bgzip wgEncodeDukeMapabilityUniqueness35bp.bed
tabix wgEncodeDukeMapabilityUniqueness35bp.bed.gz

#repeatmasker track
curl http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/rmsk.txt.gz \
| bgzip -d -c \
| cut -f 6,7,8,13 \
| bedtools sort -i stdin \
| bgzip -c > rmsk.bed.gz

 -Tで遺伝子のnogff3を指定し、-Aでリピートtrackとbedを指定する。

samplot/src/samplot.py \
-n NA12878,NA12889,NA12890 \
-b samplot/test/data/NA12878_restricted.bam,samplot/test/data/NA12889_restricted.bam,samplot/test/data/NA12890_restricted.bam \
-o 4_115928726_115931880.d100.genes_reps_map.png \
-c chr4 \
-s 115928726 \
-e 115931880 \
-t DEL \
-d 100 \
-T Homo_sapiens.GRCh37.82.sort.gff3.gz \
-A rmsk.bed.gz,wgEncodeDukeMapabilityUniqueness35bp.bed.gz 

出力

f:id:kazumaxneo:20180428160020j:plain

一番下に遺伝子、その上にリピート領域としてマスクされた部位が示されている。

 

こちらのツールでも用いられています。

https://github.com/jbelyeu/SV-plaudit

 

引用

GitHub - ryanlayer/samplot: Plot structural variant signals from many BAMs and CRAMs