C++使用htslib库读入和写出bam文件的实例-创新互联
有时候我们需要使用C++处理bam文件,比如取出read1或者read2等符合特定条件的序列,根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等,这时我们可以使用htslib库。htslib可以用来处理SAM, BAM,CRAM 和VCF文件,是samtools、bcftools的核心库。
成都创新互联公司2013年成立,是专业互联网技术服务公司,拥有项目做网站、成都网站制作网站策划,项目实施与项目整合能力。我们以让每一个梦想脱颖而出为使命,1280元铁西做网站,已为上家服务,为铁西各地企业和个人服务,联系电话:18982081108#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) { 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); }
另外有需要云服务器可以了解下创新互联建站www.cdcxhl.com,海内外云服务器15元起步,三天无理由+7*72小时售后在线,公司持有idc许可证,提供“云服务器、裸金属服务器、高防服务器、香港服务器、美国服务器、虚拟主机、免备案服务器”等云主机租用服务以及企业上云的综合解决方案,具有“安全稳定、简单易用、服务可用性高、性价比高”等特点与优势,专为企业上云打造定制,能够满足用户丰富、多元化的应用场景需求。
分享题目:C++使用htslib库读入和写出bam文件的实例-创新互联
文章来源:http://myzitong.com/article/cdshsc.html