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;

+}

+