下面小编就为大家分享一篇C++使用htslib库读入和写出bam文件的实例，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
有时候我们需要使用C++处理bam文件，比如取出read1或者read2等符合特定条件的序列，根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等，这时我们可以使用htslib库。htslib可以用来处理SAM, BAM,CRAM 和VCF文件，是samtools、bcftools的核心库。
#include #include #include 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) { 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 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文件的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持gaodaima搞代码网。
以上就是C++使用htslib库读入和写出b来源gaodai#ma#com搞*代#码网am文件的实例的详细内容，更多请关注gaodaima搞代码网其它相关文章！
搞代码网（gaodaima.com）提供的所有资源部分来自互联网，如果有侵犯您的版权或其他权益，请说明详细缘由并提供版权或权益证明然后发送到邮箱[email protected]，我们会在看到邮件的第一时间内为您处理，或直接联系QQ：872152909。本网站采用BY-NC-SA协议进行授权
转载请注明原文链接：C++使用htslib库读入和写出bam文件的实例
转载请注明原文链接：C++使用htslib库读入和写出bam文件的实例
