C++使用htslib庫讀入和寫出bam文件的實(shí)例
有時(shí)候我們需要使用C++處理bam文件，比如取出read1或者read2等符合特定條件的序列，根據(jù)cigar值對(duì)序列指定位置的堿基進(jìn)行統(tǒng)計(jì)或者對(duì)序列進(jìn)行處理并輸出等，這時(shí)我們可以使用htslib庫。htslib可以用來處理SAM, BAM,CRAM 和VCF文件，是samtools、bcftools的核心庫。
#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，通過bam_get_cigar獲得地址，aln->core.n_cigar返回cigar operation的個(gè)數(shù)
•低 4位存儲(chǔ)cigar operation；通過函數(shù)bam_cigar_op()獲得operation
•高28位存儲(chǔ)cigar值的長度；通過函數(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ù)較小的堿基，可通過bam_seqi(seq,i)遍歷。
以上這篇C++使用htslib庫讀入和寫出bam文件的實(shí)例就是小編分享給大家的全部內(nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持腳本之家。
相關(guān)文章
在 C語言中，EOF 是一個(gè)宏定義，EOF 常常用于文件的輸入輸出中，當(dāng)讀取到文件結(jié)束時(shí)，會(huì)返回 EOF，本文就詳細(xì)的介紹一下具體使用方法，感興趣的可以一起來了解一下2023-05-05
Matlab實(shí)現(xiàn)三維投影繪制的示例代碼
這篇文章系小編為大家?guī)砹艘粋€(gè)三維投影繪制函數(shù)(三視圖繪制)，函數(shù)支持三維曲線、曲面、三維多邊形、參數(shù)方程曲線、參數(shù)方程曲面的投影繪制，需要的可以參考一下2022-08-08
迭代器（iterator）是一種用于遍歷數(shù)據(jù)集合的的對(duì)象,它提供了一種訪問數(shù)據(jù)集合中元素的方式,而無需暴露數(shù)據(jù)集合內(nèi)部的細(xì)節(jié),使用迭代器,我們可以對(duì)數(shù)據(jù)集合中的每個(gè)元素進(jìn)行處理,本文介紹了C++中關(guān)于vector迭代器失效問題的原因及解決方案,需要的朋友可以參考下2024-09-09
Matlab實(shí)現(xiàn)灰色預(yù)測的示例代碼
這篇文章主要為大家詳細(xì)介紹了如何利用Matlab實(shí)現(xiàn)灰色預(yù)測，文中的示例代碼講解詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴可以了解一下2022-05-05
本文給大家分享的是一則使用C語言實(shí)現(xiàn)動(dòng)態(tài)數(shù)組的代碼，完美解決內(nèi)存溢出以及內(nèi)存回收問題，有需要的小伙伴可以參考下2021-09-09
本文主要介紹了VS2017中配置QT5.12.0的圖文教程，文中通過示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2023-06-06
linux c程序中獲取shell腳本輸出的實(shí)現(xiàn)方法
以下是對(duì)在linux下c程序中獲取shell腳本輸出的實(shí)現(xiàn)方法進(jìn)行了詳細(xì)的分析介紹，需要的朋友可以過來參考下2013-08-08
C++找出字符串中出現(xiàn)最多的字符和次數(shù),時(shí)間復(fù)雜度小于O(n^2)
今天小編就為大家分享一篇關(guān)于C++找出字符串中出現(xiàn)最多的字符和次數(shù),時(shí)間復(fù)雜度小于O(n^2)，小編覺得內(nèi)容挺不錯(cuò)的，現(xiàn)在分享給大家，具有很好的參考價(jià)值，需要的朋友一起跟隨小編來看看吧2018-12-12
C語言中find_package()的搜索路徑的實(shí)現(xiàn)
本文主要介紹了C語言中find_package()的搜索路徑的實(shí)現(xiàn)，文中通過示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2021-12-12
這篇文章主要介紹了C++14中binary literals的使用,本文通過實(shí)例代碼給大家介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或工作具有一定的參考借鑒價(jià)值，需要的朋友可以參考下2023-06-06