macでインフォマティクス

macでインフォマティクス

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

Quasi-Mappingによって高速にRNA seqの定量を行う Kallisto

2020  4/30 help更新、condaによるインストール追記

2021 2/3 タイトル修正

 

 典型的な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で導入できる。

#conda (bioconda link)
conda create -n kallisto -y
conda activate kallisto
conda install -c bioconda -y kallisto

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

 > kallisto -h

$ kallisto -h

Error: invalid command -h

kallisto 0.46.2

 

Usage: kallisto <CMD> [arguments] ..

 

Where <CMD> can be one of:

 

    index         Builds a kallisto index 

    quant         Runs the quantification algorithm 

    bus           Generate BUS files for single-cell data 

    pseudo        Runs the pseudoalignment step 

    merge         Merges several batch runs 

    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>

kallisto index -h

$ kallisto index -h

 

kallisto: invalid option -- h

Error: no FASTA files specified

Error: need to specify kallisto index name

kallisto 0.46.2

Builds a kallisto index

 

Usage: kallisto index [arguments] FASTA-files

 

Required argument:

-i, --index=STRING          Filename for the kallisto index to be constructed 

 

Optional argument:

-k, --kmer-size=INT         k-mer (odd) length (default: 31, max value: 31)

    --make-unique           Replace repeated target names with unique names

 

> kallisto quant -h

$ kallisto quant -h

kallisto: invalid option -- h

 

Error: kallisto index file missing

Error: Missing read files

Error: need to specify output directory 

 

Usage: kallisto quant [arguments] FASTQ-files

 

Required arguments:

-i, --index=STRING            Filename for the kallisto index to be used for

                              quantification

-o, --output-dir=STRING       Directory to write output to

 

Optional arguments:

    --bias                    Perform sequence based bias correction

-b, --bootstrap-samples=INT   Number of bootstrap samples (default: 0)

    --seed=INT                Seed for the bootstrap sampling (default: 42)

    --plaintext               Output plaintext instead of HDF5

    --fusion                  Search for fusions for Pizzly

    --single                  Quantify single-end reads

    --single-overhang         Include reads where unobserved rest of fragment is

                              predicted to lie outside a transcript

    --fr-stranded             Strand specific reads, first read forward

    --rf-stranded             Strand specific reads, first read reverse

-l, --fragment-length=DOUBLE  Estimated average fragment length

-s, --sd=DOUBLE               Estimated standard deviation of fragment length

                              (default: -l, -s values are estimated from paired

                               end data, but are required when using --single)

-t, --threads=INT             Number of threads to use (default: 1)

    --pseudobam               Save pseudoalignments to transcriptome to BAM file

    --genomebam               Project pseudoalignments to genome sorted BAM file

-g, --gtf                     GTF file for transcriptome information

                              (required for --genomebam)

-c, --chromosomes             Tab separated file with chromosome names and lengths

                              (optional for --genomebam, but recommended)

    --verbose                 Print out progress information every 1M proccessed reads

 

 

ラン

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

 

参考ページ 1

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

 

参考ページ2

https://pmelsted.wordpress.com/2018/01/29/bam/