# 序列演算法
## 1. Basic Algorithmic Strategies
- Huffman Codes
- Input：字元 with 頻率
- Output：把字元編碼，使得長度越小越好
- 方法：依照頻率最小的兩個 node 依序合成後建一棵樹，然後 left edge 編碼成 0，right edge 編碼成 1
- Divide and Conquer
- $T(n)=2T(n/2)+O(logn) \rightarrow O(n)$
- $T(n)=2T(n/2)+O(n) \rightarrow O(nlogn)$
- $T(n)=2T(n/2)+O(nlogn) \rightarrow O(nlog^2n)$
- Dynamic Programming
- Principle of optimality
- Overlapping subproblems
- **LCS**：目前幾乎都是 $O(mn)$ 頂多 $O(mn/logmn)$ 的方法
- **LIS**：$O(nlogn)$ 演算法
- 
- [**LCIS 2005 Yang, Huang and Chao**](https://reader.elsevier.com/reader/sd/pii/S0020019004003242?token=215D879C620B6BE779392F0891D1A9EBF8F619FA951C0A55ACA1EF25A2FE7C90FCA75664A7041E1BAB5317F60E20B091)
## 2. Alignment
- 問題： 給定兩字串 $A$、$B$，找出對齊的方式使得 cost 最小。
- 定義 match / mismatch / insertion gap / deletion gap 的 cost
- 使用 dash "-" 代表 gap ( insertion / deletion )
- Dot Matrix： 把兩字串畫成表格，相同部分標點
- Global Alignment :
- 找出對齊方式，使得兩字串全部對齊的 score 最大
- $S_{i,j}=max(~S_{i-1,j}+gap~,~S_{i,j-1}+gap~,~S_{i-1,j-1}+matchscore~)$
- 同樣的問題可以用來解決 LCS / Edit Distance
- LCS : 1/0/0/0
- ED : 0/-1/-1/-1
- Local Alignment :
- 從兩個字串個找出某一小片段，使得對齊的 score 是全部裡面最大 → 找出最相似的部分
- 定義 $S_{i,j}$ 為 $a_1a_2...a_i$ 和 $b_1b_2...b_j$ ，恰好結束在 $(i,j)$ 的 optimal local alignment 分數
- $S_{i,j}=max(~S_{i-1,j}+gap~,~S_{i,j-1}+gap~,~S_{i-1,j-1}+matchscore~,~0~)$
- Maximum Sum Query :
- Query Segment Sum in O(1) : sum(i,j) = prefix(j) - prefix(i-1)
- Method 1 : 紀錄結束在位置 i 的 Maximum Sum
- dp[i] = a[i] + max(dp[i-1],0)
- Method 2 : 先計算 prefix 及 prefixminbefore
- maxsum[i] = prefix[i] - prefixminbefore[i];
- Scoring Scheme :
- Affine Gap Penalties
- 長度 $k$ 的 gap，penalty 設定為 $x+k*y$
- $x$ opening penalty，$y$ gap penalty
- $D(i, j)$ : 最佳解 ending with deletion
- $I(i, j)$ : 最佳解 ending with insertion
- $S(i, j)$ : 最佳解 ending
- $D(i,j)=max(D(i-1,j)-y,~S(i-1,j)-x-y)$
- $I(i,j)=max(I(i,j-1)-y,~S(i,j-1)-x-y)$
- $S(i,j)=max(D(i,j),~I(i,j),~S(i-1,j-1)+w(a_i,b_j))$
- Constant Gap Penalties
- 任何長度的 gap，penalty 皆為 $x$，所以等價於 affine 把 $y$ 設成 $0$
- Restricted Affine Gap Penalties
- 長度 $k$ 的 gap，penalty 設定為 $x+f(k)*y$
- 其中 $f(k)$ 最多就只有 $c$ 而已， $f(k)=min(k,c)$
- $D'(i,j)=max(D'(i-1,j),S(i-1,j)-x-cy)$
- $I'(i,j)=max(I'(i,j-1),S(i,j-1)-x-cy)$
## 3. Space Saving Alignments
- Constrained Sequence Alignment
- 在生物序列的應用中，我們可以限制 insertion 和 deletion 的次數不要差太多，如給定 $u\le v$，$a_i$ 和 $b_j$ 要滿足 $u\le i-j\le v$ 才可以配對。
- 時間複雜度從 $O(nm)$ 降為 $O(nw)$
- 定義：
- $S^-(i,j)$ 為 $a_1a_2...a_i$ 和 $b_1b_2...b_j$ 的最好 align cost
- $S^+(i,j)$ 為 $a_{i+1}a_{i+2}...a_m$ 和 $b_{j+1}b_{j+2}...b_n$ 的最好 align cost
- 性質：找一個點 maximize $S^-(m/2,j)+S^+(m/2,j)$
- Hirschberg 演算法

    Linear_Space_Align(A,B,i1,j1,i2,j2){
        forward_SCORE(A,B,i1,j1,imid,j2)    //find the best score from (i_1,j_1)
        backward_SCORE(A,B,imid,j1,i2,j2)   //find the best score to (i_2,j_2)
        jmid = max index 使得 S_-[j]+s_+[j] 最大
        Linear_Space_Align(A,B,i1,j1,imid,jmid)
        Linear_Space_Align(A,B,imid,jmid,i2,j2)
    }

- 空間複雜度：$O(n)$
- 時間複雜度：$O(mn(1+\frac{1}{2}+\frac{1}{4}...))=O(mn)$
- Hirschberg 演算法的應用(線性空間)
- Affine Gap Penalty ： $x+k*y$
- 用線性空間來計算 $D,~I,~S$
- 用 $S$ 來切割問題
- Local alignment space saving：
- 用線性空間先找到後端點，再找到前端點
- 對於該問題用 Hirschberg 演算法
- Constrained Sequence Alignment
- naive 的方法，時間複雜度為 $O(nwlogn)$
- 需要好的 subdivide 問題的作法
- 使用中間對角線來分割，根據平行四邊形，size 必定減少一半
- 使用矩形來分割
## 4. Suboptimal Alignments ($\Delta$ Point)
- 目標：找出 $S-$ 和 $S+$ 加總的矩陣中，數值 $\ge \Delta$ 的欄位
- 滿足該性質的點稱為 **$\Delta$ Point**
- 方法一：$O(MN)$ 時間, $O(MN)$ 空間
- 直接用矩陣去存 $S-$ 和 $S+$
- 方法二：$O(M^2N)$ 時間, $O(N)$ 空間
- 用 Hirschberg 演算法
- 每一個 row，用 $O(MN)$ 時間和 $O(N)$ 空間去找 $\Delta$ Point
- 方法三：
- $S-$ 只能從上往下算，$S+$ 只能從下往上算，不能逆推
- 因此若用 Hirschberg 演算法，相鄰的 row 不可以直接轉移
- 方法四：$O(MNlogM)$ 時間, $O(NlogM)$ 空間
- 
- 有 $logM$ 層，每次只有一個 subproblem 是 active
- 每一層總時間是 $O(MN)$
- 每一層 working 空間是 $O(N)$
- 方法五：$O(MN~log\min{\{M,N\}})$ 時間, $O(M+N)$ 空間
- 
- 有 $log\min{\{M,N\}}$ 層，每次只有一個 subproblem 是 active
- 每一層總時間是 $O(MN)$
- 每一層 working 空間是 $4(M+N)\times \frac{1}{2^i}$
- 方法六：$O(MN~loglog\min{\{M,N\}})$ 時間, $O(M+N)$ 空間
- 
- 盡量切成更小的問題以降低層數，但要保證每一層 working space 成比例降低
- 有 $loglog\min{\{M,N\}}$ 層，每次只有一個 subproblem 是 active
- 每一層總時間是 $O(MN)$
- 每一層 working 空間是 $4(M+N)\times \frac{1}{2^i}$
- 方法七：$O(\frac{1}{\epsilon}MN)$ 時間, $O(\frac{1}{\epsilon}M^\epsilon N)$ 空間
- 
- 有 $1/\epsilon$ 層，每次只有一個 subproblem 是 active
- 每一層總時間是 $O(MN)$
- 每一層總空間是 $O(M^\epsilon N)$
- 方法八：$O(\frac{1}{\epsilon}MN)$ 時間, $O(\frac{1}{\epsilon}M^{1+\epsilon}+N)$ 空間
- 
- 一開始花 $O(MN)$ 時間，$O(N)$ 空間準備好 $S+$ 和 $S-$
- 有 $N/M$ 個 blocks，每一個 block 中：
- $O(\frac{1}{\epsilon}M^2)$ 時間
- $O(\frac{1}{\epsilon}M^{1+\epsilon})$ 空間
## 5.Multiple Alignments
- 目標：給定多個 sequence，找出對齊的方式使得 cost 最小
- Sum of Pairs (SP)：兩兩 sequence 的 cost 加總
- Tree Score：也有人用但比較複雜所以不提
- 解決的方法：
- $k$ 個長度 $n$ 的 sequence，時間 $O(n^k)$
- 因此是 NP-hard 問題
- 近似解：$2-\frac{l}{k}$ 對任意的 $l$。 ( $l$ 越大，近似解越好，但要花更多時間得到)
- Heuristic Approach：Progressive Alignment 沒有理論保證，但通常解尚可接受
- "Once a gap, always a gap" ， 時間複雜度 $O(k^2n^2)$
- Guide Tree
- 1.用 species 的相似度
- 2.先用 $O(k^2)\times O(n^2)$ 算兩兩的相似度
- Aligning Alignments
- 
- 在 affine gap 的應用中，不知道 gap 種類時，通常會直接把 gap 當作 open
---
---
---
## 6.Fast Alignments
- 目標：使用 $o(mn)$ 時間與空間，尋找 Exact Word Matches
- 給定兩個序列 $A、B$，事先決定好 Exact Word 長度 $k$
- 是否有長度 $k$ 的 common substring
- 方法一：Hash Tables
- 將 $A$ 每長度 $k$ 的 substring 用 hash table 紀錄
- 將 $B$ 每長度 $k$ 的 substring 去同一個 hash table 查詢
- 方法二：Suffix Tree
- 將 $A$ 的每個 suffix 用 suffix Tree 紀錄
- 將 $B$ 從左到右 (每個 suffix) 去 suffix tree 查詢
- 看是否可以有長度為 $k$ 的 prefix
- 方法三：Suffix Array
- 將 $A$ 轉換成 suffix array
- 將 $B$ 從左到右 (每個 suffix) 去 suffix array 查詢
- 用 binary search 尋找是否有 entry
- 看看這兩個 string 的 common prefix 是否為 $k$
- 方法四：BWT
- Burrows-Wheeler Transform
- 給定 T=abaaba\$，定義 BWT(T)=abba\$aa
- 給定 $T$ ，我們可以構造 $BWT(T)$
- Naive 構造：$O(n^2)$
- 
- **Suffix Array 構造：$O(n)$**
- 
- 給定 $BWT(T)$ ，我們可以回推 $T$
- Bind and Sort：$O(n^2)$
- 
- **LF mapping reverse：$O(n)$**
- 
- T ranking 性質：(第一 row 的字元，從左至右編號)
- 最左邊的 column 和最右邊的 column，同一字元的出現順序保持一致
- B ranking 性質：(最後一 column 的字元，從上至下編號)
- 根據此編號及 T ranking 性質可以設計 LF mapping 方法來 reverse
- $FM$ index backward search
- 功用：給定 T=abaaba\$，BWT(T)=abba\$aa，我們想要 query $q="aba"$ 是否為 $T$ 的 substring
- 前置處理：計算 $C、OCC$
- 
- Query 階段：
- 
- 
- FASTA
- 1.挑出最相近 region ( 用 exact word match )
- 2.用 PAM Matrix rescore ( 不同的字元 match 時有不同的權重 )
- 3.刪除掉不相近的 matches
- 4.在 band 當中做 alignment
- BLAST
- 目標在找 MSP ( Maximum segment pair)
- Match : 5，Mismatch : -4
- 對每個 diagonal 去找 maximum sum segment
- PAM Matrix
- 
- 1.對 Sequence A 建立 Hash Table
- 2.用 Sequence B 去尋找 hit 的地方
- 3.延伸 hit
- PatternHunter
- Spaced Seed
- 
- $l = 18$
- $w = 11$
- 比起 blastn，不要求連續的 match
- 若兩 sequence 有一段 region 相似度 > 61 % (過多的 1)，則 seed 1 1 \* 1 必定 hit
- 證明：如果 seed 1 1 \* 1 沒有 hit ，則 similarity <= 61 %
- 如果 seed 1 1 \* 1 沒有 hit ，每 100 個位置當中，最多只有 61 個 matches
- 開頭有一些 $0$
- $10^b,~~b\ge 1$
- $110^b,~~b\ge 2$
- $1110^b,~~b\ge 2$
- 結尾最多 $\le 3$ 個 $1$
## 7.Genome Reconstruction
- Newspaper Problem
- 
- 透過許多片面的資訊，重構出完整的內容
- 一些圖論問題：
- Eulerian Cycle Problem：走過每個 "邊" 的 cycle 是否存在
- Balanced Graphs ： indegree(v) = outdegree(v) 對於每一個點 v.
- 任意 Connected 的 Balanced Graph 都有 Eulerian Cycle
- Hamiltonian Cycle Problem：走過每個 "點" 的 cycle 是否存在
- NP-hard 問題
- Fragment Assembly
- 假設：
- 每一個 k-mer 都有出現
- 每一個 k-mer 都只出現一次
- 沒有任何 error
- 是由一個 circular 的基因序列產生
- 舉例：我們有以下的 read，我們想要重構出 circular 的基因序列
- GTG GCG GCA ATG TGG TGC GGC CGT CAA AAT
- 方法一：
- 把每個 read 建立一個節點，建邊
- 
- 在此圖上找一個 Hamiltonian cycle
- 方法二：
- 建立節點為 k-1-mer，把每個 read 建立一個邊
- 
- 在此圖上找一個 Eulerian cycle
- 放寬條件假設
- Generating (nearly) all k-mers
- Handling Errors in Reads
- Handling Repeated k-mers
- From Circular to Linear Genomes
- [De Brujin Sequence](https://en.wikipedia.org/wiki/De_Bruijn_sequence) B(n,k)
- B(3,2) -> 0 0 0 1 1 1 0 1 、 1 1 1 0 0 0 1 0
- 找出 super string，包含所有長度 n、k 種 character 的 string (共有 $k^n$ 種)
- B(n,k) 長度為 $k^n$
- B(n,k) 的種類數為 $\frac{(k!)^{k^{(n-1)}}}{k^n}$
- 其他有趣問題
- Counting number of Eularien Curcuits
- http://www.cdam.lse.ac.uk/Reports/Files/cdam-2004-12.pdf
- Shortest Superstring Problem
- http://fileadmin.cs.lth.se/cs/Personal/Andrzej_Lingas/superstring.pdf
## 8.Pattern Identification in a Haplotype Block
- 名詞定義：
- Single Nucleotide Polymorphism (**SNP**) : 單核苷酸多態性
- 序列中的某一位置，有 $\ge1\%$ 機率發生變異
- 該位置可以分成 **主等位基因、副等位基因**
- Mutation : 突變，是序列中的某一位置，有 $<1\%$ 機率發生變異
- Haplotype (pattern): 單倍型
- 一組基因序列，不同的物種、疾病會有不同的序列
- Ordered list of SNPs
- 一組基因序列可以由 recombination hotspots 切割成一些 Haplotype Blocks 區塊
- Tag SNPs : 演算法的目標，取出一部分的 SNP 來區分不同的 Haplotype patterns
- 問題定義：
- 給定一些 SNP，以及 Haplotype pattern，從 SNP 中找出最少的 Tag SNPs，使得我們可以區分任意的 haplotype sample。
- 這是一個 NP-hard 問題，又稱為 minimum test set problem
- e.g. 
- 我們可以取 Tag SNPs 為 $S_3、S_4$
- Haplotype pattern 為 $P_1\sim P_4$，每一個 pattern 用四個 SNP 描述
- 不管 haplotype sample 為何，只要看他的 $S_3、S_4$，就可以把它 match 到 $P_1\sim P_4$ 的某一個 Haplotype pattern
- 規約到別的 Np-Hard 問題
- Set Cover
- 
- 我們需要選擇適合的 tag SNPs，來區分不同的 pattern
- 相當於我們需要選適合的 set，來 cover 底下六個 element
- Greedy 演算法，可以保證被最佳解乘以 log factor 所 bound 住
- $U=\{1,2,3,4,5,6,7,8,9\}$
- $C_1=\{2,3,4,5,6,7\}$
- $C_2=\{1,2,3,4,5\}$
- $C_3=\{6,7,8,9\}$
- Greedy 解：
- $~~~~~~~~[2,3,4,5,6,7,8,9,1]$
- $S=[\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{2},\frac{1}{2},0]$
- 每一個分數 $S_i$ 不會超過最佳解的平均分數 $\frac{|C^*|}{i}$
- $|C|=\sum S_i\le\sum\frac{|C^*|}{i}=|C^*|(\frac{1}{1}+\frac{1}{2}+...\frac{1}{n})=|C^*|logn$
- Integer Programming
- 
- 為了區分 $P_1$ 和 $P_2$ $:~~\rightarrow x_3+x_4\ge 1$
- 為了區分 $P_2$ 和 $P_4$ $:~~\rightarrow x_1+x_2+x_3\ge 1$
- 可以使用 Iterative LP-relaxation Algorithm
- Missing Data
- 如果 haplotype sample 可能有些 SNP 是缺失的，仍然要能區分出來。
- 問題定義：假設 haplotype sample 最多有 $m$ 個 SNP 是缺失的，我們想要從 SNP 中找出最少的 Tag SNPs 出來，使得我們總是可以區分。
- Set Cover：
- 我們要 cover 每一個 element $m+1$ 次
- Integer Programming：
- 
## 9.Haplotype Inference
- Haplotype data 和 Genotype data 的比較
- 
- 問題定義
- 給定 $n$ 個 genotypes ，以及 $m$ 個可能的 haplotypes，求出最小的 haplotypes 集合，可以解釋所有給定的 genotypes。
- 這是一個 NP-hard 問題，又稱為 Haplotype Inference
- 轉換成整數規劃問題
- 符號定義：
- 每一個 haplotypes 對應到變數 $x_i$
- 每一個 genotypes 對應到集合 $S_j$，裏頭包含了許多 $(h_r,h_t)$ 配對
- 每一個配對可以解釋此 genotype
- 方法一：轉換成 IQP 問題，有近似演算法理論保證
- 選到為 1 不選為 -1。
- 
- 方法二：
- 選到為 1 不選為 0。
- 目標函數：$\min \sum\limits_{i=1}^mx_i$
- 限制函數：$\sum\limits_{(h_r,h_t)\in S_j}x_rx_t\ge 1$
- 方法三：
- 選到為 1 不選為 0。
- 目標函數：$\min \sum\limits_{i=1}^mx_i$
- 限制函數：$\sum\limits_{(h_r,h_t)\in S_j}\lfloor\frac{x_r+x_t}{2}\rfloor\ge 1$
- 方法四：
- $x_i\in\{0,1\},~~w_{jrt}\in\{0,1\}$
- 目標函數：$\min \sum\limits_{i=1}^mx_i$
- 限制函數：$\sum\limits_{(h_r,h_t)\in S_j}w_{jrt}=1$ 而且 $x_r\ge w_{jrt},~~x_t\ge w_{jrt}$
## 10.Heaviest Segments
- Maximum-sum segment O(n) time
- 給定一整數序列，找出子序列使得總和最大
- 方法一：停在每個位置的總和最大子序列，只要考慮是否接進上一位置
- 因此就可以知道結尾在每個位置的總和最大子序列
- 方法二：紀錄每個位置的 prefix sum 以及 prefix min
- Maximum-average segment O(n) time
- 給定一整數序列，找出子序列使得平均最大
- 方法：直接找最大值即可
- Maximum-average segment 結尾在每個位置 O(n) time
- 給定一整數序列，找出起點在每個位置的平均最大子序列
- 方法：Decreasingly right-skew decomposition
- 目標：把整數序列 partition 成遞減的 right-skew 序列
- 1.Right skew 表示 prefix 的平均小於 suffix 的平均
- 2.遞減表示各個 right-skew 的平均越來越小
- 步驟：由後往前建立 pointer 及 average
- 
- 也可以找 Increasingly left-skew decomposition
## 11.RMSQ
- **Range Maxima Query (RMQ)**
- 給定一整數序列，每次詢問時，給定一段範圍，找出此範圍的最大數字
- https://www.youtube.com/watch?v=ZBHKZF5w4YU
- https://codertw.com/%E7%A8%8B%E5%BC%8F%E8%AA%9E%E8%A8%80/561565/
- https://zh.wikipedia.org/wiki/%E8%8C%83%E5%9B%B4%E6%9C%80%E5%80%BC%E6%9F%A5%E8%AF%A2
- 演算法：

        int input[n];
        int dp[i][j];
        void preprocess(){
            for(int i=1;i<=n;i=2*i)
                for(int pos=1;(pos+i-1)<=n;pos++)
                    dp[pos][i]=min(dp[pos][i],dp[pos+i/2][i]);
        }
        int query(int a, int b){
            int k = log(b-a);
            return min(dp[a][k],dp[b-k+1][k]);
        }

- **Range Maximum-Sum Segment Query (RMSQ) 版本一**
- 給定一整數序列，每次詢問時，給定一段範圍 $[k,j]$，找出此範圍的總和最大子序列
- 1.事先算出 prefix sum，然後每個點 $i$ 紀錄兩個位置
- 1.$Left\_bound(i)$ 指向前面大於等於自己中最近的位置
- 2.$partner(i)$ 指向 $[Left\_bound(i),i]$ 中的最低位置
- 2.先用 RMQ 找出最大的 $(partner(i), i)$ 配對
- 1.沒有超出範圍，即為答案
- 2.有超出範圍，更新 $Left\_bound(i)$
- 3.答案只會是 (new_partner(i), i) 或 (partner(j), j) 其中之一
- **Range Maximum-Sum Segment Query (RMSQ) 版本二**
- 給定一整數序列，每次詢問時，給定兩段範圍，找出起點和終點分別在這兩段範圍得總和最大子序列
- 方法：首先要建立 prefix sum array
- 1.如果兩區間不重疊：直接使用 RMQ 問題的解法，找前半區間的 maxima 及後半區間的 minima
- 2.如果兩區間重疊，可以區分三種 case
- 2a.最佳解用到前半部不重疊部分：和 1. 相同
- 2b.最佳解用到後半部不重疊部分：和 1. 相同
- 2c.最佳解都在重疊部分：和 RMSQ 版本一相同
## Suffix Array Linear Time Construction
- [Suffix Array](https://www.cs.helsinki.fi/u/tpkarkka/opetus/10s/spa/lecture11.pdf)