<commit_before>#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <getopt.h>
#include <inttypes.h>

#include "bam.h"

#include "errors.h"
#include "time.h"
#include "memalloc.h"
#include "files.h"
#include "sequences.h"
#include "hashtable.h"

// This program is useful to generate information from a BAM file. The following
// information is provided as output:
// Following statistics about the alignment:
//
// QC-passed reads : 12677423
// QC-failed reads : 0
// Number of reads paired in sequencing: 0
// Number of single-fragments: 12677423
// 
// Duplicate reads : 3973834 (31.35%)
// Number of reads aligned: 6134567 (48.39%)
// Number of reads mapped: 4619967 (36.44%)
// 
// Number of read1's : 0
// Number of read2's : 0
// Number of reads properly paired: 0(nan)
// Number of reads with itself and mate aligned: 0
//
// QC-passed bases : 
// QC-failed bases : 0
// Number of bases paired in sequencing:
// Number of single-fragment bases     : 
// 
// Duplicate bases : 
// Number of bases aligned : 
// Number of bases mapped  : 
//
// Besides this we generate (if the user requests)
// a) quality value distribution with position on read.
// b) nucleotide composition for the reads. 
// c) coverage distribution of the reference  sequence.
// d) insert distance distribution.
// e) GC based coverage.

// print debug information?
int debug_flag = FALSE;

// the maximum quality value allowed
const uint maxqualityvalue = 60;

bool nowarnings;

typedef struct bamstats_t{
    // generation statistics
    uint64_t r_qcpassed;
    uint64_t r_qcfailed;
    uint64_t r_paired;
    uint64_t r_singletons;
    uint64_t b_qcpassed;
    uint64_t b_qcfailed;
    uint64_t b_paired;
    uint64_t b_singletons;
    
    // alignment statistics
    uint64_t r_duplicates;
    uint64_t r_aligned;
    uint64_t r_mapped;
    uint64_t b_duplicates;
    uint64_t b_aligned;
    uint64_t b_mapped;

    uint64_t read1s;
    uint64_t read2s;
    uint64_t properlypaired;
    uint64_t bothaligned;
}bamstats;

typedef struct chrcoverage_t
{
    uint64_t length;
    uchar* coverage;
    uchar* sequence;
} chrcoverage;

// return the number of bases in the read that are used in the alignment
static int lookatcigar(const bam1_t* const alignment)
{
    int numused = 0;

    uint32_t* cigar = bam1_cigar(alignment);

    uint i;
    for(i = 0; i < alignment->core.n_cigar; i++){
        switch(cigar[i] & 0xf){
            case BAM_CMATCH :
                numused += ((cigar[i] & 0xfffffff0) >> 4); break;
            case BAM_CINS:
                numused += ((cigar[i] & 0xfffffff0) >> 4); break;
            case BAM_CDEL:
                break;
            case BAM_CREF_SKIP:
                fatalf("Have not handled REF_SKIP in CIGAR");
            case BAM_CSOFT_CLIP:
                break;
            case BAM_CHARD_CLIP:
                break;
            case BAM_CPAD:
                break;
            default :
                fatalf("Unhandled CIGAR operation\n");
        }
    }

    return numused;
}

// print the statistics for this BAM file
static void print_stats(const bamstats* stats, FILE* fp)
{
    fprintf(fp, "Reads Summary\n");
    fprintf(fp,"-------------\n");
    fprintf(fp,"QC-passed reads : %"PRIu64"\n", stats->r_qcpassed);
    fprintf(fp,"QC-failed reads : %"PRIu64"\n", stats->r_qcfailed);
    fprintf(fp,"Number of single-fragments           : %"PRIu64"\n", stats->r_singletons);
    fprintf(fp,"Number of reads paired in sequencing : %"PRIu64"\n", stats->r_paired);
    fprintf(fp,"\n");
    fprintf(fp,"Duplicate reads         : %"PRIu64" (%2.2f%%)\n", 
            stats->r_duplicates, 
            stats->r_duplicates * 100.0 / stats->r_qcpassed);
    fprintf(fp,"Number of reads aligned : %"PRIu64" (%2.2f%%)\n", 
            stats->r_aligned,
            stats->r_aligned * 100.0 / stats->r_qcpassed);
    fprintf(fp,"Number of reads mapped  : %"PRIu64" (%2.2f%%)\n", 
            stats->r_mapped,
            stats->r_mapped * 100.0 / stats->r_qcpassed);
    fprintf(fp,"\n");
    fprintf(fp,"Number of read1's : %"PRIu64"\n", stats->read1s);
    fprintf(fp,"Number of read2's : %"PRIu64"\n", stats->read2s);
    fprintf(fp,"Number of reads properly paired              : %"PRIu64"\n", 
            stats->properlypaired);
    fprintf(fp,"Number of reads with itself and mate aligned : %"PRIu64"\n", 
            stats->bothaligned);
    fprintf(fp,"\n");
    fprintf(fp,"Bases Summary\n");
    fprintf(fp,"-------------\n");
    fprintf(fp,"QC-passed bases : %"PRIu64"\n", stats->b_qcpassed);
    fprintf(fp,"QC-failed bases : %"PRIu64"\n", stats->b_qcfailed);
    fprintf(fp,"Number of bases in single-fragments  : %"PRIu64"\n", stats->b_singletons);
    fprintf(fp,"Number of bases paired in sequencing : %"PRIu64"\n", stats->b_paired);
    fprintf(fp,"\n");
    fprintf(fp,"Duplicate bases         : %"PRIu64" (%2.2f%%)\n", 
            stats->b_duplicates, 
            stats->b_duplicates * 100.0 / stats->b_qcpassed);
    fprintf(fp,"Number of bases aligned : %"PRIu64" (%2.2f%%)\n", 
            stats->b_aligned,
            stats->b_aligned * 100.0 / stats->b_qcpassed);
    fprintf(fp,"Number of bases mapped  : %"PRIu64" (%2.2f%%)\n", 
            stats->b_mapped,
            stats->b_mapped * 100.0 / stats->b_qcpassed);
    fprintf(fp,"\n");
}

// print the read length distribution 
static void print_rlength(const char* const rlensname,
              const uint64_t* const rlens1,
              const uint64_t* const rlens2,
              const uint maxrlen,
              const bamstats* const stats)
{
    uint i;
    float perc;
    if(rlensname != NULL){
        FILE* fp = ckopen(rlensname, "w");

        for(i = 0; i < maxrlen; i++){
            perc = rlens1[i] / stats->read1s;
            fprintf(fp, "1 %d %0.3f\n", i + 1, perc);
        }
        for(i = 0; i < maxrlen; i++){
            if (stats->read2s != 0) {
                perc = rlens2[i] / stats->read2s;
                fprintf(fp, "2 %d %0.3f\n", i + 1, perc);
            }
        }

        fclose(fp);
    }
}

static void add_qualitycounts(const bam1_t* const alignment,
                              uint64_t** const qualvalues)
{
    uchar* qual = bam1_qual(alignment);
//    if((*qual) == '*'){
//        if(nowarnings == FALSE){
//            fprintf(stderr,"The quality values are not stored in this BAM for read %s\n", bam1_qname(alignment));
//        }
//        return;
//    }

    int i;
    for(i = 0;  i < alignment->core.l_qseq; i++){
        qualvalues[i][qual[i]>=maxqualityvalue ? maxqualityvalue-1 : qual[i]] += 1;
    }
}

void summary(const uint64_t* const frequencytable,    
             int* const pwmin, 
             int* const pbmin,
             int* const pmedian,
             int* const pbmax, <commit_msg>Readlength printing bug.<commit_after>#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <getopt.h>
#include <inttypes.h>

#include "bam.h"

#include "errors.h"
#include "time.h"
#include "memalloc.h"
#include "files.h"
#include "sequences.h"
#include "hashtable.h"

// This program is useful to generate information from a BAM file. The following
// information is provided as output:
// Following statistics about the alignment:
//
// QC-passed reads : 12677423
// QC-failed reads : 0
// Number of reads paired in sequencing: 0
// Number of single-fragments: 12677423
// 
// Duplicate reads : 3973834 (31.35%)
// Number of reads aligned: 6134567 (48.39%)
// Number of reads mapped: 4619967 (36.44%)
// 
// Number of read1's : 0
// Number of read2's : 0
// Number of reads properly paired: 0(nan)
// Number of reads with itself and mate aligned: 0
//
// QC-passed bases : 
// QC-failed bases : 0
// Number of bases paired in sequencing:
// Number of single-fragment bases     : 
// 
// Duplicate bases : 
// Number of bases aligned : 
// Number of bases mapped  : 
//
// Besides this we generate (if the user requests)
// a) quality value distribution with position on read.
// b) nucleotide composition for the reads. 
// c) coverage distribution of the reference  sequence.
// d) insert distance distribution.
// e) GC based coverage.

// print debug information?
int debug_flag = FALSE;

// the maximum quality value allowed
const uint maxqualityvalue = 60;

bool nowarnings;

typedef struct bamstats_t{
    // generation statistics
    uint64_t r_qcpassed;
    uint64_t r_qcfailed;
    uint64_t r_paired;
    uint64_t r_singletons;
    uint64_t b_qcpassed;
    uint64_t b_qcfailed;
    uint64_t b_paired;
    uint64_t b_singletons;
    
    // alignment statistics
    uint64_t r_duplicates;
    uint64_t r_aligned;
    uint64_t r_mapped;
    uint64_t b_duplicates;
    uint64_t b_aligned;
    uint64_t b_mapped;

    uint64_t read1s;
    uint64_t read2s;
    uint64_t properlypaired;
    uint64_t bothaligned;
}bamstats;

typedef struct chrcoverage_t
{
    uint64_t length;
    uchar* coverage;
    uchar* sequence;
} chrcoverage;

// return the number of bases in the read that are used in the alignment
static int lookatcigar(const bam1_t* const alignment)
{
    int numused = 0;

    uint32_t* cigar = bam1_cigar(alignment);

    uint i;
    for(i = 0; i < alignment->core.n_cigar; i++){
        switch(cigar[i] & 0xf){
            case BAM_CMATCH :
                numused += ((cigar[i] & 0xfffffff0) >> 4); break;
            case BAM_CINS:
                numused += ((cigar[i] & 0xfffffff0) >> 4); break;
            case BAM_CDEL:
                break;
            case BAM_CREF_SKIP:
                fatalf("Have not handled REF_SKIP in CIGAR");
            case BAM_CSOFT_CLIP:
                break;
            case BAM_CHARD_CLIP:
                break;
            case BAM_CPAD:
                break;
            default :
                fatalf("Unhandled CIGAR operation\n");
        }
    }

    return numused;
}

// print the statistics for this BAM file
static void print_stats(const bamstats* stats, FILE* fp)
{
    fprintf(fp, "Reads Summary\n");
    fprintf(fp,"-------------\n");
    fprintf(fp,"QC-passed reads : %"PRIu64"\n", stats->r_qcpassed);
    fprintf(fp,"QC-failed reads : %"PRIu64"\n", stats->r_qcfailed);
    fprintf(fp,"Number of single-fragments           : %"PRIu64"\n", stats->r_singletons);
    fprintf(fp,"Number of reads paired in sequencing : %"PRIu64"\n", stats->r_paired);
    fprintf(fp,"\n");
    fprintf(fp,"Duplicate reads         : %"PRIu64" (%2.2f%%)\n", 
            stats->r_duplicates, 
            stats->r_duplicates * 100.0 / stats->r_qcpassed);
    fprintf(fp,"Number of reads aligned : %"PRIu64" (%2.2f%%)\n", 
            stats->r_aligned,
            stats->r_aligned * 100.0 / stats->r_qcpassed);
    fprintf(fp,"Number of reads mapped  : %"PRIu64" (%2.2f%%)\n", 
            stats->r_mapped,
            stats->r_mapped * 100.0 / stats->r_qcpassed);
    fprintf(fp,"\n");
    fprintf(fp,"Number of read1's : %"PRIu64"\n", stats->read1s);
    fprintf(fp,"Number of read2's : %"PRIu64"\n", stats->read2s);
    fprintf(fp,"Number of reads properly paired              : %"PRIu64"\n", 
            stats->properlypaired);
    fprintf(fp,"Number of reads with itself and mate aligned : %"PRIu64"\n", 
            stats->bothaligned);
    fprintf(fp,"\n");
    fprintf(fp,"Bases Summary\n");
    fprintf(fp,"-------------\n");
    fprintf(fp,"QC-passed bases : %"PRIu64"\n", stats->b_qcpassed);
    fprintf(fp,"QC-failed bases : %"PRIu64"\n", stats->b_qcfailed);
    fprintf(fp,"Number of bases in single-fragments  : %"PRIu64"\n", stats->b_singletons);
    fprintf(fp,"Number of bases paired in sequencing : %"PRIu64"\n", stats->b_paired);
    fprintf(fp,"\n");
    fprintf(fp,"Duplicate bases         : %"PRIu64" (%2.2f%%)\n", 
            stats->b_duplicates, 
            stats->b_duplicates * 100.0 / stats->b_qcpassed);
    fprintf(fp,"Number of bases aligned : %"PRIu64" (%2.2f%%)\n", 
            stats->b_aligned,
            stats->b_aligned * 100.0 / stats->b_qcpassed);
    fprintf(fp,"Number of bases mapped  : %"PRIu64" (%2.2f%%)\n", 
            stats->b_mapped,
            stats->b_mapped * 100.0 / stats->b_qcpassed);
    fprintf(fp,"\n");
}

// print the read length distribution 
static void print_rlength(const char* const rlensname,
              const uint64_t* const rlens1,
              const uint64_t* const rlens2,
              const uint maxrlen,
              const bamstats* const stats)
{
    uint i;
    float perc;
    if(rlensname != NULL){
        FILE* fp = ckopen(rlensname, "w");

        for(i = 0; i < maxrlen; i++){
            perc = rlens1[i] * 1.0 / stats->read1s;
            fprintf(fp, "1 %d %0.3f\n", i + 1, perc);
        }
        for(i = 0; i < maxrlen; i++){
            if (stats->read2s != 0) {
                perc = rlens2[i] * 1.0 / stats->read2s;
                fprintf(fp, "2 %d %0.3f\n", i + 1, perc);
            }
        }

        fclose(fp);
    }
}

static void add_qualitycounts(const bam1_t* const alignment,
                              uint64_t** const qualvalues)
{
    uchar* qual = bam1_qual(alignment);
//    if((*qual) == '*'){
//        if(nowarnings == FALSE){
//            fprintf(stderr,"The quality values are not stored in this BAM for read %s\n", bam1_qname(alignment));
//        }
//        return;
//    }

    int i;
    for(i = 0;  i < alignment->core.l_qseq; i++){
        qualvalues[i][qual[i]>=maxqualityvalue ? maxqualityvalue-1 : qual[i]] += 1;
    }
}

void summary(const uint64_t* const frequencytable,    
             int* const pwmin, 
             int* const pbmin,
             int* const pmedian,
             int* const pbmax, 
{'fraction_non_alphanumeric': 0.13836431773091337, 'fraction_numerical': 0.0327724300095525, 'mean_word_length': 2.3082158483228, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 93, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6366751', 'n_tokens_mistral': 5005, 'n_tokens_neox': 4636, 'n_words': 1186}