--- # System prepended metadata title: BIGSI --- # 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，實現可被搜索的資料結構，我們需要作以下的步驟來作進行， ![](https://i.imgur.com/5BVIckS.png) ### 1.比對序列： 比對同一物種的基因組，假設相似度高(相似度是否高須利用BWA以及bowtie，則需要在可接受的時間內比對%數達一定標準以上)就返回同一對齊和比對得分。 #### 註解：BWA以及bowtie是在mapping這個環節時進行比對基因組的程式，常見的有SAM、BAM、CRAM檔 ![](https://i.imgur.com/QnUwAV9.png) 此圖為示意進行基因比對序列的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，可以進行序列比對、基因組組裝、特徵提取等分析，從而獲得序列的各種信息。 ![](https://i.imgur.com/wo0XGRh.png) ![](https://i.imgur.com/uHy3FRG.png) ### 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。 ![](https://hackmd.io/_uploads/HkyaOCYS3.png) 當要查詢（即判斷是否存在）一個元素時，同樣將其數據輸入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 #include #include #include // 求模的素數，用於將散列函數的結果映射到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