@@ -2,7 +2,7 @@ SUBDIRS = gmap gstruct

 .PHONY: all clean $(SUBDIRS)

-OBJECTS = bamreader.o bamtally.o iit.o variantsummary.o R_init_gmapR.o

+OBJECTS = bamreader.o bamtally.o iit.o variantsummary.o genome.o R_init_gmapR.o

 R_SRC_DIR = ${CURDIR}

 PREFIX = ${R_SRC_DIR}/../inst/usr

@@ -7,11 +7,14 @@

 static const R_CallMethodDef callMethods[] = {

   /* bamtally.c */

-  CALLMETHOD_DEF(R_Bamtally_iit, 17),

+  CALLMETHOD_DEF(R_Bamtally_iit, 18),

   CALLMETHOD_DEF(R_tally_iit_parse, 4),

   /* bamreader.c */

   CALLMETHOD_DEF(R_Bamread_new, 1),

+

+  /* genome.c */

+  CALLMETHOD_DEF(R_Genome_getSeq, 6),

   {NULL, NULL, 0}

 };

@@ -1,17 +1,13 @@

-#include <stdlib.h>

-#include <string.h>

-

 #include <gstruct/bool.h>

-#include <gstruct/datadir.h>

 #include <gstruct/iit-read.h>

 #include <gstruct/genome.h>

-#include <gstruct/interval.h>

 #include <gstruct/genomicpos.h>

 #include <gstruct/bamread.h>

 #include <gstruct/bamtally.h>

 #include "gmapR.h"

 #include "iit.h"

+#include "genome.h"

 SEXP

 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,

                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,

                 SEXP maximum_nhits_R,

                 SEXP need_concordant_p_R, SEXP need_unique_p_R,

-                SEXP need_primary_p_R,

+                SEXP need_primary_p_R, SEXP ignore_duplicates_p_R,

                 SEXP min_depth_R, SEXP variant_strands_R,

                 SEXP ignore_query_Ns_p_R,

                 SEXP print_indels_p_R,

@@ -38,30 +34,16 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

   bool need_concordant_p = asLogical(need_concordant_p_R);

   bool need_unique_p = asLogical(need_unique_p_R);

   bool need_primary_p = asLogical(need_primary_p_R);

+  bool ignore_duplicates_p = asLogical(ignore_duplicates_p_R);

   int min_depth = asInteger(min_depth_R);

   int variant_strands = asInteger(variant_strands_R);

   bool ignore_query_Ns_p = asLogical(ignore_query_Ns_p_R);

   bool print_indels_p = asLogical(print_indels_p_R);

   int blocksize = asInteger(blocksize_R);

   int verbosep = asLogical(verbosep_R);

-  

-  char *genomesubdir, *fileroot, *dbversion, *iitfile;

-  genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,

-                                           (char *)genome_dir, (char *)db);

-  

-  iitfile = (char *) calloc(strlen(genomesubdir) + strlen("/") +

-                            strlen(fileroot) + strlen(".chromosome.iit") + 1,

-                            sizeof(char));

-  sprintf(iitfile, "%s/%s.chromosome.iit", genomesubdir, fileroot);

-  IIT_T chromosome_iit = IIT_read(iitfile, /*name*/NULL, /*readonlyp*/true,

-                                  /*divread*/READ_ALL, /*divstring*/NULL,

-                                  /*add_iit_p*/false, /*labels_read_p*/true);

-  Genome_T genome = Genome_new(genomesubdir, fileroot, /*snps_root*/NULL,

-                               /*uncompressedp*/false, /*access*/USE_MMAP_ONLY);

-  free(iitfile);

-  free(fileroot);

-  free(dbversion);

-  free(genomesubdir);

+

+  Genome_T genome = createGenome(genome_dir, db);

+  IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);

   const char *chr = NULL;

   Genomicpos_T start = 0;

@@ -78,6 +60,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

                                  alloclength, minimum_mapq, good_unique_mapq,

                                  maximum_nhits, need_concordant_p,

                                  need_unique_p, need_primary_p,

+                                 ignore_duplicates_p,

                                  min_depth, variant_strands, ignore_query_Ns_p,

                                  print_indels_p, blocksize, verbosep);

   IIT_free(&chromosome_iit);

new file mode 100644

@@ -0,0 +1,116 @@

+#include <stdlib.h>

+#include <string.h>

+

+#include <gstruct/bool.h>

+#include <gstruct/datadir.h>

+#include <gstruct/iit-read.h>

+#include <gstruct/genome.h>

+#include <gstruct/interval.h>

+#include <gstruct/genomicpos.h>

+#include <gstruct/complement.h>

+

+#include "gmapR.h"

+#include "genome.h"

+

+IIT_T

+readChromosomeIIT (const char *genome_dir, const char *db) {

+  char *genomesubdir, *fileroot, *dbversion, *iitfile;

+  genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,

+                                           (char *)genome_dir, (char *)db);

+  

+  iitfile = (char *) calloc(strlen(genomesubdir) + strlen("/") +

+                            strlen(fileroot) + strlen(".chromosome.iit") + 1,

+                            sizeof(char));

+  sprintf(iitfile, "%s/%s.chromosome.iit", genomesubdir, fileroot);

+  IIT_T chromosome_iit = IIT_read(iitfile, /*name*/NULL, /*readonlyp*/true,

+                                  /*divread*/READ_ALL, /*divstring*/NULL,

+                                  /*add_iit_p*/false, /*labels_read_p*/true);

+  free(iitfile);

+  free(fileroot);

+  free(dbversion);

+  free(genomesubdir);

+

+  return chromosome_iit;

+}

+

+Genome_T

+createGenome (const char *genome_dir, const char *db) {

+  char *genomesubdir, *fileroot, *dbversion;

+  genomesubdir = Datadir_find_genomesubdir(&fileroot, &dbversion,

+                                           (char *)genome_dir, (char *)db);

+  

+  Genome_T genome = Genome_new(genomesubdir, fileroot, /*snps_root*/NULL,

+                               /*uncompressedp*/false, /*access*/USE_MMAP_ONLY);

+  free(fileroot);

+  free(dbversion);

+  free(genomesubdir);

+

+  return genome;

+}

+

+static char complCode[128] = COMPLEMENT_LC;

+

+static void

+make_complement_inplace (char *sequence, Genomicpos_T length) {

+  char temp;

+  int i, j;

+

+  for (i = 0, j = length-1; i < length/2; i++, j--) {

+    temp = complCode[(int) sequence[i]];

+    sequence[i] = complCode[(int) sequence[j]];

+    sequence[j] = temp;

+  }

+  if (i == j) {

+    sequence[i] = complCode[(int) sequence[i]];

+  }

+

+  return;

+}

+

+SEXP

+R_Genome_getSeq (SEXP genome_dir_R, SEXP db_R,

+                 SEXP seqnames_R, SEXP start_R, SEXP width_R, SEXP strand_R)

+{

+  const char *genome_dir =

+    genome_dir_R == R_NilValue ? NULL : CHAR(asChar(genome_dir_R));

+  const char *db = CHAR(asChar(db_R));

+  Genome_T genome = createGenome(genome_dir, db);

+  IIT_T chromosome_iit = readChromosomeIIT(genome_dir, db);

+  int *start = INTEGER(start_R);

+  int *width = INTEGER(width_R);

+

+  int max_width = 0;

+  for (int i = 0; i < length(width_R); i++) {

+    if (width[i] > max_width)

+      max_width = width[i];

+  }

+  char *buffer = (char *) R_alloc(sizeof(char), max_width + 1);

+

+  int seqname_index;

+  const char *old_seqname = NULL; // premature optimization

+  SEXP result;

+  PROTECT(result = allocVector(STRSXP, length(start_R)));

+  for (int i = 0; i < length(start_R); i++) {

+    const char *seqname = CHAR(STRING_ELT(seqnames_R, i));

+    if (old_seqname == NULL || strcmp(seqname, old_seqname)) {

+      seqname_index = IIT_find_linear(chromosome_iit, (char *)seqname);

+      if (seqname_index < 0) {

+        error("Cannot find chromosome %s in genome", seqname);

+      }

+    }

+    old_seqname = seqname;

+    int chroffset = Interval_low(IIT_interval(chromosome_iit, seqname_index));

+    Genome_fill_buffer_simple(genome, chroffset + start[i] - 1, width[i],

+                              buffer);

+    if (CHAR(STRING_ELT(strand_R, i))[0] == '-')

+      make_complement_inplace(buffer, width[i]);

+    SET_STRING_ELT(result, i, mkChar(buffer));

+  }

+

+  IIT_free(&chromosome_iit);

+  Genome_free(&genome);

+

+  UNPROTECT(1);

+  return result;

+}

+

new file mode 100644

@@ -0,0 +1,12 @@

+#ifndef GENOME_H

+#define GENOME_H

+

+#include <gstruct/iit-read.h>

+#include <gstruct/genome.h>

+

+/* Internal genome helper functions */

+

+IIT_T readChromosomeIIT (const char *genome_dir, const char *db);

+Genome_T createGenome (const char *genome_dir, const char *db);

+

+#endif

@@ -15,7 +15,7 @@ R_Bamtally_iit (SEXP bamreader_R, SEXP genome_dir_R, SEXP db_R,

                 SEXP minimum_mapq_R, SEXP good_unique_mapq_R,

                 SEXP maximum_nhits_R,

                 SEXP need_concordant_p_R, SEXP need_unique_p_R,

-                SEXP need_primary_p_R,

+                SEXP need_primary_p_R, SEXP ignore_duplicates_p_R,

                 SEXP min_depth_R, SEXP variant_strands_R,

                 SEXP ignore_query_Ns_p,

                 SEXP print_indels_p_R,

@@ -26,4 +26,8 @@ SEXP

 R_tally_iit_parse(SEXP tally_iit_R, SEXP cycle_breaks_R,

                   SEXP high_base_quality, SEXP which_R);

+SEXP

+R_Genome_getSeq (SEXP genome_dir_R, SEXP db_R,

+                 SEXP seqnames_R, SEXP start_R, SEXP width_R, SEXP strand_R);

+

 #endif