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が既存のツールを大幅に上回る。
図。3000万リードの定量にかかる時間。論文より転載。
公式HP
マニュアル
http://pachterlab.github.io/kallisto/manual.html
Kallistoをベースに作られたツール
https://pachterlab.github.io/kallisto/apps
インストール
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ディレクトリができる。出力は以下のように説明されている。
- 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
- 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.
- 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が全く発生していない)。
--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/