有时候我们供给运用C++处理bam文件,比如取出read1还是read2等符合特定条件的队列,依照cigar值对队列钦点地点的碱基进行计算恐怕对队列实行拍卖并出口等,那时大家得以应用htslib库。htslib能够用来拍卖SAM,
BAM,CRAM 和VCF文件,是samtools、bcftools的宗旨库。

  有时候大家需求运用C++处理bam文件,比如取出read1依然read2等符合特定条件的队列,依照cigar值对队列钦定地方的碱基举行总计也许对队列实行拍卖并出口等,那时我们得以行使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);
}
#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值存款和储蓄格局

32位int,通过bam_get_cigar获得地点,aln->core.n_cigar返回cigar
operation的个数

  • 低 肆人存款和储蓄cigar operation;通过函数bam_cigar_op()获得operation
    图片 1
  • 高贰13个人存款和储蓄cigar值的长短;通过函数,bam_cigar_oplen获得

cigar值存款和储蓄情势

32位int,通过bam_get_cigar获得地点,aln->core.n_cigar返回cigar
operation的个数

  • 低 4人存款和储蓄cigar operation;通过函数bam_cigar_op()获得operation
    图片 2
  • 高贰十九个人存款和储蓄cigar值的尺寸;通过函数,bam_cigar_oplen获得

seq存款和储蓄形式

八个人int,3人存款和储蓄3个碱基,1,2,4,8,1四分级代表A、C、G、T、N,高4位存款和储蓄坐标数较小的碱基,可透过bam_seqi(seq,i)遍历。

seq存款和储蓄格局

7个人int,二个人存款和储蓄1个碱基,1,2,4,8,15分别代表A、C、G、T、N,高四个人存储坐标数较小的碱基,可由此bam_seqi(seq,i)遍历。

参考资料

htslib
sam.h文件:https://github.com/samtools/htslib/blob/develop/htslib/sam.h
htslib
sam文件格式表达:https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf

参考资料

htslib
sam.h文件:https://github.com/samtools/htslib/blob/develop/htslib/sam.h
htslib
sam文件格式表达:https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf

相关文章

网站地图xml地图