Make a generic refSeq-based annotation file for CEAS program using fetchChromSizes and build_genomeBG. The resulting annotation can be used for test purposes but shouldn't be used for genuine analyses.

make_ceas_refseq_annotation.py

#!/usr/bin/env python
#
# Make generic CEAS annotation file for arbitrary genome build
#
# Needs the UCSC fetchChromSizes program plus build_genomeBG from CEAS
#
# Usage: make_ceas_refseq_annotation.py <GENOME>
#
import optparse
import os
import subprocess
import tempfile

if __name__ == "__main__":

    # Process command line
    p = optparse.OptionParser()
    options,args = p.parse_args()
    if len(args) != 1:
        p.error("Need to supply genome build name e.g. mm9")
    genome_build = args[0]
    print "Genome build: %s" % genome_build
    
    #working_dir = os.getcwd()
    working_dir = tempfile.mkdtemp()
    #genome_build = "mm10"
    chrom_sizes = os.path.join(working_dir,"%s.len" % genome_build)
    stderr_file = os.path.join(working_dir,"fetchChromSizes.stderr")
    # Fetch chrom sizes
    cmd = "fetchChromSizes %s" % genome_build
    print cmd
    proc = subprocess.Popen(args=cmd,shell=True,cwd=working_dir,
                            stdout=open(chrom_sizes,'wb'),
                            stderr=open(stderr_file,'wb'))
    proc.wait()
    # Make wig file
    chrom_wig = os.path.join(working_dir,"%s.wig" % genome_build)
    fp = open(chrom_wig,'w')
    fp.write("track type=wiggle_0\n")
    for line in open(chrom_sizes,'r'):
        chrom,size = line.strip().split()
        fp.write("fixedStep chrom=%s start=1 step=%s\n1\n" % (chrom,size))
    fp.close()
    # Generate the annotation file
    cwd = os.getcwd()
    annotation_file = os.path.join(cwd,"%s.refGene" % genome_build)
    cmd = "build_genomeBG -d %s -g refGene -w %s -o %s" % \
          (genome_build,chrom_wig,annotation_file)
    print cmd
    stderr_file = os.path.join(cwd,"build_genomeBG.stderr")
    proc = subprocess.Popen(args=cmd,shell=True,cwd=cwd,
                            stderr=open(stderr_file,'wb'))
    proc.wait()
    print "Done: %s" % annotation_file