C++使用htslib库读入和写出bam文件的实例

有时候我们需要使用C++处理bam文件，比如取出read1或者read2等符合特定条件的序列，根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等，这时我们可以使用htslib库。htslib可以用来处理SAM, BAM,CRAM 和VCF文件，是samtools、bcftools的核心库。

#include <stdio.h>
#include <stdlib.h>
#include <htslib/sam.h>

using namespace std; 

#define bam_is_read1(b) (((b)->core.flag&BAM_FREAD1) != 0)

uint8_t Base[16] = {0,65,67,0,71,0,0,0,84,0,0,0,0,0,0,78};

int main(int argc, char **argv)
{
 bam_hdr_t *header;
 bam1_t *aln = bam_init1();

 samFile *in = sam_open(argv[1], "r");
 htsFile *outR1 = hts_open(argv[2], "wb");
 header = sam_hdr_read(in);
 if (sam_hdr_write(outR1, header) < 0) {
 fprintf(stderr, "Error writing output.\n");
 exit(-1);
 }
 uint8_t *seq;
 int32_t lseq;
 uint32_t *cigar;
 char* qname;
 while (sam_read1(in, header, aln) >= 0) {
 if (bam_is_read1(aln)){
  sam_write1(outR1, header, aln);
 }
 else {
  seq = bam_get_seq(aln);
  lseq = aln->core.l_qseq;
  qname = bam_get_qname(aln);
  printf("%s\n",qname);
  cigar = bam_get_cigar(aln);
  for(int i=0; i < aln->core.n_cigar;++i){
  int icigar = cigar[i];
  printf("%d%d\n",bam_cigar_op(icigar),bam_cigar_oplen(icigar));
  }
  for(int i=0; i < lseq;++i){
  printf("%c", Base[bam_seqi(seq, i)]);
  }
  printf("\n");

 }
 }
 sam_close(in);
 sam_close(outR1);
}

cigar值存储形式

32位int，通过bam_get_cigar获得地址，aln->core.n_cigar返回cigar operation的个数

•低 4位存储cigar operation；通过函数bam_cigar_op()获得operation

•高28位存储cigar值的长度；通过函数,bam_cigar_oplen获得

seq存储形式

8位int，4位存储一个碱基，1,2,4,8,15分别代表A、C、G、T、N，高四位存储坐标数较小的碱基，可通过bam_seqi(seq,i)遍历。

以上这篇C++使用htslib库读入和写出bam文件的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持我们。

赞 (0)

相关推荐

  • C++使用htslib库读入和写出bam文件的实例

    C++使用htslib库读入和写出bam文件的实例

    有时候我们需要使用C++处理bam文件,比如取出read1或者read2等符合特定条件的序列,根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等,这时我们可以使用htslib库.htslib可以用来处理SAM, BAM,CRAM 和VCF文件,是samtools.bcftools的核心库. #include <stdio.h> #include <stdlib.h> #include <htslib/sam.h> using namespace st

  • LyScript实现计算片段Hash并写出Excel的示例代码

    本案例将学习运用LyScript计算特定程序中特定某些片段的Hash特征值,并通过xlsxwriter这个第三方模块将计算到的hash值存储成一个excel表格,本例中的知识点可以说已经具备了简单的表格输出能力,如果时间充裕完全可以实现自动化报告生成. 第一步实现计算特定片段的特征值,此类代码实现原理用户传入一个rva相对地址以及读入指令长度,并通过内置的hashlib库实现计算内存段内指令的特征,如下代码先来实现计算两段指令特征. import hashlib import zlib,bina

  • Javascript公共脚本库系列(一): 弹出层脚本

    Javascript公共脚本库系列(一): 弹出层脚本

    一.摘要 本系列文章是为了抽象通用的,跨浏览器的脚本方法. 本篇文章讲解弹出浮动层的javascript函数, 以及函数的原理和使用注意事项. 二.实现效果 用脚本弹出浮动层是我们最常用的脚本方法之一.下面是效果图:  点击图中的"航空公司"后,会在"航空公司"下面弹出浮动层. 在网上弹出框的脚本相当多, 而且还有各种第三方JS框架可供我们使用.但是其中有的脚本过于简单,仅仅粗略的实现弹出效果而忽略了灵活性,通用性和跨浏览器特性. 使用JS框架又有些杀鸡用牛刀.所以

  • 分享5个小技巧让你写出更好的 JavaScript 条件语句

    在使用 JavaScript 时,我们常常要写不少的条件语句.这里有五个小技巧,可以让你写出更干净.漂亮的条件语句. 1. 使用 Array.includes 来处理多重条件 举个栗子 : // 条件语句 function test(fruit) { if (fruit == 'apple' || fruit == 'strawberry') { console.log('red'); } } 乍一看,这么写似乎没什么大问题.然而,如果我们想要匹配更多的红色水果呢,比方说『樱桃』和『蔓越莓』?我

  • 使用Vue Composition API写出清晰、可扩展的表单实现

    表单是前端开发中最棘手的部分之一,您可能会在其中发现很多混乱的代码. 基于组件的框架,如 Vue.js,在提高前端代码的可扩展性方面做了很多工作,但是表单的问题仍然存在. 在本教程中,将向您展示新的 Vue Composition API(即将加入 Vue 3 中)如何使表单代码更清晰.更具可扩展性. 为什么表单代码经常很烂 像 Vue 这种基于组件的框架的关键设计模式是组件组合. 这种模式将应用程序的特性抽象为独立的.单一用途的组件,这些组件通信使用 props 和事件的方式. 然而,在此模式

  • 如何在面试中手写出javascript节流和防抖函数

    面试的时候我们经常会问别人是理解什么是节流和防抖,严格的可能要求你写出节流和防抖函数,这里我们抛开loadsh工具库手写节流和防抖 1.节流函数throttle // 节流方案1,每delay的时间执行一次,通过开关控制 function throttle(fn, delay, ctx) { let isAvail = true return function () { let args = arguments // 开关打开时,执行任务 if (isAvail) { fn.apply(ctx,

  • 如何写出安全的、基本功能完善的Bash脚本

    每个人或多或少总会碰到要使用并且自己完成编写一个最基础的Bash脚本的情况.真实情况是,没有人会说"哇哦,我喜欢写这些脚本".所以这也是为什么很少有人在写的时候专注在这些脚本上. 我本身也不是一个Bash脚本专家,但是我会在本文中跟你展示一个最基础最简单的安全脚本模板,会让你写的Bash脚本更加安全实用,你掌握了之后肯定会受益匪浅. 为什么要写Bash脚本 其实关于Bash脚本最好的解释如下: The opposite of "it's like riding a bike&

  • go语言中json数据的读取和写出操作

    go自带json库,在使用时需要通过 import "encoding/json"来导入该库. 在读取和写入json数据之前需要定义相关的结构体来对应被操作的json数据的格式,并且结构体中需要导出或导入的变量首字母大写. 其中,json.Marshal()用于将一个对象转换为json格式的字节数组,json.Unmarshal()用于将json格式的字节数组转换为一个对象. 具体使用示例如下所示: 首先,定义结构体: type Com struct { Name string Nod

  • 如何写出优雅的JS 代码

    变量 使用有意义和可发音的变量名 // 不好的写法 const yyyymmdstr = moment().format("YYYY/MM/DD"); // 好的写法 const currentDate = moment().format("YYYY/MM/DD"); 对同一类型的变量使用相同的词汇 // 不好的写法 getUserInfo(); getClientData(); getCustomerRecord(); // 好的写法 getUser(); 使用可

  • 浅谈JS如何写出漂亮的条件表达式

    多条件语句 多条件语句使用Array.includes 举个例子 function printAnimals(animal) { if (animal === "dog" || animal === "cat") { console.log(`I have a ${animal}`); } } console.log(printAnimals("dog")); // I have a dog 这种写法在条件比较少的情况下看起来没有问题,此时我们只

随机推荐