基因比对算法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()