Merge branch 'hotfix/v2.1.1'

Author: keiranmraine

c/bam_access.c

@@ -129,6 +129,8 @@ rg_info_t **bam_access_parse_header(bam_hdr_t *head, int *grps_size, stats_rd_t
     (*grp_stats)[j][0]->divergent= 0;
     (*grp_stats)[j][0]->mapped_bases= 0;
     (*grp_stats)[j][0]->proper= 0;
+    (*grp_stats)[j][0]->mapped_pairs= 0;
+    (*grp_stats)[j][0]->inter_chr_pairs= 0;
     (*grp_stats)[j][0]->inserts = kh_init(ins);
     (*grp_stats)[j][1] = (stats_rd_t *)malloc(sizeof(stats_rd_t));//Setup read two stats store
     check_mem((*grp_stats)[j][1]);
@@ -207,13 +209,9 @@ int bam_access_process_reads(htsFile *input, bam_hdr_t *head, rg_info_t **grps,
       uint32_t nm_val = bam_aux2i(nm);
       if(nm_val>0){
         (*grp_stats)[rg_index][read]->divergent += nm_val;
-        (*grp_stats)[rg_index][read]->mapped_bases += bam_access_get_mapped_base_count_from_cigar(b);
-      }else{
-        (*grp_stats)[rg_index][read]->mapped_bases += (bam_endpos(b) - b->core.pos) + 1;
       }
-    }else{
-      (*grp_stats)[rg_index][read]->mapped_bases += bam_access_get_mapped_base_count_from_cigar(b);
     }
+    (*grp_stats)[rg_index][read]->mapped_bases += bam_access_get_mapped_base_count_from_cigar(b);
 
     // stats that only assess read 1
     if(b->core.flag & BAM_FREAD1) {

c/c_tests/02_bam_access_tests.c

@@ -63,31 +63,31 @@ char *exp_platform_unit_2 = "5085_6";
 uint32_t exp_rd_length = 20;
 uint64_t exp_rg1_rd1_tot_count = 7;
 uint64_t exp_rg1_rd2_tot_count = 3;
-uint64_t exp_rg2_rd1_tot_count = 3;
-uint64_t exp_rg2_rd2_tot_count = 3;
+uint64_t exp_rg2_rd1_tot_count = 4;
+uint64_t exp_rg2_rd2_tot_count = 4;
 uint64_t exp_rg1_rd1_dups = 4;
 uint64_t exp_rg1_rd2_dups = 0;
 uint64_t exp_rg2_rd1_dups = 0;
 uint64_t exp_rg2_rd2_dups = 0;
 uint64_t exp_rg1_rd1_gc = 63;
 uint64_t exp_rg1_rd2_gc = 27;
-uint64_t exp_rg2_rd1_gc = 27;
-uint64_t exp_rg2_rd2_gc = 27;
+uint64_t exp_rg2_rd1_gc = 71;
+uint64_t exp_rg2_rd2_gc = 74;
 uint64_t exp_rg1_rd1_umap = 1;
 uint64_t exp_rg1_rd2_umap = 1;
 uint64_t exp_rg2_rd1_umap = 2;
 uint64_t exp_rg2_rd2_umap = 2;
 uint64_t exp_rg1_rd1_divergent = 30;
 uint64_t exp_rg1_rd2_divergent = 18;
-uint64_t exp_rg2_rd1_divergent = 11;
+uint64_t exp_rg2_rd1_divergent = 12;
 uint64_t exp_rg2_rd2_divergent = 12;
 uint64_t exp_rg1_rd1_mapped_bases = 115;
 uint64_t exp_rg1_rd2_mapped_bases = 37;
-uint64_t exp_rg2_rd1_mapped_bases = 20;
-uint64_t exp_rg2_rd2_mapped_bases = 20;
+uint64_t exp_rg2_rd1_mapped_bases = 95;
+uint64_t exp_rg2_rd2_mapped_bases = 95;
 uint64_t exp_rg1_rd1_proper = 6;
 uint64_t exp_rg1_rd2_proper = 0;
-uint64_t exp_rg2_rd1_proper = 1;
+uint64_t exp_rg2_rd1_proper = 2;
 uint64_t exp_rg2_rd2_proper = 0;
 
 char err[100];
@@ -338,22 +338,22 @@ char *test_bam_access_process_reads_rna(){ // rna flag in this method includes s
   }
   //Duplicate reads 1
   if(grp_stats[1][0]->dups!=exp_rg2_rd1_dups){
-    sprintf(err,"RG 1, read_1 duplicate count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd1_dups,grp_stats[1][0]->dups);
+    sprintf(err,"RG 2, read_1 duplicate count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd1_dups,grp_stats[1][0]->dups);
     return err;
   }
   //Duplicate reads 2
   if(grp_stats[1][1]->dups!=exp_rg2_rd2_dups){
-    sprintf(err,"RG 1, read_2 duplicate count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd2_dups,grp_stats[1][1]->dups);
+    sprintf(err,"RG 2, read_2 duplicate count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd2_dups,grp_stats[1][1]->dups);
     return err;
   }
   //GC 1
   if(grp_stats[1][0]->gc!=exp_rg2_rd1_gc){
-    sprintf(err,"RG 1, read_1 gc count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd1_gc,grp_stats[1][0]->gc);
+    sprintf(err,"RG 2, read_1 gc count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd1_gc,grp_stats[1][0]->gc);
     return err;
   }
   //GC 2
   if(grp_stats[1][1]->gc!=exp_rg2_rd2_gc){
-    sprintf(err,"RG 1, read_2 gc count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd2_gc,grp_stats[1][1]->gc);
+    sprintf(err,"RG 2, read_2 gc count incorrect. Expected %"PRIu64" but got %"PRIu64"\n",exp_rg2_rd2_gc,grp_stats[1][1]->gc);
     return err;
   }
   //Unmapped 1

lib/PCAP/Bam/Stats.pm

@@ -210,17 +210,11 @@ sub _process_reads {
     if(defined $nm) {
       if($nm > 0) {
         $rg_ref->{'total_divergent_bases_'.$read} += $nm;
-        my @increments = $a->cigar_str =~ m/([[:digit:]]+)[IMX=]/g;
-        $rg_ref->{'total_mapped_bases_'.$read} += sum0 @increments;
       }
-      else {
-        $rg_ref->{'total_mapped_bases_'.$read} += ($a->calend - $a->pos) + 1;
-      }
-    }
-    else {
-      my @increments = $a->cigar_str =~ m/([[:digit:]]+)[IMX=]/g;
-      $rg_ref->{'total_mapped_bases_'.$read} += sum0 @increments;
     }
+    my @increments = $a->cigar_str =~ m/([[:digit:]]+)[IMX=]/g;
+    $rg_ref->{'total_mapped_bases_'.$read} += sum0 @increments;
+
     # calculate mapped seq from cigar
 
     # stats that only assess read 1

t/data/Stats.bam

@@ -1,4 +1,4 @@
 bam_filename	sample	platform	platform_unit	library	readgroup	read_length_r1	read_length_r2	#_mapped_bases	#_mapped_bases_r1	#_mapped_bases_r2	#_divergent_bases	#_divergent_bases_r1	#_divergent_bases_r2	#_total_reads	#_total_reads_r1	#_total_reads_r2	#_mapped_reads	#_mapped_reads_r1	#_mapped_reads_r2	#_mapped_reads_properly_paired	#_gc_bases_r1	#_gc_bases_r2	mean_insert_size	insert_size_sd	median_insert_size	#_duplicate_reads	#_mapped_pairs	#_inter_chr_pairs
 Stats.bam	PD1234a	GAII	5178_6	PD1234a 140546_1054	29976	20	20	152	115	37	48	30	18	10	7	3	8	6	2	6	63	27	195.833	119.388	125.000	4	6	0
-Stats.bam	PD1234a	GAII	5085_6	PD1234a 140546_1054	29978	20	20	40	20	20	23	11	12	6	3	3	2	1	1	1	27	27	100.000	0.000	100.000	0	1	0
+Stats.bam	PD1234a	GAII	5085_6	PD1234a 140546_1054	29978	20	20	190	95	95	24	12	12	8	4	4	4	2	2	2	71	74	545.500	445.500	545.500	0	2	0
 Stats.bam	PD1234a	GAII	5086_6	PD1234a 140546_1054	29979	20	20	120	60	60	0	0	0	6	3	3	6	3	3	1	27	27	100.000	0.000	100.000	0	3	1

t/data/Stats.bam.bas

@@ -1,4 +1,4 @@
 bam_filename	sample	platform	platform_unit	library	readgroup	read_length_r1	read_length_r2	#_mapped_bases	#_mapped_bases_r1	#_mapped_bases_r2	#_divergent_bases	#_divergent_bases_r1	#_divergent_bases_r2	#_total_reads	#_total_reads_r1	#_total_reads_r2	#_mapped_reads	#_mapped_reads_r1	#_mapped_reads_r2	#_mapped_reads_properly_paired	#_gc_bases_r1	#_gc_bases_r2	mean_insert_size	insert_size_sd	median_insert_size	#_duplicate_reads	#_mapped_pairs	#_inter_chr_pairs
 Stats.bam	PD1234a	GAII	5178_6	PD1234a 140546_1054	29976	20	20	152	115	37	48	30	18	10	7	3	8	6	2	6	63	27	195.833	119.387	125.000	4	6	0
-Stats.bam	PD1234a	GAII	5085_6	PD1234a 140546_1054	29978	20	20	40	20	20	23	11	12	6	3	3	2	1	1	1	27	27	100.000	0.000	100.000	0	1	0
+Stats.bam	PD1234a	GAII	5085_6	PD1234a 140546_1054	29978	20	20	190	95	95	24	12	12	8	4	4	4	2	2	2	71	74	545.500	445.500	545.500	0	2	0
 Stats.bam	PD1234a	GAII	5086_6	PD1234a 140546_1054	29979	20	20	120	60	60	0	0	0	6	3	3	6	3	3	1	27	27	100.000	0.000	100.000	0	3	1

t/data/Stats.c.bam.bas

@@ -121,9 +121,9 @@ subtest 'Object funcions' => sub {
   subtest 'calc_frac_' => sub {
     my $test_obj = _create_test_object();
 
-    is(sprintf('%.4f',$test_obj->calc_frac_properly_paired()),'0.6154','calc_frac_properly_paired');
-    is(sprintf('%.4f',$test_obj->calc_frac_unmapped()),'0.2727','calc_frac_unmapped');
-    is(sprintf('%.4f',$test_obj->calc_frac_duplicate_reads()),'0.1818','calc_frac_duplicate_reads');
+    is(sprintf('%.4f',$test_obj->calc_frac_properly_paired()),'0.6429','calc_frac_properly_paired');
+    is(sprintf('%.4f',$test_obj->calc_frac_unmapped()),'0.2500','calc_frac_unmapped');
+    is(sprintf('%.4f',$test_obj->calc_frac_duplicate_reads()),'0.1667','calc_frac_duplicate_reads');
   };
 
   subtest 'med_insert_size' => sub {
@@ -135,39 +135,39 @@ subtest 'Object funcions' => sub {
   subtest 'mean_insert_size' => sub {
     my $test_obj = _create_test_object();
 
-    is(sprintf('%.4f',$test_obj->mean_insert_size()),'171.8750','mean_insert_size');
+    is(sprintf('%.4f',$test_obj->mean_insert_size()),'262.8889','mean_insert_size');
   };
 
   subtest 'properly_mapped_ratio' => sub {
     my $test_obj = _create_test_object();
-    is(sprintf('%.4f',$test_obj->properly_mapped_ratio()),'0.6154','properly_mapped_ratio');
+    is(sprintf('%.4f',$test_obj->properly_mapped_ratio()),'0.6429','properly_mapped_ratio');
   };
 
   subtest '_count_' => sub {
     my $test_obj = _create_test_object();
 
-    is($test_obj->count_total_reads(1),13,'count_total_reads read1');
-    is($test_obj->count_total_reads(2),9,'count_total_reads read2');
+    is($test_obj->count_total_reads(1),14,'count_total_reads read1');
+    is($test_obj->count_total_reads(2),10,'count_total_reads read2');
 
-    is($test_obj->count_properly_paired(),8,'count_properly_paired');
+    is($test_obj->count_properly_paired(),9,'count_properly_paired');
 
     is($test_obj->count_unmapped(1),3,'count_unmapped read1');
     is($test_obj->count_unmapped(2),3,'count_unmapped read2');
 
     is($test_obj->count_duplicate_reads(1),4,'count_duplicate_reads read1');
     is($test_obj->count_duplicate_reads(2),0,'count_duplicate_reads read2');
 
-    is($test_obj->count_total_mapped_bases(1),195,'count_total_mapped_bases read1');
-    is($test_obj->count_total_mapped_bases(2),117,'count_total_mapped_bases read2');
+    is($test_obj->count_total_mapped_bases(1),270,'count_total_mapped_bases read1');
+    is($test_obj->count_total_mapped_bases(2),192,'count_total_mapped_bases read2');
 
-    is($test_obj->count_total_divergent_bases(1),41,'count_total_divergent_bases read1');
+    is($test_obj->count_total_divergent_bases(1),42,'count_total_divergent_bases read1');
     is($test_obj->count_total_divergent_bases(2),30,'count_total_divergent_bases read2');
   };
 
   subtest 'cal_insert_size_sd' => sub {
      my $test_obj = _create_test_object();
 
-    is(sprintf('%.4f',$test_obj->cal_insert_size_sd()),'111.4096','cal_insert_size_sd');
+    is(sprintf('%.4f',$test_obj->cal_insert_size_sd()),'278.0310','cal_insert_size_sd');
   };
 
 };
@@ -285,3 +285,7 @@ IL29_5086:6:35:17751:10873	129	1	9993	0	20M	=	9993	-100	CTCTTCCGATCTTTAGGGTT	;?;
 # non-proper, inter-chr
 IL29_5086:6:35:17751:10874	65	1	9993	0	20M	X	9993	0	CTCTTCCGATCTTTAGGGTT	;?;??>>>>F<BBDEBEEFF	RG:Z:29979
 IL29_5086:6:35:17751:10874	129	X	9993	0	20M	1	9993	0	CTCTTCCGATCTTTAGGGTT	;?;??>>>>F<BBDEBEEFF	RG:Z:29979
+
+## Checks that mapped bases don't get corrupt when padding values are introduced, normally RNA issue:
+HS20_13617:3:1315:4901:63745	163	1	12670	255	52M761N23M	=	13586	991	GACGACGGCCGACTTGGATCACACTCTTGTGAGTGTCCCCAGTGTTGCAGAGGCAGCTGCACCACTGCCTGGCGC	CCCFFFFFHHHHHJJJJJJJJJJJJJJJIGIHJDHGIJJJIJGHIIIJHHHHFFFFDEEEDDDDDDDDDDDDDDB	NH:i:1	HI:i:1	NM:i:0	MD:Z:75	AS:i:146	XS:A:+	RG:Z:29978
+HS20_13617:3:1315:4901:63745	83	1	13586	255	75M	=	12670	-991	ACCCAGGAGAGTGTGGAGTCCAGAGTGATGCCAGGACCCAGGCACAGGCATTAGTGCCCGTTGGAGAAAACAGGG	EFFFFFHHHHHJIJJJJIJIJJJJJJJIJJIJJJHHIJJJJIGGIJJJJJJJJJJJJJJJJJHHHHHFFFFFCCC	NH:i:1	HI:i:1	NM:i:1	MD:Z:27T47	AS:i:146	XS:A:+	RG:Z:29978