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

最新推荐文章于 2025-03-06 06:00:17 发布

木子木木夕然

最新推荐文章于 2025-03-06 06:00:17 发布

阅读量757

点赞数

分类专栏：生信工具文章标签：生物信息学

本文链接：https://blog.csdn.net/weixin_45963508/article/details/120865605

版权

生信工具专栏收录该内容

14 篇文章

订阅专栏

该程序读取BAM文件，并根据qname（读取名称）进行处理。它检查是否已记录相同的qname，如果已存在，则根据先前随机分配的标记将read写入两个不同的输出BAM文件之一。新出现的qname随机分配到输出文件，并记录在链表中。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

#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()%10;//生成随机数 
    	    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);
	
}