将bam文件分成两个bam文件,属于同一个pe的两条reads分到一个bam文件中

#include <stdio.h>
#include <stdlib.h>
#include<time.h>
#include "htslib/sam.h"

#define bam_is_read1(b) (((b)->core.flag&BAM_FREAD1) != 0)

typedef struct Link{//用于存储qname信息
    char *qname;
    int tie;
    struct Link *next;
}Link,*LinkList;

int main(int argc,char **argv){
    samFile *in = sam_open(argv[1], "r");
    bam_hdr_t *header;
    bam1_t *aln = bam_init1();
    header = sam_hdr_read(in);  
    htsFile *outR1 = hts_open(argv[2], "wb");
    htsFile *outR2 = hts_open(argv[3], "wb");
    if (sam_hdr_write(outR1, header) < 0) {
        fprintf(stderr, "Error writing output.\n");
        exit(-1);
    }
    if (sam_hdr_write(outR2, header) < 0) {
        fprintf(stderr, "Error writing output.\n");
        exit(-1);
    }   
    Link *head = (Link *)malloc(sizeof(Link));  
    head->next = NULL; 
    Link *temp = head;
    Link *e = head;
    char *name;  
    srand((unsigned int)time(NULL));//time函数返回的时间戳作为随机数生成器的初始值
    while(sam_read1(in, header, aln) >= 0){
        
        Link *t = head;
        int sign = 0;
        name = bam_get_qname(aln);
        while(t->next != NULL){//查询是否已经录入相同qname的read 
            if(strcmp(name,t->next->qname) == 0){
                if(t->next->tie == 1){
                    int a=sam_write1(outR1, header, aln);//写入数据
                    sign = 1;
                    break;  
                }else{
                    int a=sam_write1(outR2, header, aln);
                    sign = 1;
                    break;
                }
            }
            else{
                t = t->next;
            }   
        }
        if(sign == 0){//将新出现的qname录入 
            Link *b = (Link *)malloc(sizeof(Link)); 
            b->qname =(char *)malloc(sizeof(char)*1000);
            b->qname = strcpy(b->qname,bam_get_qname(aln));         
            int i = rand(); //生成一个随机数
            if(i%2!=0){//存入一个bam文件中
                int a=sam_write1(outR1, header, aln);
                b->tie = 1; 
            }
            else{
                int a=sam_write1(outR2, header, aln);
                b->tie = 0;
            }
            b->next = NULL;
            temp->next = b;
            temp = temp->next;  
        }
    }
    
/*  int n = 0;//统计录入的qname数量 
    while(e->next != NULL){ 
        n++;
        e = e->next;
    }
    printf("%d\n",n);
*/
    sam_close(in);
    sam_close(outR1);
    sam_close(outR2);
    
}
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容