From 2e0cfa5831217afbb8af9d537c46a10777e3abe2 Mon Sep 17 00:00:00 2001 From: Nick Croucher Date: Thu, 18 Apr 2024 13:01:15 +0100 Subject: [PATCH] Restore master version to check CI --- src/branch_sequences.c | 83 ++++++++++++++++-------------------------- 1 file changed, 31 insertions(+), 52 deletions(-) diff --git a/src/branch_sequences.c b/src/branch_sequences.c index 3cc09570..8ffc88aa 100644 --- a/src/branch_sequences.c +++ b/src/branch_sequences.c @@ -86,10 +86,16 @@ int get_list_of_snp_indices_which_fall_in_downstream_recombinations(int ** curre int current_index = 0; // convert the starting coordinates of block to the nearest SNP index current_index = find_starting_index(current_block_coordinates[0][i],snp_locations,0, number_of_snps); + + //make sure that the index begins at start of block + int beginning_j = current_index; + for(beginning_j = current_index; snp_locations[beginning_j] < current_block_coordinates[0][i];beginning_j++) + { + } - // starting at the begining index of block, count all the snps until the end of the block. int j; - for(j = current_index; (j < number_of_snps && snp_locations[j] <= current_block_coordinates[1][i]); j++) + // starting at the begining index of block, count all the snps until the end of the bock. + for(j = beginning_j; (j < number_of_snps && snp_locations[j] <= current_block_coordinates[1][i]); j++) { snps_in_recombinations[num_snps_in_recombinations] = j; num_snps_in_recombinations++; @@ -329,7 +335,7 @@ char *generate_branch_sequences(newick_node *root, FILE *vcf_file_pointer,int * newick_node * child_nodes[root->childNum]; root->taxon_names = (char *) calloc(MAX_SAMPLE_NAME_SIZE*number_of_columns,sizeof(char)); - // generate pointers for each child sequence + // generate pointers for each child seuqn while (child != NULL) { @@ -738,8 +744,6 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i // create the pileup of the snps and their sphere of influence int snp_counter = 0; - int min_postion = genome_size; - int max_position = 0; for(snp_counter = 0; snp_counter < number_of_branch_snps; snp_counter++) { int j = 0; @@ -756,11 +760,6 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i snp_sliding_window_counter = 0; } - if(snp_sliding_window_counter < min_postion) - { - min_postion = snp_sliding_window_counter; - } - // Upper bound of the window around a snp int max_snp_sliding_window_counter = snp_site_coords[snp_counter]+(window_size/2); max_snp_sliding_window_counter = extend_upper_part_of_window(snp_site_coords[snp_counter] + 1, @@ -772,11 +771,6 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i { max_snp_sliding_window_counter = genome_size; } - - if(max_snp_sliding_window_counter > max_position) - { - max_position= max_snp_sliding_window_counter; - } for(j = snp_sliding_window_counter; j < max_snp_sliding_window_counter; j++) { @@ -790,7 +784,7 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i int block_lower_bound = 0; // Scan across the pileup and record where blocks are above the cutoff int i; - for(i = min_postion; i <= max_position; i++) + for(i = 0; i <= genome_size; i++) { // Just entered the start of a block if(window_count[i] > cutoff && in_block == 0) @@ -800,23 +794,20 @@ int get_blocks(int ** block_coordinates, int genome_size,int * snp_site_coords,i } // Reached end of genome - if (in_block == 1) + if(i == genome_size && in_block == 1) { - if(i == genome_size) - { - block_coordinates[0][number_of_blocks] = block_lower_bound; - block_coordinates[1][number_of_blocks] = i; - number_of_blocks++; - in_block = 0; - } - // Just left a block - else if(window_count[i] <= cutoff) - { - block_coordinates[0][number_of_blocks] = block_lower_bound; - block_coordinates[1][number_of_blocks] = i-1; - number_of_blocks++; - in_block = 0; - } + block_coordinates[0][number_of_blocks] = block_lower_bound; + block_coordinates[1][number_of_blocks] = i; + number_of_blocks++; + in_block = 0; + } + // Just left a block + else if(window_count[i] <= cutoff && in_block == 1) + { + block_coordinates[0][number_of_blocks] = block_lower_bound; + block_coordinates[1][number_of_blocks] = i-1; + number_of_blocks++; + in_block = 0; } } @@ -1224,25 +1215,9 @@ double get_block_likelihood(int branch_genome_size, int number_of_branch_snps, i int calculate_genome_length_excluding_blocks_and_gaps(char * sequence, int length_of_sequence, int ** block_coordinates, int num_blocks) { + int genome_length = length_of_sequence; - // Process list of recombinations - int j = 0; - int * filtered_start_coords; - int * filtered_end_coords; - int num_filtered_blocks = 0; - filtered_start_coords = (int *) calloc(num_blocks,sizeof(int)); - filtered_end_coords = (int *) calloc(num_blocks,sizeof(int)); - for(j = 0; j= filtered_start_coords[j] && i <= filtered_end_coords[j]) + if(block_coordinates[0][j] != -1) { - genome_length--; + if (i >= block_coordinates[0][j] && i <= block_coordinates[1][j]) + { + genome_length--; + } } } }