雷火电竞-中国电竞赛事及体育赛事平台

歡迎來(lái)到入門(mén)教程網(wǎng)!

C語(yǔ)言

當(dāng)前位置:主頁(yè) > 軟件編程 > C語(yǔ)言 >

C++使用htslib庫(kù)讀入和寫(xiě)出bam文件的實(shí)例

來(lái)源:本站原創(chuàng)|時(shí)間:2020-01-10|欄目:C語(yǔ)言|點(diǎn)擊:

有時(shí)候我們需要使用C++處理bam文件,比如取出read1或者read2等符合特定條件的序列,根據(jù)cigar值對(duì)序列指定位置的堿基進(jìn)行統(tǒng)計(jì)或者對(duì)序列進(jìn)行處理并輸出等,這時(shí)我們可以使用htslib庫(kù)。htslib可以用來(lái)處理SAM, BAM,CRAM 和VCF文件,是samtools、bcftools的核心庫(kù)。

#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值存儲(chǔ)形式

32位int,通過(guò)bam_get_cigar獲得地址,aln->core.n_cigar返回cigar operation的個(gè)數(shù)

•低 4位存儲(chǔ)cigar operation;通過(guò)函數(shù)bam_cigar_op()獲得operation

•高28位存儲(chǔ)cigar值的長(zhǎng)度;通過(guò)函數(shù),bam_cigar_oplen獲得

seq存儲(chǔ)形式

8位int,4位存儲(chǔ)一個(gè)堿基,1,2,4,8,15分別代表A、C、G、T、N,高四位存儲(chǔ)坐標(biāo)數(shù)較小的堿基,可通過(guò)bam_seqi(seq,i)遍歷。

以上這篇C++使用htslib庫(kù)讀入和寫(xiě)出bam文件的實(shí)例就是小編分享給大家的全部?jī)?nèi)容了,希望能給大家一個(gè)參考,也希望大家多多支持我們。

上一篇:基于C++ bitset常用函數(shù)及運(yùn)算符(詳解)

欄    目:C語(yǔ)言

下一篇:c++用指針交換數(shù)組的實(shí)例講解

本文標(biāo)題:C++使用htslib庫(kù)讀入和寫(xiě)出bam文件的實(shí)例

本文地址:http://www.jygsgssxh.com/a1/Cyuyan/1040.html

網(wǎng)頁(yè)制作CMS教程網(wǎng)絡(luò)編程軟件編程腳本語(yǔ)言數(shù)據(jù)庫(kù)服務(wù)器

如果侵犯了您的權(quán)利,請(qǐng)與我們聯(lián)系,我們將在24小時(shí)內(nèi)進(jìn)行處理、任何非本站因素導(dǎo)致的法律后果,本站均不負(fù)任何責(zé)任。

聯(lián)系QQ:835971066 | 郵箱:835971066#qq.com(#換成@)

Copyright © 2002-2020 腳本教程網(wǎng) 版權(quán)所有