指北 | Kallisto 使用手册 - 普通笔记本上分析转录组

写在前面

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 的功能对我来说，似乎没有意义。不折腾了