C program finding discordant reads. htslib sam bam read sv structural variant

discordant.c

/**
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;
    }