Calculate the N50 for each strain. N50 is defined as the largest contig such that at least half of the total
        genome size is contained in contigs equal to or larger than this contig
179872
def find_n50(self):
        
        for sample in self.metadata:
            # Initialise the N50 attribute in case there is no assembly, and the attribute is not created in the loop
            sample[self.analysistype].n50 = '-'
            # Initialise a variable to store a running total of contig lengths
            currentlength = 0
            for contig_length in sample[self.analysistype].contig_lengths:
                # Increment the current length with the length of the current contig
                currentlength += contig_length
                # If the current length is now greater than the total genome / 2, the current contig length is the N50
                if currentlength >= sample[self.analysistype].genome_length * 0.5:
                    # Populate the dictionary, and break the loop
                    sample[self.analysistype].n50 = contig_length
                    break
1261742
def quast_general_stats_table(self):
        

        headers = OrderedDict()
        headers['N50'] = {
            'title': 'N50 ({})'.format(self.contig_length_suffix),
            'description': 'N50 is the contig length such that using longer or equal length contigs produces half (50%) of the bases of the assembly (kilo base pairs)',
            'min': 0,
            'suffix': self.contig_length_suffix,
            'scale': 'RdYlGn',
            'modify': lambda x: x * self.contig_length_multiplier
        }
        headers['Total length'] = {
            'title': 'Length ({})'.format(self.total_length_suffix),
            'description': 'The total number of bases in the assembly (mega base pairs).',
            'min': 0,
            'suffix': self.total_length_suffix,
            'scale': 'YlGn',
            'modify': lambda x: x * self.total_length_multiplier
        }
        self.general_stats_addcols(self.quast_data, headers)
936789
def trim_N_nucleotides(prefix, suffix):
    
    if 'N' in prefix:
        # trim prefix to exclude all occurrences of N
        rightmost_index = prefix.rfind('N')
        logger.debug(
            "Trimming %d nucleotides from read prefix '%s'",
            rightmost_index + 1, prefix)
        prefix = prefix[rightmost_index + 1:]

    if 'N' in suffix:
        leftmost_index = suffix.find('N')
        logger.debug(
            "Trimming %d nucleotides from read suffix '%s'",
            len(suffix) - leftmost_index,
            suffix)
        suffix = suffix[:leftmost_index]

    return prefix, suffix
199425
public int numLines() {
    int r = header.size();
    for(MeasurementGroup m : groups) {
      r += 1 + m.measurements.size();
    }
    return r;
  }
466433
def load(self, read_tuple_name):
        
        self.prefix_width = 0
        self.read_tuple_id_width = 0
        self.genome_id_width = 0
        self.chr_id_width = 0
        self.coor_width = 0

        parts = read_tuple_name.split("__")
        self.prefix_width = len(parts[0])
        self.read_tuple_id_width = len(parts[1])

        segments = parts[2][1:-1].split("),(")
        for segment in segments:
            int_widths = list(map(len, segment.split(",")))
            self.genome_id_width = max(self.genome_id_width, int_widths[0])
            self.chr_id_width = max(self.chr_id_width, int_widths[1])
            self.coor_width = max(self.coor_width, int_widths[2], int_widths[3])
1055227
def calcTm(seq, mv_conc=50, dv_conc=0, dntp_conc=0.8, dna_conc=50,
           max_nn_length=60, tm_method='santalucia',
           salt_corrections_method='santalucia'):
    ''' 
    '''
    tm_meth = _tm_methods.get(tm_method)
    if tm_meth is None:
        raise ValueError('{} is not a valid tm calculation method'.format(
                         tm_method))
    salt_meth = _salt_corrections_methods.get(salt_corrections_method)
    if salt_meth is None:
        raise ValueError('{} is not a valid salt correction method'.format(
                         salt_corrections_method))
    # For whatever reason mv_conc and dna_conc have to be ints
    args = [pjoin(PRIMER3_HOME, 'oligotm'),
            '-mv',  str(mv_conc),
            '-dv',  str(dv_conc),
            '-n',   str(dntp_conc),
            '-d',   str(dna_conc),
            '-tp',  str(tm_meth),
            '-sc',  str(salt_meth),
            seq]
    tm = subprocess.check_output(args, stderr=DEV_NULL,
                                 env=os.environ)
    return float(tm)
316981
def main(self):
        
        self.fasta_records()
        self.fasta_stats()
        self.find_largest_contig()
        self.find_genome_length()
        self.find_num_contigs()
        self.find_n50()
        self.perform_pilon()
        self.clear_attributes()
1298531
def _process_prb_strain_genotype_view(self, limit=None):
        

        line_counter = 0
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        LOG.info("Getting genotypes for strains")
        raw = '/'.join((self.rawdir, 'prb_strain_genotype_view'))
        with open(raw, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for line in filereader:
                line_counter += 1
                if line_counter == 1:
                    continue

                (strain_key, genotype_key) = line

                if self.test_mode is True:
                    if int(genotype_key) not in self.test_keys.get('genotype') \
                            and int(strain_key) not in self.test_keys.get('strain'):
                        continue

                strain_id = self.idhash['strain'].get(strain_key)
                if strain_id is None:
                    strain_id = self._makeInternalIdentifier(
                        'strain', strain_key)
                genotype_id = self.idhash['genotype'].get(genotype_key)
                if genotype_id is None:
                    genotype_id = self._makeInternalIdentifier(
                        'genotype', genotype_key)

                if strain_id is not None and genotype_id is not None:
                    self.strain_to_genotype_map[strain_id] = genotype_id

                graph.addTriple(strain_id, self.globaltt['has_genotype'], genotype_id)
                # TODO
                # verify if this should be contingent on the exactness or not
                # if qualifier == 'Exact':
                #     gu.addTriple(
                #       graph, strain_id,
                #       self.globaltt['has_genotype'],
                #       genotype_id)
                # else:
                #     gu.addXref(graph, strain_id, genotype_id)

                if not self.test_mode and limit is not None and line_counter > limit:
                    break

        return
677102
def scan(self, seqs, nreport=100, scan_rc=True, normalize=False):
        

        if not self.threshold:
            sys.stderr.write(
                "Using default threshold of 0.95. "
                "This is likely not optimal!\n"
                )
            self.set_threshold(threshold=0.95)

        seqs = as_fasta(seqs, genome=self.genome)
           
        it = self._scan_sequences(seqs.seqs, 
                    nreport, scan_rc)
       
        if normalize:
            if len(self.meanstd) == 0:
                self.set_meanstd()
            mean_std = [self.meanstd.get(m_id) for m_id in self.motif_ids]
            means = [x[0] for x in  mean_std]
            stds = [x[1] for x in  mean_std]


        for result in it:
            if normalize:
                zresult = [] 
                for i,mrow in enumerate(result):
                    mrow = [((x[0] - means[i]) / stds[i], x[1], x[2]) for x in mrow]
                    zresult.append(mrow)
                yield zresult
            else:
                yield result
609551
def codon2weight(self, codon):
        
        length = len(codon)
        retval = int(codon)
        return retval/(10 ** (length - 1)) - 5.0
759390
def absolute_abundance(coverage, total_bases):
	
	absolute = {}
	for genome in coverage:
		absolute[genome] = []
		index = 0
		for calc in coverage[genome]:
			bases = calc[0]
			total = total_bases[index]
			absolute[genome].append((bases / total) * float(100))
			index += 1
	total_assembled = [0 for i in absolute[genome]]
	for genome in absolute:
		index = 0
		for cov in absolute[genome]:
			total_assembled[index] += cov
			index += 1
	absolute['Unassembled'] = [(100 - i) for i in total_assembled]
	return absolute
1241333
def find_contig_distribution(contig_lengths_dict):
    
    # Initialise the dictionary
    contig_len_dist_dict = dict()
    for file_name, contig_lengths in contig_lengths_dict.items():
        # Initialise integers to store the number of contigs that fall into the different bin sizes
        over_1000000 = 0
        over_500000 = 0
        over_100000 = 0
        over_50000 = 0
        over_10000 = 0
        over_5000 = 0
        other = 0
        for contig_length in contig_lengths:
            # Depending on the size of the contig, increment the appropriate integer
            if contig_length > 1000000:
                over_1000000 += 1
            elif contig_length > 500000:
                over_500000 += 1
            elif contig_length > 100000:
                over_100000 += 1
            elif contig_length > 50000:
                over_50000 += 1
            elif contig_length > 10000:
                over_10000 += 1
            elif contig_length > 5000:
                over_5000 += 1
            else:
                other += 1
        # Populate the dictionary with a tuple of each of the size range frequencies
        contig_len_dist_dict[file_name] = (over_1000000,
                                           over_500000,
                                           over_100000,
                                           over_50000,
                                           over_10000,
                                           over_5000,
                                           other)
    return contig_len_dist_dict
1150930
def N50(arr):
  
  if len(arr) == 0:
    sys.stderr.write("ERROR: no content in array to take N50\n")
    sys.exit()
  tot = sum(arr)
  half = float(tot)/float(2)
  cummulative = 0
  for l in sorted(arr):
    cummulative += l
    if float(cummulative) > half: 
      return l
  sys.stderr.write("ERROR: problem finding M50\n")
  sys.exit()
871287
def n_to_g(self, n_interval):
        """"""

        assert self.strand in [1, -1], "strand = " + str(self.strand) + "; must be 1 or -1"

        if self.strand == 1:
            frs, fre = _hgvs_coord_to_ci(n_interval.start.base, n_interval.end.base)
            start_offset, end_offset = n_interval.start.offset, n_interval.end.offset
        elif self.strand == -1:
            frs, fre = _hgvs_coord_to_ci(n_interval.start.base, n_interval.end.base)
            fre, frs = self.tgt_len - frs, self.tgt_len - fre
            start_offset, end_offset = self.strand * n_interval.end.offset, self.strand * n_interval.start.offset

        # returns the genomic range start (grs) and end (gre)
        grs, gre = self.im.map_tgt_to_ref(frs, fre, max_extent=False)
        grs, gre = grs + self.gc_offset, gre + self.gc_offset
        gs, ge = grs + start_offset, gre + end_offset
        return hgvs.location.Interval(
            start=hgvs.location.SimplePosition(
                _ci_to_hgvs_coord(gs, ge)[0], uncertain=n_interval.start.uncertain),
            end=hgvs.location.SimplePosition(
                _ci_to_hgvs_coord(gs, ge)[1], uncertain=n_interval.end.uncertain),
            uncertain=n_interval.uncertain)
194472
def mw(self):
        '''

        '''
        counter = collections.Counter(self.seq.lower())
        mw_a = counter['a'] * 313.2
        mw_t = counter['t'] * 304.2
        mw_g = counter['g'] * 289.2
        mw_c = counter['c'] * 329.2
        mw_u = counter['u'] * 306.2

        if self.material == 'dna':
            return mw_a + mw_t + mw_g + mw_c + 79.0
        else:
            return mw_a + mw_u + mw_g + mw_c + 159.0
707298
def calcHeterodimer(seq1, seq2, mv_conc=50, dv_conc=0, dntp_conc=0.8,
                     dna_conc=50, temp_c=37, max_loop=30, temp_only=False):
    ''' 
    '''
    return calcThermo(seq1, seq2, 'ANY', mv_conc, dv_conc, dntp_conc,
                      dna_conc, temp_c, max_loop, temp_only)
1213808
def combine(*rnf_profiles):
        

        for rnf_profile in rnf_profiles:
            self.prefix_width = max(self.prefix_width, rnf_profile.prefix_width)
            self.read_tuple_id_width = max(self.read_tuple_id_width, rnf_profile.read_tuple_id_width)
            self.genome_id_width = max(self.genome_id_width, rnf_profile.genome_id_width)
            self.chr_id_width = max(self.chr_id_width, rnf_profile.chr_id_width)
            self.coor_width = max(self.coor_width, rnf_profile.coor_width)