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è)參考，也希望大家多多支持腳本之家。
相關(guān)文章
C語(yǔ)言浮點(diǎn)型數(shù)據(jù)在內(nèi)存中的存儲(chǔ)方式詳解
任何數(shù)據(jù)在內(nèi)存中都是以二進(jìn)制的形式存儲(chǔ)的,例如一個(gè)short型數(shù)據(jù)1156，其二進(jìn)制表示形式為00000100 10000100,下面這篇文章主要給大家介紹了關(guān)于C語(yǔ)言浮點(diǎn)型數(shù)據(jù)在內(nèi)存中的存儲(chǔ)方式,需要的朋友可以參考下2023-03-03
C++ 內(nèi)存分配處理函數(shù)set_new_handler的使用
這篇文章主要介紹了C++ 內(nèi)存分配處理函數(shù)set_new_handler的使用，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2020-02-02
QT已有項(xiàng)目導(dǎo)入工程時(shí)注意事項(xiàng)圖文詳解
QT開(kāi)發(fā)這幾年大大小小項(xiàng)目做了不少,花了點(diǎn)時(shí)間對(duì)知識(shí)點(diǎn)總結(jié)整合了一部分,下面這篇文章主要給大家介紹了關(guān)于QT已有項(xiàng)目導(dǎo)入工程時(shí)注意事項(xiàng)的相關(guān)資料,需要的朋友可以參考下2023-11-11
這篇文章主要介紹了C語(yǔ)言入門的一些基本資源推薦和程序語(yǔ)法概覽,C語(yǔ)言是很多現(xiàn)代高級(jí)編程語(yǔ)言的基礎(chǔ),需要的朋友可以參考下2015-12-12
這篇文章介紹了，用c語(yǔ)言實(shí)現(xiàn)的一個(gè)數(shù)組指針的小例子，有需要的朋友可以參考一下2013-07-07
這篇文章主要介紹了C語(yǔ)言中二叉樹(shù)的遍歷：前序、中序、后序，認(rèn)識(shí)二叉樹(shù)結(jié)構(gòu)最簡(jiǎn)單的方式就是遍歷二叉樹(shù)，感興趣的小伙伴可以參考閱讀本文2023-03-03
基于C++實(shí)現(xiàn)kinect+opencv 獲取深度及彩色數(shù)據(jù)
本文的主要思想是Kinect SDK 讀取彩色、深度、骨骼信息并用OpenCV顯示，非常的實(shí)用，有需要的小伙伴可以參考下2015-12-12
這篇文章主要介紹了C++ 獲取進(jìn)程CPU占用率,需要的朋友可以參考下2019-04-04
C語(yǔ)言中enum關(guān)鍵字的實(shí)現(xiàn)示例
這篇文章主要介紹了C語(yǔ)言中enum關(guān)鍵字的實(shí)現(xiàn)示例，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2021-03-03
C語(yǔ)言實(shí)現(xiàn)ATM機(jī)存取款系統(tǒng)
這篇文章主要為大家詳細(xì)介紹了C語(yǔ)言實(shí)現(xiàn)ATM機(jī)存取款系統(tǒng)，文中示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2021-11-11