diff --git a/docs/Algorithm/NW.mdx b/docs/Algorithm/NW.mdx new file mode 100644 index 0000000..8a94287 --- /dev/null +++ b/docs/Algorithm/NW.mdx @@ -0,0 +1,246 @@ +--- +title: Smith-Waterman-Algorithm +description: Smith-Waterman algorithm, a local sequence alignment algorithm +tags: [algorithm, code, bioinfo] +sidebar_position: 5 +--- + +import ReferenceList from "@site/src/components/ReferenceList"; +import wikipedia from "@site/static/img/icon/wikipedia.png"; + + + +## 概念 + +Smith-Waterman 算法是一种局部序列比对算法，用于在两个序列之间寻找最佳的局部相似区域。与 Needleman-Wunsch 算法类似，Smith-Waterman 算法也是基于动态规划的思想，但是它不允许负分数的出现 (负分替换为 0)，因此可以找到局部相似区域。 + +## 理解 + +要找到两个序列的最佳比对，也就是说要找到一条路径，使得路径上的得分最大。 + +为什么要用动态规划的思想： + +比如我们要找到两个城市之间的最短距离，两个城市之间有若干个中间城市，如下图所示： + + +![distance](https://raw.githubusercontent.com/wjwei-handsome/wwjPic/main/img/20240430151241.png) + +我们要找到红色点到橙色点的最短距离，可以把问题简化成红色点分别到 A,B,C 点的最短距离 (如果 Sa, Sb, Sc 相同的话), 以此往前类推，把问题分解成子问题，最后得到最优解。 + +## 步骤解释 + +### 1. 确定置换矩阵和空位罚分 + +#### 置换矩阵 + +置换矩阵赋予每一碱基对匹配或错配的分数，相同或类似则赋予正值，不同或不相似赋予 0 分或者负分。 + +以核苷酸比对为例，若匹配 (match) 为 1 分，错配 (mismatch) 为 -1 分，则置换矩阵如下： + +| | A | T | C | G | +|---|---|---|---|---| +| A | 1 | -1 | -1 | -1 | +| T | -1 | 1 | -1 | -1 | +| C | -1 | -1 | 1 | -1 | +| G | -1 | -1 | -1 | 1 | + +也可以表示为： + +$s(a_{i},b_{j})={\begin{cases}1,&a_{i}=b_{j}\\-1,&a_{i}\neq b_{j}\end{cases}}$ + +氨基酸序列比对的置换矩阵一般更加复杂。通常性质相似的残基对具有正分数，反之，不相似的具有负分数。 + +#### 空位罚分 + +空位罚分决定了引入或延长空位的分值。 + +空位罚分决定了插入或者删除的分值。最基本的空位罚分方式为每次插入或者删除的得分相同。 + +但是从生物学意义上来讲，插入或缺失相比碱基的错配，会带来更严重的影响，比如移码突变等，因此空位罚分一般会比错配分数更高。 + +另外，单个基因突变事件可能导致一长串空位的插入。因此，一个连续的较长的空位优于多个分散的小的空位。 + +虽然多个分散的小的空位可以产生更多匹配，但一个连续的较长的空位代表这个区域只在一个序列中出现，使用后者可以避免为了得到高分而过度匹配这段序列。 + +实现该方法只需要引入**空位起始罚分**和**空位延长罚分**的概念。空位起始罚分通常高于空位延长罚分。 + +##### 空位权值恒定模型 + +该模型空位的罚分正比于空位的长度。 + +##### 空位延伸罚分模型 + +$W_{k}=u(k-1)+v\quad (u>0,v>0)$ + +该模型考虑空位起始罚分和空位延长罚分，其中$v$为开始空位的罚分，$u$为每次延长空位的罚分。例如，一个长度为 2 的空位的罚分为$u+v$, 长度为 10 的空位的罚分为$9u+v$。 + + +### 2. 初始化得分矩阵 + +得分矩阵的长度和宽度分别为两序列的长度 +1。其首行和首列所有元素均设为 0。 + +额外的首行和首列得以让一序列从另一序列的**任意位置**开始进行比对，分值为零使其不受罚分。 + + +### 3. 计算得分矩阵 + +对得分矩阵的每一元素进行从左到右、从上到下的打分，考虑匹配或错配（对角线得分），引入空位（水平或垂直得分）分别带来的结果，取最高值作为该元素的分值。 + +如果分值低于 0，则该元素分值为 0。打分的同时记录下每一个分数的来源用来回溯。 + + +举例： + +给定两条序列 TGTTACGG 和 GGTTGACTA. 并使用如下置换矩阵和空位罚分： + +- $s(a_{i},b_{j})={\begin{cases}1,&a_{i}=b_{j}\\-1,&a_{i}\neq b_{j}\end{cases}}$ + + +- $W_{k}=kW_{1}, W_{1} = 2$ + +初始化得分矩阵，然后进行打分，前三个碱基如下图所示： + +其中黄色表示正在计算的两个碱基，黑色箭头和箭头上的分数表示分数来源和对应分数，得到矩阵上的值为计算不同来源得到的最大值。 + +比如第一对碱基，有三个方向： + +- 上方的箭头为$0-2=-2$ +- 左方的箭头为$0-2=-2$ +- 左上角对角线的箭头为$0-3=-2$ + +这里最大的值为$-2$, 负数取$0$ + +再比如第三张小图，为`G-G`的碱基匹配对，也有三个方向： + +- 上方的箭头为$0-2=-2$ +- 左方的箭头为$0-2=-2$ +- 左上角对角线的箭头为$0+3=3$ + +取最大值为$3$ + + +![step](https://upload.wikimedia.org/wikipedia/commons/2/2c/Smith-Waterman-Algorithm-Example-Step1.png) + +:::tip + +横向和竖向的箭头分数都是空位罚分，对角线的箭头分数是置换矩阵所得 + +::: + +最终以此理，我们可以得到如下的得分矩阵： + +![matrix](https://upload.wikimedia.org/wikipedia/commons/2/28/Smith-Waterman-Algorithm-Example-Step2.png) + +其中红色箭头代表分数的来源用于回溯，蓝色方块表示得分矩阵的最大值。 + + +### 2. 回溯寻找最优比对 + +通过动态规划的方法，从得分矩阵的最大分值的元素开始回溯直至分数为 0 的元素。 + +具有局部最高相似性的片段在此过程中产生，如下图所示： + +![huishu](https://upload.wikimedia.org/wikipedia/commons/e/e6/Smith-Waterman-Algorithm-Example-Step3.png) + +最终得到比对结果为： + +``` +G T T - A C +| | | | | +G T T G A C +``` + +:::note + +具有第二高相似性的片段可以通过从最高相似性回溯过程之外的最高分位置开始回溯，即完成首次回溯之后，从首次回溯区域之外的最高分元素开始回溯，以得到第二个局部相似片段。 + +::: + + +## 代码实现 + +用 rust 写一个该算法的得分矩阵计算的基本实现： + +```rust + +// 1. define the score of match, mismatch and gap +const MATCH_SCORE: i64 = 3; +const MISMATCH_SCORE: i64 = -3; +const GAP_SCORE: i64 = -2; + +fn smith_waterman(seq1: &str, seq2: &str) { + + // 2. initialize the matrix + let mut matrix = vec![vec![0; seq2.len() + 1]; seq1.len() + 1]; + + // 3. calculate the matrix + for i in 1..=seq1.len() { + for j in 1..=seq2.len() { + // calculate the score from diagonal + let diag_score = if seq1.chars().nth(i - 1) == seq2.chars().nth(j - 1) { + MATCH_SCORE + } else { + MISMATCH_SCORE + }; + + // compare the three possible ways to get the current cell + matrix[i][j] = std::cmp::max( + matrix[i - 1][j - 1] + diag_score, // diagonal + std::cmp::max(matrix[i - 1][j] + GAP_SCORE, matrix[i][j - 1] + GAP_SCORE), // left and up + ); + + // if matrix_score < 0, set it to 0 + matrix[i][j] = std::cmp::max(0, matrix[i][j]); + } + } + // 4. display the matrix + display_matrix(&matrix); +} + +// pretty print the matrix +fn display_matrix(matrix: &Vec>) { + for row in matrix { + for cell in row { + print!("{:4} ", cell); + } + println!(); + } +} + +fn main() { + let seq1 = "GGTTGACTA"; + let seq2 = "TGTTACGG"; + + smith_waterman(seq1, seq2); +} + +``` + +可以得到如下结果： + +``` + 0 0 0 0 0 0 0 0 0 + 0 0 3 1 0 0 0 3 3 + 0 0 3 1 0 0 0 3 6 + 0 3 1 6 4 2 0 1 4 + 0 3 1 4 9 7 5 3 2 + 0 1 6 4 7 6 4 8 6 + 0 0 4 3 5 10 8 6 5 + 0 0 2 1 3 8 13 11 9 + 0 3 1 5 4 6 11 10 8 + 0 1 0 3 2 7 9 8 7 +``` + +和我们推理得到的一致。 + +## 参考 + + diff --git a/docs/Algorithm/_category_.json b/docs/Algorithm/_category_.json index dcdb9fb..fc2593e 100644 --- a/docs/Algorithm/_category_.json +++ b/docs/Algorithm/_category_.json @@ -5,5 +5,6 @@ "type": "generated-index", "description": "Algorithm!", "image": "/img/icon/algorithm.png" - } + }, + "collapsed": false } \ No newline at end of file diff --git a/docs/Algorithm/array_list/_category_.json b/docs/Algorithm/array_list/_category_.json new file mode 100644 index 0000000..1c4ba68 --- /dev/null +++ b/docs/Algorithm/array_list/_category_.json @@ -0,0 +1,9 @@ +{ + "label": "数组和链表", + "position": 2, + "link": { + "type": "generated-index", + "description": "基本数据结构：数组和链表" + }, + "collapsed": true +} \ No newline at end of file diff --git a/docs/Algorithm/basic/_category_.json b/docs/Algorithm/basic/_category_.json new file mode 100644 index 0000000..e02f5ba --- /dev/null +++ b/docs/Algorithm/basic/_category_.json @@ -0,0 +1,9 @@ +{ + "label": "算法基础知识", + "position": 1, + "link": { + "type": "generated-index", + "description": "算法基础知识，概念，评估，复杂度" + }, + "collapsed": true +} \ No newline at end of file diff --git a/docs/Algorithm/basic/efficiency_assess.mdx b/docs/Algorithm/basic/efficiency_assess.mdx new file mode 100644 index 0000000..e378379 --- /dev/null +++ b/docs/Algorithm/basic/efficiency_assess.mdx @@ -0,0 +1,41 @@ +--- +title: 算法效率评估 +description: 算法效率评估 +tags: [algorithm,] +sidebar_position: 1 +last_update: + date: 04/16/2024 + author: WenjieWei +--- + +## 概念 + +在正确解决问题的前提下，希望找到**高效**的算法。 + +包括以下两个维度： + +- 时间效率：算法执行所需的时间 +- 空间效率：算法占用内存的大小 + +简而言之，又快又省是目标。 + +## 评估方法 1: 实际测算 + +现有算法`A`和算法`B`，两者都能正确解决问题。找一台计算机，分别执行，记录时间和内存占用情况，这样可以真实反映算法的效率。但是存在两个较大的局限性： + +- 测算环境的干扰：不同的硬件配置\编译器\操作系统等，会对测算结果产生影响。 +- 数据规模不完整：不同的数据量表现不同，若测试不同量的数据，耗时耗力。 + +综上，因为这些局限性，经常使用理论估算来评估算法效率。 + +## 评估方法 2: 理论估算 + +也被称为复杂度分析，complexity analysis. + +定义：**它描述了随着输入数据大小的增加，算法执行所需时间和空间的增长趋势** + +重点： + +1. 划分为时间复杂度和空间复杂度 +2. 反映了效率与输入数据量的关系 +3. 关注时间空间增长的快慢 \ No newline at end of file diff --git a/docs/Algorithm/basic/space_complex.mdx b/docs/Algorithm/basic/space_complex.mdx new file mode 100644 index 0000000..69982fa --- /dev/null +++ b/docs/Algorithm/basic/space_complex.mdx @@ -0,0 +1,64 @@ +--- +title: 空间复杂度 +description: 算法占用内存空间的增长趋势 +tags: [algorithm,] +sidebar_position: 3 +last_update: + date: 04/26/2024 + author: WenjieWei +--- + +## 概念 + +空间复杂度（space complexity）用于衡量算法占用内存空间随着数据量变大时的增长趋势。这个概念与[时间复杂度](./time_complex)非常类似，只需将“运行时间”替换为“占用内存空间”。 + +## 使用的内存空间 + +算法运行过程中使用的内存空间包括： +- 输入空间：存储输入数据 +- 暂存空间：存储变量\对象\函数调用等 +- 输出空间：存储输出数据 + +其中暂存空间进一步划分为： + +- 暂存数据：保存常量\变量\对象 +- 栈帧空间：保存调用函数的上下文数据，系统在每次调用函数时都会在栈顶部创建一个栈帧，函数返回后，栈帧空间会被释放 +- 指令空间：用于保存编译后的程序指令，在实际统计中通常忽略不计 + +示例： + +```rust +use std::rc::Rc; +use std::cell::RefCell; + +/* 结构体 */ +struct Node { + val: i32, + next: Option>>, +} + +/* 创建 Node 结构体 */ +impl Node { + fn new(val: i32) -> Self { + Self { val: val, next: None } + } +} + +/* 函数 */ +fn function() -> i32 { + // 执行某些操作... + return 0; +} + +fn algorithm(n: i32) -> i32 { // 输入数据 + const a: i32 = 0; // 暂存数据（常量） + let mut b = 0; // 暂存数据（变量） + let node = Node::new(0); // 暂存数据（对象） + let c = function(); // 栈帧空间（调用函数） + return a + b + c; // 输出数据 +} +``` + +## 推算方法 + +与时间复杂度类似 \ No newline at end of file diff --git a/docs/Algorithm/basic/time_complex.mdx b/docs/Algorithm/basic/time_complex.mdx new file mode 100644 index 0000000..a38abe2 --- /dev/null +++ b/docs/Algorithm/basic/time_complex.mdx @@ -0,0 +1,211 @@ +--- +title: 时间复杂度 +description: 算法运行时间的度量 +tags: [algorithm,] +sidebar_position: 2 +last_update: + date: 04/16/2024 + author: WenjieWei +--- + +## 概念 + +时间复杂度分析不是统计算法的运行时间，而是运行时间随着数据规模增大的增长趋势。 + +举例： + +```rust +// 算法 A 的时间复杂度：常数阶 +fn algorithm_A(n: i32) { + println!("{}", 0); +} +// 算法 B 的时间复杂度：线性阶 +fn algorithm_B(n: i32) { + for _ in 0..n { + println!("{}", 0); + } +} +// 算法 C 的时间复杂度：常数阶 +fn algorithm_C(n: i32) { + for _ in 0..1000000 { + println!("{}", 0); + } +} +``` +- 算法 A: 只有打印操作，不管输入的 n 是多少，都只执行一次打印操作，所以时间复杂度是常数阶 +- 算法 B: 打印操作的次数和 n 成正比，所以时间复杂度是线性阶 +- 算法 C: 无论输入的 n 是多少，打印操作都是固定的 1000000 次，所以时间复杂度仍是常数阶 + +## 函数的渐进上界 + +给定一个函数，其接受参数的输入大小为 n: + +```rust +fn algo(n: i32) { + let mut a = 1; // +1 + a = a + 1; // +1 + a = a * 2; // +1 + + // 循环 n 次 + for _ in 0..n { // +1（每轮都执行 i ++） + println!("{}", 0); // +1 + } +} +``` + +该函数记为 $T(n) = 3 + 2n$，其中 $3$ 为常数项，$2n$ 为线性项。 + +我们将线性阶的时间复杂度记为$O(n)$，这个数学符号称为大$O$记号，表示函数的渐近上界（asymptotic upper bound）。 + +:::tip 渐进上界 + +计算渐近上界就是寻找一个函数 $f(n)$ ，使得当$n$趋向于无穷大时，$T(n)$和$f(n)$ 处于相同的增长级别，仅相差一个常数项的倍数。 + +::: + +## 推算方法 + +定义太过拗口，直接看推算方法。 + +### 1. 统计操作数量 + +针对代码，从上而下： + +1. 忽略常数项，因为对时间复杂度不产生影响 +2. 省略系数，例如循环 $2n$次和 $5n+1$次，都简化记为$n$次 +3. 循环的嵌套用乘法 + + +例如： + +```rust +fn algo(n: i32) { + let mut a = 1; // +0（技巧 1） + a = a + n; // +0（技巧 1） + + // +n（技巧 2） + for i in 0..(5 * n + 1) { + println!("{}", 0); + } + + // +n*n（技巧 3） + for i in 0..(2 * n) { + for j in 0..(n + 1) { + println!("{}", 0); + } + } +} +``` + +得到： + +$T(n) = n^2+n$ + +### 2. 判断渐进上限 + +**时间复杂度是由最高阶的项决定的**, 其他项可以忽略。 + +一些例子： + +| 函数 | 时间复杂度 | +| --- | --- | +| $T(n) = 1000$| $O(1)$ | +| $T(n) = 3 + 2n$ | $O(n)$ | +| $T(n) = 2n^2 + 3n + 1$ | $O(n^2)$ | +| $T(n) = n^3 + 8888n^2 + 1$ | $O(n^3)$ | + + +## 常见类型 + +常见的时间复杂度类型可以按如下排序： + +$O(1) < O(log\ n) < O(n) < O(n\ log\ n) < O(n^2) < O(2^n) < O(n!)$ + +常数阶 < 对数阶 < 线性阶 < 线性对数阶 < 平方阶 < 指数阶 < 阶乘阶 + +![img](https://www.hello-algo.com/chapter_computational_complexity/time_complexity.assets/time_complexity_common_types.png) + +### 1. 常数阶 + +```rust +/* 常数阶 */ +fn constant(n: i32) -> i32 { + _ = n; + let mut count = 0; + let size = 100_000; + for _ in 0..size { + count += 1; + } + count +} +``` + +### 2. 线性阶 + +```rust +/* 线性阶 */ +fn linear(n: i32) -> i32 { + let mut count = 0; + for _ in 0..n { + count += 1; + } + count +} +``` + +### 3. 平方阶 + +```rust +/* 平方阶 */ +fn quadratic(n: i32) -> i32 { + let mut count = 0; + // 循环次数与数据大小 n 成平方关系 + for _ in 0..n { + for _ in 0..n { + count += 1; + } + } + count +} +``` + +### 4. 指数阶 + +有点像细胞分裂的例子，每次分裂都会产生两倍的细胞。 + +```rust +/* 指数阶（循环实现） */ +fn exponential(n: i32) -> i32 { + let mut count = 0; + let mut base = 1; + // 细胞每轮一分为二，形成数列 1, 2, 4, 8, ..., 2^(n-1) + for _ in 0..n { + for _ in 0..base { + count += 1 + } + base *= 2; + } + // count = 1 + 2 + 4 + 8 + .. + 2^(n-1) = 2^n - 1 + count +} +``` + +计算操作的总数为： + +$2^0 + 2^1 + 2^2 +... + 2^{n-1} = 2^n - 1$ + +### 5. 对数阶 + +与指数阶相反，对数阶反映了“每轮缩减到一半”的情况。设输入数据大小为 $n$ ，由于每轮缩减到一半，因此循环次数是 $log_2(n)$. + +```rust +/* 对数阶（循环实现） */ +fn logarithmic(mut n: i32) -> i32 { + let mut count = 0; + while n > 1 { + n = n / 2; + count += 1; + } + count +} +``` \ No newline at end of file diff --git a/docs/Algorithm/bwt.mdx b/docs/Algorithm/bwt.mdx index 58761c8..917d588 100644 --- a/docs/Algorithm/bwt.mdx +++ b/docs/Algorithm/bwt.mdx @@ -2,7 +2,7 @@ title: BWT-Algorithm description: Burrus-Wheeler Transform Algorithm tags: [algorithm, code, bioinfo] -sidebar_position: 1 +sidebar_position: 3 --- ## BWT-Algorithm @@ -147,7 +147,7 @@ def bwt_encode(s: str) -> dict: raise ValueError("The parameter s must not be empty.") rotaions = all_rotations(s) - rotaions.sort() ## 排序, ASICC 码值从小到大 + rotaions.sort() ## 排序，ASICC 码值从小到大 bwt_composed_str = "".join([row[-1] for row in rotaions]) ## 取出每个元素的最后一个字符 original_index = rotaions.index(s) ## 找到原始字符串的位置 diff --git a/docs/Algorithm/graph-1.mdx b/docs/Algorithm/graph-1.mdx new file mode 100644 index 0000000..e69de29