序列比对(十三)——后验解码

最新推荐文章于 2024-08-17 21:30:16 发布
最新推荐文章于 2024-08-17 21:30:16 发布
阅读量379 收藏
点赞数
分类专栏： # 序列算法 文章标签： 序列比对 生物信息 C语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/biocity/article/details/96774163
版权

原创: hxj7

本文介绍了如何利用后验概率进行解码，可称为后验解码。

前文《序列比对（12）计算后验概率》介绍了如何计算某一位置可能状态的后验概率。那么可以据此找到某一位置最有可能的状态。即
在这里插入图片描述

从上面的公式可以看出，有两种方法求解某一位置的最可能状态。如果是依据公式(1)，先计算出后验概率，然后找到其中最大后验概率对应的状态；如果是依据公式(3)，无需计算后验概率，比较简单。

本文所用代码采用了公式(1)的方法，虽然稍嫌麻烦，但是可以在原先的代码基础上稍加增改即可。运行效果如下：
在这里插入图片描述

具体代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
//#define MIN_LOG_VALUE -15
//#define SAFE_EXP(x) ((x) < MIN_LOG_VALUE ? 0 : exp(x))

typedef char State;
typedef char Result;
State state[] = {'F', 'L'};   // 所有的可能状态
Result result[] = {'1', '2', '3', '4', '5', '6'};   // 所有的可能符号
double init[] = {0.9, 0.1};    // 初始状态的概率向量
double emission[][6] = {   // 发射矩阵：行对应着状态，列对应着符号
  1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.5
};
double trans[][2] = {   // 转移矩阵：行和列都是状态
  0.95, 0.05,
  0.1, 0.9
};
const int nstate = 2;
const int nresult = 6;

double** fscore;  // 前向算法的得分矩阵
double** bscore;  // 后向算法的得分矩阵
double* scale;   // 缩放因子向量
double logScaleSum;

int random(double* prob, const int n);
void randSeq(State* st, Result* res, const int n);
int getResultIndex(Result r);
void printState(State* st, const int n);
void printResult(Result* res, const int n);
double forward(Result* res, const int n);
double backward(Result* res, const int n);
double** getPostProb(const int n);
void postDecode(double** prob, const int n);

int main(void) {
  int i;
  int n = 5;
  State* rst;   // 一串随机状态序列
  Result* rres;  // 一串随机符号序列
  double** postProb;
  if ((rst = (State*) malloc(sizeof(State) * n)) == NULL || \
      (rres = (Result*) malloc(sizeof(Result) * n)) == NULL || \
      (scale = (double*) malloc(sizeof(double) * n)) == NULL || \
      (fscore = (double**) malloc(sizeof(double*) * nstate)) == NULL || \
      (bscore = (double**) malloc(sizeof(double*) * nstate)) == NULL) {
    fputs("Error: out of space!\n", stderr);
    exit(1);
  }
  for (i = 0; i < nstate; i++) {
    if ((fscore[i] = (double*) malloc(sizeof(double) * n)) == NULL || \
        (bscore[i] = (double*) malloc(sizeof(double) * n)) == NULL) {
      fputs("Error: out of space!\n", stderr);
      exit(1);
    }
  }
  randSeq(rst, rres, n);
  printState(rst, n);
  printResult(rres, n);
  forward(rres, n);
  backward(rres, n);
  postProb = getPostProb(n);
  postDecode(postProb, n);
  free(rst);
  free(rres);
  free(scale);
  free(fscore);
  free(bscore);
  for (i = 0; i < nstate; i++)
    free(postProb[i]);
  free(postProb);
}

// 根据一个概率向量从0到n-1随机抽取一个数
int random(double* prob, const int n) {
  int i;
  double p = rand() / 1.0 / (RAND_MAX + 1);
  for (i = 0; i < n - 1; i++) {
    if (p <= prob[i])
      break;
    p -= prob[i];
  }
  return i;
}

// 根据转移矩阵和发射矩阵生成一串随机状态和符号
void randSeq(State* st, Result* res, const int n) {
  int i, ls, lr;
  srand((unsigned int) time(NULL));
  ls = random(init, nstate);
  lr = random(emission[ls], nresult);
  st[0] = state[ls];
  res[0] = result[lr];
  for (i = 1; i < n; i++) {
    ls = random(trans[ls], nstate);
    lr = random(emission[ls], nresult);
    st[i] = state[ls];
    res[i] = result[lr];
  }
}

int getResultIndex(Result r) {
  return r - result[0];
}

// 前向算法计算P(x)
double forward(Result* res, const int n) {
  int i, l, k, idx;
  double logpx;
  // 缩放因子向量初始化
  for (i = 0; i < n; i++)
    scale[i] = 0;
  // 计算第0列分值
  idx = getResultIndex(res[0]);
  for (l = 0; l < nstate; l++) {
    fscore[l][0] = emission[l][idx] * init[l];
    scale[0] += fscore[l][0];
  }
  for (l = 0; l < nstate; l++)
    fscore[l][0] /= scale[0];
  // 计算从第1列开始的各列分值
  for (i = 1; i < n; i++) {
    idx = getResultIndex(res[i]);
    for (l = 0; l < nstate; l++) {
      fscore[l][i] = 0;
      for (k = 0; k < nstate; k++) {
        fscore[l][i] += fscore[k][i - 1] * trans[k][l];
      }
      fscore[l][i] *= emission[l][idx];
      scale[i] += fscore[l][i];
    }
    for (l = 0; l < nstate; l++)
      fscore[l][i] /= scale[i];
  }
  // P(x) = product(scale)
  // P(x)就是缩放因子向量所有元素的乘积
  logpx = 0;
  for (i = 0; i < n; i++)
    logpx += log(scale[i]);
  printf("forward: logP(x) = %f\n", logpx);
  logScaleSum = logpx;
/*
  for (l = 0; l < nstate; l++) {
    for (i = 0; i < n; i++)
      printf("%f ", fscore[l][i]);
    printf("\n");
  }
*/
  return exp(logpx);
}

// 后向算法计算P(x)
// backward算法中使用的缩放因子和forward中的一样
double backward(Result* res, const int n) {
  int i, l, k, idx;
  double tx, logpx;
  // 计算最后一列分值
  for (l = 0; l < nstate; l++)
    bscore[l][n - 1] = 1 / scale[n - 1];
  // 计算从第n - 2列开始的各列分值
  for (i = n - 2; i >= 0; i--) {
    idx = getResultIndex(res[i + 1]);
    for (k = 0; k < nstate; k++) {
      bscore[k][i] = 0;
      for (l = 0; l < nstate; l++) {
        bscore[k][i] += bscore[l][i + 1] * trans[k][l] * emission[l][idx];
      }
    }
    for (l = 0; l < nstate; l++)
      bscore[l][i] /= scale[i];
  }
  // 计算P(x)
  tx = 0;
  idx = getResultIndex(res[0]);
  for (l = 0; l < nstate; l++)
    tx += init[l] * emission[l][idx] * bscore[l][0];
  logpx = log(tx) + logScaleSum;
  printf("backward: logP(x) = %f\n", logpx);
/*
  for (l = 0; l < nstate; l++) {
    for (i = 0; i < n; i++)
      printf("%f ", bscore[l][i]);
    printf("\n");
  }
*/
  return exp(logpx);  
}

// 计算后验概率
double** getPostProb(const int n) {
  int i, k;
  double** postProb;
  //double logdiff;
  if ((postProb = (double**) malloc(sizeof(double*) * nstate)) == NULL) {
    fputs("Error: out of space!\n", stderr);
    exit(1);  
  }
  for (k = 0; k < nstate; k++) {
    if ((postProb[k] = (double*) malloc(sizeof(double) * n)) == NULL) {
      fputs("Error: out of space!\n", stderr);
      exit(1);
    }
  }
  // 计算后验概率
  for (i = 0; i < n; i++) {
    for (k = 0; k < nstate; k++) {
      postProb[k][i] = scale[i] * fscore[k][i] * bscore[k][i];
    }
  }
  printf("\n");
  printf("Posterior Probabilities:\n");
  for (k = 0; k < nstate; k++) {
    for (i = 0; i < n; i++)
      printf("%f ", postProb[k][i]);
    printf("\n");
  }
  return postProb;
}

void postDecode(double** prob, const int n) {
  int i, k;
  double maxCol;
  int idx;
  State* st;
  if ((st = (State*) malloc(sizeof(State) * n)) == NULL) {
    fputs("Error: out of space!\n", stderr);
    exit(1);    
  }
  for (i = 0; i < n; i++) {
    idx = 0;
    maxCol = prob[0][i];
    for (k = 1; k < nstate; k++)
      if (prob[k][i] > maxCol) {
        maxCol = prob[k][i];
        idx = k;
      }
    st[i] = state[idx];
  }
  printf("\n");
  printf("Posterior Decode:\n");
  printState(st, n);
  free(st);
}

void printState(State* st, const int n) {
  int i;
  for (i = 0; i < n; i++)
    printf("%c", st[i]);
  printf("\n");
}

void printResult(Result* res, const int n) {
  int i;
  for (i = 0; i < n; i++)
    printf("%c", res[i]);
  printf("\n");  
}

(公众号：生信了)

生信了（公众号同名）
关注
专栏目录
【Image captioning】AI算法说——图像描述（Image captioning）
安静到无声
06-11 1万+
COCO数据集是一个大规模常见物体检测、分割和图像描述的数据集。它包含超过330k张图像，其中包括2.5 million个目标实例，每张图像都有5个不同的图像描述。COCO数据集已成为自然语言处理和计算机视觉领域中一个非常优秀的基准数据集。2、
[网络安全自学篇] 二十三.基于机器学习的恶意请求识别及安全领域中的机器学习
杨秀璋的专栏
11-01 1万+
这是作者的系列网络安全自学教程，主要是关于网安工具和实践操作的在线笔记，特分享出来与博友共勉，希望您们喜欢，一起进步。前文分享了Web渗透的第一步工作，涉及网站信息、域名信息、端口信息、敏感信息及指纹信息收集。这篇文章换个口味，将分享机器学习在安全领域的应用，并复现一个基于机器学习（逻辑回归）的恶意请求识别。
详解基于python的全局与局部序列比对的实现(DNA)
12-16
程序能实现什么 a.完成gap值的自定义输入以及两条需比对序列的输入 b.完成得分矩阵的计算及输出 c.输出序列比对结果 d.使用matplotlib对得分矩阵路径的绘制 一、实现步骤 1.用户输入步骤 a.输入自定义的gap值 b.输入需要比对的碱基序列1（A,T,C,G）换行表示输入完成 b.输入需要比对的碱基序列2（A,T,C,G）换行表示输入完成 输入（示例）： 2.代码实现步骤 1.获取到用户输入的gap，s以及t 2.调用构建得分矩阵函数，得到得分矩阵以及方向矩阵 3.将得到的得分矩阵及方向矩阵作为参数传到回溯函数中开始回溯得到路径，路径存储使用的是全局变量，存的仍然是方向而
python序列化对比
xiaoxiaodechongzi的博客
04-08 220
要实现的功能 可以使用的api 将Python数据类型转换为（json）字符串 json.dumps() 将json字符串转换为Python数据类型 json.loads() 将Python数据类型以json形式保存到本地磁盘 json.dump() 将本地磁盘文件中的json数据转换为Python数据类型 json.load() 将Python数据类型转换为Python特定的二进制格式 pickle.dumps() 将Python特定的的二进制格式数据转换为Python...
用Python实现生信分析——序列比对（Sequence Alignment）详细讲解
最新发布
qq_41698317的博客
08-17 1212
局部比对（Local Alignment）是指在两条序列中寻找相似的局部区域进行比对。它允许我们在长度不相等或部分相似的序列中找到最相似的片段。局部比对尤其适用于在大多数序列差异较大但在某些部分具有高度相似性的情况下。
4、python序列对比
weixin_30627381的博客
01-12 213
1、列表：[ ]，杂（什么都可以放进去），有序，可变 2、元组：（），有序，不可变 3、字典：{ }，键值对，无序，可变 4、集合：{ }，不可重复，无序，可变 转载于:https://www.cnblogs.com/television/p/8276694.html...
Python中的序列比较
hou09tian的博客
09-30 3160
在Python中，具有相同数据类型的序列可以通过关系运算符进行比较。对序列进行比较大小，实际上就是对序列中的数据以“在字典中的顺序”（lexicographical ordering）进行比较，也就是出现在一本字典前面的数据要比后面的数据小。
python数据结构之序列和其他类型的比较
点点关注不迷路
05-20 354
python数据结构之序列和其他类型的比较序列和其他类型的比较 序列和其他类型的比较 序列对象可以与相同序列类型的其他对象比较。这种比较使用 字典式 顺序：首先，比较前两个对应元素，如果不相等，则可确定比较结果；如果相等，则比较之后的两个元素，以此类推，直到其中一个序列结束。如果要比较的两个元素本身是相同类型的序列，则递归地执行字典式顺序比较。如果两个序列中所有的对应元素都相等，则两个序列相等。如果一个序列是另一个的初始子序列，则较短的序列可被视为较小（较少）的序列。 对于字符串来说，字典式顺序使用 Uni
数据压缩实验5——JPEG编解码
weixin_45756474的博客
07-08 260
一.实验目的 掌握JPEG编解码系统的基本原理。初步掌握复杂的数据压缩算法实现，并能根据理论分析需要实现所对应数据的输出。 ** 二.JPEG编码过程 ** 2.1 RGB to YUV 为了减少各分量之间的相关性，减少数据的冗余，通常会把RGB颜色空间转换成YUV来进行各分量的编码。 2.2 Level Offset 零偏置 对于灰度级是 2^n 的像素，通过减去 2^n-1 ，将无符号的整数值变成有符号数。 例如：n=8，灰度级0 ~ 255，通过减去128，转化为-128 ~ 127 作用:使像素的绝
音视频-编解码-藏羚羊血红素氧合酶1基因编码区的克隆及序列分析.pdf
04-16
总的来说，尽管这个研究看似与音视频编解码技术无关，但其背后的科学方法和技术，如基因克隆、序列比对和系统进化分析，都是现代生物信息学的重要工具，而这些工具在某些方面与音视频数据处理中的编码和解码算法的...
【ESP 保姆级教程】疯狂毕设篇 —— 案例：基于ESP8266的WiFi自动打卡考勤系统
单片机菜鸟哥的博客
12-29 2274
现代社会上，基本上人手一部智能手机，自带wifi功能。只要我们打开了WiFi功能，我们就可以通过`自动捕获手机`发出的 802.11 帧 来获取到对应的手机MAC地址。当我们在后台服务器上预先配置好 MAC地址与用户信息的关联关系（比如`用户名字、用户工号、学生编号`等），并且把捕获到的MAC地址上传到后台服务器进行对比，我们就可以完成自动考勤或者无线点名功能。这整个过程都是无感知、全自动。
学学Python_23语句_条件语句4（==,序列比较,布尔运算）
10-21
文章：https://blog.csdn.net/yysyangyangyangshan/article/details/83246872
用Biopython批量比对序列（环境python）
weixin_48406854的博客
12-09 906
Biopython进行批量序列比对
Python Module之difflib-序列比较
Gavin
01-02 3664
Python的difflib模块
怎么用python实现序列比对_Biopython教程之“多序列比对”
weixin_39549852的博客
01-14 3241
原标题：Biopython教程之“多序列比对”多序列比对(Multiple Sequence Alignment, MSA)，对多个序列进行对位排列。这通常需要保证序列间的等同位点处在同一列上，并通过引进小横线(-)以保证最终的序列具有相同的长度。在生物信息分析中，我们有时需要进行多序列比对，Biopython可以帮我们实现，特别是使用linux系统的同学，biopython值得拥有。两种读取方法...
Seq2Seq模型中的序列解码策略
语言智能技术笔记簿
02-26 3822
导读：在序列生成类任务中，如机器翻译、自动摘要等，Seq2Seq是一种非常强大的模型。作为一种条件语言模型，它直接对P(y∣x)进行建模，在生成y的过程中，始终有x作为条件。当训练好了一个这样的模型后，在预测过程中，需要进行解码来找到最有可能的输出序列。这篇文章主要讲解Sequence to Sequence模型在预测阶段中使用的序列解码策略。 一些自然语言处理任务，如...
star 序列比对2020-12-25
生信小博士的博客
12-15 332
mkdir -p arab_STAR_genome #存储索引文件用 cd rnaseq/ STAR --runThreadN 10 --runMode genomeGenerate \ --genomeDir arab_STAR_genome \ --genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \ --sjdbGTFfile 00ref/Araport11_GFF3_genes_transposo...
生信软件 | STAR（测序序列与参考序列比对）
热门推荐
Baimoc
07-14 2万+
文章目录零、介绍一、安装二、使用1、建立索引2、STAR 比对三、原理聚类、拼接和评分 零、介绍 STAR （Spliced Transcripts Alignment to a Reference），用于将测序的 Read 对齐到参考基因组的比对软件，常用于 RNAseq。因其具有较高的准确率，映射速度较其他比对软件高 50 多倍，因此作为 ENCODE 项目的御用 pipeline 工具。不得不提的是，它需要占用大量内存，对计算资源有较高的要求。 一、安装 conda install -c biocon
C++实现动态规划进行蛋白质序列比对
资源摘要信息:"动态规划算法的c++实现能实现两序列的比对" 动态规划算法是一种算法设计技巧，它将一个复杂问题分解为更小的子问题，并保存这些子问题的解（通常称为子问题的解的备忘录），以避免重复计算。在计算机...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值