Skip to content

Instantly share code, notes, and snippets.

@lindenb
Created April 12, 2019 07:26
Show Gist options
  • Save lindenb/cca7f8242069022b6317d0c5de8ba4ce to your computer and use it in GitHub Desktop.
Save lindenb/cca7f8242069022b6317d0c5de8ba4ce to your computer and use it in GitHub Desktop.
C program finding discordant reads. htslib sam bam read sv structural variant
/**
Pierre Lindenbaum, 2019
gcc -o a.out -O3 -Wall -Isamtools/include discordant.c -Lsamtools/lib -lhts
*/
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <htslib/sam.h>
int main(int argc,char** argv) {
bam_hdr_t *h;
htsFile *in = NULL;
htsFile *out= NULL;
bam1_t *b = NULL;
int r;
if(!(argc==1 || argc==2)) {
fprintf(stderr,"Usage %s (in.bam|stdin)\n",argv[0]);
return EXIT_FAILURE;
}
in = hts_open(argc==1?"-":argv[1], "r");
b = bam_init1();
if(in==0) {
fprintf(stderr,"Cannot open %s. %s\n",argv[1],strerror(errno));
return EXIT_FAILURE;
}
h = sam_hdr_read(in);
b = bam_init1();
out = hts_open("-", "wb");
if(out==0) {
fprintf(stderr,"Cannot open output for writing. %s\n",strerror(errno));
return EXIT_FAILURE;
}
if(sam_hdr_write(out, h) !=0) {
fprintf(stderr,"Cannot open output for writing. %s\n",strerror(errno));
return EXIT_FAILURE;
}
while ((r = sam_read1(in, h, b)) >= 0) {
bam1_core_t* c=&(b->core);
int flag = c->flag;
if(!( flag & BAM_FPAIRED)) continue;
if( flag & (BAM_FUNMAP|BAM_FMUNMAP|BAM_FSECONDARY|BAM_FSUPPLEMENTARY|BAM_FDUP|BAM_FQCFAIL)) continue;
if( c->tid == c->mtid) continue;
if(strcmp("hs37d5",h->target_name[c->mtid])==0) continue;
if(strcmp("hs37d5",h->target_name[c->tid])==0) continue;
if( sam_write1(out, h, b) < 0 ) {
fprintf(stderr, "I/O error.\n");
return EXIT_FAILURE;
}
}
if (r < -1) {
fprintf(stderr, "Truncated file.\n");
return EXIT_FAILURE;
}
bam_destroy1(b);
hts_close(in);
hts_close(out);
bam_hdr_destroy(h);
return EXIT_SUCCESS;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment