生物信息学系列博客索引
生物信息学（1）——双序列比对之Needleman-Wunsch(NW)算法详解及C++实现
生物信息学（2）——双序列比对之Smith-Waterman(SW)算法详解
生物信息学（3）——双序列比对之BLAST算法简介
生物信息学（4）——多序列比对之CLUSTAL算法详解及C++实现
生物信息学（5）——基于CUDA的多序列比对并行算法的设计与代码实现
项目gitee地址，内附源码、论文等文档信息
1. 设计并行算法的原因
虽然 CLUSTAL 算法具有较高的精度，但是由于其构造导向树的距离，即两 两比对的过程，需要迭代调用双序列比对算法，这将使得算法的时间复杂度来到 n 的四次方，数据愈多，需要消耗的时间就指数增加，对 CPU 造成巨大压力。 CUDA 平台调度 GPU 提供的并行计算能力有效的缓解了这一问题。
2. CUDA 下 CLUSTAL 执行流程
由于 CLUSTAL 算法分为三步:两两比对、构建向导树、渐近比对，最耗时 间的为迭代调用双序列比对算法进行两两比对，且此时两两比对数据无写相关， 公共数据仅仅是读共同的序列簇，所以就可以在此处利有 CUDA 调度 GPU 提供 的并行计算能力，并行的进行两两序列比对，以此来加快 CLUSTAL 算法的执行 速度
3. 代码
有些实现的细节见文末，这是CUDA代码，安装CUDA需要先安装vs2015，然后根据电脑GPU版本去NVIDIA官网安装对应的CUDA。
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include<iostream>
#include <string>
#include<algorithm>
#include <vector>
#include <iterator>
#include<ctime>
//声明命名空间std
using namespace std;
#define MATCH 1
#define DIS_MATCH -1
#define INDEL -2
#define INDEL_CHAR '-'
class ResUnit { //一次双序列比对后的结果
public:
char str1[30]; //原始序列1
char str2[30]; //原始序列2
char res1[30]; //结果序列1
char res2[30]; //结果序列2
int str1Len;
int str2Len;
int res1Len;
int res2Len;
int score; //序列总得分，反映两个序列的相似程度
int tag; //禁止迭代多次
};
class SingleSeq { //一个序列被整合后的样子
public:
string str; //一个序列的原始序列
string res; //一个序列被整合后的样子
};
class ResUnit2 { //一次双序列比对后的结果
public:
string str1; //原始序列1
string str2; //原始序列2
string res1; //结果序列1
string res2; //结果序列2
int score; //序列总得分，反映两个序列的相似程度
int tag; //禁止迭代多次
};
struct BacktrackingUnit {
int goUp; //是否向上回溯
int goLeftUp; //是否向左上回溯
int goLeft; //是否向左回溯
int score; //得分矩阵第(i, j)这个单元的分值
};
typedef struct BacktrackingUnit *unitLine;
__device__ int max3(int a, int b, int c);
__device__ int myCompare(char a, char b);
int max32(int a, int b, int c);
int myCompare2(char a, char b);
bool complare2(const ResUnit2 &a, const ResUnit2 &b);
ResUnit2 traceback(unitLine** item, const int i, const int j, string str1, string str2, string res1, string res2, int n, ResUnit2 resUnit);
ResUnit2 NeedlemanWunch(string str1, string str2);
int ifStrInQueueFinish(string str, vector<SingleSeq> queue_finish);
vector<SingleSeq> RegularTwo(ResUnit2 tag, ResUnit2 temp, vector<SingleSeq> queue_finish);
vector<SingleSeq> RegularSeq1(ResUnit2 tag, vector<SingleSeq> queue_finish, int item);
vector<SingleSeq> RegularSeq2(ResUnit2 tag, vector<SingleSeq> queue_finish, int item);
/**
kernel函数 数据并行处理函数
*/
__global__ void kernelFunc(ResUnit *res)
{
int x = threadIdx.x;
char *str1 = res[x].str1;
char *str2 = res[x].str2;
//字符串str1,str2长度
const int m = res[x].str1Len;
const int n = res[x].str2Len;
int m1, m2, m3, mm;
BacktrackingUnit **unit;
unit = new BacktrackingUnit*[m+1];
for (int i = 0; i < m+1; i++)
{
unit[i] = new BacktrackingUnit[n+1];
for (int j = 0; j < n + 1;j++) {
unit[i][j] = BacktrackingUnit();
unit[i][j].goUp = 0;
unit[i][j].goLeftUp = 0;
unit[i][j].goLeft = 0;
}
}
unit[0][0].score = 0;
for (int i = 1; i <= m; i++) {
unit[i][0].score = -2 * i;
unit[i][0].goUp = 1;
}
for (int j = 1; j <= n; j++) {
unit[0][j].score = -2 * j;
unit[0][j].goLeft = 1;
}
// 动态规划算法计算得分矩阵每个单元的分值
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= n; j++) {
m1 = unit[i - 1][j].score + -2;
m2 = unit[i - 1][j - 1].score + myCompare(str1[i - 1], str2[j - 1]);
m3 = unit[i][j - 1].score + -2;
mm = max3(m1, m2, m3);
unit[i][j].score = mm;
//判断路径来源
if (m1 == mm) unit[i][j].goUp = 1;
if (m2 == mm) unit[i][j].goLeftUp = 1;
if (m3 == mm) unit[i][j].goLeft = 1;
}
}
//开始回溯
int i = m, int j = n;
ResUnit resUnit;
resUnit.tag = 0;
resUnit.res1Len = 0;
resUnit.res2Len = 0;
resUnit.score = unit[m][n].score;
while (resUnit.tag != 1)
{
BacktrackingUnit temp = unit[i][j];
if (!(i || j)) { // 到矩阵单元(0, 0)才算结束，这代表初始的两个字符串的每个字符都被比对到了
resUnit.tag = 1;
}
if (temp.goUp) { // 向上回溯一格
resUnit.res1[resUnit.res1Len] = str1[i - 1];
resUnit.res2[resUnit.res2Len] = '-';
++resUnit.res1Len;
++resUnit.res2Len;
--i;
}
if (temp.goLeftUp) { // 向左上回溯一格
resUnit.res1[resUnit.res1Len] = str1[i - 1];
resUnit.res2[resUnit.res2Len] = str2[j - 1];
++resUnit.res1Len;
++resUnit.res2Len;
--i;
--j;
}
if (temp.goLeft) { // 向左回溯一格
resUnit.res1[resUnit.res1Len] = '-';
resUnit.res2[resUnit.res2Len] = str2[j - 1];
++resUnit.res1Len;
++resUnit.res2Len;
--j;
}
}
delete(unit);
for (int k = 0; k < m;k++)
{
resUnit.str1[k] = str1[k];
}
for (int k = 0; k < n;k++)
{
resUnit.str2[k] = str2[k];
}
res[x] = resUnit;
}
int main()
{
const int sequence_groups = 31;
string *ss = new string[sequence_groups];
ss[0] = "AATTCCGGAAAATTCCGGGACA";
ss[1] = "AATTCGAACGAATGGACATCGA";
ss[2] = "GTTCGGAAGTTCCGGAAAAAGG";
ss[3] = "AACTCCGCGAACTCCGCGGACA";
ss[4] = "ATTCCGGGAAATTGGACAGGG";
ss[5] = "ATTCCGGACACCGGGAAACTTAA";
ss[6] = "CTCCGGGAAATTCCGGGACAG";
ss[7] = "ATTCCGGGAAGGACATCCGGGA";
ss[8] = "ATCACCGGGAAAGGGACGGACA";
ss[9] = "CCGCCGGATTAGGGGACAGCAA";
ss[10] = "CCGGACCCGTGGGACAAATT";
ss[11] = "CTCGCCCGGACATTAGGAATT";
ss[12] = "CACCGGACATTAGGACAGGGC";
ss[13] = "CGCGCCGGCCGGGGACAACGG";
ss[14] = "CCGGACAATTAAGGGACGGACA";
ss[15] = "CCGGACAGGACATTCATTAGG";
ss[16] = "CTCGATTAGGGCACATTAGGA";
ss[17] = "CACCGGACATTAGGGGGACAC";
ss[18] = "CGATTAGGGCGCCGGACGGAC";
ss[19] = "CCGGACTAGGGCGGACTGACT";
ss[20] = "CCGTAGGGCTTAAGGGACTGACT";
ss[21] = "CCGGACAATTAGGGCTGACT";
ss[22] = "CTAGGGCATTAAGGGACTGACT";
ss[23] = "CCGGACAATTTAGGGCTGACT";
ss[24] = "CCGGTAGGGCTAAGGGACTGACT";
ss[25] = "CCGGACTAGGGCGGACTGACT";
ss[26] = "CCGTAGGGCATTAAACTGACT";
ss[27] = "CCGTTAAACTAGGGACTGACT";
ss[28] = "CCGGACTTAAACTGGACTGACT";
ss[29] = "CTTAAACTTTAAGGGACTGACT";
ss[30] = "CCGGACAATTTAAACTGACT";
vector<ResUnit2> queue_initial(0); //定义等待整合的队列，是一个ResUnit对象vector
vector<SingleSeq> queue_finish(0); //定义整合完毕的队列，是一个String类型的vector
vector<ResUnit2>::iterator it; //迭代器，访问queue_initial
vector<SingleSeq>::iterator it_finish; //迭代器，访问queue_finish
clock_t startTime, endTime;
startTime = clock();//计时开始
//定义队列长度
const int queue_length = ((sequence_groups - 1)*sequence_groups) / 2;
//初始化矩阵
//res保存的是待排序的双序列组合，对象数组长度为queue_length
ResUnit res[queue_length];
int k = 0;
//初始化，将string类型的字符串转入char[]数组
for (int i = 0;i < sequence_groups;i++) {
for (int j = i + 1;j < sequence_groups;j++) {
//只遍历上三角区域
for (int m = 0; m < ss[i].length(); m++)
{
res[k].str1[m] = ss[i][m];
}
res[k].str1Len = ss[i].length();
for (int m = 0; m < ss[j].length(); m++)
{
res[k].str2[m] = ss[j][m];
}
res[k].str2Len = ss[j].length();
k++;
}
}
//定义GPU上的变量
ResUnit *resGpu;
//分配地址
cudaMalloc((void**)&resGpu, queue_length * sizeof(ResUnit));
//将主机上的数据转到设备上
cudaMemcpy(resGpu, res, queue_length * sizeof(ResUnit), cudaMemcpyHostToDevice);
//分配线程
kernelFunc << <1, queue_length >> > (resGpu);
//将设备上的数据拷贝回主机上
cudaMemcpy(res, resGpu, queue_length * sizeof(ResUnit), cudaMemcpyDeviceToHost);
//释放
cudaFree(resGpu);
//拷贝到queue_initial
for (int i = 0; i < queue_length; i++)
{
ResUnit2 resUnit2;
string str1 = res[i].str1;
string str2 = res[i].str2;
string res1 = strrev(res[i].res1);
string res2 = strrev(res[i].res2);
resUnit2.str1 = str1;
resUnit2.str2 = str2;
resUnit2.res1 = res1;
resUnit2.res2 = res2;
resUnit2.score = res[i].score;
queue_initial.push_back(resUnit2);
//cout << "str1:" << str1 << ",str2:" << str2 << ",res1:" << res1 << ",res2:" << res2 << endl;
}
//clustal算法
//排序
sort(queue_initial.begin(), queue_initial.end(), complare2);
//最多循环queue_length次
for (int i = 0;i < queue_length;i++) {
//当结果队列长度与sequence_groups相等之后，就说明全部放入了结果队列，可以跳出循环了
if (sequence_groups == queue_finish.size())
break;
//一个ResUnit对象的str1属性和str2均不在
if (ifStrInQueueFinish(queue_initial.at(i).str1, queue_finish) < 0 && ifStrInQueueFinish(queue_initial.at(i).str2, queue_finish) < 0) {
SingleSeq singleSeq1, singleSeq2;
singleSeq1.str = queue_initial.at(i).str1;
singleSeq1.res = queue_initial.at(i).res1;
singleSeq2.str = queue_initial.at(i).str2;
singleSeq2.res = queue_initial.at(i).res2;
//如果结果队列已经有元素，，且又来了俩不相干的，却很匹配的序列对
if (queue_finish.size()>0) {
//将结果队列第一个的序列和queue_initial.at(i).str1进行双序列比对
ResUnit2 temp = NeedlemanWunch(queue_finish.front().str, queue_initial.at(i).str1);
//进行规整操作
queue_finish = RegularTwo(queue_initial.at(i), temp, queue_finish);
}
else
{
queue_finish.push_back(singleSeq1);
queue_finish.push_back(singleSeq2);
}
}
//str1在，str2不在
else if (ifStrInQueueFinish(queue_initial.at(i).str1, queue_finish) > -1 && ifStrInQueueFinish(queue_initial.at(i).str2, queue_finish) < 0) {
int item = ifStrInQueueFinish(queue_initial.at(i).str1, queue_finish);
queue_finish = RegularSeq1(queue_initial.at(i), queue_finish, item);
}
//str2在，str1不在
else if (ifStrInQueueFinish(queue_initial.at(i).str2, queue_finish) > -1 && ifStrInQueueFinish(queue_initial.at(i).str1, queue_finish) < 0) {
int item = ifStrInQueueFinish(queue_initial.at(i).str2, queue_finish);
queue_finish = RegularSeq2(queue_initial.at(i), queue_finish, item);
}
}
//声明一个迭代器，来访问vector容器
cout << "--------------------------------------------------------------------" << endl;
cout << endl << "原序列" << endl;
for (it_finish = queue_finish.begin();it_finish != queue_finish.end();it_finish++) {
cout << it_finish->str << endl;
}
cout << endl << "比对后" << endl;
for (it_finish = queue_finish.begin();it_finish != queue_finish.end();it_finish++) {
cout << it_finish->res << endl;
}
endTime = clock();//计时结束
cout << endl << endl << "The run time is: " << (double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
cout << "--------------------------------------------------------------------" << endl;
system("pause");
return 0;
}
/**
比较三种路径之间谁最大
f(i-1,j-1),f(i-1,j)+indel,f(i,j-1)+indel
*/
__device__ int max3(int a, int b, int c) {
int temp = a > b ? a : b;
return temp > c ? temp : c;
}
/**
比较两个字符类型属于什么，match，dismatch，indel
*/
__device__ int myCompare(char a, char b) {
if (a == b)
return 1;
else if (a == ' ' || b == ' ')
return -2;
else
return -1;
}
/*
规整函数，规整两个序列情况
queue_finish temp tag
A1 A2 E1
B E2 F
C
D
*/
vector<SingleSeq> RegularTwo(ResUnit2 tag, ResUnit2 temp, vector<SingleSeq> queue_finish) {
string E2 = temp.res2;
string E1 = tag.res1;
string A1 = queue_finish.front().res;
string A2 = temp.res1;
string F = tag.res2;
string tempStr = "";
vector<SingleSeq>::iterator it;//声明一个迭代器，来访问vector容器
int i = 0, j = 0;
//第一步，，整合tag与temp
while (E2 != E1 && j < E1.length() && i < E2.length()) {
if (E2[i] == E1[j]) {
i++;
j++;
}
else {
if (E2[i] == '-') {
E1.insert(j, "-");
F.insert(j, "-");
}
else if (E1[j] == '-')
{
E2.insert(i, "-");
A2.insert(i, "-");
}
}
}
if (i == E2.length()) {
//E2先到头
for (int k = 0; k < E1.length() - j; k++)
{
tempStr += "-";
}
E2 += tempStr;
A2 += tempStr;
}
else if (j == E1.length()) {
//E1先到头
for (int k = 0; k < E2.length() - i; k++)
{
tempStr += "-";
}
E1 += tempStr;
F += tempStr;
}
//将tempStr置空
tempStr = "";
//第二步 融合进queue_finish
i = 0, j = 0;
while (A1 != A2 && i < A1.length() && j < A2.length()) {
if (A1[i] == A2[j]) {
i++;
j++;
}
else {
if (A1[i] == '-') {
A2.insert(j, "-");
E1.insert(j, "-");
F.insert(j, "-");
}
else if (A2[j] == '-')
{
A1.insert(i, "-");
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res.insert(i, "-");
}
}
}
}
if (i == A1.length()) {
//A1先到头
for (int k = 0; k < A2.length() - j; k++)
{
tempStr += "-";
}
A1 += tempStr;
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res + tempStr;
}
}
else if (j == A2.length()) {
//A2先到头
for (int k = 0; k < A1.length() - i; k++)
{
tempStr += "-";
}
A2 += tempStr;
E1 += tempStr;
F += tempStr;
}
//规划好之后，，将 E F 插入queue_finish尾部
SingleSeq sE, sF;
sE.res = E1;
sE.str = tag.str1;
sF.res = F;
sF.str = tag.str2;
queue_finish.push_back(sE);
queue_finish.push_back(sF);
return queue_finish;
}
/*
规整函数，规整序列1情况
queue_finish tag
A1 A2
B E
C
D
*/
vector<SingleSeq> RegularSeq1(ResUnit2 tag, vector<SingleSeq> queue_finish, int item) {
SingleSeq main_seq = queue_finish.at(item);//找到和seq1相同的序列
string A1 = main_seq.res;
string A2 = tag.res1;
string E = tag.res2;
string tempStr = "";
vector<SingleSeq>::iterator it;//声明一个迭代器，来访问vector容器
int i = 0, j = 0;
while (A1 != A2 && i < A1.length() && j < A2.length()) {
if (A1[i] == A2[j]) {
i++;
j++;
}
else {
if (A1[i] == '-') {
A2.insert(j, "-");
E.insert(j, "-");
}
else if (A2[j] == '-') {
//遍历queue_finish,给queue_finish内res洗头
A1.insert(i, "-");
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res.insert(i, "-");
}
}
}
}
if (i == A1.length()) {
//A1先到头
for (int k = 0; k < A2.length() - j; k++)
{
tempStr += "-";
}
A1 += tempStr;
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res + tempStr;
}
}
else if (j == A2.length()) {
//A2先到头
for (int k = 0; k < A1.length() - i; k++)
{
tempStr += "-";
}
A2 += tempStr;
E += tempStr;
}
//添加
SingleSeq sE;
sE.res = E;
sE.str = tag.str2;
queue_finish.push_back(sE);
return queue_finish;
}
/*
规整函数，规整序列2情况
queue_finish tag
A1 E
B A2
C
D
*/
vector<SingleSeq> RegularSeq2(ResUnit2 tag, vector<SingleSeq> queue_finish, int item) {
SingleSeq main_seq = queue_finish.at(item);//找到和seq1相同的序列
string A1 = main_seq.res;
string A2 = tag.res2;
string E = tag.res1;
string tempStr = "";
vector<SingleSeq>::iterator it;//声明一个迭代器，来访问vector容器
int i = 0, j = 0;
while (A1 != A2 && i < A1.length() && j < A2.length()) {
if (A1[i] == A2[j]) {
i++;
j++;
}
else {
if (A1[i] == '-') {
A2.insert(j, "-");
E.insert(j, "-");
}
else if (A2[j] == '-') {
//遍历queue_finish,给queue_finish内res洗头
A1.insert(i, "-");
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res.insert(i, "-");
}
}
}
}
if (i == A1.length()) {
//A1先到头
for (int k = 0; k < A2.length() - j; k++)
{
tempStr += "-";
}
A1 += tempStr;
for (it = queue_finish.begin();it != queue_finish.end();it++) {
it->res = it->res + tempStr;
}
}
else if (j == A2.length()) {
//A2先到头
for (int k = 0; k < A1.length() - i; k++)
{
tempStr += "-";
}
A2 += tempStr;
E += tempStr;
}
//添加
SingleSeq sE;
sE.res = E;
sE.str = tag.str1;
queue_finish.push_back(sE);
return queue_finish;
}
//判断一个str是否有与queue_finish数组对象内的seq相等的,没有返回-1,有就返回序号
int ifStrInQueueFinish(string str, vector<SingleSeq> queue_finish) {
int i = 0;
vector<SingleSeq>::iterator it;//声明一个迭代器，来访问vector容器
for (it = queue_finish.begin();it != queue_finish.end();it++) {
if (str == it->str)
return i;
i++;
}
return -1;
}
ResUnit2 NeedlemanWunch(string str1, string str2) {
//字符串str1,str2长度
const int m = str1.length();
const int n = str2.length();
int m1, m2, m3, mm;
unitLine **unit;
// 初始化
if ((unit = (unitLine **)malloc(sizeof(unitLine*) * (m + 1))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
for (int i = 0; i <= m; i++) {
if ((unit[i] = (unitLine *)malloc(sizeof(unitLine) * (n + 1))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
for (int j = 0; j <= n; j++) {
if ((unit[i][j] = (unitLine)malloc(sizeof(struct BacktrackingUnit))) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
unit[i][j]->goUp = 0;
unit[i][j]->goLeftUp = 0;
unit[i][j]->goLeft = 0;
}
}
unit[0][0]->score = 0;
for (int i = 1; i <= m; i++) {
unit[i][0]->score = INDEL * i;
unit[i][0]->goUp = 1;
}
for (int j = 1; j <= n; j++) {
unit[0][j]->score = INDEL * j;
unit[0][j]->goLeft = 1;
}
// 动态规划算法计算得分矩阵每个单元的分值
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= n; j++) {
m1 = unit[i - 1][j]->score + INDEL;
m2 = unit[i - 1][j - 1]->score + myCompare2(str1[i - 1], str2[j - 1]);
m3 = unit[i][j - 1]->score + INDEL;
mm = max32(m1, m2, m3);
unit[i][j]->score = mm;
//判断路径来源
if (m1 == mm) unit[i][j]->goUp = 1;
if (m2 == mm) unit[i][j]->goLeftUp = 1;
if (m3 == mm) unit[i][j]->goLeft = 1;
}
}
//开始回溯
ResUnit2 res;
res.tag = 0;
res = traceback(unit, m, n, str1, str2, "", "", 0, res);
res.score = unit[m][n]->score;
//释放内存
for (int i = 0; i <= m; i++) {
for (int j = 0; j <= n; j++) {
free(unit[i][j]);
}
free(unit[i]);
}
free(unit);
//返回值
return res;
}
ResUnit2 traceback(unitLine** item, const int i, const int j, string str1, string str2, string res1, string res2, int n, ResUnit2 resUnit) {
unitLine temp = item[i][j];
if (resUnit.tag != 1)
{
if (!(i || j)) { // 到矩阵单元(0, 0)才算结束，这代表初始的两个字符串的每个字符都被比对到了
resUnit.str1 = str1;
resUnit.str2 = str2;
resUnit.res1 = res1;
resUnit.res2 = res2;
resUnit.tag = 1;
return resUnit;
}
if (temp->goUp) { // 向上回溯一格
res1 = str1[i - 1] + res1;
res2 = INDEL_CHAR + res2;
resUnit = traceback(item, i - 1, j, str1, str2, res1, res2, n + 1, resUnit);
}
if (temp->goLeftUp) { // 向左上回溯一格
res1 = str1[i - 1] + res1;
res2 = str2[j - 1] + res2;
resUnit = traceback(item, i - 1, j - 1, str1, str2, res1, res2, n + 1, resUnit);
}
if (temp->goLeft) { // 向左回溯一格
res1 = INDEL_CHAR + res1;
res2 = str2[j - 1] + res2;
resUnit = traceback(item, i, j - 1, str1, str2, res1, res2, n + 1, resUnit);
}
return resUnit;
}
else
{
return resUnit;
}
}
bool complare2(const ResUnit2 &a, const ResUnit2 &b)
{
return a.score > b.score;
}
/**
比较三种路径之间谁最大
f(i-1,j-1),f(i-1,j)+indel,f(i,j-1)+indel
*/
int max32(int a, int b, int c) {
int temp = a > b ? a : b;
return temp > c ? temp : c;
}
/**
比较两个字符类型属于什么，match，dismatch，indel
*/
int myCompare2(char a, char b) {
if (a == b)
return MATCH;
else if (a == ' ' || b == ' ')
return INDEL;
else
return DIS_MATCH;
}
4. 结果分析
由表 5-2 a 组数据可知，当序列长度和序列数量都很小时候，CUDA 并没有 展现出他的能力，由于要将数据从主机端拷贝到设备端，且运行结束又有将设备 端数据拷贝至主机端，所以相对于 CPU 版本代码，GPU 版代码运行的效果非常 的差，甚至于来说，运行时间根本不在一个数量级。
由表 5-2 b 组数据可知，当数据量微有上升，CUP 与 GPU 运行时间虽然还 是有差距，但是已经将时间缩短至一个数量级。
表 5-2 c 组与 d 组数据放一起比对，可得当序列数量不变，序列长度增加时， GPU 对算法的加速效果虽有但甚微。
表 5-2 c 组与 d 组数据放一起比对，可得，当序列长度不变，序列数量增加 时，GPU 明显优于 CPU 版本代码执行速度，具有明显的加速效果。
GPU加速CLUSTAL算法是通过并行执行CLUSTAL算法过程中两两比对执 行双序列比对来进行加速。设有 n 条序列，则共需迭代执行 n(n-1)/2 次双序列比 对过程。当 n 数量上升时，CPU 运行时间将明显上升。GPU 并行调用 Kernel 函 数中的 NW 算法，分配 n(n-1)/2 个 block 并行执行，只需额外增加拷贝数据的时 间，便可以主观上人为的运行一次双序列比对算法时间，完成 n(n-1)/2 次双序列 比对，所以时间明显减少。
5. 代码的一些细节
在实现细节上，由于 C++与 C 未提供对字符串的设 备端处理函数，即 CUDA 中的 kernel 核函数中不能调用 strLen()等方法处理字符 串。本文通过将字符串看成字符数组进行字符串的添加操作，由于需要比对的序 列长度大致相仿，联系实际，用预设定长度的字符数组处理字符串，亦是可行之 举。
6. 总结
显然，博文的题目是篇毕业论文，当你看到这篇博文的时候，我也已经顺利的毕业，大学四年自此画上了句号。四年的本科生涯并没有学习生物信息学的我，也是这半年加班加点、临时抱佛脚学的，有些内容或者概念可能并不是理解的很透彻，如果有错误也请您不吝赐教，加以斧正。独学而无友，孤陋而寡闻，学习之途我也迷茫过，无助过，正因为此，我才知道分享的重要性。希望屏幕前的陌生人，能从我这获取一些帮助，也算是在我平凡的生活里，增添一抹荣光了。
632
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
