基因比對演算法c語言

① 全長轉錄本與參考基因組比對

長讀長測序技術的突破使得轉錄本結構鑒定、可變剪切等分析更為准確。上一期我們介紹了Iso-seq方法鑒定全長轉錄本： 全長轉錄本鑒定 。全長轉錄本鑒定之後，如何更精確地對全長轉錄本進行參考基因組比對分析，精確識別出基因異構體，是Iso-seq下游分析的關鍵。近幾年，已經開發出多款用於三代長讀長reads的比對軟體，如GMAP、minimap2、deSALT等軟體，可以較好地用於比對分析，下面我們具體看一下這些比對軟體的使用。

傳統的使用比較多的長讀長比對軟體是GMAP，05年發表公布，最開始是用來比對低通量的est序列的，後來也有進一步升級為GSNAP支持高通量的二代測序。PacBio測序技術出現後，常用於Iso-seq轉錄本的鑒定，目前仍是相關研究引用量最高的比對軟體，該軟體也一直在持續更新升級。其可以將轉錄本序列與參考基因組序列比對，輸出gff文件，比對速度稍慢，其使用方法如下：

Minimap2是生信大牛李恆18年用c語言開發的可以用於三代數據（subreads、iso-seq）比對的長序列比對軟體，與傳統的三代比對工具GMAP相比，其速度有非常顯著的提升，當然同時消耗的內存也比較大。使用方法也比較簡單，近幾年引用次數增長的也很迅速，所以大家可以試試用minimap2進行Iso-seq的比對。其使用方法如下：

deSALT是19年底哈爾濱工業大學王亞東團隊開發的一款專門針對於三代長reads測序的比對軟體，主要可以處理四種數據類型：

該軟體可以構造基於圖的比對骨架以推斷外顯子，並使用它們來生成剪接的參考序列以產生精確的比對，更好的解決了小外顯子和測序錯誤的問題。其使用方法如下：

還有兩款軟體，STAR以及BLAT也可以用於長讀長reads的比對，這兩款軟體比較早，相信大家也有所了解，雖然可用於長讀長reads的比對，但這兩款軟體並不是專門為三代測序reads比對開發的，所以大家在做Iso-seq分析時，可以更多考慮使用上文介紹的三款。

Wu T D, Watanabe C K. GMAP: a genomic mapping and alignment program for mRNA and EST sequences[J]. Bioinformatics, 2005, 21(9): 1859-1875.

Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. doi:10.1093/bioinformatics/bty191
https://github.com/lh3/minimap2

Liu, B., Liu, Y., Li, J. et al. deSALT: fast and accurate long transcriptomic read alignment with de Bruijn graph-based index. Genome Biol 20, 274 (2019). https://doi.org/10.1186/s13059-019-1895-9

② 求遺傳演算法（GA）C語言代碼

x=220;
for(i=0;i<12;i++)
{
y=202+i*16;
for(j=bits [ i][0];j<=bits [ i][1];j++)
if(g[j]==0)
g_text(x+(j-bits [ i][0])*16,y,4,"0");
else
g_text(x+(j-bits [ i][0])*16,y,4,"1");
}
}
}
void g_disp_char(x,y,x1,y1,x2,y2,v)
int x,y,x1,y1,x2,y2;
unsigned char v;
{
char c[10];
if(x>=x1&& x<=x2-8 && y>=y1 && y<=y2-10)
{
switch(v)
{
case 0: strcpy(c,"0");break;
case 1: strcpy(c,"+");break;
case 2: strcpy(c,"-");break;
case 3: strcpy(c,"x");
}
g_text(x,y,15,c);
}
}
void remove_life(n) /* 消除第n個個體 */
int n;
{
iflg[n]=0;
world[iatr[n][0]][iatr[n][1]]=0;
g_disp_unit(iatr[n][0],iatr[n][1],0);
if(food_size+1<=MAX_FOOD)
{
food_size++;
fatr[food_size-1][0]=iatr[n][0];
fatr[food_size-1][1]=iatr[n][1];
fatr[food_size-1][2]=1;
fatr[food_size-1][3]=0;
fflg[food_size-1]=1;
world[iatr[n][0]][iatr[n][1]]=5;
g_disp_unit(iatr[n][0],iatr[n][1],5);
}
}
void remove_food(n) /* 消除第n個食物 */
int n;
{
fflg[n]=0;
world[fatr[n][0]][fatr[n][1]]=0;
g_disp_unit(fatr[n][0],fatr[n][1],0);
}
void make_lives_and_foods() /* 設置虛擬環境中生物與食物 */
{
int x,y,i,j;
pop_size=0;
food_size=0;
for(y=0;y<wy;y++)
for(x=0;x<wx;x++)
{
if(world[x][y]==1||world[x][y]==2)
{
if(pop_size+1<=MAX_POP)
{
pop_size++;
/* 生成遺傳因子 */
gene[pop_size-1][0]=world[x][y]-1;
for(i=1;i<G_LENGTH;i++)
gene[pop_size-1] [ i]=random(2);
/* 設定屬性 */
iatr[pop_size-1][0]=x;
iatr[pop_size-1][1]=y;
iatr[pop_size-1][2]=70+random(30);
iatr[pop_size-1][3]=random(SL_MIN);
}
}
if(world[x][y]==3||world[x][y]==5)
{
if(food_size+1<=MAX_FOOD)
{
food_size++;
/* 設定屬性 */
fatr[food_size-1][0]=x;
fatr[food_size-1][1]=y;
if(world[x][y]==3)
fatr[food_size-1][2]=0;
else
fatr[food_size-1][2]=1;
fatr[food_size-1][3]=random(TL1-1)+1;
}
}
}
}
void find_empty(x,y) /* 尋找虛擬環境中的空處，返回坐標 */
int *x,*y;
{
int ok;
ok=0;
while(ok==0)
{
*x=random(wx);*y=random(wy);
if(world[*x][*y]==0) ok=1;
}
}
void make_world() /* 隨機設定人工環境 */
{
int i,j,k,num,x,y;
int ok,overlap;
char choice[3];
double size;
wx=0;
while(wx<10||wx>MAX_WX)
{
setcolor(15);
disp_hz16("虛擬環境長度(10-60)",10,210,20);
gscanf(300,210,4,0,3,"%s",choice);
wx=atoi(choice);
}
wy=0;
while(wy<10||wy>MAX_WY)
{
setcolor(15);
disp_hz16("虛擬環境寬度(10-32)",10,240,20);
gscanf(300,240,4,0,3,"%s",choice);
wy=atoi(choice);
}
for(i=0;i<wy;i++)
for(j=0;j<wx;j++)
if(i==0||i==wy-1||j==0||j==wx-1)
world[j] [ i]=4;
else world[j] [ i]=0;
/* 設定障礙物 */
size=(double)(wx*wy);
num=(int)(size/40.0);
if(num>MAX_POP) num=MAX_POP;
for(i=0;i<num;i++)
{
find_empty(&x,&y);
world[x][y]=4;
}
num=(int)(size/5.0);
if(num>MAX_FOOD) num=MAX_FOOD;
for(i=0;i<num;i++)
{
ok=0;
while(ok==0)
{
x=random(wx);y=random(wy);
if((world[x][y]!=4) &&
(world[x][y-1]==4 || world[x][y+1]==4 ||
world[x-1][y]==4 || world[x+1][y]==4))
{ world[x][y]=4;
ok=1;
}
}
}
for(y=0;y<wy;y++)
for(x=0;x<wx;x++)
if(world[x][y]==0)
{
num=0;
for(i=-1;i<=1;i++)
for(j=-1;j<=1;j++)
if(get_world(x+j,y+i)==4)
num++;
if(num>=6) world[x][y]=4;
}
/* 設定生物 */
num=(int)(size*R_LIFE);
for(i=0;i<num;i++)
{ find_empty(&x,&y);
world[x][y]=random(2)+1;
}
/* 設定食物 */
num=(int)(size*R_FOOD);
for(i=0;i<num;i++)
{
find_empty(&x,&y);
world[x][y]=3;
}
}
void load_world_file() /* 讀取虛擬環境數據文件設定 */
{
FILE *fopen(),*fpt;
char st[100],c;
int i,j;
if((fpt=fopen("\ga\world","r"))==NULL) exit(-1);
else
{
fscanf(fpt,"%d",&wx);
fscanf(fpt,"%d",&wy);
for(i=0;i<wy;i++)
for(j=0;j<wx;j++)
fscanf(fpt,"%d",&world[j] [ i]);
fclose(fpt);

③ 比對演算法總結(二)——基於BWT索引結構的比對演算法-Bowite1

這是美國馬里蘭大學計算機研究所、生物信息學和計算生物學中心於2009年發表在《Genome Biology》雜志的一篇經典文章，至此以後依賴於BWT索引的比對演算法成為主流。 Bowite 是一款超快速、內存佔用低的短序列比對軟體，適用於將短reads比對至大型參考基因組。採用Burrows-Wheeler 演算法建立索引的Bowite軟體可以在1 CPU時內，將2000萬條reads 比對至人參考基因組，且內存只佔有1.3Gb。於此同時Bowite 採用了新的quality-aware backtracking（質量回溯）演算法，比對過程允許錯配。

在此之前都是採用對reads (SHRiMP, Maq, RMAP,ZOOM) 或者參考基因組 (SOAP)構建哈希表的演算法進行序列比對，該演算法已在上篇文章中進行了介紹 https://www.jianshu.com/p/f5ccff73b181 。
Bowite 採用了一種完全新的索引構建策略，適用於哺乳動物重測序。根據千人基因組計劃數據，Bowite 在35bp PE 序列上的比對速度要比Maq 軟體快35 倍，比SOAP軟體快300倍。Bowite 採用 Burrows-Wheeler 演算法對 full-text minute-space (FM) 構建索引，人參考基因組佔用的內存為1.3 GB。
為了追求速度，Bowite 針對哺乳動物重測序項目進行了很多合理的折中。例如，如果一條reads有多條最優匹配，Bowite 只會輸出一條最優匹配。當輸出的最優匹配也不是完全匹配時，Bowite並不能保證在所有情況下都能輸出最高質量的匹配。在設定了較高的匹配閾值時，一小部分含有多個錯配的reads可能會比對失敗。在默認參數條件下，Bowite 的靈敏度與SOAP 相當，略低於Maq。可以在命令行手動改變參數，在犧牲更多時間的情況下，增加靈敏度，給出reads所有可能的比對結果。目前Bowite 比對的reads長度范圍為4bp - 1024bp。

Bowite 對參考基因組建立索引的方法是 Burrows-Wheeler transform (BWT) 和 FM index。Bowite 建立的人類基因組索引在硬碟上的大小為2.2GB，在比對時的內存為1.3GB。FM index 常用的精確查找方法為 Ferragina 和 Manzini 演算法。Bowite 沒有完全使用該演算法，因為該演算法不允許錯配，不能比對含有測序錯誤和變異的reads。針對這種情況，Bowite引入了新的擴展演算法：quality-aware backtracking 演算法，允許錯配並支持高質量比對；double indexing 策略，避免過度回溯；Bowite比對策略與Maq軟體相似，允許小部分的高質量reads 含有錯配，並且對所有的錯配位點的質量值設置了上限閾值。

BWT 轉換是字元串的可逆性排列，它最早應用於文本數據的壓縮，依賴BWT建立的索引，可以在較低內存下，實現大型文本的有效搜索。它被在生物信息學中有廣泛的應用，包括重復區域計數、全基因組比對、微陣列探針設計、Smith-Waterman 比對到人參考基因組。Burrows-Wheeler transform (BWT) 的轉換步驟如圖1所示：

1、輪轉排序。如圖1a 所示，（1）將字元$ 添加到文本 T （acaacg）的末尾，但需注意其中字元$ 並未實際添加到文本 T 中，且其在字母表中邏輯順序小於 T 中所有出現過的字元。（2） 然後將當前字元串的第一個字元移到最後一位，形成一個新的字元串，再將新的字元串的第一位移到最後一位形成另一個新的字元串，就這樣不斷循環這個過程，直到字元串循環完畢（即$處於第一位）,這樣就形成了一個基於原字元串的字元矩陣M（這一步原圖1a 進行了省略，見下方小圖）。（3） 然後對矩陣M的各行字元按照字典先後順序排序，獲得排序後的字元矩陣 BWM（T），矩陣的最後一列定義為 BWT（T）。 前期經過一個小復雜的過程獲得了BWT（T）列，那這一列到底有什麼用呢？其實BWT（T）列通過簡單的演算法就可以推算出原始文本T的所有信息。而經過轉換之後的BWT（T）列大量重復字元是靠近的，只儲存該列信息，可以大大提高字元壓縮比例。

2、LF-Mapping。圖1a 轉換矩陣 BWM（T）含有一種 'last first (LF) mapping' 的特性，即最後一列L中出現某字元出現的順序與第一列F某字元出現的次序時一致的。根據Supplementary1 圖中演算法1 STEPLEFT 和 演算法2 UNPERMUTE 就可以推算出BWT（T）到 T 的過程， 圖1 b記錄了整個推算過程。 詳細推算過程可參考這個博客介紹： https://blog.csdn.net/stormlovetao/article/details/7048481 。

3、reads精確匹配。使用BWT演算法的最終目的是要將短reads比對到參考基因組上，確定短reads在參考基因組上的具體位置。轉換後的BWT（T）序列，可以利用Supplementary1 圖中演算法3 EXACTMATCH 實現reads的精確匹配。圖1c 列出了 字元串 aac 比對至acaacg 的過程 。 詳細推算過程可參考這篇介紹： https://zhuanlan.hu.com/p/158901556 。

上述的BWT轉換只能用於精確的匹配，但是測序reads是含有測序錯誤和突變的，精確匹配並不適用。這里應用了 backtracking 搜索的演算法，用於允許錯配快速比對 。含有錯配的reads只是一小部分。測序reads的每個鹼基都含有唯一的測序量值，測序質量值越該位點是測序錯誤的可能越大，只有當一條read 的所有錯配的測序質量值總和小於一定閾值時可以允許錯誤匹配。
圖2顯示了精確匹配和非精確匹配的過程，backtracking 搜索過程類似於 EXACTMATCH ，首先計算連續較長的後綴矩陣。如果矩陣中沒有搜索到相應的reads，則演算法會選擇一個已經匹配的查詢位置，替換一個不同鹼基，再次進行匹配。EXACTMATCH搜索從被替換位置之後開始，這樣就可以比對就可以允許一定的錯配。backtracking 過程發生在堆棧結構的上下文中，當有替換產生時，堆棧的結構會增長；當所有結果都不匹配時，堆棧結構會收縮。
Bowite 軟體的搜索演算法是比較貪婪的，Bowite軟體會報出遇到的第一個有效比對，並不一定是在錯配數目和變異質量上的「最佳比對」。沒有查詢最優比對的原因是尋找「最佳比對」會比現有的模型慢2-3倍。而在重測序項目上，速度是更重要的因素。Bowite 也設置了可以輸出多個比對位置（-k）和所有比對位置（-a）的參數，添加這些參數後，比對速度會顯著變慢。

目前的比對軟體會有過度回溯的情況，在reads的3『端花費大量無用時間去回溯。Bowite利用『double indexing』技術減少了過度回溯的發生。簡單來說就是對正向參考基因組進行BWT轉換，稱為 『Forward index』，同時對反向（注意不是互補配對序列,是反向序列）參考基因組也進行BWT轉換，稱為『Mirror index』。 當只允許一個錯配時，比對根據reads是前半段出現錯配，還是後半段出現錯配會有兩種情況：（1）Phase1 將Forward index 載入入內存，不允許查詢reads右半段出現錯配；（2）Phase2 將Mirror index 載入如內存，不允許查詢序列的反向reads右半段（原查詢序列的左半端） 出現錯配。這樣可以避免過度回溯，提高比比對的靈敏度。 但是，如果比對軟體允許一個reads有多個錯配時，仍然會有過度回溯的現象發生，為了減少過度回溯現象的發生，這里將回溯的上限進行了限定（默認值為：125次）。

Bowite 允許使用者在高質量reads的末端（默認是28bp）設置錯配數目（默認的錯配數目是2）。高質量reads末端的28bp序列被稱為 '種子' 序列。這個『種子』序列又可分為兩等份：14bp的高質量末端稱為 『hi-half』（通常位於5『端）,14bp的低質量末端稱為『lo-half』。 如果種子序列只允許2bp 的錯配，比對會出現4 種情況：（1）種子序列中沒有錯配(case1)；（2）hi-half區域沒有錯配，lo-half區域有一個或兩個錯配（case2）;（3）lo-half區域沒有錯配，hi-half區域有一個或兩個錯配（case3）;（4）lo-half區域有一個錯配，hi-half區域有一個錯配（case4）；
在所有情況下，reads的非種子部分允許任意數目的錯配。如圖3所示，Bowite 演算法會根據上面4 種情況交替變化『Forward index』和『Mirror index』比對策略，主要會有三種比對策略。

Bowite 建立一次參考基因組索引後，後續的比對可反復使用該索引。表1和表2列出了在默認參數條件下，Bowite、SOAP、Maq軟體性能的比較。在reads比對率相近的條件下，Bowite軟體的比對速度速度相對於SOAP、Maq軟體有較大的提升。

1、將reads 比對至人參考基因組上，Bowite相對於SOAP和Maq軟體有較大的優勢。它運行的內存非常小（1.2GB），在相同靈敏度下，速度有了較大的提升。
2、Bowite 軟體建立一次參考基因組索引後，後續的比對可反復使用該索引。
3、Bowite 速度快、內存佔用小、靈敏度高主要是因為使用了BWT演算法構建索引、利用回溯演算法允許錯配、採用Double index策略避免過度回溯。
4、Bowite 軟體目前並不支持插入、缺失比對，這個是今後需要努力的方向。

[1] Langmead B . Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome biology, 2009, 10(3):R25.
[2] BWT 推算過程參考博客 https://blog.csdn.net/stormlovetao/article/details/7048481
[3] FM index 精確查匹配過程參考文章 https://zhuanlan.hu.com/p/158901556

④ C語言 基因相關性

支持命令行和數據文件

defgenerate(n)://n為長度
foriinrange(4**n):
a=[0foriinrange(n)]//a為長度為n的一個序列
num=i
z=0
while(num!=0):
a[z]=num%4
z=z+1
num=int(num/4)

forjina:

ifj==0:
print('A',end="")
elifj==1:
print('T',end="")
elifj==2:
print('G',end="")
else:
print('C',end="")
print()