C++使用htslib庫讀入和寫出bam文件的實例

2020-05-23 13:37:05

字體：大中小

來源：轉載

供稿：網友

有時候我們需要使用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

htslib庫,讀入,寫出,bam文件,C++