企业官网建设流程全解析

背景

最近看到一篇公众号推文，讲的是“转录组原始数据修改软件”，

据称能够随意调节转录组原始数据某个基因的表达，大肆吹嘘有多有用（仅供科学教研使用）。

我一看，这不就是精细一点的NGS simulator吗，十多年前搞NGS模拟数据分析的时候就有人开发过类似idea的软件了，当然不是说RNA-seq。

然后涉及到其他的测序技术，时代往后一点，搞这种人工合成仿真的也不少，起码也是五六年前起步了，

这种玩意居然还能拉出来炒作一波，还能搞个收费服务，(─.─|||。

我这里简单总结一下，从 RNA-seq 的 TPM（或者 Count）表达矩阵“逆向”生成 FASTQ 原始测序文件，在生物信息学中是一项非常成熟的技术。它不需要用到 diffusion、VAE、flow matching 等那种复杂的深度学习“生成模型（Generative Models）”，而是基于统计学抽样和确定性的序列拼接算法

从修改 TPM 值这一步开始，一步一步逆向操作到 FASTQ

⚠️ 本文中涉及到RNA-seq分析操作具体细节的，我这里模糊处理+用粗糙例子演示，毕竟好久不做组学分析了，以原理演示为准；

顺带找了一个例子软件作为参考：polyester

核心流程图（逆向工程）

目标火山图（你想长什么样） ↓ 修改 TPM 值（某些基因上调/下调） ↓ TPM → Count（逆向标准化） ↓ Count → 负二项分布抽样（加噪声） ↓ 得到"伪造"的计数矩阵 ↓ 从参考基因组切序列 → 片段化 → 加错误 → FASTQ ↓ 送去重新分析 → 得到你想要的火山图 ✓

第一步：修改 TPM 值

假设原始数据

火山图：所有点都在中间，无聊。

假设你想画成的样子

修改 TPM

对照组 TPM 不变，处理组 TPM 改： - 基因A: 对照=10, 处理=10 × 2^3 = 80 (上调8倍) - 基因B: 对照=20, 处理=20 × 2^(-2) = 5 (下调到1/4) - 基因C: 对照=5, 处理=5 (不变)

第二步：TPM → Count（最关键的一步）

为什么需要这一步？

TPM 是相对值（总和=10^6，仅作演示），Count 是绝对 reads 数。测序仪只认识绝对数。

公式

TPM i = Count i / L i ∑ j ( Count j / L j ) × 10 6 \text{TPM}_i = \frac{\text{Count}_i / L_i}{\sum_j (\text{Count}_j / L_j)} \times 10^6TPMi​=∑j​(Countj​/Lj​)Counti​/Li​​×106

逆向求解 Count：

Count i = TPM i × L i × N 10 6 \text{Count}_i = \text{TPM}_i \times L_i \times \frac{N}{10^6}Counti​=TPMi​×Li​×106N​

其中：

• L i L_iLi​= 转录本长度
• N NN= 总 reads 数（你设定的测序深度，比如 1000万条）

具体计算

假设总 readsN = 10 , 000 , 000 N = 10,000,000N=10,000,000，基因A长度L A = 3000 L_A = 3000LA​=3000bp：

验证：对照组总 Count = 700，处理组总 Count = 2575。不完美等于N NN，因为 TPM 已经归一化过了，需要缩放让总和合理。

实际代码中会更复杂（考虑有效长度、library文库大小因子），但核心就是把相对 TPM 变回绝对 Count。

第三步：Count → 带噪声的计数矩阵（NB分布登场）

为什么需要噪声？

真实数据有生物学重复间的波动。如果每个重复都恰好是 2400，太假了。

这里说的重复的意思，是指样本的重复，前面我们演示的是两个样本，1个对照，1个处理；其实真实测序之类大家应该都是做多个重复的，所以比如说3组重复，就是3个对照 vs 3个处理，那么比如说我处理组前面算出来 2400 ，这个是绝对数值的count，这个count，一般是均值意义，但是我们如果3个处理组的样本都是 2400、2400、2400就太假了。

总而言之，这里说的是从tpm求得到的count的绝对值是作为1个均值，然后需要在多个重复样本之间做采样还原分布。

负二项分布的原理我们就不讲了

NB分布抽样（gene/转录本i，样本j）

Count i j observed ∼ NB ( μ = Count i j target , size ) \text{Count}_{ij}^{\text{observed}} \sim \text{NB}(\mu = \text{Count}_{ij}^{\text{target}}, \text{size})Countijobserved​∼NB(μ=Countijtarget​,size)

对基因A处理组（目标均值 2400）：

重复1: NB(2400, size=800) → 抽得 2356 重复2: NB(2400, size=800) → 抽得 2489 重复3: NB(2400, size=800) → 抽得 2210

size 参数控制方差：

• size 大 → 方差小，重复间很接近（太假）
• size 小 → 方差大，重复间波动大（像真实实验）

前面举例说的 Polyester，默认 用 size = μ/3，即 800，属于偏小的方差（理想化数据）。

当然仅作为演示示例参考，真实数据不表。

⚠️ 再次强调，我们只是模糊演示，对于真实的 差异表达分析经典软件中的NB分布采样实现，比如说edgeR、DESeq2，它们的真实假设，是否是“同一个条件下，重复来自同一个分布，有共同的均值和离散度”，这里没有细究。

总而言之，NB 抽样就是制造符合这个假设的假数据，让下游分析工具"以为"这是真实的生物学重复。

第四步：计数矩阵 → FASTQ（物理模拟）

现在你有哪些东西呢：

比如说处理组，

对每个基因i、每个重复j，执行：

基因A，重复1，需要造 2356 条 reads ↓ 转录本A序列（3000bp，从参考基因组+GTF拼接） ↓ 循环 2356 次： 抽片段长度 F ~ N(250, 25) → 比如 F=248 抽起始位置 S ~ Uniform(1, 3000-248+1) → 比如 S=500 切片段 [500:747] 左端 read = [500:649]（150bp，正向） 右端 read = reverse_complement([598:747])（150bp） 逐碱基加错误（Bernoulli 0.5%） 赋质量值 写入 sample_01_1.fasta, sample_01_2.fasta

⚠️ 我们耳熟能详的NB分布，也就是RNA-seq 涉及到的原理部分，发挥作用的主要地方就在于从前面 相对tpm值反推绝对count值，然后从这个count的绝对值，这个均值中采样多个重复组的各自count中。

说白了，就是只用于计数矩阵的生成。

然后上面的循环中说到的另外两个采样，也就是 reads片段的采样、reads测序位置的采样，这两个是和NB分布无关的，完全是为了模拟一个物理测序过程。当然，物理测序过程也是比较复杂的，我们这里是用简单的正态分布和均匀分布来演示。

至于为什么这里还需要这些抽样来模拟，是因为真实测序不是从基因头测到尾，是随机打断后测两段的。

也就是说前面NB分布是数的问题，后面是造序列的问题。

这里稍微提一下碱基质量值如何模拟，也就是Phred质量值，

这里也不赘述质量值和错误率的数学关系了，

要如何模拟Q值，其实就是要模拟碱基测序的一个错误率分布关系，

最粗暴的方法就是均匀模型，也就是设置1个全局错误率（比如说p=0.005，也就是0.5%），

就是位置无关的均匀分布，大家错误率一致，

那么这个时候，对于每一个碱基，决定是否出错，就是一个非常简单的伯努利分布了，

如果是测错了，那么就把这个碱基替换成其他三个碱基即可，random.choice 的replace

总而言之，fastq这4行都能凑起来

第五步：送去"重新分析"

伪造的 FASTQ ↓ [质控 FastQC] → 通过（因为错误率、质量分布都是按真实模型造的） ↓ [比对 HISAT2/STAR] → 比对到参考基因组 ↓ [定量 Salmon/RSEM] → 得到新的 TPM/Count ↓ [差异分析 DESeq2/edgeR] → 火山图 ↓ 基因A: log2FC ≈ +3, p < 0.001 ✓ 基因B: log2FC ≈ -2, p < 0.001 ✓ 基因C: log2FC ≈ 0, p > 0.05 ✓

完美符合我们预设的火山图，

当然了，偷懒一点的做法就是，修改好tpm值，然后按照前面说的逆向生成fastq测序数据，做的话直接按照修改之后的目标tpm去做各种下游分析

为什么 FastQC 查不出来？

FastQC 只能查"数据质量是否像真的"，不能查"这些 reads 是否真的从某个细胞里测出来的"。

更深一步

这里简单涉及一下EM，因为话题和朋友聊的时候扯到了

隐变量结构与EM算法的适用性

RNA-seq确实存在多层隐变量结构

真实生物学状态 (θ) ← 我们真正关心的 ↓ 泊松/伽马抽样 表达强度 (λ) ← 不可直接观测 ↓ NB抽样 (size参数控制过散) 观测Count (x) ← 实际测到的 ↓ 测序深度/长度标准化 TPM/FPKM ← 相对丰度估计

经典RNA-seq工具确实在用类似EM的思想

RSEM是最典型的例子：一个读段可能比对到多个转录本，EM迭代：

• E步：计算读段属于各转录本的后验概率
• M步：更新各转录本表达量

"逆向合成"场景其实有其特殊性

正向 vs 逆向问题的隐变量角色不同

逆向合成时，隐变量被显式控制了

目标TPM (人为设定) ↓ 确定性计算 目标Count = f(TPM, L, N) ← 没有随机性，是计算不是估计 ↓ NB抽样 (μ=目标Count) 观测Count ~ NB(μ, size) ← 这里引入随机性，但μ已知 ↓ 物理模拟 FASTQ ← 完全可控的生成过程

关键点：逆向合成中，"真实值"是我们人为指定的，不需要从数据中推断。EM算法解决的是"不知道真实值，要从观测推断"的问题。

那么，EM在哪里，或者说需要在哪里？

如果要让伪造数据更"真实"，确实可以引入更深层的隐变量

当前polyester类工具的简化假设：

给定μ → NB抽样 → Count

更真实的层次模型（可用变分推断或MCMC）：

样本条件效应 β ~ Normal(0, σ²_条件) ← 批次/处理效应 基因特异性效应 α_g ~ Normal(0, σ²_g) ← 基因间差异 样本-基因交互 γ_{sg} ~ Normal(0, σ²_γ) ← 真实生物学变异 log(μ_{sg}) = 偏移 + β_s + α_g + γ_{sg} ← 对数线性模型 Count_{sg} ~ NB(μ_{sg}, size_g) ← 观测

这种深层模型下，EM/变分推断的价值

一个更深刻的视角

我们前面描述的"修改TPM→逆向生成FASTQ"本质上是确定性反演+随机前向模拟的混合：

确定性部分: TPM → Count (可精确反解) 随机部分: Count → 带噪Count → FASTQ (需要合理的生成模型)

如果要从数学原理上最严谨地伪造数据，应该：

1. 用真实数据集训练一个深度生成模型（如VAE、扩散模型、或层次贝叶斯模型）
2. 隐空间插值/条件生成得到"看似真实但表达被调控"的数据
3. 这比手工调TPM更不易被检测，因为协方差结构、基因间相关性都是真实的

"目标表达变化"是伪造者自己定的；"让变化看起来自然的噪声结构"才是需要从真实数据学习的——EM/变分推断学的是后者，不是前者。

总结：至于吗？

总而言之，根本不需要什么生成模型，

就是看你组学测序背后的统计学原理了解的有多少。

那么回过头来，怎么评价这件事：

1个简单的NGS simulator？实现原理本身不难，写code的话有了AI更是直接工程化一键走起，

就这玩意还有人吹捧，还有人打着幌子搞收费分析？

怎么搞，无非是对那些半吊子做计算生物的新手、菜鸟，忽悠一下，

比如说：“哎呀，你是不是不愿意学RNA-seq，是不是原理搞不清楚，这些通通不用管！哎呀，你跑网上的/ai给的脚本是不是拿不到你想要的差异结果？哎呀，你是不是想要这个gene高一点，那个gene低一点。没事，用我的软件，保准出图让你满意，甚至连上游fastq都能以假乱真”。

一句话：能被菜鸟割技术韭菜的，本身也差不多是菜鸟，都是不好好静下心来学原理所背技术债的锅。