C++使用htslib庫(kù)讀入和寫(xiě)出bam文件的實(shí)例
有時(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ǔ)言
本文標(biāo)題:C++使用htslib庫(kù)讀入和寫(xiě)出bam文件的實(shí)例
本文地址:http://www.jygsgssxh.com/a1/Cyuyan/1040.html
您可能感興趣的文章
- 04-02func函數(shù)+在C語(yǔ)言 func函數(shù)在c語(yǔ)言中
 - 04-02c語(yǔ)言沒(méi)有round函數(shù) round c語(yǔ)言
 - 01-10使用OpenGL實(shí)現(xiàn)3D立體顯示的程序代碼
 - 01-10深入理解C++中常見(jiàn)的關(guān)鍵字含義
 - 01-10使用C++實(shí)現(xiàn)全排列算法的方法詳解
 - 01-10深入Main函數(shù)中的參數(shù)argc,argv的使用詳解
 - 01-10c++中inline的用法分析
 - 01-10用C++實(shí)現(xiàn)DBSCAN聚類(lèi)算法
 - 01-10全排列算法的非遞歸實(shí)現(xiàn)與遞歸實(shí)現(xiàn)的方法(C++)
 - 01-10C++大數(shù)模板(推薦)
 


閱讀排行
- 1C語(yǔ)言 while語(yǔ)句的用法詳解
 - 2java 實(shí)現(xiàn)簡(jiǎn)單圣誕樹(shù)的示例代碼(圣誕
 - 3利用C語(yǔ)言實(shí)現(xiàn)“百馬百擔(dān)”問(wèn)題方法
 - 4C語(yǔ)言中計(jì)算正弦的相關(guān)函數(shù)總結(jié)
 - 5c語(yǔ)言計(jì)算三角形面積代碼
 - 6什么是 WSH(腳本宿主)的詳細(xì)解釋
 - 7C++ 中隨機(jī)函數(shù)random函數(shù)的使用方法
 - 8正則表達(dá)式匹配各種特殊字符
 - 9C語(yǔ)言十進(jìn)制轉(zhuǎn)二進(jìn)制代碼實(shí)例
 - 10C語(yǔ)言查找數(shù)組里數(shù)字重復(fù)次數(shù)的方法
 
本欄相關(guān)
- 04-02c語(yǔ)言函數(shù)調(diào)用后清空內(nèi)存 c語(yǔ)言調(diào)用
 - 04-02func函數(shù)+在C語(yǔ)言 func函數(shù)在c語(yǔ)言中
 - 04-02c語(yǔ)言的正則匹配函數(shù) c語(yǔ)言正則表達(dá)
 - 04-02c語(yǔ)言用函數(shù)寫(xiě)分段 用c語(yǔ)言表示分段
 - 04-02c語(yǔ)言中對(duì)數(shù)函數(shù)的表達(dá)式 c語(yǔ)言中對(duì)
 - 04-02c語(yǔ)言編寫(xiě)函數(shù)冒泡排序 c語(yǔ)言冒泡排
 - 04-02c語(yǔ)言沒(méi)有round函數(shù) round c語(yǔ)言
 - 04-02c語(yǔ)言分段函數(shù)怎么求 用c語(yǔ)言求分段
 - 04-02C語(yǔ)言中怎么打出三角函數(shù) c語(yǔ)言中怎
 - 04-02c語(yǔ)言調(diào)用函數(shù)求fibo C語(yǔ)言調(diào)用函數(shù)求
 
隨機(jī)閱讀
- 08-05DEDE織夢(mèng)data目錄下的sessions文件夾有什
 - 08-05織夢(mèng)dedecms什么時(shí)候用欄目交叉功能?
 - 01-10delphi制作wav文件的方法
 - 01-10SublimeText編譯C開(kāi)發(fā)環(huán)境設(shè)置
 - 04-02jquery與jsp,用jquery
 - 01-10使用C語(yǔ)言求解撲克牌的順子及n個(gè)骰子
 - 08-05dedecms(織夢(mèng))副欄目數(shù)量限制代碼修改
 - 01-11Mac OSX 打開(kāi)原生自帶讀寫(xiě)NTFS功能(圖文
 - 01-11ajax實(shí)現(xiàn)頁(yè)面的局部加載
 - 01-10C#中split用法實(shí)例總結(jié)
 


