Filter genes by fano factor and mean.

    See `filter_genes_dispersion`.

    Reference: Weinreb et al. (2017).
345984
def filter_genes_fano_deprecated(X, Ecutoff, Vcutoff):
    
    if issparse(X):
        raise ValueError('Not defined for sparse input. See `filter_genes_dispersion`.')
    mean_filter = np.mean(X, axis=0) > Ecutoff
    var_filter = np.var(X, axis=0) / (np.mean(X, axis=0) + .0001) > Vcutoff
    gene_subset = np.nonzero(np.all([mean_filter, var_filter], axis=0))[0]
    return gene_subset
76415
def zip_fit_params(data):
    
    genes, cells = data.shape
    m = data.mean(1)
    v = data.var(1)
    M = (v-m)/(m**2+v-m)
    #M = v/(v+m**2)
    #M[np.isnan(M)] = 0.0
    M = np.array([min(1.0, max(0.0, x)) for x in M])
    L = m + v/m - 1.0
    #L = (v + m**2)/m
    L[np.isnan(L)] = 0.0
    L = np.array([max(0.0, x) for x in L])
    return L, M
234561
def _filter_by_background(base_name, back_samples, gt_vcfs, data):
    
    filtname = "InBackground"
    filtdoc = "Variant also present in background samples with same genotype"
    orig_vcf = gt_vcfs[base_name]
    out_file = "%s-backfilter.vcf" % (utils.splitext_plus(orig_vcf)[0])
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with file_transaction(data, out_file) as tx_out_file:
            with utils.open_gzipsafe(orig_vcf) as in_handle:
                inp = vcf.Reader(in_handle, orig_vcf)
                inp.filters[filtname] = vcf.parser._Filter(filtname, filtdoc)
                with open(tx_out_file, "w") as out_handle:
                    outp = vcf.Writer(out_handle, inp)
                    for rec in inp:
                        if _genotype_in_background(rec, base_name, back_samples):
                            rec.add_filter(filtname)
                        outp.write_record(rec)
    if utils.file_exists(out_file + ".gz"):
        out_file = out_file + ".gz"
    gt_vcfs[base_name] = vcfutils.bgzip_and_index(out_file, data["config"])
    return gt_vcfs
588736
def run_filter(vrn_file, align_bam, ref_file, data, items):
    
    if not should_filter(items) or not vcfutils.vcf_has_variants(vrn_file):
        return data
    else:
        raw_file = "%s-damage.vcf" % utils.splitext_plus(vrn_file)[0]
        out_plot_files = ["%s%s" % (utils.splitext_plus(raw_file)[0], ext)
                          for ext in ["_seq_bias_simplified.pdf", "_pcr_bias_simplified.pdf"]]
        if not utils.file_uptodate(raw_file, vrn_file) and not utils.file_uptodate(raw_file + ".gz", vrn_file):
            with file_transaction(items[0], raw_file) as tx_out_file:
                # Does not apply --qcSummary plotting due to slow runtimes
                cmd = ["dkfzbiasfilter.py", "--filterCycles", "1", "--passOnly",
                       "--tempFolder", os.path.dirname(tx_out_file),
                       vrn_file, align_bam, ref_file, tx_out_file]
                do.run(cmd, "Filter low frequency variants for DNA damage and strand bias")
                for out_plot in out_plot_files:
                    tx_plot_file = os.path.join("%s_qcSummary" % utils.splitext_plus(tx_out_file)[0], "plots",
                                                os.path.basename(out_plot))
                    if utils.file_exists(tx_plot_file):
                        shutil.move(tx_plot_file, out_plot)
        raw_file = vcfutils.bgzip_and_index(raw_file, items[0]["config"])
        data["vrn_file"] = _filter_to_info(raw_file, items[0])
        out_plot_files = [x for x in out_plot_files if utils.file_exists(x)]
        data["damage_plots"] = out_plot_files
        return data
652340
def get_proteome_correct_percentages(prots_filtered_feathers, outpath, length_filter_pid=None,
                                     copynum_scale=False, copynum_df=None,
                                     force_rerun=False):
    """"""
    if ssbio.utils.force_rerun(flag=force_rerun, outfile=outpath):
        prot_tracker = defaultdict(int)
        big_strain_counts_df = pd.DataFrame()
        first = True
        for feather in prots_filtered_feathers:
            loaded = load_feather(protein_feather=feather, length_filter_pid=length_filter_pid,
                                  copynum_scale=copynum_scale,
                                  copynum_df=copynum_df)

            if first:
                big_strain_counts_df = pd.DataFrame(columns=loaded.columns)
                first = False
            tmp_df = pd.DataFrame(columns=loaded.columns)
            for strain in loaded.columns:
                prot_tracker[strain] += 1
                totals = list(filter(lambda x: x.endswith('total'), loaded[strain].index))
                for t in totals:
                    counts = t.rsplit('_', 1)[0]
                    aa_counts = list(
                        filter(lambda x: (x.startswith(counts) and x not in totals), loaded[strain].index))
                    for aa_count in aa_counts:
                        tmp_df.at[aa_count.replace('count', '%'), strain] = loaded[strain][aa_count] / \
                                                                            loaded[strain][t]
            big_strain_counts_df = big_strain_counts_df.add(tmp_df, fill_value=0)

        for c, total in prot_tracker.items():
            big_strain_counts_df.loc[:, c] /= total

        if len(big_strain_counts_df) > 0:
            big_strain_counts_df.astype(float).reset_index().to_feather(outpath)
        return big_strain_counts_df
    else:
        return pd.read_feather(outpath).set_index('index')
659923
def filter_butter_coeffs(filtertype, freq, samplerate, order=5):
    # type: (str, Union[float, Tuple[float, float]], int, int) -> Tuple[np.ndarray, np.ndarray]
    
    assert filtertype in ('low', 'high', 'band')
    nyq = 0.5 * samplerate
    if isinstance(freq, tuple):
        assert filtertype == 'band'
        low, high = freq
        low  /= nyq
        high /= nyq
        b, a = signal.butter(order, [low, high], btype='band')
    else:
        freq = freq / nyq
        b, a = signal.butter(order, freq, btype=filtertype)
    return b, a
53613
def fdfilter(data, *filt, **kwargs):
    
    # parse keyword args
    inplace = kwargs.pop('inplace', False)
    analog = kwargs.pop('analog', False)
    fs = kwargs.pop('sample_rate', None)
    if kwargs:
        raise TypeError("filter() got an unexpected keyword argument '%s'"
                        % list(kwargs.keys())[0])

    # parse filter
    if fs is None:
        fs = 2 * (data.shape[-1] * data.df).to('Hz').value
    form, filt = parse_filter(filt, analog=analog, sample_rate=fs)
    lti = signal.lti(*filt)

    # generate frequency response
    freqs = data.frequencies.value.copy()
    fresp = numpy.nan_to_num(abs(lti.freqresp(w=freqs)[1]))

    # apply to array
    if inplace:
        data *= fresp
        return data
    new = data * fresp
    return new
786361
def get_proteome_counts_impute_missing(prots_filtered_feathers, outpath, length_filter_pid=None,
                                       copynum_scale=False, copynum_df=None,
                                       force_rerun=False):
    """"""
    if ssbio.utils.force_rerun(flag=force_rerun, outfile=outpath):
        big_strain_counts_df = pd.DataFrame()
        first = True
        for feather in prots_filtered_feathers:
            loaded = load_feather(protein_feather=feather, length_filter_pid=length_filter_pid,
                                  copynum_scale=copynum_scale,
                                  copynum_df=copynum_df)
            if first:
                big_strain_counts_df = pd.DataFrame(index=_all_counts, columns=loaded.columns)
                first = False

            new_columns = list(set(loaded.columns.tolist()).difference(big_strain_counts_df.columns))
            if new_columns:
                for col in new_columns:
                    big_strain_counts_df[col] = big_strain_counts_df.mean(axis=1)

            not_in_loaded = list(set(big_strain_counts_df.columns).difference(loaded.columns.tolist()))
            if not_in_loaded:
                for col in not_in_loaded:
                    big_strain_counts_df[col] = big_strain_counts_df[col] + loaded.mean(axis=1)

            big_strain_counts_df = big_strain_counts_df.add(loaded, fill_value=0)

        if len(big_strain_counts_df) > 0:
            big_strain_counts_df.astype(float).reset_index().to_feather(outpath)
        return big_strain_counts_df
    else:
        return pd.read_feather(outpath).set_index('index')
667372
def find_nb_genes(data):
    
    data_means = data.mean(1)
    data_vars = data.var(1)
    nb_indices = data_means < 0.9*data_vars
    return nb_indices
707419
def random_filter(objects, reduction_factor, seed=42):
    
    assert 0 < reduction_factor <= 1, reduction_factor
    rnd = random.Random(seed)
    out = []
    for obj in objects:
        if rnd.random() <= reduction_factor:
            out.append(obj)
    return out
571195
def fit_mle(self, data, fix_mean=False):
        

        if not fix_mean:
            sigma, _, scale = stats.lognorm.fit(data, floc=0)
            return np.log(scale), sigma

        else:
            mean = np.mean(data)

            # MLE fxn to be optmimized
            mle = lambda sigma, x, mean: -1 *\
                                    np.sum(self._pdf_w_mean(x, mean, sigma))

            sigma = optim.fmin(mle, np.array([np.std(np.log(data), ddof=1)]),
                                            args=(data, mean), disp=0)[0]

            return self.translate_args(mean, sigma)
1216622
def count_missense_per_gene(lines):
    
    
    counts = {}
    for x in lines:
        x = x.split("\t")
        gene = x[0]
        consequence = x[3]
        
        if gene not in counts:
            counts[gene] = 0
        
        if consequence != "missense_variant":
            continue
        
        counts[gene] += 1
    
    return counts
864005
def _sigmainf(N, h, m, dW, Km0, Pm0):
    """"""
    M = m*(m-1)//2
    Im = broadcast_to(np.eye(m), (N, m, m))
    IM = broadcast_to(np.eye(M), (N, M, M))
    Ims0 = np.eye(m**2)
    factor1 = broadcast_to((2.0/h)*np.dot(Km0, Ims0 - Pm0), (N, M, m**2))
    factor2 = _kp2(Im, _dot(dW, _t(dW)))
    factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
    return 2*IM + _dot(_dot(factor1, factor2), factor3)
1156483
def _cnvkit_prioritize(sample, genes, allele_file, metrics_file):
    
    mdf = pd.read_table(metrics_file)
    mdf.columns = [x.lower() for x in mdf.columns]
    if len(genes) > 0:
        mdf = mdf[mdf["gene"].str.contains("|".join(genes))]
    mdf = mdf[["chromosome", "start", "end", "gene", "log2", "ci_hi", "ci_lo"]]
    adf = pd.read_table(allele_file)
    if len(genes) > 0:
        adf = adf[adf["gene"].str.contains("|".join(genes))]
    if "cn1" in adf.columns and "cn2" in adf.columns:
        adf = adf[["chromosome", "start", "end", "cn", "cn1", "cn2"]]
    else:
        adf = adf[["chromosome", "start", "end", "cn"]]
    df = pd.merge(mdf, adf, on=["chromosome", "start", "end"])
    df = df[df["cn"] != 2]
    if len(df) > 0:
        def passes(row):
            spread = abs(row["ci_hi"] - row["ci_lo"])
            return spread < 0.25
        df["passes"] = df.apply(passes, axis=1)
    df.insert(0, "sample", [sample] * len(df))
    return df
1118889
def parseNoise(rawArray): 
    '''
    
    '''
    nnGenes=[] 
    for i in range(0,(rawArray.shape[1])): #Checks all genes
        count0=np.asarray(np.where(rawArray[:,i]==0))
        count1=np.asarray(np.where(rawArray[:,i]==1))
        if ((count1.shape[1]+count0.shape[1])<rawArray.shape[0]):
            nnGenes=np.append(nnGenes,i)
    return nnGenes.astype(int)
627645
def addRandomEffect(self,K=None,covar_type='freeform',is_noise=False,normalize=True,Ks=None,offset=1e-4,rank=1,covar_K0=None):
        
        if self.P==1:	self.addSingleTraitTerm(K=K,is_noise=is_noise,normalize=normalize,Ks=Ks)
        else:			self.addMultiTraitTerm(K=K,covar_type=covar_type,is_noise=is_noise,normalize=normalize,Ks=Ks,offset=offset,rank=rank,covar_K0=covar_K0)
1242754
def savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0):
    '''
    
    '''
    # linearly interpolate missing values before filtering
    x = np.ma.masked_invalid(pd.Series(x).interpolate())

    try:
        # start filtering from the first non-zero value since these won't be addressed by
        # the interpolation above
        ind = np.isfinite(x).nonzero()[0][0]
        x[ind:] = signal.savgol_filter(x[ind:], window_length, polyorder, deriv,
                                       delta, axis, mode, cval)
    except IndexError:
        pass
    return np.ma.masked_invalid(x)