macでインフォマティクス

macでインフォマティクス

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

Kallisto

 

 典型的なRNA-seqの転写産物レベル処理ワークフローの最初の2つのステップは、トランスクリプトーム配列またはリファレンスゲノムへのアラインメントおよび転写産物存在量の推定である。これらのステップには時間がかかることがある。例えば、広く使用されているプログラムTopHat2(論文より ref.1)を使用してそれぞれを30万リードからなる20サンプルのRNA-seqデータのアライメントを行うと、20コアで28CPU時間かかり、コンパニオンプログラムであるCufflinks2で定量するとさらに14時間かかる。そのような実行時間は、サンプルが増加しシーケンシングデータが増えると禁止になる可能性が高い。アライメントされたリードの定量は、ストリーミングアルゴリズム(ref.3)またはリードのnative count(ref.4)によって高速化できるが、これらのアプローチでは定量の精度が低下する。アライメントステップを回避するために、最近の研究では、リードからk-mersを抽出し、続いてハッシュテーブルを使用してk-mersを正確にマッチングさせることによってサンプルを定量することが提案されている(ref.5)。しかし、k-merへの分けると、各k-merがリード自体より多くの転写産物とアライメントすることができるため、完全なリードに存在している貴重な情報を捨てている。その結果、精度が大幅に低下する(論文 supplementary 図1)。
 k-mersの直接使用は正確な定量には不十分だが、ハッシュベースのアプローチはRNA-seq処理を高速化するための基礎を提供している。ここでは、アライメントベースの定量精度を維持しながらk-mers情報を使えないかどうかを調べる。著者らは、正確な定量の難しさとそのために中心的に考える事項を調べる、すなわちユニークにアライメントさせることができないリードのアサインについて考える。典型的には、これらのマルチマッピングリードはRNA-seqの統計モデルを用いて説明され、そのようなリードを確率論的にアサインし、転写産物存在量の最尤推定を行う事で説明する。しかし、正確な定量だけならリードが転写産物の内部のどこから由来するのかについての情報は必要なく、リードがどの転写産物に由来するかだけわかれば良いことが示されている(ref.7)。この情報に基づいて、著者らは、リードと断片の擬似アライメントに基づいた方法を開発する。この方法は、リードが生じた可能性のある転写産物の特定にのみ焦点を当て、リードが転写産物上にどのようにアライメントされているのか特定しようとしない。

 転写産物へのリードの正確な擬似アライメントは、 k-merの高速ハッシングとtranscriptome de Bruijn graph(T-DBG)を用いて得られた。 de BruijnのgraphはDNAとRNAアセンブリで重要な役割を果たした。 Kallistoは、トランスクリプトーム(論文 図1a)に存在するk-merから構築されたde Bruijn graphであるT-DBGと、graphのパスカバーリング(グラフのすべてのエッジをカバーするパスの集合)を使用する。経路は転写産物に対応する(論文 図1b)。このT-DBGの経路被覆は、k互換性クラスと呼ばれる頂点上のマルチセットを誘導する。互換性クラスは、graph内のパスとしてそれを表現し、graph中のパスのk互換性クラスを、その構成要素のk互換性クラスの交差点として定義することによって、リードに関連付けることができる(論文 図1c)。読み取りのための等価クラスは、読み取りに関連付けられた転写物の複数のセットである。

(以下省略)

Githubより

kallistoは、RNA-Seqデータ、またはより一般的にはハイスループットシーケンシングリードを用いて転写産物の量を定量化するためのプログラムである。 kallistoは、アライメントの必要なしにターゲットとリードの互換性を迅速に決定するため擬似アライメントという新しいアイデアを使っている。 kallistoは、標準的なRNA-Seqデータのベンチマークでは、Macデスクトップコンピュータでトランスクリプトームインデックスのみを使用して3分以内にヒトの3,000万トランスクリプトームシーケンシングリードを定量できる。擬似アライメントは定量に必要な重要な情報を保持するため、kallistoは高速であるだけでなく、既存の定量ツールと同程度精度が高い。 実際、擬似アライメントの手順はリードのエラーに強く、多くのベンチマークではkallistoが既存のツールを大幅に上回る。 

 

f:id:kazumaxneo:20180714130046p:plain

図。3000万リードの定量にかかる時間。論文より転載。

 

 

公式HP

f:id:kazumaxneo:20180714130512p:plain

マニュアル

http://pachterlab.github.io/kallisto/manual.html

Kallistoをベースに作られたツール

https://pachterlab.github.io/kallisto/apps

 

インストール

Github

brewやconda(リンク)で導入できる。ここではbrewで導入する。

brew tap brewsci/science
brew tap brewsci/bio
brew install abyss samtools bwa biobamtools

 > kallisto

$ kallisto 

kallisto 0.44.0

 

Usage: kallisto <CMD> [arguments] ..

 

Where <CMD> can be one of:

 

    index         Builds a kallisto index 

    quant         Runs the quantification algorithm 

    pseudo        Runs the pseudoalignment step 

    h5dump        Converts HDF5-formatted results to plaintext

    inspect       Inspects and gives information about an index

    version       Prints version information

    cite          Prints citation information

 

Running kallisto <CMD> without arguments prints usage information for <CMD>

 

 

ラン

1、transcripts.faのindexを作成。

kallisto index -i transcripts_index transcripts.fa
  • -i    Filename for the kallisto index
  • -k   k-mer (odd) length (default: 31, max value: 31)
  • --make-unique    Replace repeated target names with unique names

ジョブが終わると transcripts_indexができる。

 

2、ペアエンドの定量。非圧縮fastqとgzip圧縮fastqに対応している。

kallisto quant -t 8 -i transcripts_index -o output pair1.fq.gz pair2.fq.gz
  • -i     Filename for the kallisto index to be used for quantification
  • -o    Directory to write output to
  • -t     Number of threads to use (default: 1)
  • --bias    learns parameters for a model of sequences specific bias and corrects the abundances accordlingly.
  • --fr-stranded    runs kallisto in strand specific mode, only fragments where the first read in the pair pseudoaligns to the forward strand of a transcript are processed. If a fragment pseudoaligns to multiple transcripts, only the transcripts that are consistent with the first read are kept.
  • --fusion    does normal quantification, but additionally looks for reads that do not pseudoalign because they are potentially from fusion genes. All output is written to the file fusion.txt in the output folder.

outputディレクトリができる。出力は以下のように説明されている。

  1. abundances.h5 is a HDF5 binary file containing run info, abundance esimates, bootstrap estimates, and transcript length information length. This file can be read in by sleuth
  2. abundances.tsv is a plaintext file of the abundance estimates. It does not contains bootstrap estimates. Please use the --plaintext mode to output plaintext abundance estimates. Alternatively, kallisto h5dump can be used to output an HDF5 file to plaintext. The first line contains a header for each column, including estimated counts, TPM, effective length.
  3. run_info.json is a json file containing information about the run

テストデータの定量結果(abundances.tsv)。

> cat test/output/abundance.tsv

$ column -t test/output/abundance.tsv 

target_id          length  eff_length  est_counts  tpm

ENST00000513300.5  1924    1746.98     102.328     11129.2

ENST00000282507.7  2355    2177.98     1592.02     138884

ENST00000504685.5  1476    1298.98     68.6528     10041.8

ENST00000243108.4  1733    1555.98     343.499     41944.9

ENST00000303450.4  1516    1338.98     664         94221.8

ENST00000243082.4  2039    1861.98     55          5612.36

ENST00000303406.4  1524    1346.98     304.189     42908.2

ENST00000303460.4  1936    1758.98     47          5076.85

ENST00000243056.4  2423    2245.98     42          3553.05

ENST00000312492.2  1805    1627.98     228         26609.9

ENST00000040584.5  1889    1711.98     4295        476675

ENST00000430889.2  1666    1488.98     623.628     79578.2

ENST00000394331.3  2943    2765.98     85.6842     5885.85

ENST00000243103.3  3335    3157.98     962         57879.3

 

シングルエンドの定量。リード長に関する情報を-l-sで指定する必要がある。

kallisto quant -t 8 -i transcripts_index -o output --single -l 200 -s 20 single.fq
  • --single    Quantify single-end reads
  • -l    Estimated average fragment length
  • -s    Estimated standard deviation of fragment length (default: -l, -s values are estimated from paired end data, but are required when using --single)

 

複数のサンプルがある場合は順番に記載する。

ペアエンド3データの定量

kallisto quant -t 8 -i transcripts_index -o output \
pairA1.fq pairA2.fq pairB1.fq pairB2.fq pairC1.fq pairC2.fq

 

Visualization

--pseudobamや--genomebamをつけると、擬似アライメント結果をbamで書き出すことができる(--genomebamはゲノムに対して擬似アライメントするのでtranscriptsのgftもランに必要)。擬似アライメントなので、一部の情報はarbitraryな固定値がプリントされる。詳細は別ページに書かれている(リンク)。意味は薄いが、ソートすればviewerで確認することはできる。

#githubのtestデータの定量結果(--pseudobamをつけて定量した)
cd output/
samtools sort -@ 12 -m 8G pseudoalignments.bam -o pseudoalignments_sorted.bam
samtools index pseudoalignments_sorted.bam

igv -g ../transcripts.fasta pair_sorted.bam

あくまで擬似アライメントであることに注意する(ミスマッチやindel、clippinが全く発生していない)。

f:id:kazumaxneo:20180714172313p:plain

--genomebamをつけてのラン例はこちら

 

 

引用

Near-optimal probabilistic RNA-seq quantification.

Bray NL, Pimentel H, Melsted P, Pachter L

Nat Biotechnol. 2016 May;34(5):525-7

 

参考ページ 

ソースからのビルド手順も説明されています。

 

k-merを使い 進化距離や相同性を高速計算する Kmer-db

Preprintより 

 何千もの異なる生物のシーケンシング解析の過程で大量のデータが生成された(100K Pathogen Genome Project(Weimer el al、2017、NCBI Pathogen Detection(https://www.ncbi.nlm.nih.gov/ pathogens) )、これは迅速な分析方法を要求する。 k-merと呼ばれる短いヌクレオチド配列のサブストリングは、ゲノムやシークエンシングリードのいずれかから抽出でき、アセンブリフリーアプローチを可能にするため、この領域で一般に使用されている。これらは、生物間の進化的距離の正確な近似を可能にするため、系統再構成(Mash(Ondov et al。、2016)(MASH紹介))、バクテリアの同定(StrainSeeker(Roosaare et al、2017)(StrainSeeker紹介))、またメタゲノム分類MetaCache(Muller et al、2017))などに利用されている。重要なことに、ゲノムが十分に近い場合、許容可能な精度を得るにはk-merの小さなサブセットで十分であり、処理時間が大幅に短縮されることである。それとは関係なく、シーケンシングされたゲノムの数と多様性が継続的に増加しているので、既存のアルゴリズムスループットは間もなくボトルネックに達するであろう。
 著者らは、シーケンシングされた大量のサンプルコレクションのk-merベースの解析を行うツールであるKmer-dbを紹介する。新規の圧縮されたk-mers表現と並列実装のおかげで、本ソフトウェアは一般的なワークステーションマシンを使って数千のバクテリアゲノムを数分で処理することができる。 

 

公式HP

f:id:kazumaxneo:20180714114310p:plain

 Kmer-dbに関するツイート。

 

インストール

ubuntu18.04に導入した。

依存

  • KMC (Optional) 

KMCはAnaconda環境ならmaclinuxともにcondaで導入可能(リンク)。またはバイナリをダウンロードしてパスの通ったディレクトリに移動する。

conda install -c bioconda kmc

 本体 Github

リリースからlinux用のバイナリをダウンロードする。

https://github.com/refresh-bio/kmer-db/releases

"kmer-db"とリネームして、パスの通ったディレクトリに移動した。

 

 

ラン 

1、データベースの構築

3つの方法に対応している。"-f"でサンプリングレートを調整できる(defaultは1)。

方法1、 ゲノムFASTAからデータベースを構築する。

kmer-db build -k 18 sample_list_file test_database
  • サンプルリストファイルは1行に1ファイルずつゲノムのFASTAファイルを記載したテキストファイル。拡張子は.fnaか.fna.gzまたは.gz。例えばゲノムファイルがsample1.fnaなら"sampke1"だけ記載すれば認識するように書かれていたが、実際はほとんどん認識せず、実行すると途中のログにfailedが出た(エラーでも普通にランは終わるので注意)。成功したのは、例えば"test.fna"ファイルならgz圧縮して"test.fna.gz"を作成し、 sample_list fileには"test.fna"と書いた時のみ workし、"test"だとfailだった(version1.0でテスト)。
  • リストファイルと.fnaファイルが同じディレクトリにないと機能しない。違う時はフルパスで記載する(リンク)。
  • ランが終わると、指定した名前のファイルができる。ここではtest_database。

 方法2、 KMC(紹介)を使ってk-merファイルからデータベースを構築する(テスト時はエラーを起こしました。方法1と同じかと思いましたが、もしかするとサンプルリスト表記方法が違うかもしれません)。

#先にtempoを作らないとエラーになる
mkdir working_directory

#k=20のカウント。fastaを使用する。
kmc -fa -k20 -m24 -ci0 input.fasta database working_directory/
#-k k-mer length (k from 1 to 256; default: 25)
#-m max amount of RAM in GB (from 1 to 1024); default: 12
#-ci exclude k-mers occurring less than <value> times (default: 2)
#-f input in FASTA format (-fa), FASTQ format (-fq) or multi FASTA (-fm); default: FASTQ


#データベース構築
kmer-db build-kmer kmc_sample_list_file database
  • 論文のsupplementary dataでは、one2allを実行する際、シーケンスエラーを考慮して"-ci5"のオプションをつけてkmcを実行している(supplementary data(リンク))。

方法3、 minhash k-merファイルから データベースを構築する。

Githubと実行時のヘルプコマンドを参照。

 

2-A、全データベースサンプルのsimilarities/distancesを調べる(all vs all)。

  • all2all - calculating number of common k-mers between all samples in the database
  • distance - calculating similarities/distances.
#all2allコマンンド実行
kmer-db all2all test_database matrix.csv

#matrix.csvができる。

#distanceコマンド実行。
kmer-db distance matrix.csv
#matrix.csv.min、matrix.csv.max、matrix.csv.mash、matrix.csv.jaccard、matrix.csv.cosineができる。

 

2-B、ある配列(ゲノムのfastaやシーケンスデータのfastq)とデータベース配列のsimilarities/distancesを調べる。

  • one2all - calculating number of common kmers between single sample and all the samples in the database,
#salmonella.fastaとtest_databaseのゲノムとの相同性/距離を調べる。

#KMCを使い、salmonellaのk-merファイルを作成
mkdir working_directory #先にtempoを作らないとエラーになる
kmc -fa -k20 -m24 -ci1 salmonella.fasta salmonella working_directory/
#salmonella.kmc_sufとsalmonella.kmc_preができる。

#one2allコマンンド実行
kmer-db one2all test_database salmonella vector.csv
#vector.csvができる

#distanceコマンド実行。
kmer-db distance vector.csv
#vector.csv.min、vector.csv.max、vector.csv.mash、vector.csv.jaccard、vector.csv.cosineができる。


#ゲノムではなくシーケンスデータとデータベース間で進化距離を調べるなら、kmc実行時のパラメータを変更してください。
kmc -fq -k20 -m24 -ci5 input.fq output_name working_directory/
kmer-db one2all test_database output_name vector.csv
kmer-db distance vector.csv

 

まだAcceptされて間もないためfull textが読めてないが、supplementaryを読む限り、論文ではNCBIのPathogen Detectionデータベース(リンク)に登録されたゲノムをダウンロードしてall2allのコマンドの結果を評価している。one2allについては、構築したpathogenデータベースを使いSRAに登録されたいくつかの病原性バクテリアのシーケンスデータを調べている。

 

引用

Kmer-db: instant evolutionary distance estimation

Sebastian Deorowicz Adam Gudyś Maciej Długosz Marek Kokot Agnieszka Danek

Bioinformatics, bty610

 https://academic.oup.com/bioinformatics/advance-article-abstract/doi/10.1093/bioinformatics/bty610/5050791?redirectedFrom=PDF

 Preprint

https://www.biorxiv.org/content/biorxiv/early/2018/02/12/263590.full.pdf

 

参考

What is the Pathogen Detection project?

FAQ - Pathogen Detection - NCBI

 

 

Fermi-lite

 

 Fermi-liteはHeng Liさん(wiki)がGithubで公開されているスタンドアロンのCライブラリ。イルミナのショートリードを100bpから1,000万bpの領域でアセンブリするコマンドラインツールとしても機能する。Fermi-liteはオーバーラップベースのアセンブリを行う。 全データを中間ファイルを生成することなくメモリにキャッシュして働く、Fermikit(Heng Li,  2015)の軽量バージョンとなっている。 パフォーマンス、比較的小さなメモリフットプリント、およびfermikitの機能を継承している。 Fermi-liteはヘテロ接合部位を保持でき、したがって、二倍体領域を組み立てバリアントコールに使用することができる。最もシンプルなインタフェースとなっている。最終出力はunitig(*1)。ここではコマンドラインの機能を紹介する。

 

Limitations(Githubより転載)

  • Fermi-lite can efficiently assemble bacterial genomes. However, it has not been carefully tuned for this type of assembly. While on a few GAGE-B data sets fermi-lite appears to work well, it may not compete with recent mainstream assemblers in general.

  • Fermi-lite does not work with genomes more than tens of megabases as a whole. It would take too much memory to stage all data in memory. For large genomes, please use fermikit instead.

  • This is the first iteration of fermi-lite. It is still immarture. In particular, I hope fermi-lite can be smart enough to automatically figure out various parameters based on input, which is very challenging given the high variability of input data.

 

 Fermi-liteに関するツイート。

 

インストール

GIthub

git clone https://github.com/lh3/fermi-lite
cd fermi-lite && make
#"fml-asm"ができる。パスの通ったディレクトリにコピーしておく。

#テストラン
./fml-asm test/MT-simu.fq.gz > MT.fq

./fml-asm

$ ./fml-asm

Usage: fml-asm [options] <in.fq>

Options:

  -e INT          k-mer length for error correction (0 for auto; -1 to disable) [0]

  -c INT1[,INT2]  range of k-mer & read count thresholds for ec and graph cleaning [4,8]

  -l INT          min overlap length during initial assembly [33]

  -r FLOAT        drop an overlap if its length is below maxOvlpLen*FLOAT [0.7]

  -t INT          number of threads (don't use multi-threading for small data sets) [1]

  -d INT          retain a bubble if one side is longer than the other side by >INT-bp [50]

  -A              discard heterozygotes (apply this to assemble bacterial genomes; override -O)

  -O              don't apply aggressive tip trimming

  -g              output the assembly graph in the GFA format

——

 

ラン

純化したバクテリアのシーケンシングデータをアセンブリする。

#ペアエンドはインターレースにしておく。ペアエンド情報は使わないのでcatを使ってもOK
seqtk mergepe pair1.fq pair2.fq > paired.fq

#de novoアセンブリ。出力はFASTQフォーマットになっている。
fml-asm -A -t 16 paired.fq > out.fq
  • -A    discard heterozygotes (apply this to assemble bacterial genomes; override -O)
  • -t     number of threads (don't use multi-threading for small data sets) [1]

 

2倍体ゲノムをアセンブリし、GFA(*2)出力する。

fml-asm -g -A -t 16 paired.fq > out.gfa
  • -g output the assembly graph in the GFA format

できた.gfaファイルをBandageに読み込ませる。

f:id:kazumaxneo:20180712172343j:plain

Bandageは以前紹介しています。

 

Cのライブラリとして使用する場合はGithubの例を参考にしてください

 

 

引用

FermiKit: assembly-based variant calling for Illumina resequencing data.
Li H.

Bioinformatics. 2015 Nov 15;31(22):3694-6.

 

参考

*1  What is a unitig? How does it differ from a contig?

https://www.quora.com/What-is-a-unitig-How-does-it-differ-from-a-contig

 

*2  ACGT Blog. The Graphical Fragment Assembly (GFA) format

http://www.acgt.me/blog/2014/10/15/the-graphical-fragment-assembly-gfa-format

 

3*  Heng Li Blog. A proposal of the Grapical Fragment Assembly format

A proposal of the Grapical Fragment Assembly format

プライマーを自動設計したり、プライマー全アニーリング部位を確認できる FastPCR

 

 ポリメラーゼ連鎖反応(PCR)は、分子生物学の基本であり、核酸増幅に使用される最も重要な実験技術である。(一部略)
 増幅反応を引き起こすための熱サイクリングに依存しない多くの等温増幅技術も開発されている。 そのような技術の1つは、鋳型DNAをオリゴヌクレオチドプライマーおよび高いストランド置換活性を有するポリメラーゼと混合し、混合物を60〜65℃の一定温度に保持するLAMP(Loop-Mediated Isothermal Amplification)[ref.3]である。他の等温DNA増幅技術も同様に、特定のDNAポリメラーゼのストランド置換活性に依存する。(一部略)
PCR反応の成功は、オリゴヌクレオチドプライマーの正確な設計に大きく依存する。順方向プライマーは第1のDNA鎖にアニールし、逆方向プライマーは第2の相補鎖にアニールする。プライマーがアニーリングするDNA部位が適切な距離で分離されている場合、PCR産物またはアンプリコンとして知られているこれらの部位間のDNA断片はDNAポリメラーゼによってコピーされ、各温度サイクルで約2倍になる。プライマーと鋳型との安定したハイブリダイゼーションはDNAポリメラーゼによるプライマー伸長の絶対的な前提条件であるため、プライマーの正しい選択と検証はPCR反応の成功に不可欠である。
 プライマー特異性は、成功するPCRの最も重要な因子の1つである。適切に設計されたプライマーは、特に複雑なゲノムDNAが鋳型として使用される場合、標的DNA配列にのみアニールする必要がある。理想的には、ターゲットセットからのすべてのDNA配列はプライマーに対して正確に相補的であり、増幅されるが、バックグラウンドセットからのDNA配列はプライマーと一致しない事である。言い換えれば、プライマーは、標的DNA配列のみを増幅するために非常に特異的でなければならない。さらに、プライマー融解温度は、同じ温度で標的DNA配列への適切なアニーリングを可能にするために非常に類似している必要がある。プライマーが反復配列(レトロトランスポゾン、DNAトランスポゾンまたはタンデムリピート(wiki))に相補的である場合に問題が生じる可能性があり[ref.16]、プライマーが逆方向反復配列に相補的である場合には代替増幅産物が生成され得る。しかしながら、逆方向反復の生成は、一般的なDNAフィンガープリント法で利用されている。多くの場合、これらのPCR反応では1つのプライマーのみが使用され、生成物の末端はプライマー配列に相補的な逆方向反復を含む。

  真核生物では、インターレトロトランスポゾン増幅多型(IRAP)などのインターリピート増幅多型技術は、レトロトランスポゾンおよびSINE様配列の長い末端反復(LTR)などの非常に豊富な分散したリピートの存在を利用している[ ref.18]。これらの配列を互いに関連付けることにより、反復配列と相同なプライマーを用いて一連のPCR産物(DNAフィンガープリント)を生成することが可能になる。 LTRから外側を指し示す1つ、2つ以上のプライマーを使用して、2つの近くのレトロトランスポゾン間のDNA配列を増幅する。プライマーは、同じレトロトランスポゾンファミリーまたは異なるファミリー由来の配列と相同であり得る。 PCR産物および対応するフィンガープリントパターンは、ゲノム中の数百から数千の標的配列の増幅から生じる[ref.19,20]。
 DNAシーケンシング技術の開発、特にハイスループットシークエンシング(次世代シーケンシング)の出現により、膨大な量のrawゲノム配列データの蓄積が進んている。現在、多くの原核生物および真核生物のゲノムがシーケンシング解析され、データベースにアノテーションが登録されている。このため、rawシーケンシングデータから有用な情報を抽出し、ツールを使用して情報を処理し、計画段階ですでに実験結果を予測するin silicoアプローチに対する需要が高まっている。そのようなアプローチの1つは、in silico PCRである。in silico PCRの通常の目標は、既に設計されたオリゴヌクレオチドプライマーを用いて、1つまたは複数のDNAテンプレートから合成されるPCR産物を予測することである。他の一般的な目的には、ゲノム中のテンプレート配列の位置、ならびにヘアピン、G-四重鎖、二次/四次構造の融解温度の予測、self-dimersおよびross-dimers in primer pairsなどのプライマー/プローブ設計および分析が含まれる。
 現在、いくつかのWebベースのin silico PCRツールが生物学者に利用可能である[ref.23-26,21,22,27-32]。そのうちの1つであるPRIMEXは、k- mer lookupテーブルを使用して、全ゲノムから短い配列のマッチする部位を検索する[31]。ウェブサーバ「Electronic PCR」(https://www.ncbi.nlm.nih.gov/tools/epcr/)[29]は最大2つのミスマッチを許容してゲノムの検索を可能にする。別のWebサーバであるUCSCの in- silico PCRhttp://genome.ucsc.edu/cgi-bin/hgPcr)は、定義されていないゲノムを検索するアルゴリズムを使用している[ref.24]。最後に、Primer-BLAST [33]は、BLAST [26]を基礎とする検索方法として使用するWebサーバーである(紹介)。
公開されたin silico PCRアルゴリズムのほとんどはヒューリスティックであり、スタンドアロンのソフトウェアパッケージとして利用できるものはほとんどない。注目すべき例外はPRIMEX [ref.31]とWindows用FastPCR [ref.34,21,35]である。さらに、in silico PCRに一般的に使用される配列類似性検索アルゴリズムの適用は完全には成功していない。例えば、BLASTは完全なプライマー配列をカバーすることができないかもしれないローカルアライメントを作成し、可変長の任意の配列によって分離されたクエリペアの検索をサポートしていない。degenerateプライマーはアライメントシードに現れないため、誤ったネガティブを避けるために特別な処理が必要となる。許容可能な一致数を見つけるのに十分に敏感なパラメータは、多数の誤検出を生成する可能性があり、有効なヒットを識別するために広範な後処理が必要となる。 Primer-BLASTはワード長が7であるため、プライマーと7-merの一致が必要になる。例えば5番目および10番目にミスマッチがある長さ15のプライマー検索を考える。最も長い同一k-merは5ntであり、BLASTは一致を見出すことができない。

 優れたin silico PCRプログラムは、5 'または3'テイルや一塩基多型(SNP)を持ったdegenerateプライマーやプローブを取り扱うことができるものである。さらに、メチル化されたシトシンと高度に分解され化学的に修飾されたDNAのみを含有する重亜硫酸塩処理DNAなど古細菌叢およびミイラサンプルのデータともコンパチブルであるべきである。最後に、プログラムは、単一の反応でmultiplex PCR、nested PCRまたはtiling PCRを処理できなければならない。これらの考慮事項は、上記のすべてのタスクを処理できるスタンドアロンソフトウェアとオンラインソフトウェアの両方を開発する動機となった。著者らの目的は、線状または環状DNAテンプレート上でin silico PCR、大規模または小規模のデータベースからの複数のプライマーまたはプローブ検索(論文 表1)から高度な配列分析ができる効率的で使いやすい計算ツールを作成することであった。このソフトウェアは、標的配列からのプライマーまたはプローブの迅速な選択、プライマーの位置、向き、アニーリング効率、および標準オリゴヌクレオチドおよびdegenerateオリゴヌクレオチドのプライマー melting temperatures計算に有用である(以下略)。

 

公式HP

https://primerdigital.com/fastpcr.html

マニュアル(右のリストからジャンプできる)

http://primerdigital.com/fastpcr/

 

In silico PCR against whole genome(s) or a list of chromosomes with jPCR or FastPCR

 

FastPCRに関するツイート


ラン

マニュアルおよびダウンロード

http://primerdigital.com/tools/pcr.html

黄色のLaunchボタンを押してJava webスタートファイル(wiki)をダウンロードする。pcr.jnlpをダブルクリックして起動する。Appストアからダウンロードしていないので起動できないと表示されたら、環境設定から例外URLを追加する。Javaをクリック。

f:id:kazumaxneo:20180712081246p:plain

Webの設定タブに切り替え、追加ボタンを押してhttp://genomeview.orgを追加した。

f:id:kazumaxneo:20180712081329p:plain

適用をクリックし、OKをクリックして閉じる。

 

pcr.jnlpをダブルクリックして起動する。

f:id:kazumaxneo:20180712081708p:plain

 

基本的な機能

FastPCRウィンドウ下部に塩基配列をペーストする。上のメニューボタンから翻訳したり、Complementary、reverse、antisenseなどの変換ができる。

f:id:kazumaxneo:20180712105355p:plain

ツールメニューからも同様に実行できる。

f:id:kazumaxneo:20180712105514p:plain

より高度な機能を見ていく。

 

 

プライマーの自動設計

1、PCR primer designは、指定した配列を対象にプライマーを自動設計できる。FastPCRウィンドウ下部にprimerを設計したい配列(遺伝子など)をペーストする。

f:id:kazumaxneo:20180712095313p:plain

 2、FastPCRウィンドウ上部のPCR primers designタブで長さ、Tm、GCなどを指定する。Special settingsの項目では、Default criteria以外に、Long distance、Quantitative、などを選ぶこともできる。指定した条件でprimerがデザインされる。

f:id:kazumaxneo:20180712095510p:plain

条件が決まったら、▶︎ボタンを押して実行する。

 

3、デザインされたprimer配列が表示される。配列が長ければかなりの数のプライマーが設計される。

f:id:kazumaxneo:20180712095957p:plain

プライマーリストの下には好ましいプライマーペアもまとめられている。

f:id:kazumaxneo:20180712100056p:plain

名前は以下のルールに従ってつけられる。f:id:kazumaxneo:20180712100133p:plain

複数配列も同時にprimer設計できる。

f:id:kazumaxneo:20180712101001p:plain

 

 

設計したPCRプライマーの全アニーリング部位を確認する

1、まずクロモソームのFASTAファイルかGenebankファイルを読み込む(複数同時読み込み可)。File => Open Genebank FASTA files(s)

f:id:kazumaxneo:20180712074728p:plain

10kb程度の配列なら下のウィンドウに直接コピー&ペーストしてOK。

 

2、FastPCRウィンドウ下部のPre-designed primerタブに前もって設計していたPCRプライマー配列をコピー&ペーストする。

f:id:kazumaxneo:20180712075326p:plain

先ほどのprimer自動設計で作った配列をそのままのフォーマットでペーストしても機能する(先にstep1の作業を行い、検索対象をプライマー設計した配列からゲノムに切り替えておく)。

f:id:kazumaxneo:20180712105940p:plain

 ↑ 自動作成したプライマーリストをそのままコピー&ペーストしている。 

 

3、FastPCRウィンドウ上部のin silico PCRタブに繰り替え、実行条件を確認する。切り替えたタブによって実行内容が変わるので必ずin silico PCRのタブに切り替えておく。show all matching for primer alignmentにチェックが付いていることを確認する。環状ゲノムなら5'末端と3'末端にまたがってアライメントされる可能性も考慮する必要があるため、Circular sequenceにチェックをつける。この指定したパラメータ設定に従ってプライマーアニーリング部位が調べられる。

f:id:kazumaxneo:20180712075435p:plain

 

4、FastPCRウィンドウ右上の緑の▶︎Runボタンをクリックして実行する。プライマーの全アニーリング部位がサーチされる。

f:id:kazumaxneo:20180712084252p:plain

 

5、バクテリアのゲノムサイズであれば1-2秒でサーチは終わる。結果はFastPCRウィンドウ下部のIn silico PCR resultに表示される。テストしたforward primerはspecificだったが、reverse primerは 2ヶ所で100%マッチし、1ヶ所で90%マッチしていた。しかし、その中で増幅可能な組み合わせは1つだけだったので、一番下のAmplicon sequenceのプライマーペアは1つだけ表示された。

f:id:kazumaxneo:20180712082247p:plain

 

Primerリストの基本情報

Primer list analysisタブでは、たくさんのprimer配列の簡単なstatisticsを分析できる。万単位のプライマーリストでも機能する。

f:id:kazumaxneo:20180712085051p:plain

 

ファミリー遺伝子それぞれの特異的なプライマーの設計

PCR primer designのGroup specific PCRにチェックをつけておけば、それぞれの配列特異的なプライマーを設計してくれる(念のため、他のツールでチェックしてから注文してください)。これはコピー遺伝子やファミリー遺伝子のspecific primerを設計するのに便利で、しかも、動画のようにマルチプルアライメント実行したファイルを直接コピー&ペーストしても機能する。

その他にもたくさんの機能があります。webマニュアルを確認してください。

 

 

  • windowsにインストールできる有料版では、よりたくさんの機能が利用できるようです(リンク)。
  • 無料のjava webスタート版は安定に動作しますが、フォーマットが多少おかしくても動いてしまいます。ラン中にフリーズしたら、りんごマーク => 強制終了からFastPCRを停止し、もう一度立ち上げて使ってください。

 

引用

Introduction on Using the FastPCR Software and the Related Java Web Tools for PCR and Oligonucleotide Assembly and Analysis.

Kalendar R, Tselykh TV, Khassenov B, Ramanculov EM

Methods Mol Biol. 2017;1620:33-64

 

FastPCR: An in silico tool for fast primer and probe design and advanced sequence analysis
Kalendar R, Khassenov B, Ramankulov Y, Samuilova O, Ivanov KI.

Genomics. 2017 Jul;109(3-4):312-319

 

FastPCR software for PCR, in silico PCR, and oligonucleotide assembly and analysis.

Kalendar R, Lee D, Schulman AH.

Methods Mol Biol. 2014;1116:271-302

最適なバーコード配列を選び出す BARCOSEL

 ハイスループットシーケンス解析プラットフォームの能力の使用を最大限にするため、いくつかのサンプルを一緒にプールすることが一般的である。たとえば、現時点では、Illumina HiSeqX1回の実行で1つのレーンから数億回シーケンシングでき、新しいNovaSeqは1回の実行で数十億回のシーケンシングできる。アプリケーションが1サンプルあたり数千万リードしか必要としない場合、1つのサンプルにレーン全体を割り当てることは無駄になる。したがって、いくつかのシークエンシングライブラリーは一緒にプールされ、シークエンシング装置内の同じレーンを用いて並行してシークエンシングされる。これは、シーケンシング後に異なるサンプルをどのように分離するかという問題を導入する。標準的な解決方法は、異なるサンプルに短いバーコードシーケンスラベルを使用することである。これらのバーコード配列は、ライブラリー調製中に断片に結合される。シーケンシング後にサンプルを混合し、次いでそれらを分離する2つのプロセスは、それぞれmultiplexing 及び demultiplexingと呼ばれる。

 正しく機能させるためには、バーコード配列は互いに十分に異なるべきである。バーコードシーケンスにおける冗長性は、誤り訂正の可能性を提供する。例えば、バーコード検出で1塩基のミスマッチを許容するためには、異なるバーコード配列は、少なくとも互いに3ヌクレオチドミスマッチ分離れていなければならない。より一般的には、m個のミスマッチを許容するために、すべてのバーコード対間の距離は少なくとも2m + 1でなければならない。シーケンシング技術は、バーコードが最適であるというさらなる制約を与えることがある。例えば、Illuminaシーケンサーでは、A / C用の赤色レーザーとG / T用の緑色レーザーの2つのレーザーを使用してヌクレオチドを検出する。最適な検出のために、これらの2つのヌクレオチドグループは、各バーコード位置のすべてのバーコード間でバランスがとれている必要がある。実験は、ヌクレオチド組成の多様性の低下がデータ損失をもたらすことを示している[論文より ref.1 PLOS ONE]。クラスター同定においてヌクレオチド多様性が重要であることに加えて、良好な塩基組成は上手くbasecallingを行うためにも重要である。

 新規バーコード設計、すなわち、バーコードのセットをfrom scratchで構築するプロセスは解決された問題であり、いくつかのツールが利用可能である[ref.2,3,4]。最初のバーコードデザインの1つは、バーコード間の非類似性を測定するためにハミング距離(wiki)が使用された[ref.5]。ハミング距離は[ref.6]で使用されている。挿入および欠失を考慮すると、レーベンシュタイン距離(wiki)(編集距離とも呼ばれる)が得られる [ref.7]。配列類似性、複雑さ、GC含量、およびセルフハイブリダイゼーションは、[ref.8]および[ref.9]において考慮され、バーコード間のヌクレオチドバランスはイルミナプラットフォームで少数サンプルのmultiplexingを行う際に特に重要になる。

 しかしながら、ユーザが既にあるバーコードのセットから最適なバーコードのセットを選択したい場合には、前述のツールのいずれも適用できない。サブセットの選択を報告する唯一のツール[ref.4]は、ユーザーが結果のバーコード数を定義することさえできず、さらにヌクレオチドバランスが考慮されていない(論文の追加ファイル1)。候補バーコードのセットは、すべてのバーコードサブセットが等しいと仮定している。最小ペアワイズ距離の基準は満たされているが、異なるサブセットは異なるヌクレオチドバランスを有するので、これは当てはまらない。シーケンシングセンターでは、バーコードの選択は実際的な毎日の問題である。個々の実験ごとにユニークなセットを注文することは無駄になる。他の極端な場合、ヌクレオチドバランスを保持するために将来のすべてのシーケンシング解析で同じセットのバーコードを再使用するとすれば、multiplexingされるサンプル数は常に一緒にせねばならず非常に制限的である。

 現時点では、イルミナはさまざまなサンプルのmultiplexingの際に、バーコード表を提供し、バーコードを選択する手順を説明している[ref.10 Google Scholar]。これらのrecommendationsではヌクレオチドバランスが考慮されている。ただし、イルミナの表ではバーコードセットは固定されている。この論文の著者らのツールでは、ユーザーが候補バーコードのセットを提供しユーザーが選択できる。ユーザが必要なバーコードの数を定義した後、ツールは、最小ペアワイズ配列距離の閾値を満たす最適セットを見出す。重要なことに、ヌクレオチドバランスが最適化されている。

 我々(著者ら)のツールを適用できる3つのタスクを論文図1(BMC Bioinformatics link)に示す。主なアプリケーションは、与えられた候補から最適なバーコードセットを選択することである(論文 図1a)。別のアプリケーションは、ユーザーが選択したバーコードのバーコード距離とヌクレオチドバランスを確認することである(論文 図1b)。たとえば、2つのライブラリがある場合、それらのバーコードに互換性があり、一緒に配列されるかどうかをチェックするために使用できる。第3のアプリケーションは、バーコードのセットを増やすことである(論文 図1c)。これは、ユーザが既存のシーケンシングライブラリに新しいサンプルを追加したい場合である。既存のものを考慮して、新しいサンプルの最適なバーコードが見つけられます。

 

f:id:kazumaxneo:20180711165406p:plain

BARCOSELには3つの計算モードがある。論文より転載。

 

マニュアル

http://ekhidna2.biocenter.helsinki.fi/barcosel/info.txt

BARCOSELに関するツイート

 

ラン

BARCOSELにアクセスする。

f:id:kazumaxneo:20180711180634p:plain

 

ここではトップページに準備されている8塩基バーコード配列が記載されたmulti-FASTAデータをダウンロードし、それを使って、もっともバーコードに適した組み合わせを選抜する流れを確認する。このデータは論文の説明に使われた配列セットで、TagGDを使って出力し(紹介)、それからPCRに適したものをFastPCR(紹介)で選抜して得た288のバーコード配列セットである。

1、multi-FASTAをアップロードまたはbox内に直接コピー&ペーストする。

f:id:kazumaxneo:20180711190240p:plain

2、Number of ~に最終的に必要なバーコード数を記載する。ここでは4とした。

3、Submitする。

結果: デフォルトでは10秒以内に結果は出力される。

f:id:kazumaxneo:20180711190325p:plain

元のリストから、最も塩基バランスの良い4つのバーコード配列が選抜された。出力された図をみると明らかなように、塩基配列のバランス(上の図)、illuminaシーケンサーのG/TとA/Cのレーザによるバランス、いずれも均等になっている(4の倍数の時に完全に均等になりうる)。

 

バラバラな配列であることは、選抜された配列を縦読みしてもすぐに分かる。

f:id:kazumaxneo:20180711190640p:plain

配列はテキスト形式でダウンロードできる。

 

 

次は選抜数を6にした。

f:id:kazumaxneo:20180711191414p:plain

4の倍数でないので各ポジションの塩基組成は均等にはならないが(上半分の図)、A/CとG/Tのレーザーのバランスは全てのポジションで50:50になっている(下半分の図)。

 

Advancedをクリックすると、いくつかのパラメータを指定できる。バーコード間の距離はデファルトはハミング距離 (number of nucleotide differences between two barcodes in the gapless alignment) で計算されるが、レーベンシュタイン距離 (number of substitutions, insertions, and deletions to convert one barcode sequence into another)に変えることができる。Minimum distanceは、defaultの3なら1塩基のミスマッチまでは完全に区別できることになる。2塩基ミスマッチまで100% 区別するためには、Minimum distanceの値は5にする必要がある。

詳細は論文と上にリンクを貼ったマニュアルを読んで確認してください。

 

TagGD

引用 

BARCOSEL: a tool for selecting an optimal barcode set for high-throughput sequencing

Somervuo P, Koskinen P, Mei P, Holm L, Auvinen P, Paulin L

BMC Bioinformatics. 2018 Jul 5;19(1):257. 

 

メタゲノムデータの平均ゲノムサイズや総カバレッジを推定する MicrobeCensus

 

 ショットガンメタゲノミクスは、人体や環境の微生物群集の機能的構成を特徴づけるためにますます使用されてきている[論文より ref.1-4]。これらの研究の共通の目標は、遺伝子ファミリー存在量を定量化し、環境、宿主の表現型、または実験条件の間で豊富さが異なる微生物の遺伝子またはパスウエイを同定することである。比較メタゲノミクスは、微生物が地球上の無数の環境にどのように機能的に適応しているか、そして微生物群集の機能的構成の変化がヒトの健康と疾病にどのように影響するかを明らかにしている[ref.5,6]。機能的変異は、同じショットガンデータから推定された分類学的変異の状況において解釈され、進化的および生態学的プロセスについての手がかりを与えている。これらの研究は、微生物群集DNAからの遺伝子量の正確な定量を必要とする。

 見過ごされがちだが、比較的メタゲノム解析を行う場合、微生物群における細胞の平均ゲノムサイズ(AGS)を考慮することが重要である。統計的な観点から、AGSは、コミュニティ間で豊富な遺伝子を比較する際に潜在的なバイアスの原因となる可能性があります。具体的には、コミュニティからある遺伝子をサンプリングする確率は、そのコミュニティのAGSに反比例する[ref.7 pubmed]。これを考慮しなければ、細胞あたり等コピー数で存在する遺伝子が異なるという結果を引き起こす(すなわち、偽陽性)。ライブラリーの大きさ(すなわちシーケンシングデプス)と遺伝子長に加えて、微生物群集間の生態学的に意味のあるゲノムの違いを同定するために、比較メタゲノム研究ではAGSを考慮する必要がある。 

 AGSは、環境内の微生物に作用する生態学的および進化的な力を理解するためにも重要である。生態学的な観点から、微生物のゲノムサイズは、環境の複雑さ、代謝生活様式、およびコミュニティのニッチを反映している可能性がある[ref.8-12]。例えば、腸内では、より大きいゲノムを有する生物はより一般的なライフスタイルに従うと考えられ、より小さなゲノムを有する生物はより特殊化されると考えられている[ref.8 pubmed]。例えばMetanobrevibacter smithii(1.9Mb)がH2の利用に特化したメタン菌であるのに対して、大きなゲノムを有するBacteroides thetaiotamicron(6.5Mb)はホストの食事およびホスト由来のグリカンの両方を利用できる代謝能を持っている。進化の観点から、ゲノムサイズの減少は、少数の集団における遺伝的浮動または大集団におけるgenome streamlining [ref.13,14](引用 suggesting that reduction in genome size is beneficial and driven by positive selection)を反映する可能性がある。遺伝的浮動は、突然変異の欠失方向にバイアスが固定されることによって微生物ゲノムのサイズを形成することができる[15]。対照的に、genome streamliningは、栄養素が増殖をリミットするような大集団において、栄養素を効率的な使用するような選択圧がかかることでゲノムサイズの減少を促進する[ref.14]。

 現時点までにメタゲノミクスデータ用に設計されたソフトウェアがないため、微生物群のAGSを迅速かつ正確に推定することは不可能であった。これは、ヒト微生物を含む多くの環境におけるAGS変異の程度および影響の両方に対する我々の理解を制限している。唯一公的に利用可能なソフトウェアツールGAAS [ref.10]は、微生物ゲノムのデータベースに対するショットガン配列のBLAST検索に基づいてAGSを推定している。しかしメタゲノムデータおよびリファレンスゲノムの規模が大きくかつ増大しており、BLASTの速度が比較的遅いためこの方法は計算上実用的ではない。さらに、微生物群集は、しばしば、培養されていないかまたはゲノムシーケンスされていない、高い割合の「新規な」微生物を含む。よく研究されたヒトの腸内微生物でさえ、現在の微生物リファレンスゲノムでは、種の豊富さ(species abundance )で43%と、richnessで58%は捕獲できないと推定されている[ref.16]。 GAASが新規分類群からなるメタゲノムについてAGSを正確に推定できるかどうかは不明である。 RaesらはBLASTXを用いて35の必須シングルコピー遺伝子セットにアサインされたリードの密度に基づいて平均ゲノムサイズを推定することでこの問題に対処する提案をした[ref.9]。この方法は著しく速くコミュニティの組成への依存度も低いが、300塩基対(bp)未満のリードについては設計されておらず、またソフトウェアはリリースされなかった。現在のシーケンシング技術は、より長いリードを生成し始めているが、膨大な量の既存のショートリードシーケンシングデータが残っており、例えば、MetaHIT [ref1,17-19]とHuman Microbiome Project(HMP)[ref.20]プロジェクトでは、ヒト共生微生物から平均96bpの長さで1300億以上のメタゲノミックリードを生成した[ref,21]。 Raes法が近代的なショートリードライブラリのAGSを推定できるかどうかは明らかではない。その他[ref.22,23]も同様の方法を記述しているが、これらは広範に検証されておらず、ソフトウェアとして利用可能になっていない。

 この問題に対処するために、著者らは、ショットガンメタゲノミックデータからAGSを迅速かつ正確に推定するMicrobeCensusを開発し、このツールを1,352のヒト共生微生物サンプルに適用した。 Raesらにも同様のアプローチを採用したが、有意な方法論的改善がなされ、50bpのショートリードを使用してAGSを正確に推定することができた。著者らは、AGSが、ヒト体内の位置でヒト共生微生物が著しく異なり、機能的および分類学的に大きな違いがあることを発見した。例えば、便のメタゲノムのAGSは、2.5〜5.8メガベース(Mb)の範囲であり、Bacteroidesの量と代謝、生合成、および二成分系に関連する遺伝子は正に相関したが、Firmicutes(フィルミクテス門)および膜輸送に関連する遺伝子はより小さいAGSを有するメタゲノムに豊富に存在していた。さらに著者らは、AGSが比較メタゲノム解析の主要なバイアスであり、正規化が異なる存在量の遺伝子検出を改善することを確認した。これらの発見はMicrobeCensusの開発なしにはできなかった。

  MicrobeCensusはメタゲノムのシーケンスデータから微生物群集の平均ゲノムサイズ( average genome size (AGS))を推定するためのツール。AGSは、ほぼすべての細胞性微生物(バクテリアアーキア、真菌)に存在する一組の普遍的な単一コピー遺伝子ファミリーにシーケンシングリードをアサインすることよって推定する。 これらの遺伝子はゲノムごとに1回だけ出現するため、微生物群集の平均ゲノムサイズは、これらの遺伝子にヒットしたリードの割合と反比例する。AGSが得られれば、サンプル中に存在する微生物ゲノムの総カバレッジを得ることが可能となる(genome equivalents = total bp sequenced / AGS in bp)。これは遺伝子の存在量を正規化するのに有用である。

f:id:kazumaxneo:20180711115509p:plain

MicrobeCensusのAGS推定ワークフロー。

 

 

MicrobeCensusに関するツイート。

 

簡単に言えば、どのような全細胞性生物でも厳密に1コピーだけ存在している遺伝子のセット(AGS)を定義し、その遺伝子のリードカウントを全リードカウントと比較してゲノムサイズを推定する方法論です。1コピーしか存在しない遺伝子なので、AGSのリードカウントはゲノムサイズによっては変化しません。そのため、ゲノムサイズが大きければ"whole ゲノムリードカウント/ AGSリードカウント"、は大きくなり、ゲノムサイズが小さければ"whole ゲノムリードカウント/ AGSリードカウント"は小さくなります。この性質を利用すれば、全リード数とAGSのリードカウントの比率からゲノムサイズは推定できます。この方法を複数ゲノムが混じったデータセットに適用すれば、得られるゲノムサイズ値は、複数ゲノムサイズの平均値になります。

 インストール 

mac os10.12のAnaconda3.4.2でテストした。

依存

  • Supported platforms: Mac OSX, Unix/Linux; Windows not currently supported
  • Python version 2 or 3
  • Numpy
  • BioPython 

本体 Github

#Anaconda環境ならcondaで導入できる(*1)。
conda install -c bioconda microbecensus

#
Anacondaを使ってなくてもpipで導入できる。
pip install microbecensus 

*1 #condaでインストール時にpymummerとコンフリクトするとの警告があったのでpymummerをconda uninstallで消してから導入した。

> run_microbe_census.py -h

$ run_microbe_census.py -h

usage: run_microbe_census.py [-h] [-v] [-V] [-r RAPSEARCH] [-n NREADS]

                             [-t THREADS] [-e]

                             [-l {50,60,70,80,90,100,110,120,130,140,150,175,200,225,250,300,350,400,450,500}]

                             [-q MIN_QUALITY] [-m MEAN_QUALITY] [-d]

                             [-u MAX_UNKNOWN]

                             SEQFILES OUTFILE

 

Estimate average genome size from metagenomic data.

 

positional arguments:

  SEQFILES              path to input metagenome(s); for paired-end

                        metagenomes use commas to specify each file (ex:

                        read_1.fq.gz,read_2.fq.gz); can be FASTQ/FASTA; can be

                        gzip (.gz) or bzip (.bz2) compressed

  OUTFILE               path to output file containing results

 

optional arguments:

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

  -v                    print program's progress to stdout (default = False)

  -V, --version         show program's version number and exit

  -r RAPSEARCH          path to external RAPsearch2 v2.15 binary; useful if

                        precompiled RAPsearch2 v2.15 binary included with

                        MicrobeCensus does not work on your system

 

Pipeline throughput (optional):

  -n NREADS             number of reads to sample from SEQFILES and use for

                        average genome size estimation. to use all reads set

                        to 100000000. (default = 2000000)

  -t THREADS            number of threads to use for database search (default

                        = 1)

  -e                    quit after average genome size is obtained and do not

                        estimate the number of genome equivalents in SEQFILES.

                        useful in combination with -n for quick tests (default

                        = False)

 

Quality control (optional):

  -l {50,60,70,80,90,100,110,120,130,140,150,175,200,225,250,300,350,400,450,500}

                        all reads trimmed to this length; reads shorter than

                        this discarded (default = median read length)

  -q MIN_QUALITY        minimum base-level PHRED quality score (default = -5;

                        no filtering)

  -m MEAN_QUALITY       minimum read-level PHRED quality score (default = -5;

                        no filtering)

  -d                    filter duplicate reads (default = False)

  -u MAX_UNKNOWN        max percent of unknown bases per read (default = 100

                        percent; no filtering)

——

 

テストラン

git clone https://github.com/snayfach/MicrobeCensus.git
cd /path/to/MicrobeCensus/test
python test_microbe_census.py

an 3 tests in 44.751s

 

OK

——

 

ラン

シーケンスデータはFASTQまたはFASTAフォーマット、gzip (.gz) か bzip (.bz2) 圧縮していても受け付ける。ペアエンドはコンマで区切って記載する。

run_microbe_census.py -t 12 -n 2000000 pair1.fq,pair2.fq output
  • -t     number of threads to use for database search (default= 1)
  • -n    number of reads to sample from SEQFILES and use for average genome size estimation. to use all reads set to 100000000. (default = 2000000)
  • -e    quit after average genome size is obtained and do not estimate the number of genome equivalents in SEQFILES. useful in combination with -n for quick tests (default = False)

コンタミネーションがかなりあったバクテリアのシーケンスデータの出力。

reads_sampled:  200000

trimmed_length: 300

min_quality:    -5

mean_quality:   -5

filter_dups:    False

max_unknown:    100

 

Results

average_genome_size:    7452243.21927

total_bases:    246786445

genome_equivalents:     33.1157260624

output (END)

 

importして他のメタゲノム解析パイプラインのモジュールとして使うこともできます。詳細はGithubで確認してください。

引用

Average genome size estimation improves comparative metagenomics and sheds light on the functional ecology of the human microbiome

Nayfach S, Pollard KS.

Genome Biol. 2015 Mar 25;16:51.

 

ONTのロングリードのマッピングからSVを検出する NanoSV

  第二世代DNAシークエンシングは、ヒトの遺伝病の研究と診断に不可欠な技術となっている。ヒトの全エキソームシーケンシングにより、メンデル疾患(ref.1)に関与する遺伝子が同定され、ヒトの全ゲノムシーケンシングにより、遺伝子内だけでなくノンコード領域の遺伝的変化によっても無数の疾患が引き起こされることが明らかになった(ref.2)。ゲノムシーケンシングにより患者のユニークな突然変異プロファイルの完全な解明ができれば個人のオーダーメイド治療が可能になるため、臨床決定の場面で急速に採用されてきている。
 SVは一塩基変異(SNV)よりもはるかに多くの塩基変化を起こす変異の重要なクラスであり、ヒトゲノムの構造変異(SV)を検出する堅牢な方法が不可欠である。さらに、SVは広範な遺伝的障害に関与している(ref.6)。
 ゲノムシーケンシングにおける特に革新的な開発は、DNA配列を直接かつリアルタイムで測定するためのナノポアタンパク質の使用である。消費者向けデバイスで成功した最初の実装は、2014年にOxford Nanopore TechnologiesのMinionによって達成された(ref,8)。 Minionは長さが数百キロベースのDNAのシーケンシングができ、すでにいくつかの生物のゲノムシーケンスにつながった(ref.9,10)。
 ナノポアシーケンシングによって生成されるロングリードの重要かつ自然な応用は、構造変化の同定である。ロングリードシークエンシングは、これまでにないスケールと深さのSVの発見の根拠となっている。最初の成功は、Pacific BioSciences SMRTロングリードシークエンシングプラットフォーム(ref.12,13)およびBioNano (ref.14)および10X Genomics (ref.15)などの代替方法を使用して達成された。 構造変化を正確に同定するために、ショートリードの次世代シーケンシングデータでは(しばしば)複数の間接的な情報源に依存しているが、ロングリードデータは構造変化を直接明らかにできる。
 この研究では、Minionシーケンサーにより二倍体ヒトゲノムの全ゲノムシーケンシングを初めて実証する。著者らは、11-16Xのデプスで2人のクロモリプリシスから生じる先天性疾患患者のゲノムを解析した。我々(著者ら)は、ナノポアリードを使い新規の複雑なSVブレークポイントを高感度に検出する可能性を実証するために、新しい計算パイプラインを採用している。 Minionのロングリードは、ゲノムワイドな遺伝的変異(SNVsおよびSVs)の効率的なフェージングを可能にし、患者におけるクロモリピスの長距離の構造を明らかにした。さらに、同じ患者ゲノムのイルミナのショートリードシーケンシングデータにおいて検出されないSVのかなりの割合を同定した。
 これらの結果は、安価な装置を使って臨床のヒト試料をリアルタイムシーケンシングできる実用性の高さを強調している。ナノポアベースのシーケンシングはほとんど投資を必要とせず、現在のデバイスのフットプリントが非常に小さいため、これらのシーケンサを主流に採用することは、現在の中央センターでのシーケンシング解析パラダイムを根本的に変える可能性がある。

 ここで開発されたNanoSVアルゴリズムhttps://github.com/mroosmalen/nanosv)は、入力としてBAMファイルを使用する(一部略)。NanoSVは、スプリットリードのクラスタリングを使用してSVブレークポイントジャンクションを識別する。第1のステップでは、スプリットされたリードのすべてのマッピングセグメントオリジナルのリード位置に基づいて順序付けられる。整列されたリードは、その整列したセグメント間のギャップ、すなわち挿入(Supplementary Figure 8 , Supplementary Figure 9 )または単に低品質配列により、リファレンスゲノム上のどこにもアライメントしない配列を含みうる(以下略)。

 


 NanoSVに関するツイート。

 インストール 

mac os10.12でテストした。lastはpython2.7環境で実行したが、nanosvはAnaconda3.5.2環境で実行した。

依存

  • samtoolsかsambambaのどちらか
  • ロングリードマッパー(last、mininmap2、graphmap、ngmlrなどいずれか)

本体 Github

#Anaconda環境ならcondaで導入できる(リンク)。
conda install -c bioconda nanosv


#
Anacondaを使ってなくてもpipで導入できる。
pip install nanosv

nanosv -h

$ nanosv -h

usage: nanosv [-h] [-t THREADS] [-s SAMBAMBA] [-c CONFIG] [-b BED] [-o OUTPUT]

              [-f SNP_FILE]

              bam

 

Put here a description.

 

positional arguments:

  bam                   /path/to/file.bam

 

optional arguments:

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

  -t THREADS, --threads THREADS

                        Number of threads [default: 4]

  -s SAMBAMBA, --sambamba SAMBAMBA

                        Give the full path to the sambamba or samtools

                        executable [default: sambamba ]

  -c CONFIG, --config CONFIG

                        Give the full path to your own ini file [ default:

                        config.ini ]

  -b BED, --bed BED     Give the full path to your own bed file, used for

                        coverage depth calculations [default: human_hg19.bed ]

  -o OUTPUT, --output OUTPUT

                        Give the full path to the output vcf file [default:

                        <stdout> ]

  -f SNP_FILE, --snp_file SNP_FILE

                        Give full path to the SNP variant file for phasing.

                        Supporting file formats: BED and VCF

——

 

ラン

1、 はじめにマッピングしてsamを作る。2018年7月現在、いくつかのアライナが利用できる。

#1 bwaを使う(紹介
bwa index ref.fa
bwa mem -x ont2d -M -t 12 ref.fa input.fq > output.sam

#2 minimap2を使う(
紹介
minimap2 -t 12 -a
ref.fa input.fq > output.sam

#3 NGMLR
を使う 紹介
ngmlr -x ont -t 12 -r
ref.fa input.fq > output.sam

$4graphmapを使う。
graphmap align -t 12 -r ref.fa -d
input.fq -o output.sam

 オーサーらは、コンピュータへの負担が大きくなるが、 LASTでマッピングした時が最もnanoSVの検出精度が高かったと述べている。試しましたが、上のマッパーと比べ非常に時間がかかります。データサイズが大きい時は注意してください。

#index
lastdb refdb ref.fa

#dictファイル作成
java -jar picard.jar CreateSequenceDictionary REFERENCE=ref.fa OUTPUT=ref.dict


#マッピングしてデータに適したマッピングパラメータ作成(トレーニング)
last-train -P12 -Q1
refdb input.fq > my-params

#マッピング。シーケンスデータがfastaなので-Qは0。
lastal -P12 -Q0 -p my-params
refdb input.fq | last-split - > output.maf

#mafからsamに変換。
maf-convert -f ref.dict sam -r ‘ID:[id] PL:[nanopore] SM:[sample]’
output.maf > output.sam
  • -P: number of parallel threads (1)
  • -Q    input format: 0=fasta, 1=fastq-sanger, 2=fastq-solexa, 3=fastq-illumina,
    4=prb, 5=PSSM (0)
  • -r     read group info for sam format

LASTを使ったロングリードのマッピングは以前の投稿も参考にして下さい(リンク)。

 

2、sambamba(紹介)を使ってsamからbamに変換。

sambamba view -t 8 -h -S -f=bam output.sam > output.bam
sambamba sort -t 8 output.bam -o output_sorted.bam
  • -f specify which format to use for output (default is SAM)
  • -S specify that input is in SAM format

 

3、NanoSVを使い、bamからSVを検出する。sambambaかsamtoolsのパスをフルパスで与える。詳細なパラメータはconfigファイルで指定する。configファイルはリンク先のように記述する。"-c"で指定しない場合、defaultのconfig条件で実行される。

NanoSV -t 12 -s /usr/local/bin/sambamba -c config.ini -b gene.bed -o output.vcf output_sorted.bam
  • -s     Give the full path to the sambamba or samtools executable [default: sambamba ]
  •  -c    Give the full path to your own ini file [ default: config.ini ]
  • -t      Maximum number of threads to use [default: 4 ]
  • -b     Give the full path to your own bed file, used for coverage depth calculations [default: human_hg19.bed ]

    -o     Give the full path to the output vcf file [default: <stdout> ]

結果はvcf形式で出力される。

 

ハプロタイプフェージングについては、現在チューニング最中と記載されている。

 

引用

Mapping and phasing of structural variation in patient genomes using nanopore sequencing

Mircea Cretu Stancu, Markus J. van Roosmalen, Ivo Renkens, Marleen M. Nieboer, Sjors Middelkamp, Joep de Ligt, Giulia Pregno, Daniela Giachino, Giorgia Mandrile, Jose Espejo Valle-Inclan, Jerome Korzelius, Ewart de Bruijn, Edwin Cuppen, Michael E. Talkowski, Tobias Marschall, Jeroen de Ridder & Wigard P. Kloosterman

Nature Communicationsvolume 8, Article number: 1326 (2017)