Organism factory method.

        Appends organisms to the organisms list.

        Args
            genome_list (list)
            genome_dir (string)
500121
def make_organisms(self, genome_list, genome_dir): 

        '''

            '''

        for genome in genome_list:

            genome_path = genome_dir + genome
            handle = open(genome_path, "rU")
            print 'Adding organism attributes for', genome

            try:
                seq_record = SeqIO.read(handle, "genbank")

                self.organisms.append(Organism(seq_record, genome_path, self))

                del(seq_record)

            except ValueError,e:
                print genome, str(e)

            except AssertionError,e:
                print genome, str(e)

            except UnboundLocalError,e:
                print genome, str(e)

            handle.close()
108274
def install(name, provider, genome_dir, localname, mask, regex, match, annotation):
    """"""
    genomepy.install_genome(
            name, provider, genome_dir=genome_dir, localname=localname, mask=mask, 
            regex=regex, invert_match=not(match), annotation=annotation)
669003
def print_genome_matrix(hits, fastas, id2desc, file_name):
    
    out = open(file_name, 'w')
    fastas = sorted(fastas)
    print('## percent identity between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                average = '-'
            else:
                average = numpy.average([hits[fasta][other][i][3] for i in hits[fasta][other]])
            line.append(str(average))
        print('\t'.join(line), file=out)
    print('', file=out)
    print('## percent of orfs that are orthologous between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                percent = '-'
            else:
                orthologs = float(len(hits[fasta][other]))
                orfs = float(len([i for i in id2desc if id2desc[i][0] == fasta]))
                percent = float(orthologs / orfs) * 100
            line.append(str(percent))
        print('\t'.join(line), file=out)
108800
def append_fresh_table(self, fresh_table):
        
        if fresh_table.name:
            elements = []
            if fresh_table.is_array:
                elements += [element_factory.create_array_of_tables_header_element(fresh_table.name)]
            else:
                elements += [element_factory.create_table_header_element(fresh_table.name)]

            elements += [fresh_table, element_factory.create_newline_element()]
            self.append_elements(elements)

        else:
            # It's an anonymous table
            self.prepend_elements([fresh_table, element_factory.create_newline_element()])
1024956
def createBamHeader(self, baseHeader):
        
        header = dict(baseHeader)
        newSequences = []
        for index, referenceInfo in enumerate(header['SQ']):
            if index < self.numChromosomes:
                referenceName = referenceInfo['SN']
                # The sequence dictionary in the BAM file has to match up
                # with the sequence ids in the data, so we must be sure
                # that these still match up.
                assert referenceName == self.chromosomes[index]
                newReferenceInfo = {
                    'AS': self.referenceSetName,
                    'SN': referenceName,
                    'LN': 0,  # FIXME
                    'UR': 'http://example.com',
                    'M5': 'dbb6e8ece0b5de29da56601613007c2a',  # FIXME
                    'SP': 'Human'
                }
                newSequences.append(newReferenceInfo)
        header['SQ'] = newSequences
        return header
385938
def append_to_list(self, source, start=None, hasIndex=False):
        '''
        '''
        nfy = Numberify()
        try:
            if start is None:
                if type(source) is str:
                    if hasIndex is True:
                        newList = self.get_from_indexedFile(source)
                else:
                    newList = nfy.numberify_data(source,
                                                 len(self.nameDict) + 1)
            else:
                newList = nfy.numberify_data(source, start)
            self.nameDict = dict(self.nameDict.items() + newList.items())
            self.totalCount = len(self.nameDict)
        except:
            print 'Unknown error:', exc_info()[0]
302659
def add_hits_to_proteins(self, hmm_hit_list):

        ''''''

        for org in self.organisms:
            print "adding SearchIO hit objects for", org.accession

            for hit in hmm_hit_list:
                hit_org_id = hit.id.split(',')[0]
                hit_prot_id = hit.id.split(',')[1]

                if org.accession == hit_org_id:
                    for prot in org.proteins:
                        if prot.accession == hit_prot_id:
                            prot.hmm_hit_list.append(hit)
310936
def snap(args):
    
    p = OptionParser(snap.__doc__)
    p.set_home("maker")
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    species, gffile, fastafile = args
    mhome = opts.maker_home
    snapdir = "snap"
    mkdir(snapdir)

    cwd = os.getcwd()
    os.chdir(snapdir)

    newgffile = "training.gff3"
    logging.debug("Construct GFF file combined with sequence ...")
    sh("cat ../{0} > {1}".format(gffile, newgffile))
    sh('echo "##FASTA" >> {0}'.format(newgffile))
    sh("cat ../{0} >> {1}".format(fastafile, newgffile))

    logging.debug("Make models ...")
    sh("{0}/src/bin/maker2zff training.gff3".format(mhome))
    sh("{0}/exe/snap/fathom -categorize 1000 genome.ann genome.dna".format(mhome))
    sh("{0}/exe/snap/fathom -export 1000 -plus uni.ann uni.dna".format(mhome))
    sh("{0}/exe/snap/forge export.ann export.dna".format(mhome))
    sh("{0}/exe/snap/hmm-assembler.pl {1} . > {1}.hmm".format(mhome, species))

    os.chdir(cwd)
    logging.debug("SNAP matrix written to `{0}/{1}.hmm`".format(snapdir, species))
966991
def main(args):

  
  sys.stderr.write("reading reference genome\n")
  ref = FASTAData(open(args.genome).read())
  sys.stderr.write("read in "+str(len(ref.keys()))+" chromosomes\n")

  """ Next make the transcriptome """
  txome = {}
  sys.stderr.write("write the transcriptome\n")
  inf = None
  if args.gpd[-3:] == '.gz':
    inf = gzip.open(args.gpd)
  else:
    inf = open(args.gpd)
  stream = GPDStream(inf)
  tof = open(args.tempdir+'/transcriptome.fa','w')
  z = 0
  for gpd in stream:
    z += 1
    if gpd.transcript_name in txome:
      sys.stderr.write("WARNING already have a transcript "+gpd.transcript_name+" ignoring line "+str(z)+" of the gpd\n")
      continue
    txome[gpd.transcript_name] = gpd.gene_name
    tof.write('>'+gpd.transcript_name+"\n"+str(gpd.get_sequence(ref))+"\n")
  tof.close()
  inf.close()
  sys.stderr.write("wrote "+str(len(txome.keys()))+" transcripts\n")

  """Build the salmon index"""
  sys.stderr.write("building a salmon index\n")
  cmd = 'salmon index -p '+str(args.numThreads)+' -t '+args.tempdir+'/transcriptome.fa -i '+args.tempdir+'/salmon_index'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished building the index\n")

  """Use the index to quanitfy"""
  sys.stderr.write("quanitfy reads\n")
  reads = ''
  if args.rU:
    reads = '-r '+args.rU
  else:
    reads = '-1 '+args.r1+' -2 '+args.r2
  cmd = 'salmon quant -p '+str(args.numThreads)+' -i '+args.tempdir+'/salmon_index -l A '+reads+' -o '+args.tempdir+'/output_quant'
  p = Popen(cmd.split())
  p.communicate()
  sys.stderr.write("finished quanitfying\n")

  """Now parse the salmon output to add gene name"""
  salmon = {}
  with open(args.tempdir+'/output_quant/quant.sf') as inf:
    header = inf.readline()
    for line in inf:
      f = line.rstrip().split("\t")
      # by each transcript name hold a data strcture of the information
      salmon[f[0]] = {'name':f[0],'length':int(f[1]),'EffectiveLength':float(f[2]),'TPM':float(f[3]),'NumReads':float(f[4])}
  genes = {}
  for name in salmon:
    gene = txome[name]
    if gene not in genes: genes[gene] = []
    genes[gene].append(salmon[name])
  genetot = {}
  for gene in genes:
    tot = sum([x['TPM'] for x in genes[gene]])
    genetot[gene] = tot
  ordered_gene_names = sorted(genetot.keys(), key=lambda x: genetot[x],reverse=True)
      
  """Collected enough information to make output"""
  sys.stderr.write("generating output\n")
  of = sys.stdout
  if args.output != '-':
    of = open(args.output,'w')
  of.write("geneName\ttranscriptName\tlength\tEffectiveLength\ttxTPM\tNumReads\tgeneTPM\n")
  for gene in ordered_gene_names:
    txs = sorted(genes[gene],key=lambda x: x['TPM'],reverse=True)
    for tx in txs:
      of.write(gene+"\t"+tx['name']+"\t"+str(tx['length'])+"\t"+str(tx['EffectiveLength'])+"\t"+str(tx['TPM'])+"\t"+str(tx['NumReads'])+"\t"+str(genetot[gene])+"\n")
  of.close()

  # Temporary working directory step 3 of 3 - Cleanup
  if not args.specific_tempdir:
    rmtree(args.tempdir)
590422
def split(*items):
    
    out = []
    for data in [x[0] for x in items]:
        dis_orgs = data["config"]["algorithm"].get("disambiguate")
        if dis_orgs:
            if not data.get("disambiguate", None):
                data["disambiguate"] = {"genome_build": data["genome_build"],
                                        "base": True}
            out.append([data])
            # handle the instance where a single organism is disambiguated
            if isinstance(dis_orgs, six.string_types):
                dis_orgs = [dis_orgs]
            for dis_org in dis_orgs:
                dis_data = copy.deepcopy(data)
                dis_data["disambiguate"] = {"genome_build": dis_org}
                dis_data["genome_build"] = dis_org
                dis_data["config"]["algorithm"]["effects"] = False
                dis_data = run_info.add_reference_resources(dis_data)
                out.append([dis_data])
        else:
            out.append([data])
    return out
875101
def add_individual(self, genotype):
        
        logger.debug("Adding genotype {0} to variant {1}".format(
            genotype, self['variant_id']))
        self['individuals'].append(genotype)
1789
def run(self, args, pipeline_command):
        

        # output that must be moved but not renamed
        consistentNaming = ['alignments/normal_dna_fix_pg_sorted.bam',
                           'alignments/normal_dna_fix_pg_sorted.bam.bai',
                           'alignments/rna_genome_sorted.bam',
                           'alignments/rna_genome_sorted.bam.bai',
                           'alignments/rna_transcriptome.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam',
                           'alignments/tumor_dna_fix_pg_sorted.bam.bai',
                           'mutations/merged/all_merged.vcf',
                           'rankboost/mhcii_rankboost_concise_results.tsv',
                           'rankboost/mhci_rankboost_concise_results.tsv',
                           ]

        # output that must be renamed as well as moved
        # map of the original name to the final name
        renamingNeeded = {'binding_predictions': 'binding_predictions.tar',
                         'expression': 'expression.tar',
                         'haplotyping': 'haplotyping.tar',
                         'peptides': 'peptides.tar',
                         'rankboost': 'rankboost.tar',
                         'reports': 'reports.tar',
                         'mutations/snpeffed/mutations.vcf': 'all_snpeffed.vcf',
                         'mutations/transgened/mutations.vcf': 'all_transgened.vcf',
                         'mutations/merged': 'merged_perchrom.tar',
                         'mutations/muse': 'muse_perchrom.tar',
                         'mutations/mutect': 'mutect_perchrom.tar',
                         'mutations/radia': 'radia_perchrom.tar',
                         'mutations/somaticsniper': 'somaticsniper_perchrom.tar',
                         'mutations/strelka/snv': 'strelka_snv_perchrom.tar',
                         'mutations/strelka/indel': 'strelka_indel_perchrom.tar'}

        def make_output(output_dir, source_dir):
            """
            :param output_dir: dir to write the output to
            :param source_dir: dir containing the directory structure to be parsed
            :return:
            """
            def make_tar(dir, tar):
                with tarfile.open(tar, "w:gz") as tar:
                    tar.add(dir)
            # the output dir is where the real output directories are written
            protect_outputs = os.listdir(source_dir)
            for protectOut in protect_outputs:
                def getName(fileName):
                    return os.path.join(os.path.join(source_dir, protectOut), fileName)
                # move individual files out
                for fileName in consistentNaming:
                    shutil.copyfile(getName(fileName), os.path.join(output_dir, os.path.basename(fileName)))
                for src, dst in renamingNeeded.iteritems():
                    if dst.endswith('.tar'):
                        make_tar(getName(src), os.path.join(output_dir, dst))
                    else:
                        shutil.copyfile(getName(src), os.path.join(output_dir, dst))
            shutil.rmtree(source_dir)


        # prepare workdir
        mount = self._prepare_mount(args)
        self._workdir = os.path.join(mount, 'Toil-' + self._name)
        # insure the pairs are in the same directory, as protect expects
        # This is made more complicated by the fact CWLTool mounts inputs into random, read-only dirs
        # to get around this we copy all inputs into their own directories that we own
        tumor_dna_dir = os.path.expanduser('~/tumorDNA')
        tumor_rna_dir = os.path.expanduser('~/tumorRNA')
        normal_dna_dir = os.path.expanduser('~/normalDNA')
        os.mkdir(tumor_dna_dir)
        os.mkdir(tumor_rna_dir)
        os.mkdir(normal_dna_dir)
        shutil.copy(args.tumor_dna, tumor_dna_dir)
        shutil.copy(args.tumor_rna, tumor_rna_dir)
        shutil.copy(args.normal_dna, normal_dna_dir)
        shutil.copy(args.tumor_dna2, tumor_dna_dir)
        shutil.copy(args.tumor_rna2, tumor_rna_dir)
        shutil.copy(args.normal_dna2, normal_dna_dir)
        args.tumor_dna = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna))
        args.tumor_dna2 = os.path.join(tumor_dna_dir, os.path.basename(args.tumor_dna2))
        args.tumor_rna = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna))
        args.tumor_rna2 = os.path.join(tumor_rna_dir, os.path.basename(args.tumor_rna2))
        args.normal_dna = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna))
        args.normal_dna2 = os.path.join(normal_dna_dir, os.path.basename(args.normal_dna2))

        # prepare config
        args_dict = vars(args)
        args_dict['output_dir'] = mount
        self._config = textwrap.dedent(self._config.format(**args_dict))
        self._sample_name = args_dict["sample_name"]
        config_path = os.path.join(self._workdir, 'config')
        command = self._make_prefix(os.path.join(self._workdir, 'jobStore'),
                                          config_path,
                                          self._workdir) + pipeline_command
        if self._resume and args.resume:
            command.append('--restart')
        self._create_workdir(args)
        with open(config_path, 'w') as f:
            f.write(self._config)

        try:
            subprocess.check_call(command)
        except subprocess.CalledProcessError as e:
            print(e, file=sys.stderr)
        finally:
            log.info('Pipeline terminated, changing ownership of output files from root to user.')
            stat = os.stat(self._mount)
            subprocess.check_call(['chown', '-R', '{}:{}'.format(stat.st_uid, stat.st_gid),
                                   self._mount])
            make_output(self._mount, os.path.join(self._mount, 'output'))
            if self._no_clean and args.no_clean:
                log.info('Flag "--no-clean" was used, therefore %s was not deleted.', self._workdir)
            else:
                log.info('Cleaning up temporary directory: %s', self._workdir)
                shutil.rmtree(self._workdir)
1298476
def deleteGenome(species, name) :
    """"""

    printf('deleting genome (%s, %s)...' % (species, name))

    conf.db.beginTransaction()
    objs = []
    allGood = True
    try :
        genome = Genome_Raba(name = name, species = species.lower())
        objs.append(genome)
        pBar = ProgressBar(label = 'preparing')
        for typ in (Chromosome_Raba, Gene_Raba, Transcript_Raba, Exon_Raba, Protein_Raba) :
            pBar.update()
            f = RabaQuery(typ, namespace = genome._raba_namespace)
            f.addFilter({'genome' : genome})
            for e in f.iterRun() :
                objs.append(e)
        pBar.close()
        
        pBar = ProgressBar(nbEpochs = len(objs), label = 'deleting objects')
        for e in objs :
            pBar.update()
            e.delete()
        pBar.close()
        
    except KeyError as e :
        #~ printf("\tWARNING, couldn't remove genome form db, maybe it's not there: ", e)
        raise KeyError("\tWARNING, couldn't remove genome form db, maybe it's not there: ", e)
        allGood = False
    printf('\tdeleting folder')
    try :
        shutil.rmtree(conf.getGenomeSequencePath(species, name))
    except OSError as e:
        #~ printf('\tWARNING, Unable to delete folder: ', e)
        OSError('\tWARNING, Unable to delete folder: ', e)
        allGood = False
        
    conf.db.endTransaction()
    return allGood
54099
def _importGenomeObjects(gtfFilePath, chroSet, genome, batchSize, verbose = 0) :
    """"""

    class Store(object) :
        
        def __init__(self, conf) :
            self.conf = conf
            self.chromosomes = {}
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

        def batch_save(self) :
            self.conf.db.beginTransaction()
            
            for c in self.genes.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
                
            for c in self.transcripts.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c.exons)
                conf.removeFromDBRegistery(c)
            
            for c in self.proteins.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
            
            self.conf.db.endTransaction()
            
            del(self.genes)
            del(self.transcripts)
            del(self.proteins)
            del(self.exons)
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

            gc.collect()

        def save_chros(self) :
            pBar = ProgressBar(nbEpochs = len(self.chromosomes))
            for c in self.chromosomes.itervalues() :
                pBar.update(label = 'Chr %s' % c.number)
                c.save()
            pBar.close()
        
    printf('Importing gene set infos from %s...' % gtfFilePath)
    
    printf('Backuping indexes...')
    indexes = conf.db.getIndexes()
    printf("Droping all your indexes, (don't worry i'll restore them later)...")
    Genome_Raba.flushIndexes()
    Chromosome_Raba.flushIndexes()
    Gene_Raba.flushIndexes()
    Transcript_Raba.flushIndexes()
    Protein_Raba.flushIndexes()
    Exon_Raba.flushIndexes()
    
    printf("Parsing gene set...")
    gtf = GTFFile(gtfFilePath, gziped = True)
    printf('Done. Importation begins!')
    
    store = Store(conf)
    chroNumber = None
    pBar = ProgressBar(nbEpochs = len(gtf))
    for line in gtf :
        chroN = line['seqname']
        pBar.update(label = "Chr %s" % chroN)
        
        if (chroN.upper() in chroSet or chroN.lower() in chroSet):
            strand = line['strand']
            gene_biotype = line['gene_biotype']
            regionType = line['feature']
            frame = line['frame']

            start = int(line['start']) - 1
            end = int(line['end'])
            if start > end :
                start, end = end, start

            chroNumber = chroN.upper()
            if chroNumber not in store.chromosomes :
                store.chromosomes[chroNumber] = Chromosome_Raba()
                store.chromosomes[chroNumber].set(genome = genome, number = chroNumber)
            
            try :
                geneId = line['gene_id']
                geneName =  line['gene_name']
            except KeyError :
                geneId = None
                geneName = None
                if verbose :
                    printf('Warning: no gene_id/name found in line %s' % gtf[i])
            
            if geneId is not None :
                if geneId not in store.genes :
                    if len(store.genes) > batchSize :
                        store.batch_save()
                    
                    if verbose > 0 :
                        printf('\tGene %s, %s...' % (geneId, geneName))
                    store.genes[geneId] = Gene_Raba()
                    store.genes[geneId].set(genome = genome, id = geneId, chromosome = store.chromosomes[chroNumber], name = geneName, strand = strand, biotype = gene_biotype)
                if start < store.genes[geneId].start or store.genes[geneId].start is None :
                        store.genes[geneId].start = start
                if end > store.genes[geneId].end or store.genes[geneId].end is None :
                    store.genes[geneId].end = end
            try :
                transId = line['transcript_id']
                transName = line['transcript_name']
                try :
                    transcript_biotype = line['transcript_biotype']
                except KeyError :
                    transcript_biotype = None
            except KeyError :
                transId = None
                transName = None
                if verbose > 2 :
                    printf('\t\tWarning: no transcript_id, name found in line %s' % gtf[i])
            
            if transId is not None :
                if transId not in store.transcripts :
                    if verbose > 1 :
                        printf('\t\tTranscript %s, %s...' % (transId, transName))
                    store.transcripts[transId] = Transcript_Raba()
                    store.transcripts[transId].set(genome = genome, id = transId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), name = transName, biotype=transcript_biotype)
                if start < store.transcripts[transId].start or store.transcripts[transId].start is None:
                    store.transcripts[transId].start = start
                if end > store.transcripts[transId].end or store.transcripts[transId].end is None:
                    store.transcripts[transId].end = end
            
                try :
                    protId = line['protein_id']
                except KeyError :
                    protId = None
                    if verbose > 2 :
                        printf('Warning: no protein_id found in line %s' % gtf[i])

                # Store selenocysteine positions in transcript
                if regionType == 'Selenocysteine':
                    store.transcripts[transId].selenocysteine.append(start)
                        
                if protId is not None and protId not in store.proteins :
                    if verbose > 1 :
                        printf('\t\tProtein %s...' % (protId))
                    store.proteins[protId] = Protein_Raba()
                    store.proteins[protId].set(genome = genome, id = protId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), name = transName)
                    store.transcripts[transId].protein = store.proteins[protId]

                try :
                    exonNumber = int(line['exon_number']) - 1
                    exonKey = (transId, exonNumber)
                except KeyError :
                    exonNumber = None
                    exonKey = None
                    if verbose > 2 :
                        printf('Warning: no exon number or id found in line %s' % gtf[i])
                
                if exonKey is not None :
                    if verbose > 3 :
                        printf('\t\t\texon %s...' % (exonId))
                    
                    if exonKey not in store.exons and regionType == 'exon' :
                        store.exons[exonKey] = Exon_Raba()
                        store.exons[exonKey].set(genome = genome, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), protein = store.proteins.get(protId, None), strand = strand, number = exonNumber, start = start, end = end)
                        store.transcripts[transId].exons.append(store.exons[exonKey])
                    
                    try :
                        store.exons[exonKey].id = line['exon_id']
                    except KeyError :
                        pass
                    
                    if regionType == 'exon' :
                        if start < store.exons[exonKey].start or store.exons[exonKey].start is None:
                            store.exons[exonKey].start = start
                        if end > store.transcripts[transId].end or store.exons[exonKey].end is None:
                            store.exons[exonKey].end = end
                    elif regionType == 'CDS' :
                        store.exons[exonKey].CDS_start = start
                        store.exons[exonKey].CDS_end = end
                        store.exons[exonKey].frame = frame
                    elif regionType == 'stop_codon' :
                        if strand == '+' :
                            if store.exons[exonKey].CDS_end != None :
                                store.exons[exonKey].CDS_end += 3
                                if store.exons[exonKey].end < store.exons[exonKey].CDS_end :
                                    store.exons[exonKey].end = store.exons[exonKey].CDS_end
                                if store.transcripts[transId].end < store.exons[exonKey].CDS_end :
                                    store.transcripts[transId].end = store.exons[exonKey].CDS_end
                                if store.genes[geneId].end < store.exons[exonKey].CDS_end :
                                    store.genes[geneId].end = store.exons[exonKey].CDS_end
                        if strand == '-' :
                            if store.exons[exonKey].CDS_start != None :
                                store.exons[exonKey].CDS_start -= 3
                                if store.exons[exonKey].start > store.exons[exonKey].CDS_start :
                                    store.exons[exonKey].start = store.exons[exonKey].CDS_start
                                if store.transcripts[transId].start > store.exons[exonKey].CDS_start :
                                    store.transcripts[transId].start = store.exons[exonKey].CDS_start
                                if store.genes[geneId].start > store.exons[exonKey].CDS_start :
                                    store.genes[geneId].start = store.exons[exonKey].CDS_start
    pBar.close()
    
    store.batch_save()
    
    conf.db.beginTransaction()
    printf('almost done saving chromosomes...')
    store.save_chros()
    
    printf('saving genome object...')
    genome.save()
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring core indexes...')
    # Genome.ensureGlobalIndex(('name', 'species'))
    # Chromosome.ensureGlobalIndex('genome')
    # Gene.ensureGlobalIndex('genome')
    # Transcript.ensureGlobalIndex('genome')
    # Protein.ensureGlobalIndex('genome')
    # Exon.ensureGlobalIndex('genome')
    Transcript.ensureGlobalIndex('exons')
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring user indexes')
    pBar = ProgressBar(label = "restoring", nbEpochs = len(indexes))
    for idx in indexes :
        pBar.update()
        conf.db.execute(idx[-1].replace('CREATE INDEX', 'CREATE INDEX IF NOT EXISTS'))
    pBar.close()
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    return store.chromosomes.values()
1187584
def infer(args):
    
    from jcvi.apps.grid import WriteJobs
    from jcvi.formats.bed import sort

    p = OptionParser(infer.__doc__)
    p.set_cpus()
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    scaffoldsf, genomef = args
    inferbed = "infer-components.bed"
    if need_update((scaffoldsf, genomef), inferbed):
        scaffolds = Fasta(scaffoldsf, lazy=True)
        genome = Fasta(genomef)
        genome = genome.tostring()
        args = [(scaffold_name, scaffold, genome) \
                for scaffold_name, scaffold in scaffolds.iteritems_ordered()]

        pool = WriteJobs(map_one_scaffold, args, inferbed, cpus=opts.cpus)
        pool.run()

    sort([inferbed, "-i"])
    bed = Bed(inferbed)
    inferagpbed = "infer.bed"
    fw = open(inferagpbed, "w")
    seen = []
    for b in bed:
        r = (b.seqid, b.start, b.end)
        if check_seen(r, seen):
            continue
        print("\t".join(str(x) for x in \
                    (b.accn, 0, b.span, b.seqid, b.score, b.strand)), file=fw)
        seen.append(r)
    fw.close()

    frombed([inferagpbed])
1176308
def list_species(model, seeds=None):
    """"""
    seeds = None or []

    from ecell4_base.core import Species

    if not isinstance(seeds, list):
        seeds = list(seeds)

    expanded = model.expand([Species(serial) for serial in seeds])
    species_list = [sp.serial() for sp in expanded.list_species()]
    species_list = sorted(set(seeds + species_list))
    return species_list
724720
def add_reference_resources(data, remote_retriever=None):
    
    aligner = data["config"]["algorithm"].get("aligner", None)
    if remote_retriever:
        data["reference"] = remote_retriever.get_refs(data["genome_build"],
                                                      alignment.get_aligner_with_aliases(aligner, data),
                                                      data["config"])
    else:
        data["reference"] = genome.get_refs(data["genome_build"], alignment.get_aligner_with_aliases(aligner, data),
                                            data["dirs"]["galaxy"], data)
        _check_ref_files(data["reference"], data)
    # back compatible `sam_ref` target
    data["sam_ref"] = utils.get_in(data, ("reference", "fasta", "base"))
    ref_loc = utils.get_in(data, ("config", "resources", "species", "dir"),
                           utils.get_in(data, ("reference", "fasta", "base")))
    if remote_retriever:
        data = remote_retriever.get_resources(data["genome_build"], ref_loc, data)
    else:
        data["genome_resources"] = genome.get_resources(data["genome_build"], ref_loc, data)
    data["genome_resources"] = genome.add_required_resources(data["genome_resources"])
    if effects.get_type(data) == "snpeff" and "snpeff" not in data["reference"]:
        data["reference"]["snpeff"] = effects.get_snpeff_files(data)
    if "genome_context" not in data["reference"]:
        data["reference"]["genome_context"] = annotation.get_context_files(data)
    if "viral" not in data["reference"]:
        data["reference"]["viral"] = viral.get_files(data)
    if not data["reference"]["viral"]:
        data["reference"]["viral"] = None
    if "versions" not in data["reference"]:
        data["reference"]["versions"] = _get_data_versions(data)

    data = _fill_validation_targets(data)
    data = _fill_prioritization_targets(data)
    data = _fill_capture_regions(data)
    # Re-enable when we have ability to re-define gemini configuration directory
    if False:
        data["reference"]["gemini"] = population.get_gemini_files(data)
    return data