Kallisto 软件的上、下游分析流程

昨天，我在答疑一个学员的问题时，写了一篇稿子：在什么情况下基因ID转换会100%失败？，曾老板看了后说这个数据提供的是Kallisto定量，需要我前去学习一下它的上下游。数据为 GSE270697：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE270697。

简单的搜了一下相关资料，发现这个软件不需要经过比对就可以直接进行转录本层面的定量，于2016年5月份发表在Nat Biotechnol（IF=33.1），标题为《Near-optimal RNA-Seq quantification》。

关于软件的用法，我们生信技能树的咖啡大师张永杰在他的菜鸟团专辑：《生物信息学百款软件》中已经介绍过，见 Kallisto — 基于伪比对的转录本定量。

首先是安装

依然是使用conda软件，简单方便的安装：

# 激活小环境
conda activate rna

# 安装
conda install kallisto -y

构建索引

物种为小鼠，去 ensemble 数据库下载最新参考基因组：https://www.ensembl.org/index.html。

注意下载的fa为转录本序列，转录组序列比较小，很快就会下载完成：

# 进入参考目录
cd /nas1/zhangj/database/genome/Mus_musculus/release113

# 下载
wget -c https://ftp.ensembl.org/pub/release-113/fasta/mus_musculus/cdna/Mus_musculus.GRCm39.cdna.all.fa.gz 

# 解压
gunzip Mus_musculus.GRCm39.cdna.all.fa.gz

建索引

• -i 指定构建索引的输出文件名
• -t 8 设置线程
• -k 指定用于索引构建的 k-mer 长度。默认值为 31，这是 k-mer 长度的最大值。k-mer 长度影响索引的精确性和内存使用量

# conda activate rna
kallisto index -k 31 -t 8 -i Mus_musculus.GRCm39.cdna.kalliso.idx  Mus_musculus.GRCm39.cdna.all.fa

速度很快，日志如下：

[build] loading fasta file Mus_musculus.GRCm39.cdna.all.fa [build] k-mer length: 31 [build] warning: clipped off poly-A tail (longer than 10) from 644 target sequences [build] warning: replaced 1 non-ACGUT characters in the input sequence with pseudorandom nucleotides KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2) KmerStream::KmerStream(): Start computing k-mer cardinality estimations (1/2) KmerStream::KmerStream(): Finished CompactedDBG::build(): Estimated number of k-mers occurring at least once: 100660025 CompactedDBG::build(): Estimated number of minimizer occurring at least once: 24681028 CompactedDBG::filter(): Processed 210946368 k-mers in 116116 reads CompactedDBG::filter(): Found 100551210 unique k-mers CompactedDBG::filter(): Number of blocks in Bloom filter is 688107 CompactedDBG::construct(): Extract approximate unitigs (1/2) CompactedDBG::construct(): Extract approximate unitigs (2/2) CompactedDBG::construct(): Closed all input files CompactedDBG::construct(): Splitting unitigs (1/2) CompactedDBG::construct(): Splitting unitigs (2/2) CompactedDBG::construct(): Before split: 680500 unitigs CompactedDBG::construct(): After split (1/1): 680500 unitigs CompactedDBG::construct(): Unitigs split: 1459 CompactedDBG::construct(): Unitigs deleted: 0 CompactedDBG::construct(): Joining unitigs CompactedDBG::construct(): After join: 630450 unitigs CompactedDBG::construct(): Joined 50326 unitigs [build] building MPHF [build] creating equivalence classes ... [build] target de Bruijn graph has k-mer length 31 and minimizer length 23 [build] target de Bruijn graph has 630450 contigs and contains 100639601 k-mers

fastq文件下载

项目编号：PRJNA1128028，使用 kingfisher 软件进行下载，万能的张永杰也在他的专辑中介绍过了：一行代码下载原始数据—Kingfisher，快去学习吧！

# 项目编号：PRJNA1128028
# fq 下载整个Bioproject
nohup kingfisher get -p PRJNA1128028 -m ena-ascp ena-ftp prefetch aws-http 1>down_PRJNA1128028.log 2>&1 &

数据比对

其实这里还差一步数据过滤，这里为了简单就先不过滤了。这个数据为单端，按照单端定量方式：

• -i 索引文件
• -o 输出目录
• -t 设定线程
• --plaintext 仅输出文本格式。默认是HDF5格式
• --single 指定输入的是单端测试数据
• -l 指定估计的平均片段长度
• -s 指定片段长度的标准偏差
• -l 和 -s 这些参数对于单端数据是必需的，因为无法从单端数据中直接估计它们

## 单端-小鼠测序数据
ref=/nas1/zhangj/database/genome/Mus_musculus/release113/Mus_musculus.GRCm39.cdna.kalliso.idx

ls *gz | while read id
do
  name=${id%.fastq.gz}
  echo "kallisto quant -i $ref -o ./${name}.quant -t 4 --single -l 100 -s 30 $id "
done >kallisto.sh

kallisto.sh 内容如下：

多任务并行神器：同时提交七个任务。ParaFly 这个软件又是张永杰介绍过的！额，我发现他没有介绍，这就去告诉他让他写个帖子！

用起来也是超简单，一个参数-c搞定！

nohup ParaFly -c kallisto.sh -CPU 7 1>kallisto.log  2>&1 &

这样就拿到了每个样本的定量结果：

tree *quant

SRR29546216.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546217.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546218.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546219.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546220.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546221.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json
SRR29546222.quant
├── abundance.h5
├── abundance.tsv
└── run_info.json

下游分析

得到每个样本的定量后，这个数据的abundance.tsv 就是 GSE270697 中提供的数据，薅一个样本看看每列：

GSE270697 这个数据对应的文章为：《Nicotinamide mononucleotide restores impaired metabolism, endothelial cell proliferation and angiogenesis in old sedentary male mice》。这个文章提供的代码在：https://github.com/Longchamp-Lab/kiesworo-et-al-cell-reports-2024，我看了下他的代码是如何处理这个数据的：https://github.com/Longchamp-Lab/Kiesworo_et_al_iScience_2024/blob/main/R%20scripts%20figure%204/Script_figPUB_KK.R，这里已经有了一个标准化后的cpm值。

问一下人工智能kimi

问：介绍 kallisto 软件定量后的这些列的含义：target_id length eff_length est_counts tpm

在使用 Kallisto 软件进行 RNA-seq 数据定量分析后，输出的 abundance.tsv 文件中包含以下列的含义：

说明

• 有效长度（eff_length）：在计算 TPM 时，有效长度用于校正转录本长度对表达量的影响。例如，如果转录本的实际长度为 2000 bp，插入片段长度为 225 bp，则有效长度为 1775 bp。
• TPM（tpm）：TPM 是一种标准化的表达量单位，通过将估计的计数（est_counts）除以有效长度并乘以 1e6 来计算。TPM 使得不同长度的转录本之间的表达量可以直接比较。
• 估计计数（est_counts）：这是经过校正后的计数，考虑了转录本长度和片段长度的影响。

再问：做差异分析应该用哪个值

在使用 Kallisto 定量结果进行差异表达分析时，通常推荐使用 TPM（Transcripts Per Million）值 或 est_counts（估计计数）。选择哪一个取决于具体的分析需求和下游工具的要求。 1. TPM（Transcripts Per Million）

• 含义：TPM 是一种标准化的表达量单位，通过将估计的计数（est_counts）除以有效长度并乘以 1e6 来计算。它校正了转录本长度的影响，使得不同长度的转录本之间的表达量可以直接比较。
• 优点：
  • 可比性：TPM 值在不同样本之间可以直接比较，不受转录本长度的影响。
  • 标准化：TPM 值经过了标准化处理，更适合用于差异表达分析。
• 缺点：
  • 非整数：TPM 值是浮点数，不适合某些需要整数计数的统计方法（如泊松分布或负二项分布）。
• 适用场景：
  • 如果你使用的是基于标准化表达量的差异表达分析工具（如 Sleuth），TPM 是更好的选择。
  • 如果你需要在不同样本之间比较表达量，TPM 更合适。

  2. est_counts（估计计数）

• 含义：est_counts 是经过校正后的计数，考虑了转录本长度和片段长度的影响。它是一个整数，表示每个转录本的预期计数。
• 优点：
  • 整数计数：est_counts 是整数，适合用于基于计数的统计方法（如 DESeq2、edgeR）。
  • 适合计数模型：这些工具通常基于泊松分布或负二项分布，整数计数更符合这些模型的假设。
• 缺点：
  • 未校正长度：est_counts 未校正转录本长度的影响，因此在比较不同长度的转录本时可能会有偏差。
• 适用场景：
  • 如果你使用的是基于计数的差异表达分析工具（如 DESeq2、edgeR），est_counts 是更好的选择。
  • 如果你需要进行下游的统计分析，如主成分分析（PCA）或聚类分析，est_counts 也更合适。

  总结

• 如果你使用 Sleuth 进行差异表达分析，推荐使用 **TPM**。
• 如果你使用 DESeq2 或 edgeR，推荐使用 **est_counts**。
• 如果你需要在不同样本之间比较表达量，**TPM** 更合适。
• 如果你需要进行基于计数的统计分析，**est_counts** 更合适。

在实际应用中，建议根据具体的分析工具和需求选择合适的值。

但是，est_counts不是一个整数，他又在这里乱说，我找到一个帖子，里面对 est_counts 使用 round 函数取了整数然后用DESeq2算法做的差异分析：

https://rpubs.com/zbengt/differential-expression-analysis-tutorial

可以看看这种方法！