Calculate Jensen-Shannon divergence between query and reference gene expression profile.
500920
def njsd(network, ref_gene_expression_dict, query_gene_expression_dict, gene_set):
    
    gene_jsd_dict = dict()
    reference_genes = ref_gene_expression_dict.keys()
    assert len(reference_genes) != 'Reference gene expression profile should have > 0 genes.'

    for gene in gene_set:
        if gene not in network.nodes:
            continue

        neighbors = find_neighbors(network, gene)

        query_expression_vec = get_neighbor_expression_vector(neighbors, query_gene_expression_dict)
        ref_expression_vec = get_neighbor_expression_vector(neighbors, ref_gene_expression_dict)
        assert len(query_expression_vec) == len(ref_expression_vec), 'Topology of reference network and query network differs. Please check.'

        # A gene which has non-expressed neighbors is ignored.
        if np.sum(query_expression_vec) == 0 and np.sum(ref_expression_vec) == 0:
            continue

        query_p_vec = exp2prob(query_expression_vec)
        ref_p_vec = exp2prob(ref_expression_vec)

        gene_jsd_dict[gene] = jsd(query_p_vec, ref_p_vec)

    return np.mean(list(gene_jsd_dict.values()))
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
339880
def calculate(bam_file, data, sv_bed):
    
    params = {"min": dd.get_coverage_depth_min(data)}
    variant_regions = dd.get_variant_regions_merged(data)
    if not variant_regions:
        variant_regions = _create_genome_regions(data)
    # Back compatible with previous pre-mosdepth callable files
    callable_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                                 dd.get_sample_name(data))),
                                 "%s-coverage.callable.bed" % (dd.get_sample_name(data)))
    if not utils.file_uptodate(callable_file, bam_file):
        vr_quantize = ("0:1:%s:" % (params["min"]), ["NO_COVERAGE", "LOW_COVERAGE", "CALLABLE"])
        to_calculate = [("variant_regions", variant_regions,
                         vr_quantize, None, "coverage_perbase" in dd.get_tools_on(data)),
                        ("sv_regions", bedutils.clean_file(sv_bed, data, prefix="svregions-"),
                         None, None, False),
                        ("coverage", bedutils.clean_file(dd.get_coverage(data), data, prefix="cov-"),
                         None, DEPTH_THRESHOLDS, False)]
        depth_files = {}
        for target_name, region_bed, quantize, thresholds, per_base in to_calculate:
            if region_bed:
                cur_depth = {}
                depth_info = run_mosdepth(data, target_name, region_bed, quantize=quantize, thresholds=thresholds,
                                          per_base=per_base)
                for attr in ("dist", "regions", "thresholds", "per_base"):
                    val = getattr(depth_info, attr, None)
                    if val:
                        cur_depth[attr] = val
                depth_files[target_name] = cur_depth
                if target_name == "variant_regions":
                    callable_file = depth_info.quantize
    else:
        depth_files = {}
    final_callable = _subset_to_variant_regions(callable_file, variant_regions, data)
    return final_callable, depth_files
692744
def _run_bcbio_variation(vrn_file, rm_file, rm_interval_file, base_dir, sample, caller, data):
    
    val_config_file = _create_validate_config_file(vrn_file, rm_file, rm_interval_file,
                                                   base_dir, data)
    work_dir = os.path.join(base_dir, "work")
    out = {"summary": os.path.join(work_dir, "validate-summary.csv"),
           "grading": os.path.join(work_dir, "validate-grading.yaml"),
           "discordant": os.path.join(work_dir, "%s-eval-ref-discordance-annotate.vcf" % sample)}
    if not utils.file_exists(out["discordant"]) or not utils.file_exists(out["grading"]):
        bcbio_variation_comparison(val_config_file, base_dir, data)
    out["concordant"] = filter(os.path.exists,
                                [os.path.join(work_dir, "%s-%s-concordance.vcf" % (sample, x))
                                 for x in ["eval-ref", "ref-eval"]])[0]
    return out
394069
def profiler(self):
        
        # Initialise variables
        profiledata = defaultdict(make_dict)
        profileset = set()
        genedict = {}
        # Find all the unique profiles to use with a set
        for sample in self.metadata.samples:
            if sample[self.analysistype].profile != 'NA':
                profileset.add(sample[self.analysistype].profile[0])
        # Extract the profiles for each set
        for sequenceprofile in profileset:
            # Clear the list of genes
            genelist = []
            for sample in self.metadata.samples:
                if sequenceprofile == sample[self.analysistype].profile[0]:
                    genelist = [os.path.split(x)[1].split('.')[0] for x in sample[self.analysistype].alleles]
            try:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Revert to standard comma separated values
            except KeyError:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Iterate through the rows
            for row in profile:
                # Iterate through the genes
                for gene in genelist:
                    # Add the sequence profile, and type, the gene name and the allele number to the dictionary
                    try:
                        profiledata[sequenceprofile][row['ST']][gene] = row[gene]
                    except KeyError:
                        pass
            # Add the gene list to a dictionary
            genedict[sequenceprofile] = sorted(genelist)
            # Add the profile data, and gene list to each sample
            for sample in self.metadata.samples:
                if sample.general.bestassemblyfile != 'NA':
                    if sequenceprofile == sample[self.analysistype].profile[0]:
                        # Populate the metadata with the profile data
                        sample[self.analysistype].profiledata = profiledata[sample[self.analysistype].profile[0]]
                        # Add the allele directory to a list of directories used in this analysis
                        self.allelefolders.add(sample[self.analysistype].alleledir)
                        dotter()
1353755
def _get_alignment_ranges(self):
    """"""
    if not self.is_aligned(): return None
    alignment_ranges = []
    cig = [x[:] for x in self.cigar_array]
    target_pos = self.entries.pos
    query_pos = 1
    while len(cig) > 0:
      c = cig.pop(0)
      if re.match('[S]$',c[1]): # hard or soft clipping
        query_pos += c[0]
      elif re.match('[ND]$',c[1]): # deleted from reference
        target_pos += c[0]
      elif re.match('[I]$',c[1]): # insertion to the reference
        query_pos += c[0]
      elif re.match('[MI=X]$',c[1]): # keep it
        t_start = target_pos
        q_start = query_pos
        target_pos += c[0]
        query_pos += c[0]
        t_end = target_pos-1
        q_end = query_pos-1
        alignment_ranges.append([GenomicRange(self.entries.rname,t_start,t_end),GenomicRange(self.entries.qname,q_start,q_end)])
    return alignment_ranges
1049117
def rnaseqc_general_stats (self):
        
        headers = OrderedDict()
        headers['Expression Profiling Efficiency'] = {
            'title': '% Expression Efficiency',
            'description': 'Expression Profiling Efficiency: Ratio of exon reads to total reads',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'YlGn',
            'modify': lambda x: float(x) * 100.0
        }
        headers['Genes Detected'] = {
            'title': '# Genes',
            'description': 'Number of genes detected with at least 5 reads.',
            'min': 0,
            'scale': 'Bu',
            'format': '{:,.0f}'
        }
        headers['rRNA rate'] = {
            'title': '% rRNA Alignment',
            'description': ' rRNA reads (non-duplicate and duplicate reads) per total reads',
            'max': 100,
            'min': 0,
            'suffix': '%',
            'scale': 'Reds',
            'modify': lambda x: float(x) * 100.0
        }

        self.general_stats_addcols(self.rna_seqc_metrics, headers)
371979
def inferred_var_seqs_plus_flanks(self, ref_seq, flank_length):
        ''''''
        flank_start = max(0, self.POS - flank_length)
        flank_end = min(len(ref_seq) - 1, self.ref_end_pos() + flank_length)
        seqs = [ref_seq[flank_start:self.POS] + self.REF + ref_seq[self.ref_end_pos() + 1: flank_end + 1]]

        for alt in self.ALT:
            seqs.append(ref_seq[flank_start:self.POS] + alt + ref_seq[self.ref_end_pos() + 1: flank_end + 1])

        return flank_start, seqs
710441
def get_kl_loss(self, means, log_vars, means_p=None, log_vars_p=None):
    """"""
    kl_loss = 0.0
    if means_p is None:
      means_p = tf.unstack(tf.zeros_like(means))
    if log_vars_p is None:
      log_vars_p = tf.unstack(tf.zeros_like(log_vars))
    enumerated_inputs = enumerate(zip(means, log_vars, means_p, log_vars_p))
    if self.is_training and self.hparams.stochastic_model:
      for i, (mean, log_var, mean_p, log_var_p) in enumerated_inputs:
        kl_loss += common_layers.kl_divergence(mean, log_var, mean_p, log_var_p)
        tf.summary.histogram("posterior_mean_%d" % i, mean)
        tf.summary.histogram("posterior_log_var_%d" % i, log_var)
        tf.summary.histogram("prior_mean_%d" % i, mean_p)
        tf.summary.histogram("prior_log_var_%d" % i, log_var_p)
      tf.summary.scalar("kl_raw", tf.reduce_mean(kl_loss))

    beta = self.get_beta(kl_loss)
    # information capacity from "Understanding disentangling in beta-VAE"
    if self.hparams.information_capacity > 0.0:
      kl_loss = tf.abs(kl_loss - self.hparams.information_capacity)
    return beta * kl_loss
1208950
def get_PSL(self,min_intron_size=68):
    
    from seqtools.format.psl import PSL
    matches = sum([x[0].length for x in self.alignment_ranges]) # 1. Matches - Number of matching bases that aren't repeats
    stats = AlignmentStats(matches,0,0,0,0,0,0,0)
    sub = self.query_sequence
    ref = self._options.reference
    if ref and sub:
      stats = self._analyze_alignment(min_intron_size=min_intron_size)
    strand = self.strand # 9. strand 
    qName = self.alignment_ranges[0][1].chr # 10. qName - Query sequence name
    qSize = self.query_sequence_length
    qStart = self.alignment_ranges[0][1].start-1
    qEnd = self.alignment_ranges[-1][1].end
    tName = self.alignment_ranges[0][0].chr
    tSize = self.target_sequence_length
    tStart = self.alignment_ranges[0][0].start-1
    tEnd = self.alignment_ranges[-1][0].end
    blockCount = len(self.alignment_ranges)
    blockSizes = ','.join([str(x[0].length) for x in self.alignment_ranges])+','
    qStarts = ','.join([str(x[1].start-1) for x in self.alignment_ranges])+','
    tStarts = ','.join([str(x[0].start-1) for x in self.alignment_ranges])+','

    psl_string = str(stats.matches)+"\t"+\
    str(stats.misMatches)+"\t"+\
    str(stats.repMatches)+"\t"+\
    str(stats.nCount)+"\t"+\
    str(stats.qNumInsert)+"\t"+\
    str(stats.qBaseInsert)+"\t"+\
    str(stats.tNumInsert)+"\t"+\
    str(stats.tBaseInsert)+"\t"+\
    strand+"\t"+\
    qName+"\t"+\
    str(qSize)+"\t"+\
    str(qStart)+"\t"+\
    str(qEnd)+"\t"+\
    tName+"\t"+\
    str(tSize)+"\t"+\
    str(tStart)+"\t"+\
    str(tEnd)+"\t"+\
    str(blockCount)+"\t"+\
    blockSizes+"\t"+\
    qStarts+"\t"+\
    tStarts
    return PSL(psl_string,PSL.Options(query_sequence=self.query_sequence,reference=self._options.reference,query_quality=self.query_quality))
432595
def validate_mutations(vcffile, bams, reffa, chr_reffa, samples, output_format, outfile,
                       mpileup_parameters=mpileup_parser.MPILEUP_DEFAULT_PARAMS):
    """"""
    output_header = "sample chrom pos ref cov A C G T * - + " \
        "val_ref val_alt val_al_type val_al_count val_maf "\
        "most_common_indel most_common_indel_count most_common_indel_maf most_common_indel_type most_common_al " \
        "most_common_al_count most_common_al_maf most_common_count most_common_maf".split()

    # for backwards compatibility
    # if bam or samples is a string, convert to list instead
    if isinstance(samples, six.string_types):
        samples = [samples]
    if isinstance(bams, six.string_types):
        bams = [bams]

    if output_format == 'vcf':
        vcf_reader = vcf.Reader(open(vcffile))
        vcf_reader.samples = samples
        vcf_reader.formats['GT'] = vcf.parser._Format(id='GT', num=1, type='String', desc="Genotype")
        vcf_reader.formats['AD'] = vcf.parser._Format(id='AD', num='R', type='Integer', desc="Allelic depth")
        vcf_reader.formats['DP'] = vcf.parser._Format(id='DP', num=1, type='Integer', desc="Depth")
        vcf_writer = vcf.Writer(outfile, vcf_reader)
    else:
        vcf_reader = open(vcffile)


    if output_format == "sufam":
        outfile.write("\t".join(output_header))
        outfile.write("\n")
    for record in vcf_reader:
        if output_format != 'vcf':
            line = record
            if line.startswith("#CHROM"):
                header = line[1:].rstrip('\n').split("\t")
                # create spoof pyvcf record if vcf_reader is not used
                _Record = namedtuple('Record', header)
            if line.startswith("#"):
                continue
            if len(header) == 0:
                raise(Exception("No header found in vcf file #CHROM not found"))
            # zip all column values, except alt (needs to be list in pyvcf)
            record_args = dict(zip(header, line.rstrip('\n').split("\t")))
            record_args['ALT'] = [record_args['ALT']]
            record = _Record(**record_args)

        # determine type of mutation
        record_type = "snv"
        if len(record.ALT) > 1:
            warnings.warn("Multiple ALT in one record is not implemented - using first")
        if len(record.REF) > len(record.ALT[0]):
            record_type = "deletion"
        elif len(record.ALT[0]) > len(record.REF):
            record_type = "insertion"

        # no coverage results
        no_cov = pd.Series({
            "chrom": str(record.CHROM), "pos": str(record.POS),
            "ref": str(record.REF),
            "cov": 0, "A": 0, "C": 0, "G": 0, "T": 0,
            "val_ref": str(record.REF), "val_alt": str(record.ALT[0]),
            "val_al_type": record_type, "val_al_count": 0, "val_maf": 0})

        # collect mpileup baseparser results per bam
        bps = []
        for i, bam in enumerate(bams):
            sample = samples[i]
            no_cov['sample'] = sample
            bp_lines = mpileup_parser.run_and_parse(bam, str(record.CHROM), str(record.POS), str(record.POS), reffa,
                                                    chr_reffa, mpileup_parameters)
            bpdf = get_baseparser_extended_df(sample, bp_lines, str(record.REF), str(record.ALT[0]))
            if bpdf is None:
                bp = no_cov
            else:
                bp = bpdf.ix[0, :]
            bps += [bp]

        # output call
        if output_format == "vcf":
            _write_bp_vcf(outfile, bps, vcf_writer, record)
        else:
            # only one bam file supported for outputs other than vcf
            _write_bp(outfile, bps[0], output_header, output_format)
1074404
def align_transcriptome(fastq_file, pair_file, ref_file, data):
    
    work_bam = dd.get_work_bam(data)
    base, ext = os.path.splitext(work_bam)
    out_file = base + ".transcriptome" + ext
    if utils.file_exists(out_file):
        data = dd.set_transcriptome_bam(data, out_file)
        return data
    # bwa mem needs phred+33 quality, so convert if it is Illumina
    if dd.get_quality_format(data).lower() == "illumina":
        logger.info("bwa mem does not support the phred+64 quality format, "
                    "converting %s and %s to phred+33.")
        fastq_file = fastq.groom(fastq_file, data, in_qual="fastq-illumina")
        if pair_file:
            pair_file = fastq.groom(pair_file, data, in_qual="fastq-illumina")
    bwa = config_utils.get_program("bwa", data["config"])
    gtf_file = dd.get_gtf_file(data)
    gtf_fasta = index_transcriptome(gtf_file, ref_file, data)
    args = " ".join(_bwa_args_from_config(data["config"]))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    samtools = config_utils.get_program("samtools", data["config"])
    cmd = ("{bwa} mem {args} -a -t {num_cores} {gtf_fasta} {fastq_file} "
           "{pair_file} ")
    with file_transaction(data, out_file) as tx_out_file:
        message = "Aligning %s and %s to the transcriptome." % (fastq_file, pair_file)
        cmd += "| " + postalign.sam_to_sortbam_cl(data, tx_out_file, name_sort=True)
        do.run(cmd.format(**locals()), message)
    data = dd.set_transcriptome_bam(data, out_file)
    return data
876452
def find_features(seqs, locus_tag="all", utr_len=200):
    """"""
    found_features = []

    for seq_i in seqs:
        for feature in seq_i.features:
            if feature.type == "CDS" and (locus_tag == "all" or \
                    ('locus_tag' in feature.qualifiers and \
                    feature.qualifiers['locus_tag'][0] == locus_tag)):
                start = max(0, feature.location.nofuzzy_start - utr_len)
                stop  = max(0, feature.location.nofuzzy_end + utr_len)
                feature_seq = seq_i.seq[start:stop]
                f_match = FeatureMatch(feature, feature_seq, feature.strand,
                                       utr_len)
                found_features.append(f_match)

    return found_features
667708
def sequencetyper(self):
        """"""
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Initialise variables
                    header = 0
                    # Iterate through the genomes
                    # for sample in self.metadata:
                    genome = sample.name
                    # Initialise self.bestmatch[genome] with an int that will eventually be replaced by the # of matches
                    self.bestmatch[genome] = defaultdict(int)
                    if sample[self.analysistype].profile != 'NA':
                        # Create the profiledata variable to avoid writing self.profiledata[self.analysistype]
                        profiledata = sample[self.analysistype].profiledata
                        # For each gene in plusdict[genome]
                        for gene in sample[self.analysistype].allelenames:
                            # Clear the appropriate count and lists
                            multiallele = []
                            multipercent = []
                            # Go through the alleles in plusdict
                            for allele in self.plusdict[genome][gene]:
                                percentid = list(self.plusdict[genome][gene][allele].keys())[0]
                                # "N" alleles screw up the allele splitter function
                                if allele != "N":
                                    # Use the alleleSplitter function to get the allele number
                                    # allelenumber, alleleprenumber = allelesplitter(allele)
                                    # Append as appropriate - alleleNumber is treated as an integer for proper sorting
                                    multiallele.append(int(allele))
                                    multipercent.append(percentid)
                                # If the allele is "N"
                                else:
                                    # Append "N" and a percent identity of 0
                                    multiallele.append("N")
                                    multipercent.append(0)
                                if not multiallele:
                                    multiallele.append("N")
                                    multipercent.append(0)
                            # if self.analysistype == 'rmlst':
                            #     # For whatever reason, the rMLST profile scheme treat multiple allele hits as 'N's.
                            #     multiallele = multiallele if len(multiallele) == 1 else ['N']
                            #     if multipercent:
                            #         multipercent = multipercent if len(multiallele) == 1 else [0, 0]
                            #     else:
                            #         multipercent = [0]
                            # Populate self.bestdict with genome, gene, alleles joined with a space (this was made like
                            # this because allele is a list generated by the .iteritems() above
                            self.bestdict[genome][gene][" ".join(str(allele)
                                                                 for allele in sorted(multiallele))] = multipercent[0]
                            # Find the profile with the most alleles in common with the query genome
                            for sequencetype in profiledata:
                                # The number of genes in the analysis
                                header = len(profiledata[sequencetype])
                                # refallele is the allele number of the sequence type
                                refallele = profiledata[sequencetype][gene]
                                # If there are multiple allele matches for a gene in the reference profile e.g. 10 692
                                if len(refallele.split(" ")) > 1:
                                    # Map the split (on a space) alleles as integers - if they are treated as integers,
                                    # the alleles will sort properly
                                    intrefallele = map(int, refallele.split(" "))
                                    # Create a string of the joined, sorted alleles
                                    sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                else:
                                    # Use the reference allele as the sortedRefAllele
                                    sortedrefallele = refallele
                                for allele, percentid in self.bestdict[genome][gene].items():
                                    # If the allele in the query genome matches the allele in the reference profile, add
                                    # the result to the bestmatch dictionary. Genes with multiple alleles were sorted
                                    # the same, strings with multiple alleles will match: 10 692 will never be 692 10
                                    if allele == sortedrefallele and float(percentid) == 100.00:
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                                    # Special handling of BACT000060 and BACT000065 genes. When the reference profile
                                    # has an allele of 'N', and the query allele doesn't, set the allele to 'N', and
                                    # count it as a match
                                    elif gene == 'BACT000060' or gene == 'BACT000065':
                                        if sortedrefallele == 'N' and allele != 'N':
                                            # Increment the number of matches to each profile
                                            self.bestmatch[genome][sequencetype] += 1
                                    elif allele == sortedrefallele and sortedrefallele == 'N':
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                        # Get the best number of matches
                        # From: https://stackoverflow.com/questions/613183/sort-a-python-dictionary-by-value
                        try:
                            sortedmatches = sorted(self.bestmatch[genome].items(), key=operator.itemgetter(1),
                                                   reverse=True)[0][1]
                        # If there are no matches, set :sortedmatches to zero
                        except IndexError:
                            sortedmatches = 0
                        # Otherwise, the query profile matches the reference profile
                        if int(sortedmatches) == header:
                            # Iterate through best match
                            for sequencetype, matches in self.bestmatch[genome].items():
                                if matches == sortedmatches:
                                    for gene in profiledata[sequencetype]:
                                        # Populate resultProfile with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                            list(self.bestdict[genome][gene]
                                                .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                    sample[self.analysistype].sequencetype = sequencetype
                                    sample[self.analysistype].matchestosequencetype = matches
                        # If there are fewer matches than the total number of genes in the typing scheme
                        elif 0 < int(sortedmatches) < header:
                            mismatches = []
                            # Iterate through the sequence types and the number of matches in bestDict for each genome
                            for sequencetype, matches in self.bestmatch[genome].items():
                                # If the number of matches for a profile matches the best number of matches
                                if matches == sortedmatches:
                                    # Iterate through the gene in the analysis
                                    for gene in profiledata[sequencetype]:
                                        # Get the reference allele as above
                                        refallele = profiledata[sequencetype][gene]
                                        # As above get the reference allele split and ordered as necessary
                                        if len(refallele.split(" ")) > 1:
                                            intrefallele = map(int, refallele.split(" "))
                                            sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                        else:
                                            sortedrefallele = refallele
                                        # Populate self.mlstseqtype with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        if self.updateprofile:
                                            self.mlstseqtype[genome][sequencetype][sortedmatches][gene][
                                                str(list(self.bestdict[genome][gene]
                                                    .keys())[0])][sortedrefallele] = str(list(self.bestdict[genome][gene]
                                                                                        .values())[0])
                                        else:
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene].keys())[0]] \
                                                = str(list(self.bestdict[genome][gene])[0])
                                                #= str(list(self.bestdict[genome][gene].values())[0])
                                            if sortedrefallele != list(self.bestdict[sample.name][gene].keys())[0]:
                                                mismatches.append(
                                                    ({gene: ('{} ({})'.format(list(self.bestdict[sample.name][gene]
                                                                              .keys())[0], sortedrefallele))}))
                                        if not self.updateprofile or self.analysistype == 'mlst':
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene]
                                                    .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                            sample[self.analysistype].mismatchestosequencetype = mismatches
                                            sample[self.analysistype].sequencetype = sequencetype
                                            sample[self.analysistype].matchestosequencetype = matches
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile and self.analysistype != 'mlst':
                                self.reprofiler(int(header), genome, sample)
                        elif sortedmatches == 0:
                            for gene in sample[self.analysistype].allelenames:
                                # Populate the results profile with negative values for sequence type and sorted matches
                                self.resultprofile[genome]['NA'][sortedmatches][gene]['NA'] = 0
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile:
                                self.reprofiler(int(header), genome, sample)
                            sample[self.analysistype].sequencetype = 'NA'
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                        else:
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                            sample[self.analysistype].sequencetype = 'NA'
                        dotter()
                else:
                    sample[self.analysistype].matchestosequencetype = 'NA'
                    sample[self.analysistype].mismatchestosequencetype = 'NA'
                    sample[self.analysistype].sequencetype = 'NA'

            else:
                sample[self.analysistype].matchestosequencetype = 'NA'
                sample[self.analysistype].mismatchestosequencetype = 'NA'
                sample[self.analysistype].sequencetype = 'NA'
317284
def qualNorm_filter_genes(data, qualitative, pval_threshold=0.05, smoothing=1e-5, eps=1e-5):
    
    genes, cells = data.shape
    clusters = qualitative.shape[1]
    output = np.zeros((genes, clusters))
    missing_indices = []
    genes_included = []
    qual_indices = []
    thresholds = qualitative.min(1) + (qualitative.max(1) - qualitative.min(1))/2.0
    pvals = np.zeros(genes)
    for i in range(genes):
        if qualitative[i,:].max() == -1 and qualitative[i,:].min() == -1:
            missing_indices.append(i)
            continue
        qual_indices.append(i)
        threshold = thresholds[i]
        data_i = data[i,:]
        if sparse.issparse(data):
            data_i = data_i.toarray().flatten()
        assignments, means = poisson_cluster(data_i.reshape((1, cells)), 2)
        means = means.flatten()
        high_i = 1
        low_i = 0
        if means[0]>means[1]:
            high_i = 0
            low_i = 1
        # do a p-value test
        p_val = poisson_test(data_i[assignments==low_i], data_i[assignments==high_i], smoothing=smoothing)
        pvals[i] = p_val
        if p_val <= pval_threshold:
            genes_included.append(i)
        else:
            continue
        high_mean = np.median(data_i[assignments==high_i])
        low_mean = np.median(data_i[assignments==low_i]) + eps
        for k in range(clusters):
            if qualitative[i,k]>threshold:
                output[i,k] = high_mean
            else:
                output[i,k] = low_mean
    output = output[genes_included,:]
    pvals = pvals[genes_included]
    return output, pvals, genes_included
936848
def variation(self,parents):
        """"""
        # downselect to features that are important
        if (self.valid(parents) and
            self.ml_type != 'SVC' and self.ml_type != 'SVR'):
            # this is needed because svm has a bug that throws valueerror on
            # attribute check
            if hasattr(self.pipeline.named_steps['ml'],'coef_'):
                # for l1 regularization, filter individuals with 0 coefficients
                if self.weight_parents:
                    weights = self.pipeline.named_steps['ml'].coef_
                    if len(weights.shape)>1: # handle multi-coefficient models
                        weights = [np.mean(abs(c)) for c in weights.transpose()]
                    # softmax transformation of the weights
                    weights = np.exp(weights)/np.sum(np.exp(weights))
                    offspring = copy.deepcopy(
                        list(self.random_state.choice(self.valid(parents),
                                              self.population_size, p=weights)))
                else:
                    offspring = copy.deepcopy(list(
                        x for i,x in zip(self.pipeline.named_steps['ml'].coef_,
                                         self.valid(parents)) if  (i != 0).any()))
            elif hasattr(self.pipeline.named_steps['ml'],'feature_importances_'):
                # for tree methods, filter our individuals with 0 feature importance
                if self.weight_parents:
                    weights = self.pipeline.named_steps['ml'].feature_importances_
                    # softmax transformation of the weights
                    weights = np.exp(weights)/np.sum(np.exp(weights))
                    offspring = copy.deepcopy(list(
                        self.random_state.choice(self.valid(parents),
                                         self.population_size, p=weights)))
                else:
                    offspring = copy.deepcopy(list(
                        x for i,x in zip(self.pipeline.named_steps['ml'].feature_importances_,
                                         self.valid(parents)) if  i != 0))
            else:
                offspring = copy.deepcopy(self.valid(parents))
        else:
            offspring = copy.deepcopy(self.valid(parents))

        if self.elitism: # keep a copy of the elite individual
            elite_index = np.argmin([x.fitness for x in parents])
            elite = copy.deepcopy(parents[elite_index])

        # Apply crossover and mutation on the offspring
        if self.verbosity > 2: print("variation...")
        for child1, child2 in it.zip_longest(offspring[::2], offspring[1::2],
                                             fillvalue=None):

            if self.random_state.rand() < self.crossover_rate and child2 != None:
            # crossover
                self.cross(child1.stack, child2.stack, self.max_depth)
                # update ids
                child1.parentid = [child1.id,child2.id]
                child1.id = uuid.uuid4()
                child2.parentid = [child1.id,child2.id]
                child2.id = uuid.uuid4()
                # set default fitness
                child1.fitness = -1
                child2.fitness = -1
            elif child2 == None:
            # single mutation
                self.mutate(child1.stack,self.func_set,self.term_set)
                # update ids
                child1.parentid = [child1.id]
                child1.id = uuid.uuid4()
                # set default fitness
                child1.fitness = -1
            else:
            #double mutation
                self.mutate(child1.stack,self.func_set,self.term_set)
                self.mutate(child2.stack,self.func_set,self.term_set)
                # update ids
                child1.parentid = [child1.id]
                child1.id = uuid.uuid4()
                child2.parentid = [child2.id]
                child2.id = uuid.uuid4()
                # set default fitness
                child1.fitness = -1
                child2.fitness = -1

        while len(offspring) < self.population_size:
            #make new offspring to replace the invalid ones
            offspring.append(Ind())
            self.make_program(offspring[-1].stack,self.func_set,self.term_set,
                              self.random_state.randint(self.min_depth,
                                                        self.max_depth+1),
                              self.otype)
            offspring[-1].stack = list(reversed(offspring[-1].stack))
        
        return offspring,elite,elite_index
142000
def profiler(self):
        """"""
        # Initialise variables
        profiledata = defaultdict(make_dict)
        profileset = set()
        # supplementalset = ''
        genedict = {}
        # Find all the unique profiles to use with a set
        for sample in self.metadata:
            if sample[self.analysistype].profile != 'NA':
                profileset.add(sample[self.analysistype].profile[0])
                # if self.analysistype == 'rmlst':
                #     supplementalset = sample[self.analysistype].supplementalprofile
        # Extract the profiles for each set
        for sequenceprofile in profileset:
            # Clear the list of genes
            genelist = []
            for sample in self.metadata:
                if sequenceprofile == sample[self.analysistype].profile[0]:
                    # genelist = [os.path.split(x)[1].split('.')[0] for x in sample[self.analysistype].alleles]
                    genelist = sample[self.analysistype].allelenames
            try:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile), dialect='excel-tab')
            # Revert to standard comma separated values
            except KeyError:
                # Open the sequence profile file as a dictionary
                profile = DictReader(open(sequenceprofile))
            # Iterate through the rows
            for row in profile:
                # Iterate through the genes
                for gene in genelist:
                    # Add the sequence profile, and type, the gene name and the allele number to the dictionary
                    try:
                        profiledata[sequenceprofile][row['ST']][gene] = row[gene]
                    except KeyError:
                        try:
                            profiledata[sequenceprofile][row['rST']][gene] = row[gene]
                        except KeyError:
                            raise

            # # Load the supplemental profile definitions
            # if self.analysistype == 'rmlst':
            #     supplementalprofile = DictReader(open(supplementalset), dialect='excel-tab')
            #     # Do the same with the supplemental profile
            #     for row in supplementalprofile:
            #         # Iterate through the genes
            #         for gene in genelist:
            #             # Add the sequence profile, and type, the gene name and the allele number to the dictionary
            #             profiledata[sequenceprofile][row['rST']][gene] = row[gene]
            # Add the gene list to a dictionary
            genedict[sequenceprofile] = sorted(genelist)
            # Add the profile data, and gene list to each sample
            for sample in self.metadata:
                if sample.general.bestassemblyfile != 'NA':
                    if sequenceprofile == sample[self.analysistype].profile[0]:
                        # Populate the metadata with the profile data
                        sample[self.analysistype].profiledata = profiledata[sample[self.analysistype].profile[0]]
                        # Add the allele directory to a list of directories used in this analysis
                        self.allelefolders.add(sample[self.analysistype].alleledir)
                        dotter()