nuget
Gibbs Sampler 是一种基于马尔可夫链蒙特卡洛（MCMC）思想的随机抽样算法，被广泛用于从一组序列中发现未知的保守 motif。其基本流程为： 假设每条输入序列中均包含一个长度为 W 的 motif 实例，算法首先在每条序列中随机选取一个长度为 W 的窗口作为 motif 位点的初始状态； 在随后每一轮迭代中，随机抽出一条序列并暂时移除，用其余序列中已定位的 motif 片段构建位置特异性评分矩阵（PWM），并结合背景模型计算 被移除序列中每个候选窗口作为 motif 位点的后验概率，再依据该概率分布随机抽样新的位点并放回；如此反复迭代更新各序列的 motif 位置， 直至全体位点趋于收敛。由于每一步均以概率采样（而非贪心选取最优）的方式更新，Gibbs Sampler 能有效规避 EM 类算法（如 MEME）易陷入 局部最优的问题；实践中通常设置伪计数和平滑项、从多个随机初始状态独立运行并保留总得分最高的比对结果，以提高发现真实 motif 的准确性与稳健性。
GibbsSampler
01 Syntax
SMRUCC.genomics.Analysis.SequenceTools.SequencePatterns.Motif.GibbsSampler
02 Methods
| Name | Overloads | Summary |
|---|---|---|
| .ctor | 1 | |
| find | 1 | Runs numSamples gibbsSamples to find a prediction on the sites and motifs with the highest information content in the sequences |
| findTopN | 1 | 迭代式（find -> mask -> resample）的多 motif 发现： 1. |
| DefaultRestarts | 1 | 在未显式指定重启次数时，依据序列规模推算一个合理的重启次数 |
| RunRestarts | 1 | 以 restarts 个互不相同的随机初始状态并行运行吉布斯采样， 各重启之间共享的最优结果汇总在 sampler 之中。 |
| MaskSites | 1 | 把本轮发现的位点窗口（两侧各外扩 pad 个碱基）就地屏蔽为 N。 屏蔽是「等长替换」：只改写字符而不改变序列长度，因为 GibbsSampler.m_sequenceLength、随机初始化以及候选起点的计算全部依赖于长度不变。 被屏蔽的位置 Utils.indexOfBase()) 会返回 -1， 计数矩阵与似然比计算都会自动跳过它们。 |
| BuildMotifResult | 1 | 依据一轮多重启采样的结果构造 MSAMotif； 当所有重启都没有产出有效结果时（例如信息含量恒为 NaN）返回 Nothing。 |
| SiteOdds | 1 | 计算某一条 motif 实例在 PWM 模型之下的概率 q 与在背景模型之下的概率 p。 两者都以「逐列几何平均」的形式给出，以避免较长 motif 的连乘下溢， 因此 MSAMotif.score（= q / p）是逐列的优势比(odds ratio)： 大于 1 表示该位点相比于随机背景更像是 motif 实例。 |
| EstimateEvalue | 1 | 对 motif 的显著性做启发式的保守估计： 1. |
| gibbsSample | 1 | Implements the Gibbs Sampling algorithm found in the lawrence93.pdf |
| predictiveUpdateStep | 1 | One of the sequenceLength sequences, z, is chosen either at random The pattern description q_{i,j} frequency is then calculated from the current positions a_k in all sequences… |
| getMotifStrings | 1 | Grabs the motif strings of length motifLength from each sequence and site |
| samplingStep | 1 | Every possible segment of width motifLength within sequence z is considered as a possible instance of the pattern. |
| calculateMotifProbability | 1 | calculates the log probability of a character appearing at a specific index in a motif |
| weightedChooseIndex | 1 | 依据 log 概率分布做轮盘赌随机采样（而非贪心地选取最大值）， 这正是吉布斯采样能够跳出局部最优的关键所在。 |
| getRandomSites | 1 | Creates a list of sequenceLength random numbers using the random object supplied the numbers are from 0 to sequenceLength-motifLength-1 inclusive |
03 Properties
| Name | Overloads | Summary |
|---|---|---|
| Sequences | 1 | populate all fasta FastaSeq.SequenceData in upper case. |
| SequenceCount | 1 |
04 Fields
| Name | Overloads | Summary |
|---|---|---|
| MIN_BACKGROUND | 1 | 背景概率的下限值：输入序列中完全没有出现过的碱基其背景频率为 0， 会在计算似然比 q/p 时产生除零(+Inf)，需要钳制到一个极小的概率。 |
| MAX_ICPC | 1 | 信息含量的理论上限：log2(4) = 2 bits/column（DNA 四个碱基） |
| ICPC_NOISE_FACTOR | 1 | 自适应信息含量阈值所使用的倍数：随机对齐的期望信息含量约为 Probability.E())（即有限样本带来的偏置）， 默认为其 5 倍，也就是要求 motif 的信息含量至少是随机噪声水平的 5 倍。 |
| MASKED_LOG_WEIGHT | 1 | 落在屏蔽区（窗口内含有 N）之上的候选起点所使用的 log 权重。 这里取一个足够小的常数而不是负无穷：若某条序列的候选起点在若干轮屏蔽之后 全部落入屏蔽区，采样分布依然是可归一化的（此时退化为均匀分布）。 |
| m_sequenceLength | 1 | sequence length of the input sequence collection |
05 Members
#ctor(
IEnumerable(Of FastaSeq), Int32)Parameters
| Name | Type | Description |
|---|---|---|
fastaFile | IEnumerable(Of FastaSeq) | the un-aligned raw sequence data, any sequence object with length less than the required motifLength will be ignored. |
motifLength | Int32 | recommended use value of 0.8 multiply of the average length of the fasta sequence |
find(
Int32)Runs numSamples gibbsSamples to find a prediction on the sites and motifs with the highest information content in the sequences
Parameters
| Name | Type | Description |
|---|---|---|
maxIterations | Int32 | maximum number of times to iterate in a Gibbs Sample |
Returns
信息含量最高的那一个 motif；当所有重启都没有产出有效结果时返回 Nothing。 若需要一次发现多个 motif，请改用 GibbsSampler.findTopN()。
findTopN(
Int32, Int32, Int32, Double, Double, Double)迭代式（find -> mask -> resample）的多 motif 发现：
- 在当前的序列集合上运行一轮标准的吉布斯采样，得到 motif Mk 及其在每条序列上的位点；
- 若 Mk 的信息含量或者 E-value 达不到阈值要求，则终止整个发现过程；
- 把 Mk 的所有位点窗口（两侧各外扩 maskPadding 倍 motif 宽度） 屏蔽为 N，然后重新随机初始化位点，在屏蔽之后的序列上继续发现下一个 motif。
由于每一轮都会把上一轮的发现屏蔽掉，同一个 motif 不会被重复发现， 各轮结果之间的位点窗口也互不重叠。发现的数量允许少于 topN： 一旦新的 motif 达不到阈值要求就立即终止，绝不会返回低质量的结果。
返回值按照「发现顺序」排列。需要注意每一轮面对的是不同的搜索空间 （越往后屏蔽区越大），因此各轮的信息含量并不保证严格单调递减： 某一轮若陷入了较差的局部最优，其后一轮完全可能找到信息含量更高的 motif。
Parameters
| Name | Type | Description |
|---|---|---|
topN | Int32 | 期望发现的 motif 数量上限 |
maxIterations | Int32 | maximum number of times to iterate in a Gibbs Sample |
restarts | Int32 | 每一轮发现所使用的随机重启次数；小于等于 0 时按照序列规模自动推算。 注意：一旦某个重启达到了信息含量的理论上限，其余重启会提前空转退出。 |
maskPadding | Double | 位点窗口两侧的屏蔽外扩量，单位为 motif 宽度的倍数（默认 ±w/2） |
icpcCutoff | Double | 单位列信息含量(bits/column)的下限，其理论上限为 log2(4) = 2.0。 信息含量的绝对值会随序列条数的增加而系统性下降（随机对齐的期望信息含量约为 |
evalueCutoff | Double | E-value 的上限，默认是 |
Returns
按发现顺序排列的 motif 数组，其长度允许小于 topN。
DefaultRestarts
在未显式指定重启次数时，依据序列规模推算一个合理的重启次数
以 restarts 个互不相同的随机初始状态并行运行吉布斯采样， 各重启之间共享的最优结果汇总在 sampler 之中。
MaskSites(
String(), List(Of Int32), Int32)把本轮发现的位点窗口（两侧各外扩 pad 个碱基）就地屏蔽为 N。
屏蔽是「等长替换」：只改写字符而不改变序列长度，因为 GibbsSampler.m_sequenceLength、随机初始化以及候选起点的计算全部依赖于长度不变。 被屏蔽的位置 Utils.indexOfBase() 会返回 -1， 计数矩阵与似然比计算都会自动跳过它们。
BuildMotifResult(RunSample,
Int32)依据一轮多重启采样的结果构造 MSAMotif； 当所有重启都没有产出有效结果时（例如信息含量恒为 NaN）返回 Nothing。
SiteOdds(SequenceMatrix,
String)计算某一条 motif 实例在 PWM 模型之下的概率 q 与在背景模型之下的概率 p。
两者都以「逐列几何平均」的形式给出，以避免较长 motif 的连乘下溢， 因此 MSAMotif.score（= q / p）是逐列的优势比(odds ratio)： 大于 1 表示该位点相比于随机背景更像是 motif 实例。
EstimateEvalue(SequenceMatrix)
对 motif 的显著性做启发式的保守估计：
- 逐列计算信息含量 bits（沿用 MSAMotif.CreateMotif() 的口径， 其中包含
Probability.E()的小样本校正）； - 以整个 motif 的 bits 之和作为得分（它是 PWM 之下的期望位点对数似然比，单位为 bit）， 再按 E = 候选位点总数 × 2^(-score) 折算。
对于以 bit 计量的对数似然比得分，P(随机位点得分 ≥ score) ≤ 2^(-score) 成立， 因此上面给出的是 E-value 的 Chernoff 上界，而不是严格意义上的 Karlin-Altschul E-value； 作为上界它偏保守（同样的 motif 会给出偏大的 E 值）。
gibbsSample(
Int32, List(Of String))Implements the Gibbs Sampling algorithm found in the lawrence93.pdf
Parameters
| Name | Type | Description |
|---|---|---|
maxIterations | Int32 | , maximum number of iterations sampling may take |
Returns
Sets of int predicting the position motifs located in each sequence
predictiveUpdateStep(
List(Of String), List(Of Int32))One of the sequenceLength sequences, z, is chosen either at random The pattern description q_{i,j} frequency is then calculated from the current positions a_k in all sequences excluding z
Parameters
| Name | Type | Description |
|---|---|---|
S | List(Of String) | , the sequences other than z |
A | List(Of Int32) | , the sites for the sequences other than z |
getMotifStrings(
IEnumerable(Of String), List(Of Int32))Grabs the motif strings of length motifLength from each sequence and site
Parameters
| Name | Type | Description |
|---|---|---|
S | IEnumerable(Of String) | , sequences |
A | List(Of Int32) | , sites |
Returns
sequenceCount motif strings
Every possible segment of width motifLength within sequence z is considered as a possible instance of the pattern. The probabilities Q_x of generating each segment x according to the current pattern probabilities q_{i,j} are calculated The weight A_x = Q_x/P_x is assigned to segment x, and with each segment so weighted, a random one is selected. Its position then becomes the new a_z.
Parameters
| Name | Type | Description |
|---|---|---|
z | SequenceMatrix | , sequence we are iterating through |
calculates the log probability of a character appearing at a specific index in a motif
Parameters
| Name | Type | Description |
|---|---|---|
q_ij | SequenceMatrix | , motif weight matrix |
z | String | , string of characters |
x | Int32 | , index of site in z |
P | List(Of Double) | , background frequencies |
Returns
log probability
weightedChooseIndex(
List(Of Double))依据 log 概率分布做轮盘赌随机采样（而非贪心地选取最大值）， 这正是吉布斯采样能够跳出局部最优的关键所在。
Parameters
| Name | Type | Description |
|---|---|---|
weightDistribution | List(Of Double) | , 每一个候选起点上的 log 似然比 |
Returns
new index of the site
getRandomSites
Creates a list of sequenceLength random numbers using the random object supplied the numbers are from 0 to sequenceLength-motifLength-1 inclusive
Returns
sequenceLength random ints
Sequences
populate all fasta FastaSeq.SequenceData in upper case.
SequenceCount
Returns
the size of the list sequences
MIN_BACKGROUND
背景概率的下限值：输入序列中完全没有出现过的碱基其背景频率为 0， 会在计算似然比 q/p 时产生除零(+Inf)，需要钳制到一个极小的概率。
MAX_ICPC
信息含量的理论上限：log2(4) = 2 bits/column（DNA 四个碱基）
ICPC_NOISE_FACTOR
自适应信息含量阈值所使用的倍数：随机对齐的期望信息含量约为 Probability.E()（即有限样本带来的偏置）， 默认为其 5 倍，也就是要求 motif 的信息含量至少是随机噪声水平的 5 倍。
MASKED_LOG_WEIGHT
落在屏蔽区（窗口内含有 N）之上的候选起点所使用的 log 权重。
这里取一个足够小的常数而不是负无穷：若某条序列的候选起点在若干轮屏蔽之后 全部落入屏蔽区，采样分布依然是可归一化的（此时退化为均匀分布）。
m_sequenceLength
sequence length of the input sequence collection