macでインフォマティクス

macでインフォマティクス

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

cBarでプラスミド配列を区別する

 

cBarでFASTAは(論文発表当時では)大規模なトレーニングデータを用いて学習されたメタゲノムなどのデータ(FASTA)中の プラスミドゲノムを区別する方法論。5量体頻度(pentamer frequencies)を元に判定を行う。入力はFASTAは配列。

 

インストール

macOSXではビルドでエラーを吐いたのでcent OSに導入。

 

cBar公式HP (ダウンロードリンクあり)

http://csbl.bmb.uga.edu/~ffzhou/cBar/

 

上記ページからダウンロードし、解凍してビルドする。

unzip cBar.1.2.zip
cd cBar.1.2/
make

テストラン 

mkdir temp #必要
./cBar.pl Synechococcus_elongatus_PCC_7942.fna output.txt

 出力。

user$ cat output.txt 

#SeqID Length Prediction

gi|81230333|ref|NC_007595.1| 46366 Plasmid

gi|81298811|ref|NC_007604.1| 2695903 Chromosome

プラスミド1つとクロモソーム1つを持つラン藻ゲノムデータを入力として、plasmidとクロモソームの配列が予測された。

 

 

ラン

 モデルシアノバクテリア(S.6803)でテスト。

#SeqID  Length  Prediction

chr     3573470 Chromosome

plasmid_pSYSA   103307  Plasmid

plasmid_pSYSG   44343   Plasmid

plasmid_pSYSM   119895  Plasmid

plasmid_pSYSX   106004  Plasmid

PCC5.2  5214    Chromosome

pCA2.4  2378    Chromosome

pCB2.4  2345    Plasmid

 5.2 kbと2.3 kbのスモールプラスミド2つはchromosomeと予測されてしまったが、他は正しい。

 

 

 

 

 

 

 

引用

cBar: a computer program to distinguish plasmid-derived from chromosome-derived sequence fragments in metagenomics data.

Zhou F, Xu Y.

Bioinformatics. 2010 Aug 15;26(16):2051-2. doi: 10.1093/bioinformatics/btq299. Epub 2010 Jun 10.

 

異なるk-merの割合を計算し、エラー率推定やゲノムサイズ推定に使える KmerStream

 

KmerStreamは異なるk-merの数を計算する方法論。シーケンス業者のクオリティに依存せず純粋にk-merの頻度からエラー率を見積もることができるため、上手く使えばシーケンスの品質管理などに使用することができる。サンプリングを行うためメモリ使用量が少ないのも特徴となっている。

 

公式HP 

https://bioweb.pasteur.fr/packages/pack@KmerStream@1.1

 

インストール

Github

https://github.com/pmelsted/KmerStream

 

brewでインストールできる。またはgitからソースをダウンロードしてmakeする。

brew install Kmerstream

KmerStream #動作確認

user $ KmerStream

Need to specify files for input

KmerStream 1.1

 

Estimates occurrences of k-mers in fastq or fasta files and saves results

 

Usage: KmerStream [options] ... FASTQ files

 

-k, --kmer-size=INT      Size of k-mers, either a single value or comma separated list

-q, --quality-cutoff=INT Comma separated list, keep k-mers with bases above quality threshold in PHRED (default 0)

-o, --output=STRING      Filename for output

-e, --error-rate=FLOAT   Error rate guaranteed (default value 0.01)

-t, --threads=INT        SNumber of threads to use (default value 1)

-s, --seed=INT           Seed value for the randomness (default value 0, use time based randomness)

-b, --bam                Input is in BAM format (default false)

    --binary             Output is written in binary format (default false)

    --tsv                Output is written in TSV format (default false)

    --verbose            Print lots of messages during run

    --online             Prints out estimates every 100K reads

    --q64                set if PHRED+64 scores are used (@...h) default used PHRED+33

 

ラン

k=33,127

KmerStream -k 33,127 -t 8 -o output

出力

user $ cat output 

Q = 0, k = 33

F0 = 6251288

f1 = 2380496

F1 = 69409025

Q = 0, k = 127

F0 = 7626737

f1 = 3738549

F1 = 44803640

 

 

付属のスクリプトを使いゲノムサイズを推定

KmerStream -k 33,127 -t 8 --tsv -o output.tsv
KmerStreamEstimate.py output.tsv #出力されたtsvファイルを指定

 出力

KmerStreamEstimate.py output 

Q k F0 f1 F1 F0-f1 G e lambda

0 33 6233464 2364673 69409025 3868791.0 3861418.52635 0.0342813550606 17.9750069893

0 127 7630471 3741247 44803640 3889224.0 3885642.33901 0.0836912695204 11.5305620258

F0、F1、の他に、33merと99merから推定されたゲノムサイズがプリントされている。ここではバクテリアのシーケンスデータ(fastq)を使用した。リピート遺伝子を100程度持つ3.95Mの生き物なので、だいたいあってる。

 

 

phread quality score≥30以上のk-merだけ使い計算。

KmerStream -k 33 -t 8 -o output -q 30 input.fq

 

 

 

引用

KmerStream: streaming algorithms for k-mer abundance estimation.

Melsted P, Halldórsson BV.

Bioinformatics. 2014 Dec 15;30(24):3541-7. doi: 10.1093/bioinformatics/btu713. Epub 2014 Oct 28.

 

参考資料

http://tech.plaid.co.jp/streaming-algorithm/

 

k-mer出現頻度を高速計算するntCard

DSK、KmerStream、Khmer、kmerGenieなどより高速に動作するk-merカウントの方法論。原理は大きく異なるが、論文中での上記ツールとの比較では、kmerGenieより100倍以上高速に処理できている。

 

インストール

brewでインストールできる。

brew install ntcard
ntcard --hep #動作確認。

user $ ntcard --help

Usage: ntcard [OPTION]... FILES...

Estimates the number of k-mers in FILES(>=1).

Accepatble file formats: fastq, fasta, sam, bam, gz, bz, zip.

 

 Options:

 

  -t, --threads=N use N parallel threads [1]

  -k, --kmer=N the length of kmer 

  -c, --cov=N the maximum coverage of kmer in output [64]

  -p, --pref=STRING the prefix for output file name [freq]

      --help display this help and exit

      --version output version information and exit

 

ラン

入力はfastq、fasta、sam、bamなど。

 

4スレッド使い32、64、96、128-merの頻度を計算。

ntcard -k32,64,96,128 input.fastq -t4
  • -k the length of kmer 

 -pを省略すると頻度ごとに別ファイル出力される。freq~というファイル名になる。

f:id:kazumaxneo:20171217114554j:plain

 -cで整数の指定がなければ、頻度が64までプリントされる。

F1はトータル、F0はユニークなk-merの合計と定義される。f1は頻度が1。 詳細はこの辺りを参考に

https://github.com/bcgsc/ntCard/issues/1

 

 

6スレッド使い複数フォーマットの同時計算。

ntcard -k64 -c100 -t6 file_1.fq.gz file_2.fa file_3.sam file_4.bam file_5.fq
  • -c the maximum coverage of kmer in output [64]
  • -t use N parallel threads [1]

 複数ファイルの結果は合算して出力される。

 

引用

ntCard: a streaming algorithm for cardinality estimation in genomics data.

Mohamadi H, Khan H, Birol I

Bioinformatics. 2017 May 1;33(9):1324-1330. doi: 10.1093/bioinformatics/btw832.

 

高速なbam/samの解析ツール Sambamba

 

Sambambaはsam、bam、cramの処理ツール。D言語で構築されている。フォーマットを変えたり、フィルタリングすることができる。SAMToolsやPicard-toolsの一部機能と重複するが、Sambambaは並列化に対応しており、SAMToolsより高速に動作するとされる。特にmpileupなどは並列化の恩恵が大きいと思われる。分割して平行処理しているだけであるが、今まで数日かかっていた処理が数時間で終わる。

 

公式サイト

http://lomereiter.github.io/sambamba/

SAMtoolsとの比較

https://github.com/biod/sambamba/wiki/Comparison-with-samtools

 Github

https://github.com/biod/sambamba 

コマンド

https://github.com/biod/sambamba/wiki/Command-line-tools

フィルタリングの仕方

https://github.com/biod/sambamba/wiki/%5Bsambamba-view%5D-Filter-expression-syntax

 

インストール

依存

  • samtools (pileupで使う)

本体

Githubでバイナリが配布されている。brewでもインストールできる。

brew install sambamba
sambanba view #動作確認。例えばviewのヘルプを見る。

 

 

ラン

viewリンク 特定のクロモソームや特定の領域にアライメントされたリードだけ抽出したり、数を数える。

bamの基本情報

sambamba view --reference-info input.bam
  •  -I, --reference-info  output to stdout only reference names and lengths in JSON 

samの基本情報

sambamba view -S --reference-info input.bam
  • -S specify that input is in SAM format

bamのマッピングクオリティ≥50のリード数

sambamba view -c -F "mapping_quality >= 50" input.bam
  • -c output to stdout only count of matching records, hHI are ignored  
  • -F set custom filter for alignments 

bamのリファレンス"chr19"に適切にアライメントされるリードだけsam出力。(カウントなら-cをつけ、-o以下を消す)

sambamba view -F "proper_pair" input.bam chr19 -t 8 -f sam -o output.sam
  • -o specify output filename
  • -t maximum number of threads to use
  • -f specify which format to use for output (default is SAM)

samのリファレンス"chr1"に適切にアライメントされるリードだけbam出力。

sambamba view -S -F "proper_pair" input.sam chr1 -t 8 -f bam -h -p -o output.bam
  • -h print header before reads (always done for BAM output)
  • -p show progressbar in STDERR

"proper_pair"の割合を知りたければ、下のflagstatが便利です。 特定の領域にオーバーラップしたリードを数えるなら、sliceの方が高速と書かれています(sliceリンク)。

 

 

sortリンク ソート。

指定がなければindexが自動でつくcoordinate-sorted。

sambamba sort input.bam -o sorted.bam -t 20 -p -l 6
  • -o output file name; if not provided, the result is written to a file with .sorted.bam extension
  • -t use specified number of threads
  • -p show progressbar in STDERR
  • -l level of compression for sorted BAM, from 0 to 9(指定がなければ"6"くらいのサイズになる)

名前でソート。(indexは付けられない)

sambamba sort input.bam -o sorted.bam -t 20 -p -l 5 
  • -o output file name; if not provided, the result is written to

メモリが潤沢にあるなら、-mでメモリの最大使用量をあげると、I/Oが改善され高速化されるようです。でHDD使ってるサーバー等なら効果が大きそうです(未検討)。

 

 

flagstatリンク bamの情報表示。

bamの情報表示。

sambamba flagstat input.bam -t 8 -p
  • -t use NTHREADS for decompression
  • -p show progressbar in STDERR

sambamba flagstat input.bam

509280 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

538 + 0 supplementary

0 + 0 duplicates

456058 + 0 mapped (89.55%:N/A)

508742 + 0 paired in sequencing

254371 + 0 read1

254371 + 0 read2

451458 + 0 properly paired (88.74%:N/A)

454670 + 0 with itself and mate mapped

850 + 0 singletons (0.17%:N/A)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

表示された情報については上のflagstatのリンク先で確認してください。bamのBGZFフォーマットしか対応していないのでsamはviewで変換してから使ってください。 

 

 

indexリンク bamのindexファイル作成

sambamba index input.bam -t 8 -p
  • -t use NTHREADS for decompression
  • -p show progressbar in STDERR

 input.bam.baiができる。

 

 

mergeリンク samtoolsのmpileupのparallelバージョン。

12スレッド使い、圧縮率最大で3つのbamを1つにマージ。

sambamba merge merge.bam input1.bam input2.bam nput3.bam -t 12 -p -l 9
  • -t Specify number of threads to use for compression/decompression tasks.
  • -l Specify compression level of output file as a number from 0 to 9
  • -p  Show progressbar in STDERR.

 

 

mpileupリンク 

20スレッド使い、pileupを実行。samtoolsの条件は--samtoolsをつけて書く。

sambamba mpileup input.bam -t 20 -p samtools -gu -S -D -d 1000 -L 1000 -m 3 -F 0.0002
  • -t Specify number of threads to use.
  • -p  Show progressbar in STDERR.
  • -o specify output filename

 

 

depthコマンドも現在準備中とのことです。 

 

 

引用

Sambamba: fast processing of NGS alignment formats.

Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P.

Bioinformatics. 2015 Jun 15;31(12):2032-4. doi: 10.1093/bioinformatics/btv098. Epub 2015 Feb 19.

 

 

SVPVで構造変化の予測結果を可視化して比較する

 

SVPVは 構造変化の検出結果のvcfファイルを読み込んで、異なる構造変化の検出ツールでの解析結果を比較できるツール。特定の条件でフィルタリングすることが可能になっている。GUIのアプリでも提供されている。

 

 

オーサーらが作成したSVの解説wiki

https://github.com/VCCRI/SVPV/wiki

 

 

インストール

依存

GUIモード

 GraphicsMagickはbrewでインストール可能。

 

 

本体

https://github.com/VCCRI/SVPV

git clone https://github.com/VCCRI/SVPV.git SVPV

Ubuntuかcent OSならUbuntu_set_up.shかCentOS_set_up.shを実行すればSAMtoolsやbcftoolsなどの依存も含め、まとめてインストールできる(sudoで実行するので、インストール前に必ず実行内容を確認)。 

簡単にインストールできそうだったので、Ubuntu14.04にインストールした。

sudo sh ./SVPV/set_up/Ubuntu_set_up.sh

下に書いたテストデータでのランを行なって動作確認。

 

ラン

現在のところ、SVタイプとして(DEL),  (DUP),  (CNV),  (INV), (INS)  ('BND'). Delly2-style translocations (TRA) がサポートされている。

 

python ./SVPV/SVPV -example #動作確認
python ./SVPV/SVPV -example -gui #GUIモード動作確認

exampleを実行するとsvpv_output/ができて、その中にSVタイプごとにサブフォルダができる。DEL/に移動する。 

f:id:kazumaxneo:20171214104739j:plain

 chr1の2つの領域がまとめられている(chr1上でコールされたDELはこの2箇所)。どちらかのサブフォルダに移動する。

f:id:kazumaxneo:20171214110053j:plain

3つのbamのフォルダと比較部位のpdfファイルがある。

f:id:kazumaxneo:20171214110208j:plain

3つのbamが上下に並べて比較される。もともとexampleには3つのbamファイルとそのSVが入っているので、画像出力も縦に3つの結果が並べて表示されている(下は2つだけ表示)。

f:id:kazumaxneo:20171214104751j:plain

CNVnatorとDelly、LumpyがこのDELを予測していた。

 

GUIモードで起動する。

python ./SVPV/SVPV -example -gui

ウィンドウが出現する。

f:id:kazumaxneo:20171214110454j:plain

 

左上のbamを3つとも選択(selectボタン)、1/1のホモのSVだけ選択。

f:id:kazumaxneo:20171214110631j:plain

フィルターやサイズの制限は設けずに、apply  filterをクリック。

f:id:kazumaxneo:20171214111238j:plain

コール部位はゼロだった。

 

今度はGenotypeやフィルターやサイズの制限は設けずに、chr13のDELをクリック。

f:id:kazumaxneo:20171214110821j:plain

Get Genotypesを押すと、一番下のウィンドウに選択中のSVの情報が出る。Plot all SVsを押すと、CUIと同じようにウィンドウ中の全データが描画される。Plot Selected SVVを押してこのDELだけ描画。

f:id:kazumaxneo:20171214111632j:plain

下の2データだけヘテロの欠失が検出されている。 画面の見方については、wikiを参照してください。

 

このようにGUIモードは選んだbamの特定のSVだけフィルターをかけながら描画できるので、GUIモードの方が使いやすいと思います。

 

 

引用

SVPV: a structural variant prediction viewer for paired-end sequencing datasets.

Munro JE, Dunwoodie SL, Giannoulatou

Bioinformatics. 2017 Jul 1;33(13):2032-2033. doi: 10.1093/bioinformatics/btx117.

 

バクテリアの保存されたgene clusterを探し、結果をビジュアル表示する Gecko3

Gecko3は複数ゲノムを比較して、保存された遺伝子クラスターを検出する方法論。ユーザーが指定した特定の遺伝子群について関連のある遺伝子や遺伝子クラスターを検索することができるSTRINGなどのデータベースと異なり、Gecko3は調べたい生物群の全遺伝子を入力値として、推定遺伝子クラスターをp-value付きで全て検出する。Gecko3はGecko2の改良版で、Gecko2よりメモリ使用量が大きく削減されている。GUI版とコマンドライン版がある。ここではGUI版を紹介する。

 

テストデータ: クラスター探索済みのデータと探索前のデータがある(リンク先の中盤にある "Sample data~")

Gecko3 | Lehrstuhl Bioinformatik Jena

Synechocystis sp. PCC6803のクラスター探索済みデータがあるので、S6803やそれに近いラン藻の遺伝子クラスターを探したいなら、このデータをそのまま使って解析をスタートすることができます(S6803とSTRINGデータベースから関連が見いだされた677ゲノムとの比較 => リンク先のSynechocystis.Default.zip

 

インストール

依存

  • java7

本体のダウンロード

https://bio.informatik.uni-jena.de/software/gecko3/ 

 

 

ラン

./Gecko3

実行するとwindowが立ち上がる。

 

f:id:kazumaxneo:20171211211154j:plain

ダウンロードしたテストデータstring_5genomes_default.gckを開く。

f:id:kazumaxneo:20171211211252j:plain

 

5ゲノム分のgeneが並んで表示されている。

f:id:kazumaxneo:20171211211404j:plain

虫眼鏡のアイコンで拡大縮小ができる。

f:id:kazumaxneo:20171211211443j:plain

 

中央の青い矢印アイコンをクリックしてcluster探索を開始する。

f:id:kazumaxneo:20171211211617j:plain

デフォルトでは以下のようなパラメータ設定になっている。クラスターと定義する最小遺伝子数(トータル)、どの程度近い遺伝子をクラスターと見なすか、などを設定することができる。

f:id:kazumaxneo:20171211211557j:plain

OKで比較開始。ゲノム数が少ないので、10秒程度で計算は終わる。

 

5つのクラスターが検出された。

f:id:kazumaxneo:20171211212153j:plain

 

5遺伝子見つかるクラスターを選択 (ID番号をダブルクリック)。 

f:id:kazumaxneo:20171211212257j:plain

ATP合成酵素のF0F1のサブユニットをコードする遺伝子群だった。

ゲノムが増えるとそれなりの時間がかかります。

 

入力ファイルのフォーマットがやや特殊です。新規に検索セットを作成する場合はREADMEを読んでみてください。

 

オペロンを探すなら、DOORなどが利用できます(予測精度90%以上らしい)。

=> DOOR検索

 

引用

Finding approximate gene clusters with GECKO 3

Sascha Winter, Katharina Jahn, Stefanie Wehner, Leon Kuchenbecker, Manja Marz, Jens Stoye and Sebastian Bo ̈cker

9600–9610 Nucleic Acids Research, 2016, Vol. 44, No. 20 Published online 26 September 2016 doi: 10.1093/nar/gkw843

 

Efficient Computation of Approximate Gene Clusters Based on Reference Occurrences

Katharina Jahn.

Journal of Computational Biology. September 2011, 18(9): 1255-1274. https://doi.org/10.1089/cmb.2011.0132

 

Gecko2

http://bibiserv2.cebitec.uni-bielefeld.de/gecko

 

javaの複数バージョンの使い分けについて

複数バージョンのJavaを切り替えて使用できるようにする。 - Qiita

 

 

Mugsyでゲノムのマルチプルアライメントを行う

Mugsyはnucmerを内部で動かし、all against allのペアワイズアライメントを行い、ゲノムサイズのマルチプルアライメントを可能にする方法論。論文では31のバクテリアゲノムを2時間以内に解析できたと記載されている。

 

公式サイト

http://mugsy.sourceforge.net

 

ダウンロード

macでは動作しなかったのでcent OS6に導入した。 

sorceforge

https://sourceforge.net/projects/mugsy/files/

mummerも内部に含まれている。 

中のmugsyenv.shを開き、2行目のパスをmugsyがあるパスに変更し、sourceする。

source mugsyenv.sh
mugsy -h #動作確認

 

ラン

アライメント

mugsy --directory ./output -p mygenomes genome1.fa genome2.fa genome3.fa genome4.fa
  • --directory directory used to store output and temporary files. Must be a absolute path
  • -p prefix for output files

 

maf形式で出力される。

> mygenomes.maf

##maf version=1 scoring=mugsy

a score=1239890 label=1 mult=5

s 1788.1788             873135 1116720 + 1989855 TCAAG--AAT-ATGACAATACAGGGAGTGGAAATTTATAACCTGAAAAGTGGAATGAATAATAGAAACGAAAAAGGCAAGGAGTTTGCCATGATTGGAAAGAACATAAAATCCTTACGTAAAACACATGACTTAACACAACACGAATTTGCACGGATTGTAGGTATTTCACGAAATAGTCTGAGTCGTTATGAAAATGGAACGAGT---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ATTGTGAAAGAAAGAGGTGCTAATCTATTATCTCGACTCTATCGCTATCAAGATAGTCAGGGAATTAGCATTGATGATGAATCTAATCCTTGGATT----------------------------------------T-TAATGAGTGATGATCTTTCTGATTTGATTCATACGA-A-AAT-------------------------------------CTGAAAAGCGGATGGTAGCTTAATGGAAATCCAAGATTATACTGATAGTGAATTCAAACATGCTTTAGCAAGGAATCTTCGTTCACTGACAAGAGGAAAAAAGTCCAGTAAGCAACCTATAGCGATTTTGCTTGGAGGTCAAAGTGGTGCCGGTAAGACTA

mafは例えばmaftoolsなどでパースすることができる(maftools)。

 

テストデータとして、harvest suitsで提供されているデータなどが使えると思います。

https://www.cbcb.umd.edu/software/harvest

 

引用

Mugsy: fast multiple alignment of closely related whole genomes 

Samuel V. Angiuoli Steven L. Salzberg Author Notes

Bioinformatics, Volume 27, Issue 3, 1 February 2011, Pages 334–342