寫(xiě)在前面
Kallisto 是一款 “Alignmnet-Free” 的轉(zhuǎn)錄組表達(dá)量估計(jì)的工具。其提出了一種 “pseudoalignment ” 的概念，實(shí)則基于 Kmer 索引，然后 Hash 計(jì)數(shù)和分析。經(jīng)過(guò)幾年開(kāi)發(fā)，其實(shí)能做的事情已經(jīng)遠(yuǎn)超過(guò)前述我推文中介紹的內(nèi)容，比如早已支持鏈特異文庫(kù)，也支持單細(xì)胞測(cè)序數(shù)據(jù)。
就我個(gè)人接下來(lái)不長(zhǎng)時(shí)間的使用場(chǎng)景，應(yīng)是不會(huì)涉及單細(xì)胞，故不做單細(xì)胞相關(guān)記錄，僅看看 Kallisto 在普通轉(zhuǎn)錄組分析中的應(yīng)用和效果。
關(guān)注《生信札記》的朋友應(yīng)該清楚，我多次表達(dá)了過(guò)一個(gè)觀點(diǎn)（當(dāng)然更多是在小群里或者是語(yǔ)音上），即新手才最喜歡翻譯軟件手冊(cè)。我自認(rèn)為這個(gè)觀點(diǎn)很 solid，所以對(duì)自己的師弟師妹，往往建議看文獻(xiàn)，自己看官網(wǎng)手冊(cè)，而不要去翻網(wǎng)絡(luò)上其他人翻譯的手冊(cè)，因?yàn)榇蟛糠制鋵?shí)存在問(wèn)題，比如我現(xiàn)在在寫(xiě)的這份。
事實(shí)上，并非別人翻的，別人寫(xiě)的就不好，而只是說(shuō)，這些往往是別人消化過(guò)得，對(duì)與錯(cuò)并無(wú)法保證。所以我個(gè)人建議的最好方式一直是：
- 看軟件文獻(xiàn)，多少了解下軟件實(shí)現(xiàn)邏輯與大體算法（生物背景出身，新手一般要開(kāi)很多遍，掌握不了就大體了解，這有助于后續(xù)明白參數(shù)設(shè)置）
- 看軟件官方手冊(cè)，事實(shí)上，這往往是最好的參考，畢竟開(kāi)發(fā)者自己寫(xiě)的，幾乎不存在可能的原則性錯(cuò)誤
- 如果實(shí)在看不懂，或者覺(jué)得有一些理解不清晰，那么才去看別人寫(xiě)的手冊(cè)（比如你現(xiàn)在在看的這一篇）
好，廢話說(shuō)太多，還是動(dòng)手干活。
下載 kallisto
除非你使用的平臺(tái)不是 Windows MacOS 或者 Linux 操作系統(tǒng)，否則，直接下載編譯好的二進(jìn)制文件最為合適，具體鏈接 https://github.com/pachterlab/kallisto/releases 獲得鏈接后，下載即可
wget https://github.com/pachterlab/kallisto/releases/download/v0.46.2/kallisto_linux-v0.46.2.tar.gz
# 解壓即可獲得 可執(zhí)行程序
tar -zxvf kallisto_linux-v0.46.2.tar.gz
雖然這里寫(xiě)的是在 Linux 上，其實(shí)我是在 Windows 的 CMD 下執(zhí)行的，也不是 WSL...
查看可用命令
# 直接運(yùn)行即可
kallisto
kallisto 0.46.1
Usage: kallisto <CMD> [arguments] ..
Where <CMD> can be one of:
index 建立索引 - Builds a kallisto index
quant 轉(zhuǎn)錄本定量 - Runs the quantification algorithm
bus 單細(xì)胞測(cè)序，不看 - Generate BUS files for single-cell data
pseudo 單細(xì)胞測(cè)序，不看 - Runs the pseudoalignment step 生成BAM？
merge 單細(xì)胞測(cè)序，不看 - 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>
于是，只是做個(gè)轉(zhuǎn)錄本定量，值得我們關(guān)注的只有兩個(gè)子程序 index 和 quant，即 建立索引 和 表達(dá)定量
建立索引
對(duì)于每個(gè)子程序，一般我們直接
kallisto index
即可看到所有的子程序參數(shù)說(shuō)明
kallisto 0.46.2
Builds a kallisto index
Usage: kallisto index [arguments] FASTA-files
Required argument:
-i, --index=STRING 輸出的索引文件名（輸出基因一個(gè)文件，可包含路徑） Filename for the kallisto index to be constructed
Optional argument:
-k, --kmer-size=INT Kmer的長(zhǎng)度，越長(zhǎng)，多匹配就越少，但也容易錯(cuò)漏一些Kmer匹配，一般這類默認(rèn)參數(shù)，除非必要，也懶得去調(diào)整 k-mer (odd) length (default: 31, max value: 31)
--make-unique 對(duì)于序列名字相同的Fasta記錄，自動(dòng)標(biāo)識(shí)，應(yīng)是大多數(shù)時(shí)候用不上 Replace repeated target names with unique names
Emmm, 發(fā)現(xiàn)要寫(xiě)完這個(gè)，得下載一些數(shù)據(jù)....我先去寫(xiě)個(gè) TBtools 下載數(shù)據(jù)的插件再說(shuō)吧
搞了大半天，寫(xiě)完數(shù)據(jù)下載插件了，具體見(jiàn) 插件 | 人人-點(diǎn)點(diǎn)點(diǎn)-光速下載 NCBI/ENA NGS原始數(shù)據(jù)
暫時(shí)沒(méi)時(shí)間寫(xiě)完這個(gè)，回頭再說(shuō)吧。
于是，對(duì)于 Kallisto 的索引構(gòu)建常用命令為
kallisto index -k 31 -i PathToIndexOutFile InTranscriptome.fa
轉(zhuǎn)錄本定量
kallisto quant
參數(shù)略多，可以看看
kallisto 0.46.2
Computes equivalence classes for reads and quantifies abundances
Usage: kallisto quant [arguments] FASTQ-files
Required arguments:
-i, --index=STRING 輸入的索引文件，F(xiàn)ilename for the kallisto index to be used for
quantification
-o, --output-dir=STRING 輸出文件目錄，注意到輸出文件較多，所以指定的是輸出目錄，Directory to write output to
Optional arguments:
--bias 開(kāi)啟序列偏好矯正，按理說(shuō)給上這個(gè)參數(shù)合適，Perform sequence based bias correction
-b, --bootstrap-samples=INT BootStrap重抽樣測(cè)試結(jié)果穩(wěn)定性，如果只是要Counts字，不需要開(kāi)；除非用的配套下游分析軟件需要這個(gè)穩(wěn)定性參數(shù)，Number of bootstrap samples (default: 0)
--seed=INT 隨機(jī)數(shù)產(chǎn)生種子，不管，Seed for the bootstrap sampling (default: 42)
--plaintext 輸出文本文件，而不輸出HDF5，給上，用不到 Output plaintext instead of HDF5
--fusion 檢測(cè)融合基因，用不上 Search for fusions for Pizzly
--single 輸入的數(shù)據(jù)是單端數(shù)據(jù) Quantify single-end reads
--single-overhang 對(duì)于雙端數(shù)據(jù)，納入僅一端被覆蓋的誒情況，按理說(shuō)影響不大，畢竟并不會(huì)用上所有 Kmer，不用感覺(jué)也很好，畢竟邊緣的不多 Include reads where unobserved rest of fragment is
predicted to lie outside a transcript
--fr-stranded 輸入數(shù)據(jù)為鏈特異的一種，目前較少使用 Strand specific reads, first read forward
--rf-stranded 輸入數(shù)據(jù)為鏈特異的一種，目前最常用的dUTP，NSR等建庫(kù)，使用這一參數(shù) Strand specific reads, first read reverse
-l, --fragment-length=DOUBLE 單端數(shù)據(jù)指定建庫(kù)時(shí)插入片段長(zhǎng)度，如果不知道，常見(jiàn) RNAseq 約 200，具體看實(shí)驗(yàn) Estimated average fragment length
-s, --sd=DOUBLE 單端數(shù)據(jù)指定建庫(kù)時(shí)插入片段長(zhǎng)度標(biāo)準(zhǔn)差，如果不知道，常見(jiàn) RNAseq 約 30 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 設(shè)定運(yùn)行時(shí)使用的線程數(shù) Number of threads to use (default: 1)
--pseudobam 輸出比對(duì)到轉(zhuǎn)錄組的BAM文件 Save pseudoalignments to transcriptome to BAM file
--genomebam 【目前，加上這個(gè)參數(shù)會(huì)產(chǎn)生段錯(cuò)誤，且似乎四五年過(guò)去了，尚未解決，退回上一版本即可 https://github.com/pachterlab/kallisto/issues/254】輸出比對(duì)到基因組的BAM文件（需要其他輸入） Project pseudoalignments to genome sorted BAM file
-g, --gtf 設(shè)定基因結(jié)構(gòu)注釋信息，GTF格式（--genomebam 時(shí)必須帶上） GTF file for transcriptome information
(required for --genomebam)
-c, --chromosomes 設(shè)定染色體長(zhǎng)度信息，制表符分隔（染色體ID 染色體長(zhǎng)度）【這與TBtools一些功能類似，主要是避免要求用戶保存染色體序列，如果不給這個(gè)參數(shù)，相信是直接從 gtf文件中獲取染色體長(zhǎng)度，不過(guò)這往往影響不大，只有細(xì)微偏差；但建議還是給上】Tab separated file with chromosome names and lengths
(optional for --genomebam, but recommended)
--verbose Print out progress information every 1M proccessed reads
于是，文檔翻完了，可以確定大體命令。
最基礎(chǔ)的命令
雙端數(shù)據(jù)
kallisto quant -i PathToKallistoIndex -o PathToOutDirectory R1.fq R2.fq
單端數(shù)據(jù) （如果確實(shí)不知道建庫(kù)時(shí)用的插入片段大?。?/p>
kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --single -l 200 -s 30 R1.fq
保守并推薦的命令
非鏈特異測(cè)序數(shù)據(jù)（以雙端為例，單端數(shù)據(jù)，自行調(diào)整）
kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --bias -t NumberOfThreads R1.fq R2.fq
鏈特異測(cè)序數(shù)據(jù)（以目前最常見(jiàn)的 dUTP 建庫(kù)方式為準(zhǔn)）
kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --bias -t NumberOfThreads --rf-stranded R1.fq R2.fq
最全面的命令
非鏈特異測(cè)序數(shù)據(jù)（鏈特異的自行調(diào)整）
kallisto quant \
-i PathToKallistoIndex \ # 索引文件
-o PathToOutDirectory \ # 輸出目錄
--bias \ # 測(cè)序偏好矯正
-t 1 \ # 指定線程數(shù)，此處用 1
-b 0 \ # 認(rèn)為幾乎只有做方法學(xué)的才會(huì)用到這個(gè) bootstrap
--single-overhang \ # 納入雙端數(shù)據(jù)中僅有一端匹配的情況
--pseudobam \ # 輸出轉(zhuǎn)錄組的 BAM 文件
--genomebam \ # 輸出基因組比對(duì)的BAM文件
-g GeneStructureAnno.gtf \ # 基因結(jié)構(gòu)注釋信息，用于 --genomebam
-c ChrLen.tab.txt \ # 染色體長(zhǎng)度信息文件
R1.fq R2.fq
寫(xiě)在最后
Emmm，搞個(gè)使用手冊(cè)，想想確實(shí)沒(méi)啥用。后面就把他打包了，弄成 TBtools 插件吧。
中間寫(xiě)了 Aspera 打包工具，感覺(jué)還不錯(cuò)。
似乎應(yīng)該補(bǔ)充一個(gè)鏈特異數(shù)據(jù)判斷功能....再說(shuō)吧
補(bǔ)充
Emmm，我嘗試了下生成基因組比對(duì)的 BAM，不過(guò)似乎沒(méi)啥用
為此我還專門(mén)優(yōu)化了 TBtools 的一個(gè)功能，但還是沒(méi)生成符合要求的 GTF 格式，想想找不到問(wèn)題，還是算了。
java -cp /tools/TBtools_JRE1.6.jar biocjava.bioDoer.GXFUtils.RecallmRNAFeature --in Lchinesis_genome.Chr.gtf --out Lchinesis_genome.Chr.mRNA.gff3
sed 1d Lchinesis_genome.Chr.mRNA.gff3|sort -k1,1 -k4,4n -k5,5nr|perl -lane 'print join qq{\t},@F[0..7],(join qq{ },map {s/ID/transcript_id/r} map {s/=(\S+)/ "$1";/r} map {split /;/} $F[8])' |perl -F'\t' -lape 'if($F[2] eq "mRNA"){$F[2]="gene"}print join qq{\t},@F'|perl -F'\t' -lane '$F[2] eq "mRNA" and $F[2] = "transcript";if($F[2] eq "exon"){$F[8]=~s/Parent/transcript_id/}print join qq{\t},@F'|perl -lane 'print unless $F[2] eq "CDS"' > TheBest.gtf
PS: 轉(zhuǎn)念一想，之前的想法有點(diǎn)問(wèn)題，這個(gè) genomeBam 的功能對(duì)我來(lái)說(shuō)，似乎沒(méi)有意義。不折騰了
