get genome lengths
762248
def genome_lengths(fastas, info):
    
    if info is False:
        info = {}
    for genome in fastas:
        name = genome.rsplit('.', 1)[0].rsplit('/', 1)[-1].rsplit('.contigs')[0]
        if name in info:
            continue
        length = 0
        fragments = 0
        for seq in parse_fasta(genome):
            length += len(seq[1])
            fragments += 1
        info[name] = {'genome size (bp)':length, '# contigs':fragments}
    return info
336082
def calculate_content_length(self):
        """"""
        try:
            self._ensure_sequence()
        except RuntimeError:
            return None
        return sum(len(x) for x in self.iter_encoded())
458316
def _get_biodata(base_file, args):
    
    with open(base_file) as in_handle:
        config = yaml.safe_load(in_handle)
    config["install_liftover"] = False
    config["genome_indexes"] = args.aligners
    ann_groups = config.pop("annotation_groups", {})
    config["genomes"] = [_setup_genome_annotations(g, args, ann_groups)
                         for g in config["genomes"] if g["dbkey"] in args.genomes]
    return config
270281
def getGenomeList() :
	""""""
	import rabaDB.filters as rfilt
	f = rfilt.RabaQuery(Genome_Raba)
	names = []
	for g in f.iterRun() :
		names.append(g.name)
	return names
321299
def init_genome_coordinates(self) :
        """"""
        self.seqpos2genome = {}

        # record genome positions for each sequence position
        seq_pos = 0
        for estart, eend in self.exons:
            for genome_pos in range(estart, eend):
                if self.strand == '+':
                    self.seqpos2genome[seq_pos] = genome_pos
                elif self.strand == '-':
                    tmp = self.cds_len - seq_pos - 1
                    self.seqpos2genome[tmp] = genome_pos
                seq_pos += 1

        # recode 5' splice site locations
        for i in range(0, self.five_ss_len):
            seq_pos = self.cds_len + i
            ss_ix = i // 2  # the ss_ix'th 5'ss starting from upstream tx
            pos_in_ss = i % 2  # whether first/second nuc in splice site

            # determine genome coordinates for 5' splice site
            if self.strand == '+':
                self.seqpos2genome[seq_pos] = self.exons[ss_ix][1] + pos_in_ss
            else:
                exon_pos = -1 - ss_ix
                self.seqpos2genome[seq_pos] = self.exons[exon_pos][0] - pos_in_ss - 1

        # recode 3' splice site locations
        for i in range(0, self.three_ss_len):
            seq_pos = self.cds_len + self.five_ss_len + i
            ss_ix = i // 2  # the ss_ix'th 3'ss starting from upstream tx
            pos_in_ss = i % 2  # whether first/second nuc in splice site

            # determine genome coordinates for 3' splice site
            if self.strand == '+':
                self.seqpos2genome[seq_pos] = self.exons[ss_ix+1][0] - 2 + pos_in_ss
            else:
                exon_pos = -1 - ss_ix
                self.seqpos2genome[seq_pos] = self.exons[exon_pos-1][1] + 1 - pos_in_ss
1174249
def six_frame(genome, table, minimum = 10):
    
    for seq in parse_fasta(genome):
        dna = Seq(seq[1].upper().replace('U', 'T'), IUPAC.ambiguous_dna)
        counter = 0
        for sequence in ['f', dna], ['rc', dna.reverse_complement()]:
            direction, sequence = sequence
            for frame in range(0, 3):
                for prot in \
                            sequence[frame:].\
                            translate(table = table, to_stop = False).split('*'):
                    if len(prot) < minimum:
                        continue
                    counter += 1
                    header = '%s_%s table=%s frame=%s-%s %s' % \
                                (seq[0].split()[0], counter, table, frame+1, \
                                direction, ' '.join(seq[0].split()[1:]))
                    yield [header, prot]
291260
def _guess_seqtype_from_file(handle):
    ""
    if isinstance(handle, basestring):
        handle = StringIO(handle)

    for line in handle:
        if not line.strip():
            continue
        if line.lstrip().split()[0] in ('LOCUS', 'FEATURES', 'source', 'CDS',
                                        'gene'):
            return 'genbank'
        if len(line) > 2 and line[:3] in ('ID ', 'FT '):
            return 'embl'
        if line.startswith('>'):
            return 'fasta'
    handle.seek(0)
    import string
    from Bio.Data import IUPACData as iupac
    all_input_letters = set(handle.read().lower())
    all_valid = set(string.digits)
    all_valid.update(set(iupac.protein_letters.lower()))
    all_valid.update(set(iupac.unambiguous_dna_letters.lower()))
    all_valid.update(set('- \n'))
    if all_valid.issuperset(all_input_letters):
        return 'fasta'

    raise ValueError("Failed to guess format for input")
26359
def calculate_new_length(gene_split, gene_results, hit):
   ''' 
   '''
   # Looping over splitted hits and calculate new length
   first = 1
   for split in gene_split[hit['sbjct_header']]:
      new_start = int(gene_results[split]['sbjct_start'])
      new_end = int(gene_results[split]['sbjct_end'])
      
      # Get the frist HSP
      if first == 1:
         new_length = int(gene_results[split]['HSP_length'])
         old_start = new_start
         old_end = new_end
         first = 0
         continue
      if new_start < old_start:
         new_length = new_length + (old_start - new_start)
         old_start = new_start
      
      if new_end > old_end:
         new_length = new_length + (new_end - old_end)
         old_end = new_end
   
   return(new_length)
136930
def parse_sampleinfo(data: dict) -> dict:
    

    genome_build = data['human_genome_build']
    genome_build_str = f"{genome_build['source']}{genome_build['version']}"
    if 'svdb' in data['program']:
        svdb_outpath = (f"{data['program']['svdb']['path']}")
    else:
        svdb_outpath = ''
    outdata = {
        'date': data['analysis_date'],
        'family': data['family'],
        'genome_build': genome_build_str,
        'rank_model_version': data['program']['genmod']['rank_model']['version'],
        'is_finished': True if data['analysisrunstatus'] == 'finished' else False,
        'pedigree_path': data['pedigree_minimal'],
        'peddy': {
            'ped': (data['program']['peddy']['peddy']['path'] if
                    'peddy' in data['program'] else None),
            'ped_check': (data['program']['peddy']['ped_check']['path'] if
                          'peddy' in data['program'] else None),
            'sex_check': (data['program']['peddy']['sex_check']['path'] if
                          'peddy' in data['program'] else None),
        },
        'qcmetrics_path': data['program']['qccollect']['path'],
        'samples': [],
        'snv': {
            'bcf': data['most_complete_bcf']['path'],
            'clinical_vcf': data['vcf_binary_file']['clinical']['path'],
            'gbcf': data['gbcf_file']['path'],
            'research_vcf': data['vcf_binary_file']['research']['path'],
        },
        'svdb_outpath': svdb_outpath,
        'sv': {
            'bcf': data.get('sv_bcf_file', {}).get('path'),
            'clinical_vcf': (data['sv_vcf_binary_file']['clinical']['path'] if
                             'sv_vcf_binary_file' in data else None),
            'merged': svdb_outpath,
            'research_vcf': (data['sv_vcf_binary_file']['research']['path'] if
                             'sv_vcf_binary_file' in data else None),
        },
        'version': data['mip_version'],
    }

    for sample_id, sample_data in data['sample'].items():
        sample = {
            'id': sample_id,
            'bam': sample_data['most_complete_bam']['path'],
            'sambamba': list(sample_data['program']['sambamba_depth'].values())[0]['path'],
            'sex': sample_data['sex'],
            # subsample mt is only for wgs data
            'subsample_mt': (list(sample_data['program']['samtools_subsample_mt'].values())[0]['path'] if
                             'samtools_subsample_mt' in sample_data['program'] else None),
            'vcf2cytosure': list(sample_data['program']['vcf2cytosure'].values())[0]['path'],
        }
        chanjo_sexcheck = list(sample_data['program']['chanjo_sexcheck'].values())[0]
        sample['chanjo_sexcheck'] = chanjo_sexcheck['path']
        outdata['samples'].append(sample)

    return outdata
964434
def _get_all_lengths(self):
        
        utt_lengths = {}

        for utt_idx in self.utt_ids:
            per_container = [c._file[utt_idx].shape[0] for c in self.containers]
            utt_lengths[utt_idx] = tuple(per_container)

        return utt_lengths
806742
def estimate_genome_size(self):
        
        logging.info('Estimating genome size using kmercountexact')
        for sample in self.metadata:
            # Initialise the name of the output file
            sample[self.analysistype].peaksfile = os.path.join(sample[self.analysistype].outputdir, 'peaks.txt')
            # Run the kmer counting command
            out, err, cmd = bbtools.kmercountexact(forward_in=sorted(sample.general.fastqfiles)[0],
                                                   peaks=sample[self.analysistype].peaksfile,
                                                   returncmd=True,
                                                   threads=self.cpus)
            # Set the command in the object
            sample[self.analysistype].kmercountexactcmd = cmd
            # Extract the genome size from the peaks file
            sample[self.analysistype].genomesize = bbtools.genome_size(sample[self.analysistype].peaksfile)
            write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
448978
def n_to_g(self, n_interval):
        """"""

        frs, fre = n_interval.start.base - 1, n_interval.end.base - 1
        start_offset, end_offset = n_interval.start.offset, n_interval.end.offset

        if self.strand == -1:
            fre, frs = self.tgt_len - frs - 1, self.tgt_len - fre - 1
            start_offset, end_offset = -end_offset, -start_offset

        # returns the genomic range start (grs) and end (gre)
        grs, _, grs_cigar = self._map(
            from_pos=self.tgt_pos, to_pos=self.ref_pos, pos=frs, base="start")
        gre, _, gre_cigar = self._map(
            from_pos=self.tgt_pos, to_pos=self.ref_pos, pos=fre, base="end")
        grs, gre = grs + self.gc_offset + 1, gre + self.gc_offset + 1
        gs, ge = grs + start_offset, gre + end_offset

        # The returned interval would be uncertain when locating at alignment gaps
        return hgvs.location.Interval(
            start=hgvs.location.SimplePosition(gs, uncertain=n_interval.start.uncertain),
            end=hgvs.location.SimplePosition(ge, uncertain=n_interval.end.uncertain),
            uncertain=grs_cigar in 'DI' or gre_cigar in 'DI')
293669
def find_genome_length(self):
        
        for sample in self.metadata:
            # Use the sum() method to add all the contig lengths in the list
            sample[self.analysistype].genome_length = sum(sample[self.analysistype].contig_lengths)
1251410
def generate_sizes(name, genome_dir):
    """"""
    fa = os.path.join(genome_dir, name, "{}.fa".format(name))
    sizes = fa + ".sizes"
    g = Fasta(fa)
    with open(sizes, "w") as f:
        for seqname in g.keys():
            f.write("{}\t{}\n".format(seqname, len(g[seqname])))
364255
def getUI(prog_name, args):
  """"""
  longDescription = "Given a set of BED intervals, compute a profile of " +\
                    "conservation by averaging over all intervals using a " +\
                    "whole genome alignment to a set of relevent species." +\
                    "\n\n" +\
                    "Usage: " + prog_name + " [options] regions.bed " +\
                    "genome-alig species" +\
                    "\n\n" +\
                    "genome-alig can be either a single MAF file, or a " +\
                    "directory of MAF files. In the latter case, the " +\
                    "directory may also optionally contain index files for " +\
                    "the alignment files."
  shortDescription = longDescription

  ui = CLI(prog_name, shortDescription, longDescription)
  # gotta have two args -- MAF dir/file and BED regions.
  # Input by stdin not allowed
  ui.minArgs = 3
  ui.maxArgs = 4
  ui.addOption(Option(short="o", long="output", argName="filename",
                      description="output to given file, else stdout",
                      required=False, type=str))
  ui.addOption(Option(short="w", long="window", argName="size",
                      description="size of window to compute around each " +
                                  "interval; " +
                                  str(DEFAULT_WINDOW_SIZE) +
                                  " to use whole interval. " +
                                  "Default " + str(DEFAULT_WINDOW_SIZE),
                      required=False, type=int))
  ui.addOption(Option(short="e", long="extensions", argName="extension",
                      description="if genome-alig specifies a directory, " +
                                  "treat files with this extension as " +
                                  "alignment files.", required=False,
                                  type=str))
  ui.addOption(Option(short="i", long="index-extensions", argName="extension",
                      description="if genome-alig specifies a directory, " +
                                  "treat files with this extension as " +
                                  "index files for alignments.",
                                  required=False, type=str))
  ui.addOption(Option(short="f", long="fail-no-index",
                      description="fail if an alignment file without an " +
                                  "index is found; otherwise index-less " +
                                  "alignment files are loaded whole (which " +
                                  "might be slow if they're large, and " +
                                  "might require a lot of memory)",
                      default=False, required=False))
  ui.addOption(Option(short="m", long="missing", argName="strategy",
                      description="how to treat missing sequences in " +
                                  "blocks. Options are " +
                                  ", ".join([str(x.name) for x in
                                             MissingSequenceHandler]),
                                  required=False, type=str))
  ui.addOption(Option(short="s", long="species", argName="species",
                      description="consider only these species. Default is " +
                                  "all.", required=False, type=str))
  ui.addOption(Option(short="v", long="verbose",
                      description="output additional messages to stderr " +
                                  "about run (default: " +
                                  str(DEFAULT_VERBOSITY) + ")",
                      default=DEFAULT_VERBOSITY, required=False))
  ui.addOption(Option(short="h", long="help",
                      description="show this help message ", special=True))
  ui.addOption(Option(short="u", long="test",
                      description="run unit tests ", special=True))

  ui.parseCommandLine(args)
  return ui
266464
def length_by_id(self):
        """"""
        hashmap = dict((seq.id, len(seq)) for seq in self)
        tmp = hashmap.copy()
        hashmap.update(tmp)
        return hashmap
1177028
def get_encoded_seq_len(self, seq_len: int) -> int:
        
        for encoder in self.encoders:
            seq_len = encoder.get_encoded_seq_len(seq_len)
        return seq_len