文章阅读目录大纲
零分布（null distribution）是指在假设零假设（null hypothesis）成立的情况下，某个统计量随机取值的概率分布。在统计假设检验中，我们通常提出一个零假设（例如“两组数据没有显著差异”或“观察到的模式仅由随机因素造成”），然后根据观测数据计算一个检验统计量。零分布描述了这个统计量在零假设为真时的分布情况。通过将实际观测到的统计量与零分布进行比较，我们可以计算出P-value：即在零假设下，出现等于或更极端观测结果的概率。如果P-value很小（例如低于预设的显著性水平α），我们就认为零假设不太可能成立，从而拒绝零假设，认为观测结果是统计显著的。
构建零分布的方法有多种，取决于问题的复杂性和可用的模型：
- 解析模型（参数模型）：当统计量的分布已知且可推导时，可以使用理论公式计算P-value。例如，二项分布、泊松分布、正态分布等都可用于描述特定零假设下的统计量分布。在这种情况下，零分布是解析已知的，可以直接计算P-value。
- 置换检验（Permutation Test）：当没有现成的理论分布或数据不满足参数假设时，可以通过随机置换来构建经验零分布。具体做法是将样本标签随机打乱，重新计算统计量，重复成千上万次，从而得到该统计量在“无真实效应”情况下的分布。然后根据实际观测统计量在这个经验分布中的位置计算P-value。置换检验不依赖数据分布假设，因此非常灵活，但计算成本较高。
- 蒙特卡洛模拟：与置换检验类似，但通过生成符合零假设的随机数据来构建零分布。例如，在序列分析中，可以生成大量与背景组成相同的随机序列，计算某模式出现的频率，从而得到该模式在随机情况下的分布。
无论采用哪种方法，零分布的核心思想都是量化“如果零假设成立，我们观测到的结果有多极端”。P-value就是衡量这种极端程度的指标。需要注意的是，P-value本身并不衡量零假设为真的概率，而是衡量数据与零假设的兼容程度。P-value越小，表示在零假设下观测到当前结果或更极端结果的概率越低，从而提供证据反对零假设。基于上面所描述的情景和定义，我们在进行编写代码进行零分布的Pvalue计算的公式非常的简单，如下所示：
在上面的公式中，n 是随机序列得分优于或等于观测得分的次数。分母是蒙特卡洛实验置换检验的次数。从上面的公式中不难看出，所想要表达的意思就是假设我们进行了10000次蒙特卡洛实验，然后假设出现了n等于若干次（例如n=1）的得分会高于我们的目标得分的结果。那计算出来的概率就是万分之一，这就是一个极罕见的事件。即随机生成的东西他的得分高于目标观测到的得分这样子的事件是一个极罕见的事件，也就是说，我们的目标观测结果是一种极显著的结果。
在生物信息学中，零分布和P-value计算被广泛应用于评估各种分析结果的统计显著性，例如motif发现、序列比对、通路富集分析以及宏基因组reads的k-mer分类等。
VisualBasic中实现一个零分布检验框架
在下面所示的一个VB.NET代码中，我们实现了一个基于蒙特卡洛模拟或置换检验的通用框架，用于计算零假设检验中的 P 值。在这段代码框架中，他的核心原理是通过计算机模拟生成大量的“零分布”数据，以此来量化观测到的统计量出现的概率。按照上面的零分布pvalue计算的原理，我们可以理解下面的代码实现逻辑为：
首先，我们看到的这个计算框架通过抽象方法 ZeroSet 和 Score 将算法逻辑与具体的业务场景解耦。ZeroSet 负责生成在零假设成立情况下的模拟数据集（即零分布样本），而 Score 则负责定义如何计算检验统计量。
在 Pvalue 函数的计算过程中，代码利用了 AsParallel 实现了并行计算，以高效处理成千上万次模拟运算。核心步骤是对模拟生成的统计量进行“极端性计数”：程序会遍历所有模拟样本，根据传入的备择假设类型——包括 Greater（右侧检验）、Less（左侧检验）或 TwoSided（双侧检验）——统计模拟统计量比实际观测统计量（score）更为极端的次数（变量 n）。例如，在双侧检验中，代码会比较统计量的绝对值，以衡量偏离零假设的程度。
在最后，代码框架使用了 (n + 1) / (Permutation + 1) 这一公式计算最终的 P 值。这里分子分母同时加 1 是为了进行连续性修正，确保在模拟次数有限时，P 值不会因为完全没有出现更极端情况而变为 0，从而保证了统计推断的稳健性；对于双侧检验，代码则直接将该单侧概率乘以 2。
Public MustInherit Class NullHypothesis(Of T)
Public ReadOnly Property Permutation As Integer
Sub New(Optional permutation As Integer = 1000)
_Permutation = permutation
End Sub
''' <summary>
''' generates the random set with <see cref="permutation"/> elements.
''' </summary>
''' <returns></returns>
Public MustOverride Function ZeroSet() As IEnumerable(Of T)
Public MustOverride Function Score(x As T) As Double
Public Function Pvalue(score As Double, Optional alternative As Hypothesis = Hypothesis.Greater) As Double
Dim zero As T() = ZeroSet.ToArray
Dim n As Integer
Select Case alternative
Case Hypothesis.Greater
' mu > mu0
n = Aggregate x As T
In zero.AsParallel
Let per_score As Double = Me.Score(x)
Where per_score >= score
Into Count
Case Hypothesis.Less
' mu < mu0
n = Aggregate x As T
In zero.AsParallel
Let per_score As Double = Me.Score(x)
Where per_score <= score
Into Count
Case Hypothesis.TwoSided
' 计算观测统计量的绝对值，用于双侧检验判断极端性
Dim abs_score As Double = std.Abs(score)
' 双侧检验：置换统计量的绝对值 >= 观测统计量的绝对值
n = Aggregate x As T
In zero.AsParallel
Let per_score As Double = Me.Score(x)
Let abs_per_score As Double = std.Abs(per_score)
Where abs_per_score >= abs_score
Into Count
Case Else
Throw New InvalidProgramException($"unknown alternative hypothesis: {alternative}!")
End Select
If alternative = Hypothesis.TwoSided Then
Return 2 * (n + 1) / (Permutation + 1)
Else
Return (n + 1) / (Permutation + 1)
End If
End Function
End Class零分布在生物信息学中的应用例子
基于上面所展示的一个零分布检验的框架，在这里主要讲解两个我经常会遇到的生物信息学工作作为例子，来讲解我么可以如何基于具体的问题来建立起对应的零分布检验模型。
Motif扫描显著性检验
现在假设我们通过一些算法建立起了一个基因转录调控结合位点的Motif的PWM矩阵模型，现在我们想要利用这个PWM矩阵模型从我们的一个目标基因组序列上找出所有的和这个PWM模型匹配的序列片段，用来做后面的调控网络的研究。我们可以基于下面的算法流程模块来实现这样子的一个基于PWM矩阵的扫描：首先基于 Smith-Waterman 局部比对算法进行匹配区域的查找。在这里我们将PWM 模型视为“查询序列”，将基因组目标序列（target）视为“数据库”。这样子，基于比对算法就可以输出了一个在基因组上找到与 PWM 模型局部相似度最高的若干区域的候选列表。当找到一个候选位点后，代码需要计算该位点序列与 PWM 模型的具体匹配得分。在 PWM 模型中，矩阵的每一列代表在该位置上出现某种碱基（A/T/C/G）的权重（通常是Log-odds值）。匹配得分就是序列上每个位置的碱基权重之和。在这里对于每一个候选的 Match 对象 m，代码提取出对应的 PWM 片段（motifSlice）和基因组上的序列片段（site），然后通过上述逻辑计算出 total，即该位点的观测得分（Observed Score）。最后，在得到了具体的计算得分之后，我们就可以利用置换检验或蒙特卡洛模拟的方法来计算Motif匹配的统计学显著性（P-value）。
在这里我们可以基于下面的代码函数来完成基于上面所描述的计算流程的PWM扫描的工作：
' 基于PWM的概率匹配，从目标基因组序列中找到Motif位点用于建立调控关系
Public Iterator Function ScanSites(PWM As Residue(), target As FastaSeq,
Optional cutoff# = 0.6,
Optional minW% = 6,
Optional pvalue_cut As Double = 0.05,
Optional identities As Double = 0.8,
Optional n As Integer = 500,
Optional top As Integer = 9) As IEnumerable(Of MotifMatch)
Dim chars As Char() = target.SequenceData.ToCharArray
Dim subject As Residue() = target.ToResidues
' define a general smith-waterman symbol comparer
' mapping the char in target sequence as the PWM column object
' example as: A mapping as [1 0 0 0]
' T mapping as [0 1 0 0]
' G mapping as [0 0 1 0]
' C mapping as [0 0 0 1]
Dim symbol As New GenericSymbol(Of Residue)(
equals:=Function(a, b) Compare(a, b) >= 0.85,
similarity:=AddressOf Compare,
toChar:=AddressOf Residue.Max,
empty:=AddressOf Residue.GetEmpty
)
' create a general smith-waterman(GSW) alignment algorithm
Dim core As New GSW(Of Residue)(PWM, subject, symbol)
' use PWM seuqnece to search target sequence
' find all best local alignment HSP as motif site candidates
Dim result As Match() = core.BuildMatrix.GetMatches(cutoff * core.MaxScore) _
.Where(Function(m) (m.toB - m.fromB) >= minW) _
.OrderByDescending(Function(a) a.score) _
.Take(top) _
.ToArray
Dim seqTitle As String = target.Title
Dim background As New ZERO(background:=NT _
.ToDictionary(Function(b) b,
Function(b)
Return target.SequenceData.Count(b) / target.Length
End Function))
' evaluate each local best alignment HSP result
For Each m As Match In result
Dim len As Integer = std.Min(m.toA - m.fromA, m.toB - m.fromB)
Dim motifSpan As New Span(Of Residue)(PWM, m.fromA, len)
Dim motifSlice As Residue() = motifSpan.SpanView
Dim site As New Span(Of Char)(chars, m.fromB, motifSpan.Length)
Dim motifStr As String = motifSpan.SpanCopy.JoinBy("")
Dim v As Double() = New Double(motifStr.Length - 1) {}
Dim total As Double = 0
Dim one As Vector = Vector.Ones(v.Length)
For i As Integer = 0 To motifStr.Length - 1
v(i) = motifSpan(i)(site(i))
total += v(i)
Next
Dim null As New NullTest(background, motifSlice, v.Length, permutation:=n)
Dim pvalue As Double = null.Pvalue(total, Hypothesis.Greater)
If pvalue >= pvalue_cut Then
Continue For
End If
Dim score2 As Double = one.SSM(v.AsVector)
If score2 < identities Then
Continue For
End If
' found a significant motif site match on target sequence
' as the TFBS candidates, this site will be used for
' re-construct of the TRN network
Yield New MotifMatch With {
.identities = score2,
.segment = site.SpanCopy.CharString,
.motif = motifStr,
.score1 = m.score,
.score2 = total,
.title = seqTitle,
.start = m.fromB,
.ends = m.toB,
.pvalue = pvalue
}
Next
End Function在上面的代码中，零分布检验的实现主要是基于NullTest对象中的置换检验来完成的。这是代码中最核心的统计学部分。代码通过模拟“零分布”来计算 P-value，判断这个高分是真实的生物学信号，还是随机碰巧产生的。在上面的代码中，使用的是基于背景模型的随机模拟。在下面的一段代码中，通过枚举所有的字符来统计目标序列中 A、T、C、G 的频率（target.SequenceData.Count(b) / target.Length）。然后基于这个生成的背景频率，ZERO 类利用这些频率构建了一个累积概率分布，这样子 NextSequence 方法就可以根据这个分布生成一条长度指定、完全随机的“虚假DNA序列”。
Public Class ZERO
ReadOnly nucleotides As Char()
ReadOnly cumulativeProbs As IReadOnlyCollection(Of Double)
Sub New(background As Dictionary(Of Char, Double))
Dim cumulativeProbs As New List(Of Double)()
Dim nucleotides As Char() = background.Keys.ToArray
' 构建累积概率分布
Dim cumulative As Double = 0
For Each NT As Char In background.Keys
cumulative += background(NT)
cumulativeProbs.Add(cumulative)
Next
Me.nucleotides = nucleotides
Me.cumulativeProbs = cumulativeProbs
End Sub
Public Function NextSequence(length As Integer) As String
' 生成随机序列
Dim sequence As Char() = New Char(length - 1) {}
For i As Integer = 1 To length
Dim rndValue As Double = randf.NextDouble()
' 选择核苷酸
For j As Integer = 0 To nucleotides.Length - 1
If rndValue <= cumulativeProbs(j) Then
sequence(i - 1) = nucleotides(j)
Exit For
End If
Next
Next
Return New String(sequence)
End Function
End Class
Dim background As New ZERO(background:=NT _
.ToDictionary(Function(b) b,
Function(b)
Return target.SequenceData.Count(b) / target.Length
End Function))ZERO 类代表背景分布模型。它告诉我们，如果在这个基因组上随机生成一段序列，每个位置出现各个碱基的概率是多少。这通常是基于0阶马尔科夫模型（即各碱基独立分布）。这样子，基于上面的代码所生成的随机序列，每一条序列我们都可以按照与真实的序列扫描结果一样的计算方法同样也计算出一个评分。这样子就得到了一个由设定的置换检验次数个随机得分组成的分布，这就是零分布。基于此，我们就可以做出下面的检验假设：
- 零假设 (H0)：观察到的基因组序列是随机的，其碱基出现频率符合该基因组的背景分布。在这个假设下，任何与 PWM 的匹配都只是随机噪声。
- 备择假设 (H1)：观察到的匹配得分显著高于随机背景，说明存在真实的 Motif 位点。
至于，如何基于这个零分布计算出Pvalue，在上面的代码中，将真实序列的观测得分与零分布进行比较。然后统计在零分布中，随机得分 大于等于 观测得分的次数。然后通过比较分析发现，如果 P-value 很小（例如 < 0.05），意味着在随机背景下生成这么高分序列的概率极低，从而拒绝原假设，认为该Motif匹配是显著的。
Public Class NullTest : Inherits NullHypothesis(Of String)
ReadOnly length As Integer
ReadOnly zero As ZERO
ReadOnly motifSlice As Residue()
Sub New(zero As ZERO, motifSlice As Residue(), length As Integer, Optional permutation As Integer = 1000)
Call MyBase.New(permutation)
Me.motifSlice = motifSlice
Me.zero = zero
Me.length = length
End Sub
Public Overrides Iterator Function ZeroSet() As IEnumerable(Of String)
For i As Integer = 1 To Permutation
Yield zero.NextSequence(length)
Next
End Function
<MethodImpl(MethodImplOptions.AggressiveInlining)>
Public Overrides Function Score(x As String) As Double
Return score(x.ToCharArray, motifSlice)
End Function
Friend Overloads Shared Function score(seq As String, PWM As IReadOnlyCollection(Of Residue)) As Double
Dim total As Double = 0
For i As Integer = 0 To seq.Length - 1
total += PWM(i)(seq(i))
Next
Return total
End Function
End Class按照上面的算法步骤描述，我们可以总结出如下的如何使用零分布做Motif扫描结果的显著性的检验方法过程：
- 建模：根据实际基因组序列计算碱基频率，建立背景概率模型 (ZERO)。
- 模拟：基于背景模型，生成大量（如500条）随机序列 (NullTest.ZeroSet)。
- 评分：用PWM模型给这些随机序列打分，构成一个背景得分的分布。
- 检验：看真实位点的得分在这个背景分布中是否处于极端位置（右侧尾部）。如果是，则认为该位点生物学显著。
GCModeller中做基于PWM模型的motif位点扫描例子
require(GCModeller);
imports "bioseq.patterns" from "seqtoolkit";
imports "bioseq.fasta" from "seqtoolkit";
imports "GenBank" from "seqtoolkit";
imports "genomics_context" from "seqtoolkit";
# read sequence data and found motif via gibbs scan method
let raw = read.fasta("Staphylococcaceae_LexA___Staphylococcaceae.fasta");
let motif = gibbs_scan(raw, width = 18);
# draw sequence logo of the generated motif
bitmap(file = "LexA.png") {
plot(motif, title = "LexA");
}
let gb_asm = GenBank::read.genbank("Escherichia coli str. K-12 substr. MG1655.gbff");
let nt = gb_asm |> TSS_upstream( upstream_len =150);
let motif_sites = motif.find_sites(motif, nt, parallel = TRUE, motif_name="LexA");
# cast motif site data result as dataframe and export to table file
motif_sites = as.data.frame(motif_sites);
motif_sites = motif_sites[order(motif_sites$identities, decreasing=TRUE),];
print(motif_sites, max.print = 13);
write.csv(motif_sites, file = "LexA_sites.csv");
GSEA富集计算分析中的零分布检验
在比较两个样本（例如 癌组织 vs 正常组织，A vs B）时，我们通常会得到成千上万个基因的差异表达列表。如果只关注单个基因（例如基因P53上调了），往往会陷入细节，且受实验噪声影响大。我们利用富集分析可以将具有相似功能的基因打包成一个集合（例如一个KEGG pathway，比如“细胞周期通路”、“DNA修复通路”）。我们通过分析这一组基因是否整体倾向于排在列表的前面，然后就可以基于此推断出一些生物学结论：“在这个生物学过程中，发生了一些系统性的协调变化。”
针对基因表达结果，我们可以建立如下的数据结构用来进行建模。在下面代码中的GeneExpressionRank对象定义中，rank 通常代表了基因的变化程度（如 Log2 FoldChange）或统计显著性（如 -log10 P-value）。Rank 越高，变化越剧烈。那基于这个rank值，我们可以对基因做出一些排序，然后基于排序的结果做出一些计算得分，然后通过得分来得到一些生物学解释。这样子，这里所说的得到生物学解释的得分计算过程就是一个富集计算的过程。
Public Class GeneExpressionRank
Implements INamedValue
Public Property gene_id As String Implements INamedValue.Key
''' <summary>
''' gene id tagged with the expression ranking value, example as -log10 of t-test pvalue
''' </summary>
''' <returns></returns>
Public Property rank As Double
End Class从上面的Motif注释的代码学习中我们已经很清楚的了解到，在统计学假设检验中，“零分布”是指在零假设成立的情况下，统计量（这里是富集得分）的分布情况。在这里我们所想要描述的GSEA富集分析问题中，零假设（H0）是：基因的表达排序（rank）与该通路中的基因成员身份无关。换句话说，该通路中的基因在排序列表中是随机分布的，没有聚集在顶部（高表达）或底部（低表达）。为了验证这一假设，GSEA分析使用了置换检验。
基于上面所定义的基因表达结果，我们可以实现如下所示的一个基于置换检验的基因集富集分析算法：
Public Module GSEACalculate
<Extension>
Public Iterator Function Enrichment(background As Background,
geneExpression As GeneExpressionRank(),
Optional permutations As Integer = 1000) As IEnumerable(Of EnrichmentResult)
Dim geneInputs As String() = geneExpression.Keys
For Each pathway As Cluster In TqdmWrapper.Wrap(background.clusters, wrap_console:=App.EnableTqdm)
Dim geneSet As Index(Of String) = pathway.memberIds
Dim enrich = geneExpression.Enrich(geneSet, permutations)
Dim intersect = geneInputs.Intersect(geneSet.Objects).ToArray
Yield New EnrichmentResult With {
.cluster = geneSet.Count,
.description = pathway.description,
.enriched = intersect.Length,
.IDs = intersect,
.name = pathway.names,
.pvalue = enrich.pvalue,
.score = enrich.score,
.term = pathway.ID
}
Next
End Function
''' <summary>
''' Enrichment for one specific pathway
''' </summary>
''' <param name="geneExpression">
''' all gene mean expression value, usually be the foldchange value of the a vs b comparision result.
''' </param>
''' <param name="geneSet">
''' gene set in current pathway
''' </param>
''' <param name="permutations"></param>
''' <returns></returns>
<Extension>
Public Function Enrich(geneExpression As GeneExpressionRank(),
geneSet As Index(Of String),
Optional permutations As Integer = 1000) As (score As Double, pvalue As Double)
Dim score As Double = geneExpression.enrich_score(geneSet)
Dim nulltest As New PermutationTest(geneSet, geneExpression, permutations)
Dim pval As Double = nulltest.Pvalue(score, Hypothesis.Greater)
Return (score, pval)
End Function
<Extension>
Friend Function enrich_score(geneExpression As GeneExpressionRank(), geneSet As Index(Of String)) As Double
Dim sortedGenes = geneExpression.OrderByDescending(Function(a) a.rank).ToArray
Dim enrichmentScore As Double = 0
Dim maxEnrichmentScore As Double = Double.MinValue
' --- 计算分数所需的参数 ---
Dim totalGenes As Integer = sortedGenes.Length
Dim hitsCount As Integer = geneSet.Count
' 预先计算“命中”和“未命中”时的分数步长，避免在循环中重复除法
Dim hitStep As Double = 1.0 / hitsCount
Dim missStep As Double = 1.0 / (totalGenes - hitsCount)
For Each gene As GeneExpressionRank In sortedGenes
If gene.gene_id Like geneSet Then
' 使用归一化的分数 (1 / N_hits)
enrichmentScore += hitStep
If enrichmentScore > maxEnrichmentScore Then
maxEnrichmentScore = enrichmentScore
End If
Else
' 使用归一化的分数 (1 / (N_total - N_hits))
enrichmentScore -= missStep
End If
Next
Return maxEnrichmentScore
End Function
End Module在上面所实现的富集得分计算函数 enrich_score 函数中，这个函数会按照基因的 rank（通常是变化倍数或显著性指标）从高到低遍历，实现了一个针对基因按照表达排序后的fractions得分统计量的计算。如果基因属于目标通路 geneSet，分数增加；如果不属于，分数减少。这个过程会记录一个最大的偏离值 maxEnrichmentScore，这就是我们要检验的观测统计量。那现在我们已经得到了一个具体的观测值了，下面就可以通过构建置换检验来生成零分布数据，和我们现在得到的这个观测值做比较，完成统计检验计算。在下面的置换检验代码模块中，ZeroSet函数是 PermutationTest 类的核心工作。为了判断观测到的得分是“真实显著”还是“运气好碰巧得到的”，我们需要在ZeroSet函数中模拟一个“没有生物学关联”的世界。例如可以随机打乱基因表达的rank标签。
Public Class PermutationTest : Inherits NullHypothesis(Of GeneExpressionRank())
ReadOnly geneExpression As GeneExpressionRank()
ReadOnly geneSet As Index(Of String)
Sub New(geneSet As Index(Of String), geneExpression As GeneExpressionRank(), permutations As Integer)
MyBase.New(permutations)
Me.geneExpression = geneExpression
Me.geneSet = geneSet
End Sub
Public Overrides Iterator Function ZeroSet() As IEnumerable(Of GeneExpressionRank())
For n As Integer = 1 To Permutation
' make copy of the raw array
Dim permutedGeneExpression = geneExpression.ToArray
For i As Integer = 0 To permutedGeneExpression.Length - 1
Dim k As Integer = randf.NextInteger(permutedGeneExpression.Length)
Dim swapRank = permutedGeneExpression(k).rank
Dim temp = permutedGeneExpression(i).rank
permutedGeneExpression(i) = New GeneExpressionRank With {
.gene_id = permutedGeneExpression(i).gene_id,
.rank = swapRank
}
permutedGeneExpression(k) = New GeneExpressionRank With {
.gene_id = permutedGeneExpression(k).gene_id,
.rank = temp
}
Next
Yield permutedGeneExpression
Next
End Function
Public Overrides Function Score(x() As GeneExpressionRank) As Double
Return x.enrich_score(geneSet)
End Function
End Class在上面的置换检验代码中，我们保留了基因ID的位置，但是将对应的 rank 值进行随机交换。通过交换 rank，我们破坏了基因与其表达变化之间的真实联系。在打乱后的数据中，某个基因是否属于某个通路，跟它现在持有的 rank 值大小毫无关系。这个过程重复 Permutations 次（在上面的代码中，代码默认为1000次），每次都计算出一个在“随机情况下的富集得分”。这1000个得分就构成了零分布。最后，代码将真实的观测得分放到零分布中去比较。
GCModeller中做富集计算分析
在下面的一段R#脚本代码展示了我们如何针对某一个基因表达矩阵做limma差异检验，然后针对差异检验结果做通路富集计算分析的流程：
require(GCModeller);
imports ["geneExpression", "sampleInfo"] from "phenotype_kit";
imports ["background", "GSEA"] from "gseakit";
let expr = load.expr("expression_matrix.csv")
|> filterZeroGenes()
|> impute_missing(by.features=TRUE)
|> totalSumNorm()
|> geneExpression::relative(median = TRUE)
;
let sample_info = read.csv("sample_metadata.csv", row.names = NULL,
check.names = FALSE);
sample_info = sampleInfo(
ID = sample_info$Sample,
sample_info = sample_info$Group
);
cat("\n");
message("inspect of the sample group information:");
print(as.data.frame(sample_info), max.print = 6);
let design = sample_info |> make.analysis(
control = "Strain_Y203",
treatment = "High_Osmolarity");
let deg = limma(expr, design);
let i = ([deg]::class == "up") || ([deg]::class == "down");
cat("\n");
message("view of the DEG analysis result which is measured via the limma:");
print(as.data.frame(deg), max.print = 6);
deg = deg[i];
let geneSet = JSON::json_decode(readText("gene_clusters.json"));
let kb = background::fromList(geneSet);
let result = kb |> enrichment([deg]::id, expression = -log10([deg]::adj_P_Val),
permutations = 3000);
result = as.data.frame(result);
result[,"class"] =NULL;
result[,"category"]=NULL;
cat("\n");
message("get GSEA analysis output:");
print(result);
write.csv(result, file = "Strain_Y203_vs_High_Osmolarity.csv");
- 建立KEGG的KO序列数据库 - 2026年1月4日
- 环境微生物群落GEMs建模综述 - 2025年12月29日
- 【虚拟细胞】翻译事件建模 - 2025年12月21日
No responses yet