fixed gc background

Author: simonvh

gimmemotifs/background.py

@@ -177,8 +177,16 @@ def _weighted_random(self, l):
 
 def matched_gc_bedfile(bedfile, matchfile, genome, number):
     N_FRACTION = 0.1
-    genome_size = "/data/genomes/{0}/bwa/{0}.fa.sizes".format(genome)
-    genome_fa = "/data/genomes/{0}/bwa/{0}.fa".format(genome)
+    
+    config = MotifConfig()
+    index = os.path.join(config.get_index_dir(), genome)
+    
+    genome_size = os.path.join(index, "genome.size")
+    genome_fa = os.path.join(index, "genome.fa")
+
+    if not os.path.exists(genome_size) or not os.path.exists(genome_fa):
+        raise RuntimeError, "genome files not found, please re-index {} "  \
+                "with a recent version of gimme index".format(genome)
 
     try:
         fa = Fasta(matchfile)
@@ -222,7 +230,7 @@ def matched_gc_bedfile(bedfile, matchfile, genome, number):
     #sys.stderr.write("{}\n".format(number))
 
     r = rnd.random(l=length, n=number * 15, g=genome_size).nucleotide_content(fi=genome_fa)
-    #sys.stderr.write("Retrieving\n")
+   #sys.stderr.write("Retrieving\n")
     features = [f[:3] + [float(f[7])] for f in r if float(f[12]) <= length * N_FRACTION]
     gc = [f[3] for f in features]
 

gimmemotifs/commands/index.py

@@ -1,23 +1,49 @@
 #!/usr/bin/env python
-# Copyright (c) 2009-2013 Simon van Heeringen <s.vanheeringen@ncmls.ru.nl>
+# Copyright (c) 2009-2015 Simon van Heeringen <s.vanheeringen@ncmls.ru.nl>
 #
 # This module is free software. You can redistribute it and/or modify it under 
 # the terms of the MIT License, see the file COPYING included with this 
 # distribution.
 
 import sys
 import os
+from tempfile import NamedTemporaryFile
+
+import pybedtools
 
 from gimmemotifs.genome_index import GenomeIndex
+from gimmemotifs.shutils import find_by_ext
+from gimmemotifs.config import FASTA_EXT
+
+
 
 def index(args):
 
     if not os.path.exists(args.indexdir):
-    	print "Index_dir %s does not exist!" % (args.indexdir)
-    	sys.exit(1)
+        print "Index_dir %s does not exist!" % (args.indexdir)
+        sys.exit(1)
 
     fasta_dir = args.fastadir
     index_dir = os.path.join(args.indexdir, args.indexname)
 
     g = GenomeIndex()
-    g = g.create_index(fasta_dir, index_dir)
+    g.create_index(fasta_dir, index_dir)
+
+    # Create genome FASTA file for use with bedtools
+    with open(os.path.join(index_dir, "genome.fa"), 'w') as out:
+        for f in find_by_ext(fasta_dir, FASTA_EXT):
+            for line in open(f):
+                out.write(line)
+
+    test_chr = g.get_chromosomes()[0]
+    tmp = NamedTemporaryFile()
+    tmp.write("{}\t1\t2\n".format(test_chr))
+    tmp.flush()
+    
+    b = pybedtools.BedTool(tmp.name)
+    try:
+        b.nucleotide_content(fi=os.path.join(index_dir, "genome.fa"))
+    except pybedtools.helpers.BEDToolsError as e:
+        if str(e).find("generating") == -1:
+            raise
+

gimmemotifs/config.py

@@ -15,6 +15,8 @@
 FA_VALID_BGS = ["random", "promoter", "gc", "user"]
 BED_VALID_BGS = ["random", "genomic", "gc", "promoter", "user"]
 BG_RANK = {"user":1, "promoter":2, "gc":3, "random":4, "genomic":5}
+FASTA_EXT = [".fasta", ".fa", ".fsa"]
+
 
 class MotifConfig:
     __shared_state = {}