写在前面Kallisto 是一款 “Alignmnet-Free” 的转录组表达量估计的工具。其提出了一种 “pseudoalignment ” 的概念，实则基于 Kmer 索引，然后 Hash 计数和分析。经过几年开发，其实能做的事情已经远超过前述我推文中介绍的内容，比如早已支持链特异文库，也支持单细胞测序数据。
就我个人接下来不长时间的使用场景，应是不会涉及单细胞，故不做单细胞相关记录，仅看看 Kallisto 在普通转录组分析中的应用和效果。
关注《生信札记》的朋友应该清楚，我多次表达了过一个观点（当然更多是在小群里或者是语音上），即新手才最喜欢翻译软件手册。我自认为这个观点很 solid，所以对自己的师弟师妹，往往建议看文献，自己看官网手册，而不要去翻网络上其他人翻译的手册，因为大部分其实存在问题，比如我现在在写的这份。
事实上，并非别人翻的，别人写的就不好，而只是说，这些往往是别人消化过得，对与错并无法保证。所以我个人建议的最好方式一直是： 看软件文献，多少了解下软件实现逻辑与大体算法（生物背景出身，新手一般要开很多遍，掌握不了就大体了解，这有助于后续明白参数设置） 看软件官方手册，事实上，这往往是最好的参考，毕竟开发者自己写的，几乎不存在可能的原则性错误 如果实在看不懂，或者觉得有一些理解不清晰，那么才去看别人写的手册（比如你现在在看的这一篇）
好，废话说太多，还是动手干活。 下载 kallisto除非你使用的平台不是 Windows MacOS 或者 Linux 操作系统，否则，直接下载编译好的二进制文件最为合适，具体链接 https://github.com/pachterlab/kallisto/releases 获得链接后，下载即可 wget https://github.com/pachterlab/kallisto/releases/download/v0.46.2/kallisto_linux-v0.46.2.tar.gz
# 解压即可获得 可执行程序
tar -zxvf kallisto_linux-v0.46.2.tar.gz
虽然这里写的是在 Linux 上，其实我是在 Windows 的 CMD 下执行的，也不是 WSL... 查看可用命令# 直接运行即可
kallisto
kallisto 0.46.1
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 生成BAM？
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>
于是，只是做个转录本定量，值得我们关注的只有两个子程序 index 和 quant，即 建立索引 和 表达定量 建立索引对于每个子程序，一般我们直接 kallisto index
即可看到所有的子程序参数说明 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 Kmer的长度，越长，多匹配就越少，但也容易错漏一些Kmer匹配，一般这类默认参数，除非必要，也懒得去调整 k-mer (odd) length (default: 31, max value: 31)
--make-unique 对于序列名字相同的Fasta记录，自动标识，应是大多数时候用不上 Replace repeated target names with unique names
Emmm, 发现要写完这个，得下载一些数据....我先去写个 TBtools 下载数据的插件再说吧
搞了大半天，写完数据下载插件了，具体见 插件 | 人人-点点点-光速下载 NCBI/ENA NGS原始数据
暂时没时间写完这个，回头再说吧。
于是，对于 Kallisto 的索引构建常用命令为 kallisto index -k 31 -i PathToIndexOutFile InTranscriptome.fa
转录本定量kallisto quant
参数略多，可以看看 kallisto 0.46.2
Computes equivalence classes for reads and quantifies abundances
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 BootStrap重抽样测试结果稳定性，如果只是要Counts字，不需要开；除非用的配套下游分析软件需要这个稳定性参数，Number of bootstrap samples (default: 0)
--seed=INT 随机数产生种子，不管，Seed for the bootstrap sampling (default: 42)
--plaintext 输出文本文件，而不输出HDF5，给上，用不到 Output plaintext instead of HDF5
--fusion 检测融合基因，用不上 Search for fusions for Pizzly
--single 输入的数据是单端数据 Quantify single-end reads
--single-overhang 对于双端数据，纳入仅一端被覆盖的诶情况，按理说影响不大，毕竟并不会用上所有 Kmer，不用感觉也很好，毕竟边缘的不多 Include reads where unobserved rest of fragment is
predicted to lie outside a transcript
--fr-stranded 输入数据为链特异的一种，目前较少使用 Strand specific reads, first read forward
--rf-stranded 输入数据为链特异的一种，目前最常用的dUTP，NSR等建库，使用这一参数 Strand specific reads, first read reverse
-l, --fragment-length=DOUBLE 单端数据指定建库时插入片段长度，如果不知道，常见 RNAseq 约 200，具体看实验 Estimated average fragment length
-s, --sd=DOUBLE 单端数据指定建库时插入片段长度标准差，如果不知道，常见 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 设定运行时使用的线程数 Number of threads to use (default: 1)
--pseudobam 输出比对到转录组的BAM文件 Save pseudoalignments to transcriptome to BAM file
--genomebam 【目前，加上这个参数会产生段错误，且似乎四五年过去了，尚未解决，退回上一版本即可 https://github.com/pachterlab/kallisto/issues/254】输出比对到基因组的BAM文件（需要其他输入） Project pseudoalignments to genome sorted BAM file
-g, --gtf 设定基因结构注释信息，GTF格式（--genomebam 时必须带上） GTF file for transcriptome information
(required for --genomebam)
-c, --chromosomes 设定染色体长度信息，制表符分隔（染色体ID 染色体长度）【这与TBtools一些功能类似，主要是避免要求用户保存染色体序列，如果不给这个参数，相信是直接从 gtf文件中获取染色体长度，不过这往往影响不大，只有细微偏差；但建议还是给上】Tab separated file with chromosome names and lengths
(optional for --genomebam, but recommended)
--verbose Print out progress information every 1M proccessed reads
于是，文档翻完了，可以确定大体命令。 最基础的命令双端数据 kallisto quant -i PathToKallistoIndex -o PathToOutDirectory R1.fq R2.fq
单端数据 （如果确实不知道建库时用的插入片段大小） kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --single -l 200 -s 30 R1.fq
保守并推荐的命令非链特异测序数据（以双端为例，单端数据，自行调整） kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --bias -t NumberOfThreads R1.fq R2.fq
链特异测序数据（以目前最常见的 dUTP 建库方式为准） kallisto quant -i PathToKallistoIndex -o PathToOutDirectory --bias -t NumberOfThreads --rf-stranded R1.fq R2.fq
最全面的命令非链特异测序数据（链特异的自行调整） kallisto quant \
-i PathToKallistoIndex \ # 索引文件
-o PathToOutDirectory \ # 输出目录
--bias \ # 测序偏好矫正
-t 1 \ # 指定线程数，此处用 1
-b 0 \ # 认为几乎只有做方法学的才会用到这个 bootstrap
--single-overhang \ # 纳入双端数据中仅有一端匹配的情况
--pseudobam \ # 输出转录组的 BAM 文件
--genomebam \ # 输出基因组比对的BAM文件
-g GeneStructureAnno.gtf \ # 基因结构注释信息，用于 --genomebam
-c ChrLen.tab.txt \ # 染色体长度信息文件
R1.fq R2.fq
写在最后Emmm，搞个使用手册，想想确实没啥用。后面就把他打包了，弄成 TBtools 插件吧。
中间写了 Aspera 打包工具，感觉还不错。
似乎应该补充一个链特异数据判断功能....再说吧 补充Emmm，我尝试了下生成基因组比对的 BAM，不过似乎没啥用
为此我还专门优化了 TBtools 的一个功能，但还是没生成符合要求的 GTF 格式，想想找不到问题，还是算了。 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: 转念一想，之前的想法有点问题，这个 genomeBam 的功能对我来说，似乎没有意义。不折腾了
|
