parses and plots GC content and genome GC content files
559135
def parse_gc_content(self):
        """"""
        # Find and parse GC content:
        for f in self.find_log_files('homer/GCcontent', filehandles=True):
            # Get the s_name from the parent directory
            s_name = os.path.basename(f['root'])
            s_name = self.clean_s_name(s_name, f['root'])
            parsed_data = self.parse_twoCol_file(f)
            if parsed_data is not None:
                if s_name in self.tagdir_data['GCcontent']:
                    log.debug("Duplicate GCcontent sample log found! Overwriting: {}".format(s_name))

                self.add_data_source(f, s_name, section='GCcontent')
                self.tagdir_data['GCcontent'][s_name] = parsed_data

        ## get esimated genome content distribution:
        for f in self.find_log_files('homer/genomeGCcontent', filehandles=True):
            parsed_data = self.parse_twoCol_file(f)
            if parsed_data is not None:
                if s_name + "_genome" in self.tagdir_data['GCcontent']:
                    log.debug("Duplicate genome GCcontent sample log found! Overwriting: {}".format(s_name+ "_genome"))

                self.add_data_source(f, s_name + "_genome", section='GCcontent')
                self.tagdir_data['GCcontent'][s_name + "_genome"] = parsed_data

        self.tagdir_data['GCcontent'] = self.ignore_samples(self.tagdir_data['GCcontent'])

        if len(self.tagdir_data['GCcontent']) > 0:
            self.add_section (
                name = 'Per Sequence GC Content',
                anchor = 'homer_per_sequence_gc_content',
                description = 'This plot shows the distribution of GC content.',
                plot = self.GCcontent_plot()
            )
56338
def parse_rRNA(insertion, seq, gff):
    
    offset = insertion['offset']
    strand = insertion['strand']
    for rRNA in parse_masked(seq, 0)[0]:
        rRNA = ''.join(rRNA)
        Start = seq[1].find(rRNA) + 1
        End = Start + len(rRNA) - 1
        if strand == '-':
            Start, End = End - 2, Start - 2
        pos = (abs(Start + offset) - 1, abs(End + offset) - 1)
        Start, End = min(pos), max(pos)
        source = insertion['source']
        annot = '%s rRNA' % (source.split('from', 1)[0])
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(source)
        gff['feature'].append('rRNA')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('Name=%s' % (annot))
    return gff
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()
1291224
def fixpartials(args):
    
    p = OptionParser(fixpartials.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gffile, gfasta, partials, = args

    gff = make_index(gffile)
    genome = Fasta(gfasta, index=True)
    partials = LineFile(partials, load=True).lines

    # all_transcripts = [f.id for f in gff.features_of_type("mRNA", \
    #     order_by=("seqid", "start"))]
    seen = set()
    fw = must_open(opts.outfile, "w")
    for gene in gff.features_of_type('gene',  order_by=('seqid', 'start')):
        children = AutoVivification()
        cflag = False
        transcripts = list(gff.children(gene, level=1, order_by=('start')))
        for transcript in transcripts:
            trid, seqid, strand = transcript.id, transcript.seqid, transcript.strand

            for child in gff.children(transcript, order_by=('start')):
                ftype = child.featuretype
                if ftype not in children[trid]: children[trid][ftype] = []
                children[trid][ftype].append(child)

            five_prime, three_prime = True, True
            nstart, nstop = (None, None), (None, None)
            cds_span = [children[trid]['CDS'][0].start, \
                children[trid]['CDS'][-1].stop]
            new_cds_span = [x for x in cds_span]

            start_codon = (cds_span[0], cds_span[0] + 2)
            stop_codon = (cds_span[1] - 2, cds_span[1])
            if strand == '-': start_codon, stop_codon = stop_codon, start_codon

            if trid in partials:
                seen.add(trid)
                start_codon_fasta = _fasta_slice(genome, seqid, \
                    start_codon[0], start_codon[1], strand)
                stop_codon_fasta = _fasta_slice(genome, seqid, \
                    stop_codon[0], stop_codon[1], strand)

                if not is_valid_codon(start_codon_fasta, type='start'):
                    five_prime = False
                    nstart = scan_for_valid_codon(start_codon, strand, \
                        seqid, genome, type='start')

                if not is_valid_codon(stop_codon_fasta, type='stop'):
                    three_prime = False
                    nstop = scan_for_valid_codon(stop_codon, strand, \
                        seqid, genome, type='stop')

                logging.debug("feature={0} ({1})".format(trid, strand) + \
                ", 5'={0}, 3'={1}".format(five_prime, three_prime) + \
                ", {0} <== {1} ==> {2}".format(nstart if strand == '+' else nstop, \
                    cds_span, nstop if strand == '+' else nstart))

            if not five_prime or not three_prime:
                if nstart != (None, None) and (start_codon != nstart):
                    i = 0 if strand == '+' else 1
                    new_cds_span[i] = nstart[i]
                if nstop != (None, None) and (stop_codon != nstop):
                    i = 1 if strand == '+' else 0
                    new_cds_span[i] = nstop[i]
                new_cds_span.sort()

            if set(cds_span) != set(new_cds_span):
                cflag = True
                # if CDS has been extended, appropriately adjust all relevent
                # child feature (CDS, exon, UTR) coordinates
                for ftype in children[trid]:
                    for idx in range(len(children[trid][ftype])):
                        child_span = (children[trid][ftype][idx].start, \
                            children[trid][ftype][idx].stop)
                        if ftype in ('exon', 'CDS'):
                            # if exons/CDSs, adjust start and stop according to
                            # new CDS start and stop, respectively
                            if child_span[0] == cds_span[0]:
                                children[trid][ftype][idx].start = new_cds_span[0]
                            if child_span[1] == cds_span[1]:
                                children[trid][ftype][idx].stop = new_cds_span[1]
                        elif ftype.endswith('UTR'):
                            # if *_prime_UTR, adjust stop according to new CDS start and
                            #                 adjust start according to new CDS stop
                            if child_span[1] == cds_span[0]:
                                children[trid][ftype][idx].stop = new_cds_span[0]
                            if child_span[0] == cds_span[1]:
                                children[trid][ftype][idx].start = new_cds_span[1]

                transcript.start, transcript.stop = \
                    children[trid]['exon'][0].start, children[trid]['exon'][-1].stop

        if cflag:
            _gene_span = range_minmax([(tr.start, tr.stop) for tr in transcripts])
            gene.start, gene.stop = _gene_span[0], _gene_span[1]

        # print gff file
        print(gene, file=fw)
        for transcript in transcripts:
            trid = transcript.id
            print(transcript, file=fw)
            for cftype in children[trid]:
                for child in children[trid][cftype]:
                    print(child, file=fw)

    fw.close()
483341
def get_genome_ref(genome_build, aligner, galaxy_base):
    
    ref_files = dict(
            bowtie = "bowtie_indices.loc",
            bwa = "bwa_index.loc",
            samtools = "sam_fa_indices.loc",
            maq = "bowtie_indices.loc")
    remap_fns = dict(
            maq = _remap_to_maq
            )
    out_info = []
    for ref_get in [aligner, "samtools"]:
        ref_file = os.path.join(galaxy_base, "tool-data", ref_files[ref_get])
        with open(ref_file) as in_handle:
            for line in in_handle:
                if not line.startswith("#"):
                    parts = line.strip().split()
                    if parts[0] == "index":
                        parts = parts[1:]
                    if parts[0] == genome_build:
                        out_info.append(parts[-1])
                        break
        try:
            out_info[-1] = remap_fns[ref_get](out_info[-1])
        except KeyError:
            pass
        except IndexError:
            raise IndexError("Genome %s not found in %s" % (genome_build,
                ref_file))

    if len(out_info) != 2:
        raise ValueError("Did not find genome reference for %s %s" %
                (genome_build, aligner))
    else:
        return tuple(out_info)
945038
def handle_synonyms(in_file, ref_file, genome_build, work_dir, data):
    
    if genome_build in GMAP and ref_file:
        mappings = GMAP[genome_build]
        contigs = set([c.name for c in ref.file_contigs(ref_file)])
        out_file = os.path.join(work_dir, "%s-fixed_contigs%s" % utils.splitext_plus(os.path.basename(in_file)))
        if not utils.file_exists(out_file):
            if out_file.endswith(".gz"):
                out_file = out_file.replace(".gz", "")
                needs_bgzip = True
            else:
                needs_bgzip = False
            checked_file = "%s.checked" % utils.splitext_plus(out_file)[0]
            if not _matches_contigs(in_file, contigs, checked_file):
                with file_transaction(data, out_file) as tx_out_file:
                    _write_newname_file(in_file, tx_out_file, mappings)
                if needs_bgzip:
                    out_file = vcfutils.bgzip_and_index(out_file, data["config"])
                return out_file
    return in_file
1313511
def parse_kegg_gene_metadata(infile):
    
    metadata = defaultdict(str)

    with open(infile) as mf:
        kegg_parsed = bs_kegg.parse(mf.read())

    # TODO: additional fields can be parsed

    if 'DBLINKS' in kegg_parsed.keys():
        if 'UniProt' in kegg_parsed['DBLINKS']:
            unis = str(kegg_parsed['DBLINKS']['UniProt']).split(' ')
            # TODO: losing other uniprot ids by doing this
            if isinstance(unis, list):
                metadata['uniprot'] = unis[0]
            else:
                metadata['uniprot'] = unis
        if 'NCBI-ProteinID' in kegg_parsed['DBLINKS']:
            metadata['refseq'] = str(kegg_parsed['DBLINKS']['NCBI-ProteinID'])
    if 'STRUCTURE' in kegg_parsed.keys():
        metadata['pdbs'] = str(kegg_parsed['STRUCTURE']['PDB']).split(' ')
    else:
        metadata['pdbs'] = None
    if 'ORGANISM' in kegg_parsed.keys():
        metadata['taxonomy'] = str(kegg_parsed['ORGANISM'])

    return metadata
466547
def parse_reports(self):
    

    # Set up vars
    self.read_gc = dict()
    self.read_gc_pct = dict()

    # Go through files and parse data
    for f in self.find_log_files('rseqc/read_gc'):

        if f['f'].startswith('GC%	read_count'):
            gc = list()
            counts = list()
            for l in f['f'].splitlines():
                s = l.split()
                try:
                    gc.append(float(s[0]))
                    counts.append(float(s[1]))
                except:
                    pass
            if len(gc) > 0:
                sorted_gc_keys = sorted(range(len(gc)), key=lambda k: gc[k])
                total = sum(counts)
                if f['s_name'] in self.read_gc:
                    log.debug("Duplicate sample name found! Overwriting: {}".format(f['s_name']))
                self.add_data_source(f, section='read_GC')
                self.read_gc[f['s_name']] = OrderedDict()
                self.read_gc_pct[f['s_name']] = OrderedDict()
                for i in sorted_gc_keys:
                    self.read_gc[f['s_name']][gc[i]] = counts[i]
                    self.read_gc_pct[f['s_name']][gc[i]] = (counts[i]/total)*100

    # Filter to strip out ignored sample names
    self.read_gc = self.ignore_samples(self.read_gc)

    if len(self.read_gc) > 0:

        # Add line graph to section
        pconfig = {
            'id': 'rseqc_read_gc_plot',
            'title': 'RSeQC: Read GC Content',
            'ylab': 'Number of Reads',
            'xlab': "GC content (%)",
            'xmin': 0,
            'xmax': 100,
            'tt_label': "<strong>{point.x}% GC</strong>: {point.y:.2f}",
            'data_labels': [
                {'name': 'Counts', 'ylab': 'Number of Reads'},
                {'name': 'Percentages', 'ylab': 'Percentage of Reads'}
            ]
        }
        self.add_section (
            name = 'Read GC Content',
            anchor = 'rseqc-read_gc',
            description = '<a href="http://rseqc.sourceforge.net/#read-gc-py" target="_blank">read_GC</a>' \
                " calculates a histogram of read GC content.</p>",
            plot = linegraph.plot([self.read_gc, self.read_gc_pct], pconfig)
        )

    # Return number of samples found
    return len(self.read_gc)
1044224
def parse_cov(cov_table, scaffold2genome):
    
    size   = {} # size[genome] = genome size
    mapped = {} # mapped[genome][sample] = mapped bases
    # parse coverage files
    for line in open(cov_table):
        line = line.strip().split('\t')
        if line[0].startswith('#'):
            samples = line[1:]
            samples = [i.rsplit('/', 1)[-1].split('.', 1)[0] for i in samples]
            continue
        scaffold, length = line[0].split(': ')
        length = float(length)
        covs  = [float(i) for i in line[1:]]
        bases = [c * length for c in covs]
        if scaffold not in scaffold2genome:
            continue
        genome = scaffold2genome[scaffold]
        if genome not in size:
            size[genome] = 0
            mapped[genome] = {sample:0 for sample in samples}
        # keep track of genome size
        size[genome] += length
        # keep track of number of mapped bases
        for sample, count in zip(samples, bases):
            mapped[genome][sample] += count
    # calculate coverage from base counts and genome size
    coverage = {'genome':[], 'genome size (bp)':[], 'sample':[], 'coverage':[]}
    for genome, length in size.items():
        for sample in samples:
            cov = mapped[genome][sample] / length
            coverage['genome'].append(genome)
            coverage['genome size (bp)'].append(length)
            coverage['sample'].append(sample)
            coverage['coverage'].append(cov)
    return pd.DataFrame(coverage)
483873
def make_pizzly_gtf(gtf_file, out_file, data):
    
    if file_exists(out_file):
        return out_file
    db = gtf.get_gtf_db(gtf_file)
    with file_transaction(data, out_file) as tx_out_file:
        with open(tx_out_file, "w") as out_handle:
            for gene in db.features_of_type("gene"):
                children = [x for x in db.children(id=gene)]
                for child in children:
                    if child.attributes.get("gene_biotype", None):
                        gene_biotype = child.attributes.get("gene_biotype")
                gene.attributes['gene_biotype'] = gene_biotype
                gene.source = gene_biotype[0]
                print(gene, file=out_handle)
                for child in children:
                    child.source = gene_biotype[0]
                    # gffread produces a version-less FASTA file
                    child.attributes.pop("transcript_version", None)
                    print(child, file=out_handle)
    return out_file
1348160
def wgsim(args):
    
    p = OptionParser(wgsim.__doc__)
    p.add_option("--erate", default=.01, type="float",
                 help="Base error rate of the read [default: %default]")
    p.add_option("--noerrors", default=False, action="store_true",
                 help="Simulate reads with no errors [default: %default]")
    p.add_option("--genomesize", type="int",
                 help="Genome size in Mb [default: estimate from data]")
    add_sim_options(p)
    opts, args = p.parse_args(args)

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

    fastafile, = args
    pf = op.basename(fastafile).split(".")[0]

    genomesize = opts.genomesize
    size = genomesize * 1000000 if genomesize else Fasta(fastafile).totalsize
    depth = opts.depth
    readlen = opts.readlen
    readnum = int(math.ceil(size * depth / (2 * readlen)))

    distance = opts.distance
    stdev = distance / 10

    outpf = opts.outfile or "{0}.{1}bp.{2}x".format(pf, distance, depth)

    logging.debug("Total genome size: {0} bp".format(size))
    logging.debug("Target depth: {0}x".format(depth))
    logging.debug("Number of read pairs (2x{0}): {1}".format(readlen, readnum))

    if opts.noerrors:
        opts.erate = 0

    cmd = "dwgsim -e {0} -E {0}".format(opts.erate)
    if opts.noerrors:
        cmd += " -r 0 -R 0 -X 0 -y 0"

    cmd += " -d {0} -s {1}".format(distance, stdev)
    cmd += " -N {0} -1 {1} -2 {1}".format(readnum, readlen)
    cmd += " {0} {1}".format(fastafile, outpf)
    sh(cmd)
146830
def _get_G2P(data):
    
    G2P_file = os.path.realpath(tz.get_in(("genome_resources", "variation", "genotype2phenotype"), data))
    args = ["--plugin", "G2P,file:%s" % (G2P_file)]
    if G2P_file:
        return args
    else:
        return []
51291
def shuffle_genome(genome, cat, fraction = float(100), plot = True, \
        alpha = 0.1, beta = 100000, \
        min_length = 1000, max_length = 200000):
    
    header = '>randomized_%s' % (genome.name)
    sequence = list(''.join([i[1] for i in parse_fasta(genome)]))
    length = len(sequence)
    shuffled = []
    # break genome into pieces
    while sequence is not False:
        s = int(random.gammavariate(alpha, beta))
        if s <= min_length or s >= max_length:
            continue
        if len(sequence) < s:
            seq = sequence[0:]
        else:
            seq = sequence[0:s]
        sequence = sequence[s:]
#        if bool(random.getrandbits(1)) is True:
#            seq = rev_c(seq)
#            print('fragment length: %s reverse complement: True' % ('{:,}'.format(s)), file=sys.stderr)
#        else:
#            print('fragment length: %s reverse complement: False' % ('{:,}'.format(s)), file=sys.stderr)
        shuffled.append(''.join(seq))
        if sequence == []:
            break
    # shuffle pieces
    random.shuffle(shuffled)
    # subset fragments
    if fraction == float(100):
        subset = shuffled
    else:
        max_pieces = int(length * fraction/100)
        subset, total = [], 0
        for fragment in shuffled:
            length = len(fragment)
            if total + length <= max_pieces:
                subset.append(fragment)
                total += length
            else:
                diff = max_pieces - total
                subset.append(fragment[0:diff])
                break
    # combine sequences, if requested
    if cat is True:
        yield [header, ''.join(subset)]
    else:
        for i, seq in enumerate(subset):
            yield ['%s fragment:%s' % (header, i), seq]
682287
def parse_gc_dist(self, f):
    
    # Get the sample name from the parent parent directory
    # Typical path: <sample name>/raw_data_qualimapReport/mapped_reads_gc-content_distribution.txt
    s_name = self.get_s_name(f)

    d = dict()
    reference_species = None
    reference_d = dict()
    avg_gc = 0
    for l in f['f']:
        if l.startswith('#'):
            sections = l.strip("\n").split("\t", 3)
            if len(sections) > 2:
                reference_species = sections[2]
            continue
        sections = l.strip("\n").split("\t", 3)
        gc = int(round(float(sections[0])))
        content = float(sections[1])
        avg_gc += gc * content
        d[gc] = content
        if len(sections) > 2:
            reference_content = float(sections[2])
            reference_d[gc] = reference_content

    # Add average GC to the general stats table
    self.general_stats_data[s_name]['avg_gc'] = avg_gc

    # Save results
    if s_name in self.qualimap_bamqc_gc_content_dist:
        log.debug("Duplicate Mapped Reads GC content distribution sample name found! Overwriting: {}".format(s_name))
    self.qualimap_bamqc_gc_content_dist[s_name] = d
    if reference_species and reference_species not in self.qualimap_bamqc_gc_by_species:
        self.qualimap_bamqc_gc_by_species[reference_species] = reference_d
    self.add_data_source(f, s_name=s_name, section='mapped_gc_distribution')
1048826
def parse_re_dist(self):
        """"""
        # Find and parse homer restriction distribution reports
        for f in self.find_log_files('homer/RestrictionDistribution', filehandles=True):
            s_name = os.path.basename(f['root'])
            s_name = self.clean_s_name(s_name, f['root'])
            parsed_data = self.parse_restriction_dist(f)
            if parsed_data is not None:
                if s_name in self.tagdir_data['restriction']:
                    log.debug("Duplicate Restriction Distribution sample log found! Overwriting: {}".format(s_name))

                self.add_data_source(f, s_name, section='restriction')
                self.tagdir_data['restriction'][s_name] = parsed_data
                self.tagdir_data['restriction_norm'][s_name] = self.normalize(parsed_data)

        self.tagdir_data['restriction'] = self.ignore_samples(self.tagdir_data['restriction'])
        self.tagdir_data['restriction_norm'] = self.ignore_samples(self.tagdir_data['restriction_norm'])

        if len(self.tagdir_data['restriction']) > 0:
            self.add_section (
                name = 'Restriction Site Tag Dist',
                anchor = 'homer-restrictionDist',
                description = 'This plot shows the distribution of tags around restriction enzyme cut sites.',
                helptext = '''Hi-C data often involves the digestion of DNA using a restriction enzyme.
                    A good quality control for the experiment is the centering of reads around the
                    restriction enzyme cut site.''',
                plot = self.restriction_dist_chart()
            )
517725
def gcs_files(prefix_filter=None):
  """"""
  top_level_xml_str = download_gcs_file("", prefix_filter=prefix_filter)
  xml_root = ElementTree.fromstring(top_level_xml_str)
  filenames = [el[0].text for el in xml_root if el.tag.endswith("Contents")]
  return filenames
500017
def _parse_glob_spec(self, spec):
        '''
        
        '''
        some_globs = glob.glob(spec['glob'])
        if not some_globs:
            return []
        el_globs = []
        for g in some_globs:
            _spec = copy.copy(spec)
            _spec['file'] = g
            el_globs.append(_spec)
        return el_globs