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);

									}