# BIGSI
## 前言簡介、BIGSI的功用

筆者前言：本篇BIGSI(中國翻譯：超高速搜索現存細菌和病毒基因組)，是在2019年才被探索出來並發表在nature、ICCC等期刊上發表，小弟翻了多篇研究，還有與以前的做比對比較，未縮減前全文超過分析兩萬五千字，後續發現延伸太多出去，偏離主題內容，故進行刪減。

用你我一般人非這方面領域學者、專家能理解的話來解釋：就是現在全球細菌以及病毒基因組呈指數成長，而這個BIGSI結合了基因組學以及網路搜索的計算方法來檢索基因，進一步我們可以檢索、儲存新病毒的基因組更快速。
知道基因組之後，研究人員可以更好的去研究基因組的功能、與疾病的關係、演化的過程。
以下我都會用方便理解的白話文進行解說。

##### 5/6 與教授討論過後，希望著重在資料結構以及算法的部分，生物基因應用則不再更新。

如有興趣了解，歡迎來信討論。

Bitsliced Genomic Signature Index(BIGSI)是個可被搜索的資料結構，並源於一個對於細菌以及病毒的基因組序列進行索引BIGSI的搜索功能，其中利用到了Bloom filter的概念
可以快速找到抗性基因並確定宿主範圍，量化存檔資料集中的抗生素耐藥性，將新的序列資料儲存而逐漸成長，並且可以成長到數百萬個資料集。

以下是針對BIGSI的資料結構部分進行的整理

| K-mer | 序列中長度為K的連續子序列，如在DNA中的ATCG，RNA中的AUCG。在BIGSI中，每個基因組都被分解為 k-mer（長度為 k 的 DNA 序列，把它當字串就好）。這種資料結構使得基因組的索引和搜索變得更有效。 |
| ------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
| Bloom Filter | 在BIGSI 中，每個基因組的 k-mer 都被hash過後，變成一組數字，存儲在一個 Bloom Filter 中。 |
| BIGSI Matrix | 所有基因組的 Bloom Filter 被合併成一個巨大的矩陣，稱為 BIGSI 矩陣。每一行對應一個 k-mer，每一列對應一個基因組。如果一個 k-mer 存在於一個基因組中，則該矩陣的對應位置為 1，否則為 0。 |
| Bit Vector | BIGSI 矩陣的每一行都可被視為一個 bit vector，代表了 k-mer 在所有基因組中的存在情況。這種資料結構可以使查詢特定 k-mer 變得非常快速。 |

## 實作方法
為了要實作BIGSI，實現可被搜索的資料結構，我們需要作以下的步驟來作進行，
### 1.比對序列：
比對同一物種的基因組，假設相似度高(相似度是否高須利用BWA以及bowtie，則需要在可接受的時間內比對%數達一定標準以上)就返回同一對齊和比對得分。

#### 註解：BWA以及bowtie是在mapping這個環節時進行比對基因組的程式，常見的有SAM、BAM、CRAM檔
此圖為示意進行基因比對序列的Bwa、bowtie2結果，可看到有多樣匹配數據。
### 2.BLAST的應用：
BLAST將一個查詢字符串與一個系統發育的參考基因組資料庫進行比較。BLAST從查詢中獲取k-mer，對於每個k-mer，在固定的編輯距離內創建一個k-mer的"neighborhood(鄰域)"並在參考基因組搜索，比對只能透過特定點擴展完成。
BLAST主要用於核甘酸跟蛋白質的搜索，並能找到近距離同源匹配。

#### 註解：BLAST是在生物領域常被使用的程式，用於比對DNA、RNA、蛋白質序列的排序，並可用於比較查詢序列與已知得序列的大型資料庫，尋找查詢序列與資料庫中的序列之間的相似區域或對齊。

#### 註解：核甘酸(對，高中生物學過的那個DNA、RNA)是遺傳訊息的主要因子，其中會有AGCT四種碱基，K-mer字串搜尋配對就是為了要找這些的組成。

#### 註解：k-mer是一種生物訊息學中的概念，代表的是序列中長度為K的連續子序列，如在DNA中的ATCG，RNA中的AUCG。以下是簡單的範例實驗k-mer，運用python實作

def generate_kmers(sequence, k):
    kmers = set()
    for i in range(len(sequence) - k + 1):
        kmer = sequence[i:i+k]
        kmers.add(kmer)
    return kmers

sequence = "ATGCGATCGA"
k = 3
kmers = generate_kmers(sequence, k)
print(kmers)

這個函數接受一個序列的和k的值，然後生成一個所有K-mer的集合，例如：對於序列"ATGCGATCGA"和k=3，生成的結果會是{"ATG", "TGC", "GCG", "CGA", "GAT", "ATC"}。

#### 註解：實作K-mer的基本過程可以分為以下步驟：

1.獲取序列
首先需要從基因組或轉錄組等序列數據庫中獲取序列。

2.確定K值
決定所需的K-mer值，K值越大，K-mer越具體，能夠提供更多的信息，但是也會增加算法的時間和空間複雜度。通常，K值在3到10之間。

3.分割序列
根據K值將序列分割成一系列K-mer，如果序列長度小於K值，則可能需要使用padding進行填充。

4.處理K-mer
在處理K-mer時，需要注意避免出現重複的K-mer，可以使用哈希表等資料結構進行去重。在某些情況下，需要計算K-mer在序列中的頻率，以便進一步分析序列特性。

5.分析K-mer
根據分割和處理後的K-mer，可以進行序列比對、基因組組裝、特徵提取等分析，從而獲得序列的各種信息。


### 3.MASH(MinHash)
在資料庫中存儲每個參考資料的微小指紋，通過對組裝序列集的查詢，將組裝序列的指紋與RefSeq的指紋進行比較，以找到最接近的參考序列(圖中使用RefSeq)。

#### 註解：MinHash是為一種算法，主要估計兩個集合之間的相似度，用於比較和分類大規模生物序列，而在使用MASH進行比較時，需要將每個參考資料的微小指紋M用於比較和分類大規模生物序列資料集。再運用MASH比對時，需要將每個參考資料的微小指紋

#### 註解：RefSeq可供公開核苷酸序列及其蛋白質產物，由美國國家生物訊息中心建立(NCBI)

### 4.配對 1 個低多樣性物種的序列
利用序列開花樹(Sequence Bloom Tree)來索引資料中的k-mers，壓縮索引來搜索原始未組裝的序列集（未組裝的序列集顯示為“堆(piles)”的序列（短線），所有這些序列的顏色都相同，表示相同的種類）。設計用於人類資訊，SBT可以用來尋找哪些RNA測序資料集包含指定的轉錄本。

#### 註解：SBT 是一種資料結構，專門用在快速、高效地搜索多個序列資料集，多用於大規模基因組資料。它是Bloom filter的概念擴展， Bloom filter是一種概率資料結構，可以快速確定元素是否在集合中，並可以接受較小的誤差率。

更進一步來說，SBT其實就是把Bloom filter組織成一個二叉樹。每個節點（包括內部節點和葉節點）都有一個相應的 Bloom Filter。對於內部節點，其 Bloom Filter 是它的兩個子節點的 Bloom Filter 的位元邏輯 OR。

from bloom_filter import BloomFilter
import math

# 定義 SequenceBloomTree 類別
class SequenceBloomTree:
    # 初始化 SequenceBloomTree
    def __init__(self, seq_list, kmer_size, bf_capacity, bf_error_rate):
        # 設定 k-mer 的大小
        self.kmer_size = kmer_size
        # Bloom Filter 的容量
        self.bf_capacity = bf_capacity
        # Bloom Filter 的錯誤率
        self.bf_error_rate = bf_error_rate
        # 建立樹的根節點
        self.root = self._build_tree(seq_list)

    # 建立樹
    def _build_tree(self, seq_list):
        # 如果只有一個序列，創建一個 Bloom Filter
        if len(seq_list) == 1:
            return self._create_bloom_filter(seq_list[0])

        # 如果有多個序列，則分割序列列表並遞迴地建立子樹
        middle = len(seq_list) // 2
        left_child = self._build_tree(seq_list[:middle])
        right_child = self._build_tree(seq_list[middle:])
        return Node(bf=None, left=left_child, right=right_child)

    # 創建 Bloom Filter
    def _create_bloom_filter(self, sequence):
        bloom_filter = BloomFilter(capacity=self.bf_capacity, error_rate=self.bf_error_rate)
        # 從序列中提取 k-mer
        kmers = [sequence[i:i + self.kmer_size] for i in range(len(sequence) - self.kmer_size + 1)]
        # 將每個 k-mer 添加到 Bloom Filter
        for kmer in kmers:
            bloom_filter.add(kmer)
        return Node(bf=bloom_filter)

    # 查詢某個序列是否在樹中
    def query(self, sequence):
        # 從序列中提取 k-mer
        kmers = [sequence[i:i + self.kmer_size] for i in range(len(sequence) - self.kmer_size + 1)]
        return self._query_helper(self.root, kmers)

    # 查詢輔助函數
    def _query_helper(self, node, kmers):
        # 如果節點為空，返回 False
        if node is None:
            return False

        # 如果節點有 Bloom Filter，檢查每個 k-mer 是否都在 Bloom Filter 中
        if node.bf is not None:
            return all(kmer in node.bf for kmer in kmers)

        # 如果節點是內部節點，則在左子樹和右子樹中查詢
        return self._query_helper(node.left, kmers) or self._query_helper(node.right, kmers)

# 定義節點類別
class Node:
    # 初始化節點
    def __init__(self, bf=None, left=None, right=None):
        # bf 是這個節點的 Bloom Filter
        self.bf = bf
        # left 是這個節點的左子節點
        self.left = left
        # right 是這個節點的右子節點
        self.right = right

以上程式碼為SBT的實作，以上本程式使用ChatGPT-4實作案例。

# 假設我們有四個序列
sequences = [
    "ATGCACTAGCATCGACGACTAGCTAGC",
    "TGCATGCAAGCTAGCTAGCTAGCTAGC",
    "AGCTAGCTAGCTAGCTATGCACTGCTA",
    "CGTACGTACGTAGCTAGCTAGCTAGCT"
]

創建一個SequenceBloomTree
kmer_size = 4
bf_capacity = 1000
bf_error_rate = 0.1
sbt = SequenceBloomTree(sequences, kmer_size, bf_capacity, bf_error_rate)

現在，我們可以使用SBT查詢序列是否可能存在於我們的數據集中
query_seq1 = "ATGCACTAGCATCGAC" # 存在於第一個序列
query_seq2 = "TGCATGCAAGCTAGCT" # 存在於第二個序列
query_seq3 = "ACGTACGTACGTA" # 不在任何序列中

print(sbt.query(query_seq1)) # 應該返回True
print(sbt.query(query_seq2)) # 應該返回True
print(sbt.query(query_seq3)) # 可能返回False（根據bloom fliter的誤報率）

依照教授要求，著重在資料結構以及算法上進行更多介紹
因此，將著重介紹Bloom fliter、k-mer的部分。

## Bloom fliter
由於BIGSI，是在SBT基於Bloom fliter的基礎之上做發展。因此，筆者得先介紹何為Bloom Fliter。

Bloom Fliter，Bloom Filter是一種以空間效率為主的資料結構，用於檢查一個元素是否在一個集合中。而這個資料結構允許一定程度上的誤判（即可能產生該元素表示存在但實際上不存在），但不會出現反例（表示不存在但實際上存在）。簡單來說，它可以幫助我們節省存儲空間，但會有一定程度上誤差值的產生。

也就是說，當一個元素被加入集合時，通過K個hash函數將這個元素映射成一個位數組中的K個點，把它們置為1。檢索時，我們只要看看這些點是不是都是1就（大約）知道集合中有沒有它了：如果這些點有任何一個0，則被檢元素一定不在；如果都是1，則被檢元素很可能在。這就是布隆過濾器的想法。

Bloom fliter是由一個長度為m的bit的位數組（bit array）與k個hash function組成的資料結構。位數組均初始化為0，所有hash function都可以分別把輸入數據盡量均勻地散列。

當要插入一個元素時，將其數據分別輸入k個hash function，產生k個hash值。以哈希值作為位數組中的下標，將所有k個對應的bit置為1。

當要查詢（即判斷是否存在）一個元素時，同樣將其數據輸入hash function，然後檢查對應的k個比特。如果有任意一個bit為0，表明該元素一定不在集合中。如果所有bit均為1，表明該集合有（較大的）可能性在集合中。為什麼不是一定在集合中呢？因為一個bit被置為1有可能會受到其他元素的影響( hash碰撞)，這就是所謂的“假陽性”（false positive）。相對地，“假陰性”（false negative）在布隆過濾器中是絕不會出現的。

總的來說，Bloom fliter是一種提供存儲和查詢效率的資料結構，特別適合用於大數據和分布式系統中，可以幫助節省存儲空間和提高查詢速度。但同時，它也存在一定的假陽性錯誤率，這是由於hash碰撞所導致的，也就是說，在Bloom fliter中，一個位置可能被多個元素映射到，所以當我們查詢某個元素是否存在時，如果該元素映射到的位置都是1，我們只能說該元素可能存在，不能確定一定存在。

####

時間複雜度：主要取決於有多少個hash function的數量，對於每個hash function，需要透過k個hash function的數量進行hash，並在陣列中檢查相應的位置，所以時間複雜度為O(k)，k為hash function的數量。

空間複雜度：Bloom Filter 的空間複雜度與陣列的大小有關，即 O(m)。在實際應用中，根據集合的元素數量 n 和期望的誤判率 p 來選擇合適的陣列大小 m 和hash function的數量 k。通常，根據最佳實踐，可以使用以下公式計算 m 和 k：

m = - (n * log(p)) / (log(2)^2)

k = (m / n) * log(2)

其中 log 表示自然對數。

### hash function

Hash Function是一種將任意大小的輸入數據映射到固定大小的輸出空間的函數。換句話說，

Hash Function接收一個輸入(key)，然後通過函式的轉換算法將其轉換成固定範圍內的一個數值被稱做hash value。應具有以下特性：
單向性、確定性、均勻分布

單向性：對於給定的輸入，計算出相應的Hash value；但對於給定的Hash value，推算出原始輸入應該困難(不可逆向)。

確定性：相同的輸入應該產生相同的Hash Value。

均勻分布：Hash Function應該盡可能地將輸入數據均勻分布在輸出空間中，以減少不同輸入產生相同Hash Value的情況

有許多常用的Hash Function，如MD5、SHA-1、SHA-256 以及 MurmurHash 等。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>

// 求模的素數，用於將散列函數的結果映射到array中
#define BIT_ARRAY_SIZE 1000

// 位數組，用來存儲元素的狀態
bool bit_array[BIT_ARRAY_SIZE] = {false};

// hash function 1
unsigned int hash1(const char *str) {
    unsigned int hash = 5381;
    int c;

    while ((c = *str++))
        hash = ((hash << 5) + hash) + c; // hash * 33 + c

    return hash % BIT_ARRAY_SIZE;
}

// hash function 2
unsigned int hash2(const char *str) {
    unsigned int hash = 0;
    int c;

    while ((c = *str++))
        hash = c + (hash << 6) + (hash << 16) - hash;

    return hash % BIT_ARRAY_SIZE;
}

// 將元素添加到Bloom Filter中
void add(const char *str) {
    bit_array[hash1(str)] = true;
    bit_array[hash2(str)] = true;
}

// 檢查元素是否可能在Bloom Filter中
bool contains(const char *str) {
    return bit_array[hash1(str)] && bit_array[hash2(str)];
}

int main() {
    // 向Bloom Filter中添加元素
    add("apple");
    add("banana");
    add("cherry");

    // 檢查元素是否在Bloom Filter中
    printf("apple is %spresent.\n", contains("apple") ? "" : "not ");
    printf("banana is %spresent.\n", contains("banana") ? "" : "not ");
    printf("cherry is %spresent.\n", contains("cherry") ? "" : "not ");
    printf("durian is %spresent.\n", contains("durian") ? "" : "not ");

    return 0;
}

這也是為什麼bloom fliter會出現偽陽性的狀況，偽陽性發生的原因是因為不同的元素可能被hash到同一個位置。
舉例來說，假設我們有兩個元素 A 和 B，它們分別被hash到bit array的第 1 和第 2 位。現在，如果我們檢查一個元素 C 是否在 Bloom Filter 中，並且 C 剛好被hash到第 1 和第 2 位，那麼 Bloom Filter 會誤判 C 是在集合中的，即產生偽陽性(表示存在於集合中，但實際上不存在)

也就是說，我們將k-mer的基因

## BIGSI 矩陣

BIGSI Matrix 是由多個 Bloom Filter 構成的。

在 BIGSI 中，每個基因組都被分解為 k-mer（長度為 k 的 DNA 序列），並且每個基因組的 k-mer 都被存儲在一個 Bloom Filter 中。然後，這些 Bloom Filter 被合併成一個大的矩陣。

在 BIGSI Matrix 中，每一行對應一個 k-mer，每一列對應一個基因組。如果一個 k-mer 存在於一個基因組中，則該矩陣的對應位置為 1，否則為 0。

因此，BIGSI Matrix 實際上是一個二維的bool array，其中的每一個元素都是 0 或 1

5/23補充
一、資料結構被破壞：使用二維矩陣的很大一個原因是因為行代表數據集，列代表其中的k-mer，這樣的結構可以快速查找我們所需要的資訊。如果我們將這個二維矩陣轉化為一維陣列，我們就失去了這種快速的查找能力。

在一維陣列中，我們喪失了二維矩陣提供的結構信息。在二維矩陣中，我們可以獨立地查看每一個Hash
Function（行）在每個數據集（Bloom
Filter，也就是列）中的結果。這種結構讓我們可以在對一個特定的k-mer進行查詢時，迅速地找到包含該k-mer的所有數據集。這是因為我們可以將Hash
Function應用到該k-mer上，然後查詢該Hash Function對應的行，看看哪些列（也就是哪些Bloom
Filter）在這一行中的值為1。這種查詢方式的時間複雜度為O(1)，因為行索引和行位向量都被存儲在一個Hash Table中存儲。

如果我們將這個二維矩陣轉化為一維陣列，我們就無法再像之前那樣快速地查找。在一維陣列中，每一個元素都代表一個Bloom
Filter。當我們對一個特定的k-mer進行查詢時，我們必須遍歷整個陣列，對每一個元素（也就是每一個Bloom
Filter）應用所有的Hash
Function，然後檢查是否所有的哈希函數都返回1。這種查詢方式的時間複雜度為O(n)，其中n是陣列的長度，即我們的數據集的數量。這顯然比二維矩陣的查詢方式效率要低得多。

二、查詢效率降低：我們在BIGSI中使用二維陣列時，時間複雜度是O(1)，因為行索引以及行位向量都被儲存在一個Hash
Table裡，但在一維陣列裡我們無法這樣做。
今天在一維陣列的話，每一個元素都代表一個Bloom Fliter。如果說我們想要查詢其中一個K-mer，必須得將每一個Hash
Function應用到K-mer上，然後再回來一維矩陣中找對應的元素
這樣的時間複雜度會是O(n)，n是一維陣列的長度。這樣的效率會比使用二維來說低很多，所以這方面才是一個創新以及革新。

三、添加數據的成本：如果我們想在二維陣列中添加入數據的話，我們只須將新的Bloom Fliter作為新的列加入到現有的Bit
Array就好，但在一維裡，我們需要調整整個陣列的結構，會導致大量的計算以及儲存資訊成本。

說到這裡，不曉得看正在看這篇文章的你至今有無明白這三者的關係。
再次整理一次，k-mer是用於整理ATCG並進行排序比較，k-mer實際上就是字串處理，然後在bloom fliter裡hash完之後儲存，最後這些bloom fliter在放進BIGSI Martrix經歷bool 得到一個二維大矩陣。

## Bit Vector
所以，綜上所述 BIGSI 矩陣的每一行都可以被視為一個 bit vector，表 k-mer 在所有基因組中的存在情況。這種資料結構使得查詢特定 k-mer 變得非常方便快速。

來源：ISCA2021 Sieve: Scalable In-situ DRAM-based AcceleratorDesigns for Massively Parallel k-mer Matching
https://ieeexplore.ieee.org/abstract/document/9499937

Kraken: ultrafast metagenomic sequence classification using exact alignments
https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-3-r46

Reducing storage requirements for biologicalsequence comparison
https://academic.oup.com/bioinformatics/article/20/18/3363/202143

Ultra-fast search of all deposited bacterial and viral genomic data
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6420049/

nature子刊
https://www.nature.com/articles/s41587-018-0010-1

Improved representation of sequence bloom trees(2020/3)
https://academic.oup.com/bioinformatics/article/36/3/721/5553093

https://blog.csdn.net/qq_39473674/article/details/124510538
Bloomfilter 解說
https://blog.csdn.net/jiaomeng/article/details/1495500