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

发布时间：2020-12-15 03:31:10 所属栏目：C语言来源：网络整理

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

有时候我们需要使用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,71,84,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,aln);
 }
 else {
  seq = bam_get_seq(aln);
  lseq = aln->core.l_qseq;
  qname = bam_get_qname(aln);
  printf("%sn",qname);
  cigar = bam_get_cigar(aln);
  for(int i=0; i < aln->core.n_cigar;++i){
  int icigar = cigar[i];
  printf("%d%dn",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文件的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持编程小技巧。

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!