minor fix in bam generation

grev-uv · Nov 29, 2018 · 4df3838 · 4df3838
1 parent 4f2e9a4
commit 4df3838
Show file tree

Hide file tree

Showing 8 changed files with 120 additions and 6 deletions.
diff --git a/lib/c/src/bioformats/bam/alignment.c b/lib/c/src/bioformats/bam/alignment.c
@@ -535,8 +535,11 @@ bam1_t* convert_to_bam(alignment_t* alignment_p, int base_quality) {
     char* optional_fields;
     size_t optional_fields_length;
 
+
     optional_fields = bam_tag_list_to_string(alignment_p->optional_tags, &optional_fields_length);
 
+
+
     sequence_length = strlen(alignment_p->sequence);
 
     //data length is the sum of lengths if five codified fields (cigar, query name, sequence, quality and optional info)

diff --git a/lib/c/src/commons/workflow_scheduler.c b/lib/c/src/commons/workflow_scheduler.c
@@ -568,7 +568,7 @@ void workflow_run_with(int num_threads, void *input, workflow_t *wf) {
   pthread_t threads[num_threads];
   pthread_attr_t attr;
 
-  int num_cpus = 64;
+  int num_cpus = NUM_MAX_CPUS;
   int cpuArray[num_cpus];
 
   for (int i = 0; i < num_cpus; i++) {

diff --git a/lib/c/src/commons/workflow_scheduler.h b/lib/c/src/commons/workflow_scheduler.h
@@ -21,6 +21,8 @@ extern "C" {
 #define WORKFLOW_STATUS_RUNNING  1
 #define WORKFLOW_STATUS_FINISHED 2
 
+#define NUM_MAX_CPUS 128
+
 //----------------------------------------------------------------------------------------
 // work_item
 //----------------------------------------------------------------------------------------

diff --git a/man/manual.pdf b/man/manual.pdf
diff --git a/src/bs/bs_writer.c b/src/bs/bs_writer.c
@@ -19,8 +19,8 @@ int bs_writer(void *data) {
   mapping_batch_t *mapping_batch = (mapping_batch_t *) batch->mapping_batch;
   bs_context_t *bs_context = (mapping_batch->bs_context);
 
-  // Set the sequences of the mapping to the original //RICARDO why revert? (if restore, cigar problem with hard clipping)
-  //revert_mappings_seqs(mapping_batch->mapping_lists, mapping_batch->mapping_lists2, mapping_batch->fq_batch);
+  // Set the sequences of the mapping to the original //RICARDO adapted to avoid cigar problem with hard clipping
+  revert_mappings_seqs(mapping_batch->mapping_lists, mapping_batch->mapping_lists2, mapping_batch->fq_batch);
 
   batch_writer_input_t *writer_input = batch->writer_input;
   bam_file_t *bam_file = writer_input->bam_file;     
@@ -443,6 +443,9 @@ void write_mapped_read(array_list_t *array_list, bam_file_t *bam_file) {
     LOG_DEBUG("writting bam..\n");
 
     if (alig != NULL) {
+      //Ricardo - quitar
+    	//alig->quality[0] = "\0";
+
       bam1 = convert_to_bam(alig, 33);
       bam_fwrite(bam1, bam_file);
       bam_destroy1(bam1);

diff --git a/src/bs/methylation.c b/src/bs/methylation.c
@@ -135,7 +135,58 @@ void revert_mappings_seqs(array_list_t **src1, array_list_t **src2, array_list_t
         free(align_tmp->sequence);
       }
 
-      align_tmp->sequence = strdup(fastq_orig->sequence);
+      if (align_tmp->num_cigar_operations > 1)	//Check for Hard clipping in order to remove bases and Soft clipping in order to restore quality
+           {
+         	  //Find hard clipping
+
+         	  int i=0;
+         	  char aux[10000], aux2[10000];
+
+         	  char operation;
+         	  int bases;
+         	  int inicio=0;
+         	  int fin = strlen(fastq_orig->sequence);
+         	  bool restore_quality = false;
+
+         	  strcpy(aux, align_tmp->cigar);
+         	  for (i=0; i<align_tmp->num_cigar_operations;i++)
+         	  {
+         		  sscanf(aux,"%d%c%s", &bases, &operation, aux2);
+         		  if (operation == 'H')	//Hard Clipping
+         		  {
+         			  if (i==0) //Hard clipping in the start of the read
+         			  {
+         				inicio = bases;
+
+         			  }
+         			  else	//Hard clipping in the end of the read
+         			  {
+         				 fin = fin - bases;
+         			  }
+         		  }
+         		  else  if (operation == 'S')	//Soft Clipping
+         		  {
+         			  restore_quality = true;
+         		  }
+
+         		  strcpy(aux,aux2);
+
+
+         	  }
+
+         	  align_tmp->sequence = strndup(fastq_orig->sequence+inicio,fin-inicio);
+         	  if (restore_quality)
+         	  {
+         		 if (align_tmp->quality != NULL) {
+         		         free(align_tmp->quality);
+         		       }
+         		 align_tmp->quality = strndup(fastq_orig->quality+inicio,fin-inicio);
+         	  }
+           }
+
+      else{
+    	        align_tmp->sequence = strdup(fastq_orig->sequence);
+      }
     }
 
     // Go over all the alignments in list 2
@@ -149,7 +200,57 @@ void revert_mappings_seqs(array_list_t **src1, array_list_t **src2, array_list_t
         free(align_tmp->sequence);
       }
 
-      align_tmp->sequence = strdup(fastq_orig->sequence);
+      if (align_tmp->num_cigar_operations > 1)	//Check for Hard clipping in order to remove bases and Soft clipping in order to restore quality
+      {
+    	  //Find hard clipping
+    	  int i=0;
+    	  char aux[10000], aux2[10000];
+
+    	  char operation;
+    	  int bases;
+    	  int inicio=0;
+    	  int fin = strlen(fastq_orig->sequence);
+    	  bool restore_quality = false;
+
+    	  strcpy(aux, align_tmp->cigar);
+    	  for (i=0; i<align_tmp->num_cigar_operations;i++)
+    	  {
+    		  sscanf(aux,"%d%c%s", &bases, &operation, aux2);
+    		  if (operation == 'H')	//Hard Clipping
+    		  {
+    			  if (i==0) //Hard clipping in the start of the read
+    			  {
+    				inicio = bases;
+
+    			  }
+    			  else	//Hard clipping in the end of the read
+    			  {
+    				 fin = fin - bases;
+    			  }
+    		  }
+    		  else  if (operation == 'S')	//Soft Clipping
+    		  {
+    			  restore_quality = true;
+    		  }
+
+    		  strcpy(aux,aux2);
+
+
+    	  }
+
+    	  align_tmp->sequence = strndup(fastq_orig->sequence+inicio,fin-inicio);
+    	  if (restore_quality)
+    	  {
+    	  	 if (align_tmp->quality != NULL) {
+    	         free(align_tmp->quality);
+    	     }
+    	  	 align_tmp->quality =  strndup(fastq_orig->quality+inicio,fin-inicio);
+    	  }
+      }
+      else
+      {
+    	  align_tmp->sequence = strdup(fastq_orig->sequence);
+      }
     }
   }
 }

diff --git a/src/cal_seeker_bs.c b/src/cal_seeker_bs.c
@@ -218,7 +218,10 @@ void fill_gaps_bs(mapping_batch_t *mapping_batch, sw_optarg_t *sw_optarg,
 						gap_read_len = gap_read_end - gap_read_start + 1;
 						gap_genome_len = gap_genome_end - gap_genome_start + 1;
 
-						cal->start = gap_genome_start;
+						if ((int)gap_genome_start >0)	//cal->start must be higher than 1
+								cal->start = gap_genome_start;
+						else
+							;//RICARDO
 
 						assert(gap_read_len != 0);
 						assert(gap_genome_len != 0);

diff --git a/src/options.h b/src/options.h
@@ -29,6 +29,8 @@
 #define DEFAULT_REGION_THREADS		     1
 #define DEFAULT_NUM_SW_THREADS		     1
 
+
+
 // BEGIN: Mariano (4/12/2014)
 #define DEFAULT_NUM_SEEDS		         10
 // END: Mariano (4/12/2014)