1. gmapR
  2. src
  3. bamtally.c
Browse code

support new version of bam_tally

Michael Lawrence authored on 12/06/2019 22:48:54
Showing 1 changed files
src/bamtally.c
History View file @ c3db7e3
... ... 
@@ -19,11 +19,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
19 19 
                 SEXP maximum_nhits_R,
20 20 
                 SEXP need_concordant_p_R, SEXP need_unique_p_R,
21 21 
                 SEXP need_primary_p_R, SEXP ignore_duplicates_p_R,
22 
-                SEXP min_depth_R, SEXP variant_strands_R,
22 
+                SEXP min_depth_R, SEXP variant_strands_R, SEXP variant_pct_R,
23 23 
                 SEXP ignore_query_Ns_p_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 
-                SEXP verbosep_R, SEXP max_softclip_R,
26 
+                SEXP verbosep_R, SEXP min_softclip_R, SEXP max_softclip_R,
27 27 
                 SEXP exon_iit_file_R,
28 28 
                 SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 29 
                 SEXP minimum_quality_score_R, SEXP noncovered_R,
... ... 
@@ -34,6 +34,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
34 34 
     genome_dir_R == R_NilValue ? NULL : CHAR(asChar(genome_dir_R));
35 35 
   const char *db = CHAR(asChar(db_R));
36 36 
   int alloclength = asInteger(alloclength_R);
37 
+  int pastlength = alloclength;
37 38 
   const char *desired_read_group =
38 39 
     desired_read_group_R == R_NilValue ? NULL :
39 40 
     CHAR(asChar(desired_read_group_R));
... ... 
@@ -46,10 +47,12 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
46 47 
   bool ignore_duplicates_p = asLogical(ignore_duplicates_p_R);
47 48 
   int min_depth = asInteger(min_depth_R);
48 49 
   int variant_strands = asInteger(variant_strands_R);
50 
+  double variant_pct = asReal(variant_pct_R);
49 51 
   bool ignore_query_Ns_p = asLogical(ignore_query_Ns_p_R);
50 52 
   bool print_indels_p = asLogical(print_indels_p_R);
51 53 
   int blocksize = asInteger(blocksize_R);
52 54 
   int verbosep = asLogical(verbosep_R);
55 
+  int min_softclip = asInteger(min_softclip_R);
53 56 
   int max_softclip = asInteger(max_softclip_R);
54 57 
   bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
55 58 
   bool print_cycles_p = asLogical(print_cycles_p_R);
... ... 
@@ -79,7 +82,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
79 82 
                                  genome,
80 83 
                                  chromosome_iit,
81 84 
                                  map_iit,
82 
-                                 alloclength,
85 
+                                 alloclength, pastlength,
83 86 
                                  (char *)desired_read_group,
84 87 
                                  minimum_mapq, good_unique_mapq,
85 88 
                                  minimum_quality_score,
... ... 
@@ -87,10 +90,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
87 90 
                                  need_unique_p, need_primary_p,
88 91 
                                  ignore_duplicates_p,
89 92 
                                  min_depth, variant_strands,
90 
-				 /* TODO: variant_pct */ 0,
93 
+				 variant_pct,
91 94 
 				 ignore_query_Ns_p,
92 95 
                                  print_indels_p, blocksize, verbosep,
93 
-                                 /*readlevel_p*/false, max_softclip,
96 
+                                 /*readlevel_p*/false,
97 
+				 min_softclip, max_softclip,
94 98 
                                  print_cycles_p,
95 99 
                                  print_nm_scores_p,
96 100 
                                  print_xs_scores_p, noncovered_p);
Browse code

update to latest GSTRUCT

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@123972 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 14/11/2016 21:54:30
Showing 1 changed files
src/bamtally.c
History View file @ 55d77a1
... ... 
@@ -86,7 +86,9 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
86 86 
                                  maximum_nhits, need_concordant_p,
87 87 
                                  need_unique_p, need_primary_p,
88 88 
                                  ignore_duplicates_p,
89 
-                                 min_depth, variant_strands, ignore_query_Ns_p,
89 
+                                 min_depth, variant_strands,
90 
+				 /* TODO: variant_pct */ 0,
91 
+				 ignore_query_Ns_p,
90 92 
                                  print_indels_p, blocksize, verbosep,
91 93 
                                  /*readlevel_p*/false, max_softclip,
92 94 
                                  print_cycles_p,
Browse code

fixes for bam_tally update

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@110596 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 13/11/2015 22:40:27
Showing 1 changed files
src/bamtally.c
History View file @ 38518dd
... ... 
@@ -24,7 +24,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 26 
                 SEXP verbosep_R, SEXP max_softclip_R,
27 
-                SEXP genome_iit_file_R,
27 
+                SEXP exon_iit_file_R,
28 28 
                 SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 29 
                 SEXP minimum_quality_score_R, SEXP noncovered_R,
30 30 
                 SEXP print_nm_scores_p_R)
... ... 
@@ -59,8 +59,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
59 59
60 60 
   Genome_T genome = createGenome(genome_dir, db);
61 61 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
62 
-  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL :
63 
-    IIT_read((char*)CHAR(asChar(genome_iit_file_R)), /*name*/ NULL,
62 
+  IIT_T map_iit = exon_iit_file_R == R_NilValue ? NULL :
63 
+    IIT_read((char*)CHAR(asChar(exon_iit_file_R)), /*name*/ NULL,
64 64 
              /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL,
65 65 
              /*add_iit_p*/false, /*labels_read_p*/true);
66 66 
   const char *chr = NULL;
Browse code

work towards supporting new bam_tally; drops quality info (now cutoff based), adds xs/nm/del and codon-level counts

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@110593 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 13/11/2015 21:50:11
Showing 1 changed files
src/bamtally.c
History View file @ 685d60e
... ... 
@@ -26,7 +26,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
26 26 
                 SEXP verbosep_R, SEXP max_softclip_R,
27 27 
                 SEXP genome_iit_file_R,
28 28 
                 SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 
-                SEXP minimum_quality_score_R, SEXP noncovered_R)
29 
+                SEXP minimum_quality_score_R, SEXP noncovered_R,
30 
+                SEXP print_nm_scores_p_R)
30 31 
 {
31 32 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
32 33 
   const char *genome_dir =
... ... 
@@ -54,6 +55,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
54 55 
   bool print_cycles_p = asLogical(print_cycles_p_R);
55 56 
   int minimum_quality_score = asInteger(minimum_quality_score_R);
56 57 
   bool noncovered_p = asLogical(noncovered_R);
58 
+  bool print_nm_scores_p = asLogical(print_nm_scores_p_R);
57 59
58 60 
   Genome_T genome = createGenome(genome_dir, db);
59 61 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
... ... 
@@ -88,7 +90,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
88 90 
                                  print_indels_p, blocksize, verbosep,
89 91 
                                  /*readlevel_p*/false, max_softclip,
90 92 
                                  print_cycles_p,
91 
-                                 /*print_nm_scores_p*/ false,
93 
+                                 print_nm_scores_p,
92 94 
                                  print_xs_scores_p, noncovered_p);
93 95 
   IIT_free(&chromosome_iit);
94 96 
   Genome_free(&genome);
Browse code

resurrect improvements that we reverted before the April release, expect breakage for a while

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@110039 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 28/10/2015 18:51:00
Showing 1 changed files
src/bamtally.c
History View file @ 8418591
... ... 
@@ -24,7 +24,9 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 26 
                 SEXP verbosep_R, SEXP max_softclip_R,
27 
-		SEXP genome_iit_file_R)
27 
+                SEXP genome_iit_file_R,
28 
+                SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 
+                SEXP minimum_quality_score_R, SEXP noncovered_R)
28 30 
 {
29 31 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
30 32 
   const char *genome_dir =
... ... 
@@ -48,11 +50,17 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
48 50 
   int blocksize = asInteger(blocksize_R);
49 51 
   int verbosep = asLogical(verbosep_R);
50 52 
   int max_softclip = asInteger(max_softclip_R);
53 
+  bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
54 
+  bool print_cycles_p = asLogical(print_cycles_p_R);
55 
+  int minimum_quality_score = asInteger(minimum_quality_score_R);
56 
+  bool noncovered_p = asLogical(noncovered_R);
51 57
52 58 
   Genome_T genome = createGenome(genome_dir, db);
53 59 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
54 
-  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL : IIT_read(CHAR(asChar(genome_iit_file_R)), /*name*/ NULL, /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL, /*add_iit_p*/false, /*labels_read_p*/true);
55 
-//  IIT_T map_iit = NULL;
60 
+  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL :
61 
+    IIT_read((char*)CHAR(asChar(genome_iit_file_R)), /*name*/ NULL,
62 
+             /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL,
63 
+             /*add_iit_p*/false, /*labels_read_p*/true);
56 64 
   const char *chr = NULL;
57 65 
   Genomicpos_T start = 0;
58 66 
   Genomicpos_T end = 0;
... ... 
@@ -67,19 +75,21 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
67 75 
                                  /* TODO: bam_lacks_chr */ NULL,
68 76 
                                  start, end,
69 77 
                                  genome,
70 
-				 /*chromosome_iit */ chromosome_iit,
71 
-				 /* map_iit*/ map_iit,
72 
-				 alloclength,
73 
-				 (char *)desired_read_group,
78 
+                                 chromosome_iit,
79 
+                                 map_iit,
80 
+                                 alloclength,
81 
+                                 (char *)desired_read_group,
74 82 
                                  minimum_mapq, good_unique_mapq,
83 
+                                 minimum_quality_score,
75 84 
                                  maximum_nhits, need_concordant_p,
76 85 
                                  need_unique_p, need_primary_p,
77 86 
                                  ignore_duplicates_p,
78 87 
                                  min_depth, variant_strands, ignore_query_Ns_p,
79 88 
                                  print_indels_p, blocksize, verbosep,
80 89 
                                  /*readlevel_p*/false, max_softclip,
81 
-				 /* print_xs_scores_p ??? */ false,
82 
-				 /* print_noncovered_p ??? */ false);
90 
+                                 print_cycles_p,
91 
+                                 /*print_nm_scores_p*/ false,
92 
+                                 print_xs_scores_p, noncovered_p);
83 93 
   IIT_free(&chromosome_iit);
84 94 
   Genome_free(&genome);
85 95 
   if(map_iit != NULL)
Browse code

Reverted to r101133, along with NAMESPACE fixes

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@102487 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 14/04/2015 21:40:44
Showing 1 changed files
src/bamtally.c
History View file @ d5ed124
... ... 
@@ -24,9 +24,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 26 
                 SEXP verbosep_R, SEXP max_softclip_R,
27 
-                SEXP genome_iit_file_R,
28 
-                SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 
-                SEXP minimum_quality_score_R, SEXP noncovered_R)
27 
+		SEXP genome_iit_file_R)
30 28 
 {
31 29 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
32 30 
   const char *genome_dir =
... ... 
@@ -50,17 +48,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
50 48 
   int blocksize = asInteger(blocksize_R);
51 49 
   int verbosep = asLogical(verbosep_R);
52 50 
   int max_softclip = asInteger(max_softclip_R);
53 
-  bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
54 
-  bool print_cycles_p = asLogical(print_cycles_p_R);
55 
-  int minimum_quality_score = asInteger(minimum_quality_score_R);
56 
-  bool noncovered_p = asLogical(noncovered_R);
57 51
58 52 
   Genome_T genome = createGenome(genome_dir, db);
59 53 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
60 
-  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL :
61 
-    IIT_read((char*)CHAR(asChar(genome_iit_file_R)), /*name*/ NULL,
62 
-             /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL,
63 
-             /*add_iit_p*/false, /*labels_read_p*/true);
54 
+  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL : IIT_read(CHAR(asChar(genome_iit_file_R)), /*name*/ NULL, /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL, /*add_iit_p*/false, /*labels_read_p*/true);
55 
+//  IIT_T map_iit = NULL;
64 56 
   const char *chr = NULL;
65 57 
   Genomicpos_T start = 0;
66 58 
   Genomicpos_T end = 0;
... ... 
@@ -75,21 +67,19 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
75 67 
                                  /* TODO: bam_lacks_chr */ NULL,
76 68 
                                  start, end,
77 69 
                                  genome,
78 
-                                 chromosome_iit,
79 
-                                 map_iit,
80 
-                                 alloclength,
81 
-                                 (char *)desired_read_group,
70 
+				 /*chromosome_iit */ chromosome_iit,
71 
+				 /* map_iit*/ map_iit,
72 
+				 alloclength,
73 
+				 (char *)desired_read_group,
82 74 
                                  minimum_mapq, good_unique_mapq,
83 
-                                 minimum_quality_score,
84 75 
                                  maximum_nhits, need_concordant_p,
85 76 
                                  need_unique_p, need_primary_p,
86 77 
                                  ignore_duplicates_p,
87 78 
                                  min_depth, variant_strands, ignore_query_Ns_p,
88 79 
                                  print_indels_p, blocksize, verbosep,
89 80 
                                  /*readlevel_p*/false, max_softclip,
90 
-                                 print_cycles_p,
91 
-                                 /*print_nm_scores_p*/ false,
92 
-                                 print_xs_scores_p, noncovered_p);
81 
+				 /* print_xs_scores_p ??? */ false,
82 
+				 /* print_noncovered_p ??? */ false);
93 83 
   IIT_free(&chromosome_iit);
94 84 
   Genome_free(&genome);
95 85 
   if(map_iit != NULL)
Browse code

Recent work towards supporting the new features... will revert back to a stable version immediately...

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@102485 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 14/04/2015 21:32:10
Showing 1 changed files
src/bamtally.c
History View file @ c8fd96e
... ... 
@@ -24,7 +24,9 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 26 
                 SEXP verbosep_R, SEXP max_softclip_R,
27 
-		SEXP genome_iit_file_R)
27 
+                SEXP genome_iit_file_R,
28 
+                SEXP print_xs_scores_p_R, SEXP print_cycles_p_R,
29 
+                SEXP minimum_quality_score_R, SEXP noncovered_R)
28 30 
 {
29 31 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
30 32 
   const char *genome_dir =
... ... 
@@ -48,6 +50,10 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
48 50 
   int blocksize = asInteger(blocksize_R);
49 51 
   int verbosep = asLogical(verbosep_R);
50 52 
   int max_softclip = asInteger(max_softclip_R);
53 
+  bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
54 
+  bool print_cycles_p = asLogical(print_cycles_p_R);
55 
+  int minimum_quality_score = asInteger(minimum_quality_score_R);
56 
+  bool noncovered_p = asLogical(noncovered_R);
51 57
52 58 
   Genome_T genome = createGenome(genome_dir, db);
53 59 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
... ... 
@@ -69,19 +75,21 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
69 75 
                                  /* TODO: bam_lacks_chr */ NULL,
70 76 
                                  start, end,
71 77 
                                  genome,
72 
-				 chromosome_iit,
73 
-				 map_iit,
74 
-				 alloclength,
75 
-				 (char *)desired_read_group,
78 
+                                 chromosome_iit,
79 
+                                 map_iit,
80 
+                                 alloclength,
81 
+                                 (char *)desired_read_group,
76 82 
                                  minimum_mapq, good_unique_mapq,
83 
+                                 minimum_quality_score,
77 84 
                                  maximum_nhits, need_concordant_p,
78 85 
                                  need_unique_p, need_primary_p,
79 86 
                                  ignore_duplicates_p,
80 87 
                                  min_depth, variant_strands, ignore_query_Ns_p,
81 88 
                                  print_indels_p, blocksize, verbosep,
82 89 
                                  /*readlevel_p*/false, max_softclip,
83 
-				 /* print_xs_scores_p */ false,
84 
-				 /* print_noncovered_p */ false);
90 
+                                 print_cycles_p,
91 
+                                 /*print_nm_scores_p*/ false,
92 
+                                 print_xs_scores_p, noncovered_p);
85 93 
   IIT_free(&chromosome_iit);
86 94 
   Genome_free(&genome);
87 95 
   if(map_iit != NULL)
Browse code

update gstruct/bamtally

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@101176 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 25/03/2015 20:25:05
Showing 1 changed files
src/bamtally.c
History View file @ f58d4e9
... ... 
@@ -51,8 +51,10 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
51 51
52 52 
   Genome_T genome = createGenome(genome_dir, db);
53 53 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
54 
-  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL : IIT_read(CHAR(asChar(genome_iit_file_R)), /*name*/ NULL, /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL, /*add_iit_p*/false, /*labels_read_p*/true);
55 
-//  IIT_T map_iit = NULL;
54 
+  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL :
55 
+    IIT_read((char*)CHAR(asChar(genome_iit_file_R)), /*name*/ NULL,
56 
+             /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL,
57 
+             /*add_iit_p*/false, /*labels_read_p*/true);
56 58 
   const char *chr = NULL;
57 59 
   Genomicpos_T start = 0;
58 60 
   Genomicpos_T end = 0;
... ... 
@@ -67,8 +69,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
67 69 
                                  /* TODO: bam_lacks_chr */ NULL,
68 70 
                                  start, end,
69 71 
                                  genome,
70 
-				 /*chromosome_iit */ chromosome_iit,
71 
-				 /* map_iit*/ map_iit,
72 
+				 chromosome_iit,
73 
+				 map_iit,
72 74 
 				 alloclength,
73 75 
 				 (char *)desired_read_group,
74 76 
                                  minimum_mapq, good_unique_mapq,
... ... 
@@ -78,8 +80,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
78 80 
                                  min_depth, variant_strands, ignore_query_Ns_p,
79 81 
                                  print_indels_p, blocksize, verbosep,
80 82 
                                  /*readlevel_p*/false, max_softclip,
81 
-				 /* print_xs_scores_p ??? */ false,
82 
-				 /* print_noncovered_p ??? */ false);
83 
+				 /* print_xs_scores_p */ false,
84 
+				 /* print_noncovered_p */ false);
83 85 
   IIT_free(&chromosome_iit);
84 86 
   Genome_free(&genome);
85 87 
   if(map_iit != NULL)
Browse code

some fiddling, disable maintainer mode, vbump

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@93390 bc3139a8-67e5-0310-9ffc-ced21a209358

Gabriel Becker authored on 14/08/2014 20:39:03
Showing 1 changed files
src/bamtally.c
History View file @ cfb01d1
... ... 
@@ -79,11 +79,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
79 79 
                                  print_indels_p, blocksize, verbosep,
80 80 
                                  /*readlevel_p*/false, max_softclip,
81 81 
 				 /* print_xs_scores_p ??? */ false,
82 
-//				 /* print_noncovered_p ??? */ false);
83 
-				 /* print_noncovered_p ??? */ true);
82 
+				 /* print_noncovered_p ??? */ false);
84 83 
   IIT_free(&chromosome_iit);
85 84 
   Genome_free(&genome);
86 
-  IIT_free(&map_iit);
85 
+  if(map_iit != NULL)
86 
+    IIT_free(&map_iit);
87 87 
   if (tally_iit == NULL) {
88 88 
     error("Could not create tally\n");
89 89 
   }
Browse code

add codon tally support. No vbump yet

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@93316 bc3139a8-67e5-0310-9ffc-ced21a209358

Gabriel Becker authored on 11/08/2014 23:42:48
Showing 1 changed files
src/bamtally.c
History View file @ c78be69
... ... 
@@ -23,8 +23,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
23 23 
                 SEXP ignore_query_Ns_p_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 
-                SEXP verbosep_R, SEXP max_softclip_R, SEXP print_xs_scores_p_R,
27 
-                SEXP noncovered_R)
26 
+                SEXP verbosep_R, SEXP max_softclip_R,
27 
+		SEXP genome_iit_file_R)
28 28 
 {
29 29 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
30 30 
   const char *genome_dir =
... ... 
@@ -48,12 +48,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
48 48 
   int blocksize = asInteger(blocksize_R);
49 49 
   int verbosep = asLogical(verbosep_R);
50 50 
   int max_softclip = asInteger(max_softclip_R);
51 
-  bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
52 
-  int noncovered = asLogical(noncovered_R);
53 51
54 52 
   Genome_T genome = createGenome(genome_dir, db);
55 53 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
56 
-
54 
+  IIT_T map_iit = genome_iit_file_R == R_NilValue ? NULL : IIT_read(CHAR(asChar(genome_iit_file_R)), /*name*/ NULL, /*readonlyp*/true, /*divread*/READ_ALL, /*divstring*/NULL, /*add_iit_p*/false, /*labels_read_p*/true);
55 
+//  IIT_T map_iit = NULL;
57 56 
   const char *chr = NULL;
58 57 
   Genomicpos_T start = 0;
59 58 
   Genomicpos_T end = 0;
... ... 
@@ -67,8 +66,11 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
67 66 
   IIT_T tally_iit = Bamtally_iit(bamreader, (char *)chr,
68 67 
                                  /* TODO: bam_lacks_chr */ NULL,
69 68 
                                  start, end,
70 
-                                 genome, chromosome_iit, alloclength,
71 
-                                 (char *)desired_read_group,
69 
+                                 genome,
70 
+				 /*chromosome_iit */ chromosome_iit,
71 
+				 /* map_iit*/ map_iit,
72 
+				 alloclength,
73 
+				 (char *)desired_read_group,
72 74 
                                  minimum_mapq, good_unique_mapq,
73 75 
                                  maximum_nhits, need_concordant_p,
74 76 
                                  need_unique_p, need_primary_p,
... ... 
@@ -76,10 +78,12 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
76 78 
                                  min_depth, variant_strands, ignore_query_Ns_p,
77 79 
                                  print_indels_p, blocksize, verbosep,
78 80 
                                  /*readlevel_p*/false, max_softclip,
79 
-                                 print_xs_scores_p, noncovered);
81 
+				 /* print_xs_scores_p ??? */ false,
82 
+//				 /* print_noncovered_p ??? */ false);
83 
+				 /* print_noncovered_p ??? */ true);
80 84 
   IIT_free(&chromosome_iit);
81 85 
   Genome_free(&genome);
82 
-
86 
+  IIT_free(&map_iit);
83 87 
   if (tally_iit == NULL) {
84 88 
     error("Could not create tally\n");
85 89 
   }
Browse code

update to new bam_tally with support for XS counting, which we now support via BamTallyParam@count_xs.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@90225 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 13/05/2014 02:04:05
Showing 1 changed files
src/bamtally.c
History View file @ f5e63d3
... ... 
@@ -23,7 +23,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
23 23 
                 SEXP ignore_query_Ns_p_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 
-                SEXP verbosep_R, SEXP max_softclip_R, SEXP noncovered_R)
26 
+                SEXP verbosep_R, SEXP max_softclip_R, SEXP print_xs_scores_p_R,
27 
+                SEXP noncovered_R)
27 28 
 {
28 29 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
29 30 
   const char *genome_dir =
... ... 
@@ -47,6 +48,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
47 48 
   int blocksize = asInteger(blocksize_R);
48 49 
   int verbosep = asLogical(verbosep_R);
49 50 
   int max_softclip = asInteger(max_softclip_R);
51 
+  bool print_xs_scores_p = asLogical(print_xs_scores_p_R);
50 52 
   int noncovered = asLogical(noncovered_R);
51 53
52 54 
   Genome_T genome = createGenome(genome_dir, db);
... ... 
@@ -73,7 +75,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
73 75 
                                  ignore_duplicates_p,
74 76 
                                  min_depth, variant_strands, ignore_query_Ns_p,
75 77 
                                  print_indels_p, blocksize, verbosep,
76 
-                                 /*readlevel_p*/false, max_softclip, noncovered);
78 
+                                 /*readlevel_p*/false, max_softclip,
79 
+                                 print_xs_scores_p, noncovered);
77 80 
   IIT_free(&chromosome_iit);
78 81 
   Genome_free(&genome);
79 82
Browse code

add support for the noncovered argument; tells bam_tally to report even those positions that have no coverage

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@89724 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 30/04/2014 17:05:27
Showing 1 changed files
src/bamtally.c
History View file @ 6d5d75f
... ... 
@@ -23,7 +23,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
23 23 
                 SEXP ignore_query_Ns_p_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 
-                SEXP verbosep_R, SEXP max_softclip_R)
26 
+                SEXP verbosep_R, SEXP max_softclip_R, SEXP noncovered_R)
27 27 
 {
28 28 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
29 29 
   const char *genome_dir =
... ... 
@@ -47,6 +47,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
47 47 
   int blocksize = asInteger(blocksize_R);
48 48 
   int verbosep = asLogical(verbosep_R);
49 49 
   int max_softclip = asInteger(max_softclip_R);
50 
+  int noncovered = asLogical(noncovered_R);
50 51
51 52 
   Genome_T genome = createGenome(genome_dir, db);
52 53 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
... ... 
@@ -72,7 +73,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
72 73 
                                  ignore_duplicates_p,
73 74 
                                  min_depth, variant_strands, ignore_query_Ns_p,
74 75 
                                  print_indels_p, blocksize, verbosep,
75 
-                                 /*readlevel_p*/false, max_softclip);
76 
+                                 /*readlevel_p*/false, max_softclip, noncovered);
76 77 
   IIT_free(&chromosome_iit);
77 78 
   Genome_free(&genome);
78 79
Browse code

update GSTRUCT (bam_tally); add include_soft_clip parameter for counting over soft clips of a given max length (more accurate allele frequency)

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@88343 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 02/04/2014 22:03:40
Showing 1 changed files
src/bamtally.c
History View file @ 17db5e7
... ... 
@@ -23,7 +23,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
23 23 
                 SEXP ignore_query_Ns_p_R,
24 24 
                 SEXP print_indels_p_R,
25 25 
                 SEXP blocksize_R,
26 
-                SEXP verbosep_R)
26 
+                SEXP verbosep_R, SEXP max_softclip_R)
27 27 
 {
28 28 
   Bamreader_T bamreader = (Bamreader_T) R_ExternalPtrAddr(bamreader_R);
29 29 
   const char *genome_dir =
... ... 
@@ -46,6 +46,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
46 46 
   bool print_indels_p = asLogical(print_indels_p_R);
47 47 
   int blocksize = asInteger(blocksize_R);
48 48 
   int verbosep = asLogical(verbosep_R);
49 
+  int max_softclip = asInteger(max_softclip_R);
49 50
50 51 
   Genome_T genome = createGenome(genome_dir, db);
51 52 
   IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
... ... 
@@ -71,7 +72,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
71 72 
                                  ignore_duplicates_p,
72 73 
                                  min_depth, variant_strands, ignore_query_Ns_p,
73 74 
                                  print_indels_p, blocksize, verbosep,
74 
-                                 /*readlevel_p*/false);
75 
+                                 /*readlevel_p*/false, max_softclip);
75 76 
   IIT_free(&chromosome_iit);
76 77 
   Genome_free(&genome);
77 78
Browse code

update bam_tally to fix counting issue with some range-restricted queries

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@85510 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 14/01/2014 19:15:57
Showing 1 changed files
src/bamtally.c
History View file @ 3a231c3
... ... 
@@ -70,7 +70,8 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
70 70 
                                  need_unique_p, need_primary_p,
71 71 
                                  ignore_duplicates_p,
72 72 
                                  min_depth, variant_strands, ignore_query_Ns_p,
73 
-                                 print_indels_p, blocksize, verbosep);
73 
+                                 print_indels_p, blocksize, verbosep,
74 
+                                 /*readlevel_p*/false);
74 75 
   IIT_free(&chromosome_iit);
75 76 
   Genome_free(&genome);
76 77
Browse code

update to latest gstruct; brings faster bam_tally (for high coverage regions) and read-group filtering support in bam_tally

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@80611 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 20/09/2013 02:05:42
Showing 1 changed files
src/bamtally.c
History View file @ 4d5cc80
... ... 
@@ -13,7 +13,7 @@
13 13
14 14 
 SEXP
15 15 
 R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
16 
-                SEXP which_R,
16 
+                SEXP which_R, SEXP desired_read_group_R,
17 17 
                 SEXP alloclength_R,
18 18 
                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,
19 19 
                 SEXP maximum_nhits_R,
... ... 
@@ -30,6 +30,9 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
30 30 
     genome_dir_R == R_NilValue ? NULL : CHAR(asChar(genome_dir_R));
31 31 
   const char *db = CHAR(asChar(db_R));
32 32 
   int alloclength = asInteger(alloclength_R);
33 
+  const char *desired_read_group =
34 
+    desired_read_group_R == R_NilValue ? NULL :
35 
+    CHAR(asChar(desired_read_group_R));
33 36 
   int minimum_mapq = asInteger(minimum_mapq_R);
34 37 
   int good_unique_mapq = asInteger(good_unique_mapq_R);
35 38 
   int maximum_nhits = asInteger(maximum_nhits_R);
... ... 
@@ -57,9 +60,12 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
57 60 
     end = asInteger(VECTOR_ELT(which_R, 2));
58 61 
   }
59 62
60 
-  IIT_T tally_iit = Bamtally_iit(bamreader, (char *)chr, start, end,
61 
-                                 genome, chromosome_iit,
62 
-                                 alloclength, minimum_mapq, good_unique_mapq,
63 
+  IIT_T tally_iit = Bamtally_iit(bamreader, (char *)chr,
64 
+                                 /* TODO: bam_lacks_chr */ NULL,
65 
+                                 start, end,
66 
+                                 genome, chromosome_iit, alloclength,
67 
+                                 (char *)desired_read_group,
68 
+                                 minimum_mapq, good_unique_mapq,
63 69 
                                  maximum_nhits, need_concordant_p,
64 70 
                                  need_unique_p, need_primary_p,
65 71 
                                  ignore_duplicates_p,
Browse code

Add <stdlib.h> include to solve undefined off_t issue.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@74948 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 27/03/2013 19:58:06
Showing 1 changed files
src/bamtally.c
History View file @ a2a29ed
... ... 
@@ -1,3 +1,5 @@
1 
+#include <stdlib.h>
2 
+
1 3 
 #include <gstruct/bool.h>
2 4 
 #include <gstruct/iit-read.h>
3 5 
 #include <gstruct/genome.h>
Browse code

Add R_Genome_getSeq for efficiently retrieving sequence from a GMAP genome index. This is super fast.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@74871 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 26/03/2013 21:38:39
Showing 1 changed files
src/bamtally.c
History View file @ 86afb16
... ... 
@@ -1,17 +1,13 @@
1 
-#include <stdlib.h>
2 
-#include <string.h>
3 
-
4 1 
 #include <gstruct/bool.h>
5 
-#include <gstruct/datadir.h>
6 2 
 #include <gstruct/iit-read.h>
7 3 
 #include <gstruct/genome.h>
8 
-#include <gstruct/interval.h>
9 4 
 #include <gstruct/genomicpos.h>
10 5 
 #include <gstruct/bamread.h>
11 6 
 #include <gstruct/bamtally.h>
12 7
13 8 
 #include "gmapR.h"
14 9 
 #include "iit.h"
10 
+#include "genome.h"
15 11
16 12 
 SEXP
17 13 
 R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
... ... 
@@ -20,7 +16,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
20 16 
                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,
21 17 
                 SEXP maximum_nhits_R,
22 18 
                 SEXP need_concordant_p_R, SEXP need_unique_p_R,
23 
-                SEXP need_primary_p_R,
19 
+                SEXP need_primary_p_R, SEXP ignore_duplicates_p_R,
24 20 
                 SEXP min_depth_R, SEXP variant_strands_R,
25 21 
                 SEXP ignore_query_Ns_p_R,
26 22 
                 SEXP print_indels_p_R,
... ... 
@@ -38,30 +34,16 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
38 34 
   bool need_concordant_p = asLogical(need_concordant_p_R);
39 35 
   bool need_unique_p = asLogical(need_unique_p_R);
40 36 
   bool need_primary_p = asLogical(need_primary_p_R);
37 
+  bool ignore_duplicates_p = asLogical(ignore_duplicates_p_R);
41 38 
   int min_depth = asInteger(min_depth_R);
42 39 
   int variant_strands = asInteger(variant_strands_R);
43 40 
   bool ignore_query_Ns_p = asLogical(ignore_query_Ns_p_R);
44 41 
   bool print_indels_p = asLogical(print_indels_p_R);
45 42 
   int blocksize = asInteger(blocksize_R);
46 43 
   int verbosep = asLogical(verbosep_R);
47 
-
48 
-  char *genomesubdir, *fileroot, *dbversion, *iitfile;
49 
-  genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,
50 
-                                           (char *)genome_dir, (char *)db);
51 
-
52 
-  iitfile = (char *) calloc(strlen(genomesubdir) + strlen("/") +
53 
-                            strlen(fileroot) + strlen(".chromosome.iit") + 1,
54 
-                            sizeof(char));
55 
-  sprintf(iitfile, "%s/%s.chromosome.iit", genomesubdir, fileroot);
56 
-  IIT_T chromosome_iit = IIT_read(iitfile, /*name*/NULL, /*readonlyp*/true,
57 
-                                  /*divread*/READ_ALL, /*divstring*/NULL,
58 
-                                  /*add_iit_p*/false, /*labels_read_p*/true);
59 
-  Genome_T genome = Genome_new(genomesubdir, fileroot, /*snps_root*/NULL,
60 
-                               /*uncompressedp*/false, /*access*/USE_MMAP_ONLY);
61 
-  free(iitfile);
62 
-  free(fileroot);
63 
-  free(dbversion);
64 
-  free(genomesubdir);
44 
+
45 
+  Genome_T genome = createGenome(genome_dir, db);
46 
+  IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);
65 47
66 48 
   const char *chr = NULL;
67 49 
   Genomicpos_T start = 0;
... ... 
@@ -78,6 +60,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
78 60 
                                  alloclength, minimum_mapq, good_unique_mapq,
79 61 
                                  maximum_nhits, need_concordant_p,
80 62 
                                  need_unique_p, need_primary_p,
63 
+                                 ignore_duplicates_p,
81 64 
                                  min_depth, variant_strands, ignore_query_Ns_p,
82 65 
                                  print_indels_p, blocksize, verbosep);
83 66 
   IIT_free(&chromosome_iit);
Browse code

Decoupling the bam tallying from summarization. Soon, we will return an IIT, which can be summarized in various ways. Currently, the only way is the variant summary.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@69806 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 25/09/2012 16:00:17
Showing 1 changed files
src/bamtally.c
History View file @ 3f89995
... ... 
@@ -11,233 +11,11 @@
11 11 
 #include <gstruct/bamtally.h>
12 12
13 13 
 #include "gmapR.h"
14 
-
15 
-enum { SEQNAMES, POS, REF, READ, N_CYCLES, N_CYCLES_REF, COUNT, COUNT_REF,
16 
-       COUNT_TOTAL, HIGH_QUALITY, HIGH_QUALITY_REF, HIGH_QUALITY_TOTAL,
17 
-       MEAN_QUALITY, MEAN_QUALITY_REF, COUNT_PLUS, COUNT_PLUS_REF, COUNT_MINUS,
18 
-       COUNT_MINUS_REF, N_BASE_COLS };
19 
-
20 
-typedef struct TallyTable {
21 
-  SEXP seqnames_R;
22 
-  int *pos;
23 
-  SEXP ref_R;
24 
-  SEXP read_R;
25 
-  int *n_cycles;
26 
-  int *n_cycles_ref;
27 
-  int *count;
28 
-  int *count_ref;
29 
-  int *count_total;
30 
-  int *high_quality;
31 
-  int *high_quality_ref;
32 
-  int *high_quality_total;
33 
-  double *mean_quality;
34 
-  double *mean_quality_ref;
35 
-  int *count_plus;
36 
-  int *count_plus_ref;
37 
-  int *count_minus;
38 
-  int *count_minus_ref;
39 
-  SEXP cycle_bins_R;
40 
-} TallyTable;
41 
-
42 
-static inline int
43 
-read_int (unsigned char **bytes) {
44 
-  int x = ((int *)*bytes)[0];
45 
-  *bytes += sizeof(int);
46 
-  return x;
47 
-}
48 
-
49 
-static inline char
50 
-read_char (unsigned char **bytes) {
51 
-  char x = (*bytes)[0];
52 
-  (*bytes)++;
53 
-  return x;
54 
-}
55 
-
56 
-static inline const char *
57 
-read_string (unsigned char **bytes) {
58 
-    const char *string = *bytes;
59 
-    (*bytes) += strlen(string) + 1;
60 
-    return string;
61 
-}
62 
-
63 
-static int
64 
-parse_indel_count(unsigned char *bytes) {
65 
-  int count = read_int(&bytes);
66 
-  return count;
67 
-}
68 
-
69 
-static int
70 
-parse_allele_count(unsigned char *bytes) {
71 
-  int n_alleles = 1; /* always have a reference */
72 
-  bytes += sizeof(int) * 4 + 1; /* skip total and reference */
73 
-  while(bytes[0] != '\0') {
74 
-    bytes += sizeof(int) * 2 + 1;
75 
-    n_alleles++;
76 
-  }
77 
-  return n_alleles;
78 
-}
79 
-
80 
-static void
81 
-read_total_counts(unsigned char **bytes, int row, int *count_total) {
82 
-  int count_total_plus = read_int(bytes);
83 
-  int count_total_minus = read_int(bytes);
84 
-  count_total[row] = count_total_plus + count_total_minus;
85 
-}
86 
-
87 
-static void
88 
-read_cycle_counts(unsigned char **bytes, int row, int *n_cycles,
89 
-                  SEXP cycle_bins_R, SEXP cycle_breaks_R)
90 
-{
91 
-  n_cycles[row] = read_int(bytes);
92 
-  if (cycle_breaks_R != R_NilValue) {
93 
-    for (int index = 0; index < n_cycles[row]; index++) {
94 
-      int cycle = abs(read_int(bytes));
95 
-      int count = read_int(bytes);
96 
-      int bin = 0;
97 
-      while(length(cycle_breaks_R) > bin &&
98 
-            cycle > INTEGER(cycle_breaks_R)[bin])
99 
-        bin++;
100 
-      if (bin > 0 && bin < length(cycle_breaks_R)) {
101 
-        SEXP bin_vector = VECTOR_ELT(cycle_bins_R, bin - 1);
102 
-        INTEGER(bin_vector)[row] += count;
103 
-      }
104 
-    }
105 
-  } else {
106 
-    (*bytes) += n_cycles[row] * 2 * sizeof(int);
107 
-  }
108 
-}
109 
-
110 
-static int
111 
-parse_indels(unsigned char *bytes, int row, SEXP cycle_breaks_R,
112 
-             int high_base_quality, TallyTable tally, bool insertion)
113 
-{
114 
-  int indel_count = read_int(&bytes);
115 
-  for (int indel = 0; indel < indel_count; indel++) {
116 
-    tally.count_plus[row] = read_int(&bytes);
117 
-    tally.count_minus[row] = read_int(&bytes);
118 
-    tally.count[row] = tally.count_plus[row] + tally.count_minus[row];
119 
-    SEXP seq_R = mkChar(read_string(&bytes));
120 
-    if (insertion)
121 
-      SET_STRING_ELT(tally.read_R, row, seq_R);
122 
-    else SET_STRING_ELT(tally.ref_R, row, seq_R);
123 
-    read_cycle_counts(&bytes, row, tally.n_cycles, tally.cycle_bins_R,
124 
-                      cycle_breaks_R);
125 
-    /* stuff not relevant for indels */
126 
-    tally.mean_quality[row] = NA_REAL;
127 
-    tally.high_quality[row] = NA_INTEGER;
128 
-  }
129 
-  return indel_count;
130 
-}
131 
-
132 
-static void
133 
-read_quality_counts(unsigned char **bytes, int row, int *high_quality,
134 
-                    double *mean_quality, int high_base_quality)
135 
-{
136 
-  int n_qualities = read_int(bytes);
137 
-  int total_quality = 0;
138 
-  int total_quality_weight = 0;
139 
-  high_quality[row] = 0;
140 
-  for (int index = 0; index < n_qualities; index++) {
141 
-    int quality = read_int(bytes);
142 
-    int count = read_int(bytes);
143 
-    if (quality > high_base_quality) {
144 
-      total_quality += quality * count;
145 
-      total_quality_weight += count;
146 
-      high_quality[row] += count;
147 
-    }
148 
-  }
149 
-  mean_quality[row] = total_quality_weight > 0 ?
150 
-    (double)total_quality / total_quality_weight : R_NaN;
151 
-}
152 
-
153 
-static int
154 
-read_allele_counts(unsigned char **bytes, int row, SEXP read_R,
155 
-                   int *count_plus, int *count_minus, int *count)
156 
-{
157 
-  int n_alleles = 0;
158 
-  char allele;
159 
-  while((allele = read_char(bytes)) != '\0') {
160 
-    SET_STRING_ELT(read_R, row, mkCharLen(&allele, 1));
161 
-    count_plus[row] = read_int(bytes);
162 
-    count_minus[row] = read_int(bytes);
163 
-    count[row] = count_plus[row] + count_minus[row];
164 
-    row++;
165 
-    n_alleles++;
166 
-  }
167 
-  return n_alleles;
168 
-}
169 
-
170 
-static int
171 
-parse_alleles(unsigned char *bytes, int row, int ref_row, SEXP cycle_breaks_R,
172 
-              int high_base_quality, TallyTable tally)
173 
-{
174 
-  read_total_counts(&bytes, row, tally.count_total);
175 
-  int n_alleles = read_allele_counts(&bytes, row, tally.read_R,
176 
-                                     tally.count_plus, tally.count_minus,
177 
-                                     tally.count);
178 
-  for (int allele = 0; allele < n_alleles; allele++, row++) {
179 
-    tally.n_cycles[row] = 0;
180 
-    for (int b = 0; b < length(cycle_breaks_R) - 1; b++) {
181 
-      INTEGER(VECTOR_ELT(tally.cycle_bins_R, b))[row] = 0;
182 
-    }
183 
-    tally.high_quality[row] = 0;
184 
-    tally.mean_quality[row] = R_NaN;
185 
-    if (tally.count[row] > 0) {
186 
-      read_cycle_counts(&bytes, row, tally.n_cycles, tally.cycle_bins_R,
187 
-                        cycle_breaks_R);
188 
-      read_quality_counts(&bytes, row, tally.high_quality, tally.mean_quality,
189 
-                          high_base_quality);
190 
-    }
191 
-  }
192 
-  int high_quality_total = 0;
193 
-  for (int r = ref_row; r < row; r++) {
194 
-    high_quality_total += tally.high_quality[r];
195 
-  }
196 
-  for (int r = ref_row; r < row; r++) {
197 
-    tally.n_cycles_ref[r] = tally.n_cycles[ref_row];
198 
-    tally.count_total[r] = tally.count_total[ref_row];
199 
-    tally.mean_quality_ref[r] = tally.mean_quality[ref_row];
200 
-    tally.high_quality_ref[r] = tally.high_quality[ref_row];
201 
-    tally.high_quality_total[r] = high_quality_total;
202 
-    tally.count_plus_ref[r] = tally.count_plus[ref_row];
203 
-    tally.count_minus_ref[r] = tally.count_minus[ref_row];
204 
-    tally.count_ref[r] = tally.count_plus_ref[r] + tally.count_minus_ref[r];
205 
-    SET_STRING_ELT(tally.ref_R, r, STRING_ELT(tally.read_R, ref_row));
206 
-  }
207 
-  bool have_ref_row = row > ref_row;
208 
-  if (have_ref_row) {
209 
-    /* clear the 'alt' columns for the 'ref' row with NAs */
210 
-    SET_STRING_ELT(tally.read_R, ref_row, NA_STRING);
211 
-    tally.n_cycles[ref_row] = NA_INTEGER;
212 
-    tally.mean_quality[ref_row] = NA_REAL;
213 
-    tally.high_quality[ref_row] = NA_REAL;
214 
-    tally.count_plus[ref_row] = NA_INTEGER;
215 
-    tally.count_minus[ref_row] = NA_INTEGER;
216 
-    tally.count[ref_row] = NA_INTEGER;
217 
-  }
218 
-  return n_alleles;
219 
-}
220 
-
221 
-
222 
-/* FORMAT
223 
-
224 
-   block header:
225 
-   [0][ins][del][mm]
226 
-
227 
-   mismatches for each position:
228 
-   [t+][t-][rn][r+][r-]([an][a+][a-])*[\0]([c#]([cv][cc])*)*([q#]([qv][qc])*)*
229 
-
230 
-   insertions:
231 
-   [i#]([seq][t+][t-][c#]([cv][cc])*)*
232 
-
233 
-   deletions:
234 
-   [d#]([seq][t+][t-][c#]([cv][cc])*)*
235 
-*/
14 
+#include "iit.h"
236 15
237 16 
 SEXP
238 17 
 R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
239 18 
                 SEXP which_R,
240 
-                SEXP cycle_breaks_R, SEXP high_base_quality_R,
241 19 
                 SEXP alloclength_R,
242 20 
                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,
243 21 
                 SEXP maximum_nhits_R,
... ... 
@@ -266,7 +44,6 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
266 44 
   bool print_indels_p = asLogical(print_indels_p_R);
267 45 
   int blocksize = asInteger(blocksize_R);
268 46 
   int verbosep = asLogical(verbosep_R);
269 
-  int high_base_quality = asInteger(high_base_quality_R);
270 47
271 48 
   char *genomesubdir, *fileroot, *dbversion, *iitfile;
272 49 
   genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,
... ... 
@@ -310,119 +87,5 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,
310 87 
     error("Could not create tally\n");
311 88 
   }
312 89
313 
-  int n_cycle_bins =
314 
-    length(cycle_breaks_R) == 0 ? 0 : length(cycle_breaks_R) - 1;
315 
-
316 
-  int n_rows = 0;
317 
-  /* loop over the IIT, getting the total number of rows
318 
-     this is (num alts + 1) for every position. */
319 
-  for (int index = 1; index <= IIT_total_nintervals(tally_iit); index++) {
320 
-    unsigned char *bytes = IIT_data(tally_iit, index);
321 
-    int width = IIT_length(tally_iit, index);
322 
-    unsigned char *base = bytes + (3 * width + 1) * sizeof(int);
323 
-    for (int pos = 0; pos < width; pos++) {
324 
-      int insertion_offset = read_int(&bytes);
325 
-      int deletion_offset = read_int(&bytes);
326 
-      int allele_offset = read_int(&bytes);
327 
-      int next_offset = read_int(&bytes);
328 
-      bytes -= 4; /* rewind from read-ahead */
329 
-      if (deletion_offset - insertion_offset > 0) {
330 
-        n_rows += parse_indel_count(base + insertion_offset);
331 
-      }
332 
-      if (allele_offset - deletion_offset > 0) {
333 
-        n_rows += parse_indel_count(base + deletion_offset);
334 
-      }
335 
-      if (next_offset - allele_offset > 0) {
336 
-        n_rows += parse_allele_count(base + allele_offset);
337 
-      }
338 
-    }
339 
-  }
340 
-
341 
-  SEXP tally_R; /* the result list */
342 
-  PROTECT(tally_R =
343 
-          allocVector(VECSXP, N_BASE_COLS + n_cycle_bins));
344 
-
345 
-  TallyTable tally;
346 
-
347 
-  SET_VECTOR_ELT(tally_R, SEQNAMES, allocVector(STRSXP, n_rows));
348 
-  tally.seqnames_R = VECTOR_ELT(tally_R, SEQNAMES);
349 
-  SET_VECTOR_ELT(tally_R, POS, allocVector(INTSXP, n_rows));
350 
-  tally.pos = INTEGER(VECTOR_ELT(tally_R, POS));
351 
-  SET_VECTOR_ELT(tally_R, REF, allocVector(STRSXP, n_rows));
352 
-  tally.ref_R = VECTOR_ELT(tally_R, REF);
353 
-  SET_VECTOR_ELT(tally_R, READ, allocVector(STRSXP, n_rows));
354 
-  tally.read_R = VECTOR_ELT(tally_R, READ);
355 
-  SET_VECTOR_ELT(tally_R, N_CYCLES, allocVector(INTSXP, n_rows));
356 
-  tally.n_cycles = INTEGER(VECTOR_ELT(tally_R, N_CYCLES));
357 
-  SET_VECTOR_ELT(tally_R, N_CYCLES_REF, allocVector(INTSXP, n_rows));
358 
-  tally.n_cycles_ref = INTEGER(VECTOR_ELT(tally_R, N_CYCLES_REF));
359 
-  SET_VECTOR_ELT(tally_R, COUNT, allocVector(INTSXP, n_rows));
360 
-  tally.count = INTEGER(VECTOR_ELT(tally_R, COUNT));
361 
-  SET_VECTOR_ELT(tally_R, COUNT_REF, allocVector(INTSXP, n_rows));
362 
-  tally.count_ref = INTEGER(VECTOR_ELT(tally_R, COUNT_REF));