<commit_before>/**   LICENSE
* Copyright (c) 2016 Genome Research Ltd.
*
* Author: Cancer Genome Project <EMAIL>
*
* This file is part of cgpBigWig.
*
* cgpBigWig is free software: you can redistribute it and/or modify it under
* the terms of the GNU Affero General Public License as published by the Free
* Software Foundation; either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
* details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*    1. The usage of a range of years within a copyright statement contained within
*    this distribution should be interpreted as being equivalent to a list of years
*    including the first and last year specified and all consecutive years between
*    them. For example, a copyright statement that reads ‘Copyright (c) 2005, 2007-
*    2009, 2011-2012’ should be interpreted as being identical to a statement that
*    reads ‘Copyright (c) 2005, 2007, 2008, 2009, 2011, 2012’ and a copyright
*    statement that reads ‘Copyright (c) 2005-2012’ should be interpreted as being
*    identical to a statement that reads ‘Copyright (c) 2005, 2006, 2007, 2008,
*    2009, 2010, 2011, 2012’."
*
*/

#include <getopt.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <inttypes.h>
#include <libgen.h>
#include "bam_access.h"
#include "utils.h"

char *out_file = NULL;
char *input_file = NULL;
char *region_store = NULL;
char **our_region_list = NULL;
char *reference = NULL;
uint32_t region_list_count = 0;
int is_regions_file = 0;
uint8_t is_base = 0;
int filter = 4;
char base = 0;
uint32_t single = 1;

void print_usage (int exit_code){
	printf("Usage: bam2bwbases -i input.[b|cr]am -o output.bw\n");
	printf("bam2bwbases can be used to generate four bw files of per base proportions.\n\n");
	printf("-i  --input [file]                                Path to the input [b|cr]am file.\n");
	printf("-F  --filter [int]                                SAM flags to filter. [default: %d]\n",filter);
	printf("-o  --outfile [file]                              Path to the output .bw file produced. Per base results wiillbe output as four bw files [ACGT].outputname.bw [default:'%s']\n\n",out_file);
	printf("Optional: \n");
	printf("-c  --region [file]                               A samtools style region (contig:start-stop) or a bed file of regions over which to produce the bigwig file\n");
	printf("-r  --reference [file]                            Path to reference genome.fa file (required for cram if ref_path cannot be resolved)\n\n");
  printf ("Other:\n");
	printf("-h  --help                                        Display this usage information.\n");
	printf("-v  --version                                     Prints the version number.\n\n");

  exit(exit_code);
}

int check_region_string(char *region){
  if(check_exist(region)==1){
    is_regions_file = 1;
    return 1;
  }else{
    //Parse this as a region
    int beg = 0;
    int end = 0;
    const char *res = hts_parse_reg(region, &beg, &end);
    if(res) return 1;
  }
  return -1;
}

int get_int_length(int input){
  return (input == 0 ? 1 : (int)(log10(input)+1));
}

void setup_options(int argc, char *argv[]){
	const struct option long_opts[] =
	{
             	{"input", required_argument, 0, 'i'},
             	{"filter", required_argument, 0, 'F'},
             	{"outfile",required_argument, 0, 'o'},
             	{"region",required_argument, 0, 'c'},
             	{"reference",required_argument, 0, 'r'},
             	{"help", no_argument, 0, 'h'},
             	{"version", no_argument, 0, 'v'},
             	{ NULL, 0, NULL, 0}

   }; //End of declaring opts

   int index = 0;
   int iarg = 0;

   //Iterate through options
   while((iarg = getopt_long(argc, argv, "F:i:o:c:r:hv",long_opts, &index)) != -1){
    switch(iarg){
      case 'F':
        if(sscanf(optarg, "%i", &filter) != 1){
      		fprintf(stderr,"Error parsing -F|--filter argument '%s'. Should be an integer > 0",optarg);
      		print_usage(1);
      	}
        break;
   		case 'i':
				input_file = optarg;
				if(check_exist(input_file) != 1){
          fprintf(stderr,"Input bam file %s does not appear to exist.\n",input_file);
          print_usage(1);
        }
   			break;
   		case 'o':
				out_file = optarg;
   			break;
			case 'h':
				print_usage (0);
				break;
			case 'v':
				print_version (0);
				break;
			case 'c':
			  region_store = optarg;
			  //First check for a region format
			  int res = check_region_string(region_store);
        if(res<0){
          fprintf(stderr,"Region %s is not in correct format or not an existing bed file.\n",region_store);
          print_usage(1);
        }
			  break;
			case 'r':
			  reference = optarg;
			  break;
			case '?':
        print_usage (1);
        break;
      default:
      	print_usage (0);
   	}; // End of args switch statement

   }//End of iteration through options

  if(input_file==NULL){
    fprintf(stderr,"Required option -i|--input-bam not defined.\n");
    print_usage(1);
  }

  if(out_file == NULL){
    out_file = "output.bam.bw";
  }
  return;
}

// callback for bam_plbuf_init()
static int perbase_pileup_func(uint32_t tid, uint32_t position, int n, const bam_pileup1_t *pl, void *data){
  tmpstruct_t *tmp = (tmpstruct_t*)data;
  int pos              = (int)position;
  int coverage         = n;
  int base_coverage    = 0;
  int i;
  for (i=0;i<n;i++){
    if (pl[i].is_del){
      coverage--;
    }else{
      if(bam_seqi(bam_get_seq(pl[i].b), pl[i].qpos) == tmp->base_bit) base_coverage++;
    }
  }
  float result = 0;
  if(base_coverage>0) result = (float)base_coverage / (float) coverage;
  if(tmp->ltid != tid || tmp->lbaseprop != result || pos > tmp->lpos+1){
    //if(tmp->lpos > 0){
      uint32_t start =  tmp->lstart;
      uint32_t stop = (tmp->lpos +1);
      float res = tmp->lbaseprop;
      if(start >= tmp->reg_start-1 && stop <= tmp->reg_stop){
        int chck = bwAddIntervals(tmp->bwout,&tmp->head->target_name[tmp->ltid],&start,&stop,&res,single);
        if(chck!=0){
          fprintf(stderr,"Error adding region to bw '%s:%"PRIu32"-%"PRIu32"\t%f'. Error code: %d.\n",tmp->head->target_name[tmp->ltid],start,stop,res,chck);
          exit(1);
        }
      }
    //}
    tmp->ltid          = tid;
    tmp->lstart        = pos;
    tmp->lbaseprop     = result;
  }
  tmp->lpos = pos;
  return 0;
}

bigWigFile_t *initialise_bw_output(char *out_file, chromList_t *chromList){
  //Open output file
  bigWigFile_t *fp = bwOpen(out_file, NULL, "w");
  check(fp, "Error opening %s for writing.",out_file);
  //10 zoom levels
  int res = bwCreateHdr(fp, 10);
  check(!res,"Received error %d in bwCreateHeader for %s",res,out_file);
  //Create the chromosome lists
  fp->cl = bwCreateChromList(chromList->chrom, chromList->len, chromList->nKeys);
  check(fp->cl != NULL, "Error generating chromlist as bw for %s.",out_file);<commit_msg>Correction to  c/bam2bwbases.c positions<commit_after>/**   LICENSE
* Copyright (c) 2016 Genome Research Ltd.
*
* Author: Cancer Genome Project <EMAIL>
*
* This file is part of cgpBigWig.
*
* cgpBigWig is free software: you can redistribute it and/or modify it under
* the terms of the GNU Affero General Public License as published by the Free
* Software Foundation; either version 3 of the License, or (at your option) any
* later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
* details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*    1. The usage of a range of years within a copyright statement contained within
*    this distribution should be interpreted as being equivalent to a list of years
*    including the first and last year specified and all consecutive years between
*    them. For example, a copyright statement that reads ‘Copyright (c) 2005, 2007-
*    2009, 2011-2012’ should be interpreted as being identical to a statement that
*    reads ‘Copyright (c) 2005, 2007, 2008, 2009, 2011, 2012’ and a copyright
*    statement that reads ‘Copyright (c) 2005-2012’ should be interpreted as being
*    identical to a statement that reads ‘Copyright (c) 2005, 2006, 2007, 2008,
*    2009, 2010, 2011, 2012’."
*
*/

#include <getopt.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <inttypes.h>
#include <libgen.h>
#include "bam_access.h"
#include "utils.h"

char *out_file = NULL;
char *input_file = NULL;
char *region_store = NULL;
char **our_region_list = NULL;
char *reference = NULL;
uint32_t region_list_count = 0;
int is_regions_file = 0;
uint8_t is_base = 0;
int filter = 4;
char base = 0;
uint32_t single = 1;

void print_usage (int exit_code){
	printf("Usage: bam2bwbases -i input.[b|cr]am -o output.bw\n");
	printf("bam2bwbases can be used to generate four bw files of per base proportions.\n\n");
	printf("-i  --input [file]                                Path to the input [b|cr]am file.\n");
	printf("-F  --filter [int]                                SAM flags to filter. [default: %d]\n",filter);
	printf("-o  --outfile [file]                              Path to the output .bw file produced. Per base results wiillbe output as four bw files [ACGT].outputname.bw [default:'%s']\n\n",out_file);
	printf("Optional: \n");
	printf("-c  --region [file]                               A samtools style region (contig:start-stop) or a bed file of regions over which to produce the bigwig file\n");
	printf("-r  --reference [file]                            Path to reference genome.fa file (required for cram if ref_path cannot be resolved)\n\n");
  printf ("Other:\n");
	printf("-h  --help                                        Display this usage information.\n");
	printf("-v  --version                                     Prints the version number.\n\n");

  exit(exit_code);
}

int check_region_string(char *region){
  if(check_exist(region)==1){
    is_regions_file = 1;
    return 1;
  }else{
    //Parse this as a region
    int beg = 0;
    int end = 0;
    const char *res = hts_parse_reg(region, &beg, &end);
    if(res) return 1;
  }
  return -1;
}

int get_int_length(int input){
  return (input == 0 ? 1 : (int)(log10(input)+1));
}

void setup_options(int argc, char *argv[]){
	const struct option long_opts[] =
	{
             	{"input", required_argument, 0, 'i'},
             	{"filter", required_argument, 0, 'F'},
             	{"outfile",required_argument, 0, 'o'},
             	{"region",required_argument, 0, 'c'},
             	{"reference",required_argument, 0, 'r'},
             	{"help", no_argument, 0, 'h'},
             	{"version", no_argument, 0, 'v'},
             	{ NULL, 0, NULL, 0}

   }; //End of declaring opts

   int index = 0;
   int iarg = 0;

   //Iterate through options
   while((iarg = getopt_long(argc, argv, "F:i:o:c:r:hv",long_opts, &index)) != -1){
    switch(iarg){
      case 'F':
        if(sscanf(optarg, "%i", &filter) != 1){
      		fprintf(stderr,"Error parsing -F|--filter argument '%s'. Should be an integer > 0",optarg);
      		print_usage(1);
      	}
        break;
   		case 'i':
				input_file = optarg;
				if(check_exist(input_file) != 1){
          fprintf(stderr,"Input bam file %s does not appear to exist.\n",input_file);
          print_usage(1);
        }
   			break;
   		case 'o':
				out_file = optarg;
   			break;
			case 'h':
				print_usage (0);
				break;
			case 'v':
				print_version (0);
				break;
			case 'c':
			  region_store = optarg;
			  //First check for a region format
			  int res = check_region_string(region_store);
        if(res<0){
          fprintf(stderr,"Region %s is not in correct format or not an existing bed file.\n",region_store);
          print_usage(1);
        }
			  break;
			case 'r':
			  reference = optarg;
			  break;
			case '?':
        print_usage (1);
        break;
      default:
      	print_usage (0);
   	}; // End of args switch statement

   }//End of iteration through options

  if(input_file==NULL){
    fprintf(stderr,"Required option -i|--input-bam not defined.\n");
    print_usage(1);
  }

  if(out_file == NULL){
    out_file = "output.bam.bw";
  }
  return;
}

// callback for bam_plbuf_init()
static int perbase_pileup_func(uint32_t tid, uint32_t position, int n, const bam_pileup1_t *pl, void *data){
  tmpstruct_t *tmp = (tmpstruct_t*)data;
  int pos              = (int)position;
  int coverage         = n;
  int base_coverage    = 0;
  int i;
  for (i=0;i<n;i++){
    if (pl[i].is_del){
      coverage--;
    }else{
      if(bam_seqi(bam_get_seq(pl[i].b), pl[i].qpos) == tmp->base_bit) base_coverage++;
    }
  }
  float result = 0;
  if(base_coverage>0) result = (float)base_coverage / (float) coverage;
  if(tmp->ltid != tid || tmp->lbaseprop != result || pos > tmp->lpos+1){
    //if(tmp->lpos > 0){
      uint32_t start =  tmp->lstart;
      uint32_t stop = pos;
      float res = tmp->lbaseprop;
      if(start >= tmp->reg_start-1 && stop <= tmp->reg_stop){
        int chck = bwAddIntervals(tmp->bwout,&tmp->head->target_name[tmp->ltid],&start,&stop,&res,single);
        if(chck!=0){
          fprintf(stderr,"Error adding region to bw '%s:%"PRIu32"-%"PRIu32"\t%f'. Error code: %d.\n",tmp->head->target_name[tmp->ltid],start,stop,res,chck);
          exit(1);
        }
      }
    //}
    tmp->ltid          = tid;
    tmp->lstart        = pos;
    tmp->lbaseprop     = result;
  }
  tmp->lpos = pos;
  return 0;
}

bigWigFile_t *initialise_bw_output(char *out_file, chromList_t *chromList){
  //Open output file
  bigWigFile_t *fp = bwOpen(out_file, NULL, "w");
  check(fp, "Error opening %s for writing.",out_file);
  //10 zoom levels
  int res = bwCreateHdr(fp, 10);
  check(!res,"Received error %d in bwCreateHeader for %s",res,out_file);
  //Create the chromosome lists
  fp->cl = bwCreateChromList(chromList->chrom, chromList->len, chromList->nKeys);
  check(fp->cl != NULL, "Error generating chromlist as bw for %s.",out_file);
{'fraction_non_alphanumeric': 0.1344589937973811, 'fraction_numerical': 0.023018607856650587, 'mean_word_length': 2.7160051216389243, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '455676', 'n_tokens_mistral': 5183, 'n_tokens_neox': 4794, 'n_words': 1489}