diff options
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/lmm.py | 107 |
1 files changed, 64 insertions, 43 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py index e5978933..3344b024 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py @@ -59,68 +59,89 @@ def run_human(pheno_vector, # Buffers for pvalues and t-stats p_values = [] t_stats = [] - + plink_input.getSNPIterator() total_snps = plink_input.numSNPs with Bench("snp iterator loop"): count = 0 for snp, this_id in plink_input: - if count > 5000: + if count > 10000: break count += 1 - percent_complete = (float(count) / total_snps) * 100 #print("percent_complete: ", percent_complete) temp_data.store("percent_complete", percent_complete) + + ps, ts = human_association(snp, + n, + keep, + lmm_ob, + pheno_vector, + covariate_matrix, + kinship_matrix, + refit) - x = snp[keep].reshape((n,1)) - #x[[1,50,100,200,3000],:] = np.nan - v = np.isnan(x).reshape((-1,)) - - # Check SNPs for missing values - if v.sum(): - keeps = True - v - xs = x[keeps,:] - # If no variation at this snp or all genotypes missing - if keeps.sum() <= 1 or xs.var() <= 1e-6: - p_values.append(np.nan) - t_stats.append(np.nan) - continue - - # Its ok to center the genotype - I used options.normalizeGenotype to - # force the removal of missing genotypes as opposed to replacing them with MAF. - - #if not options.normalizeGenotype: - # xs = (xs - xs.mean()) / np.sqrt(xs.var()) - - filtered_pheno = pheno_vector[keeps] - filtered_covariate_matrix = covariate_matrix[keeps,:] - filtered_kinship_matrix = kinship_matrix[keeps,:][:,keeps] - filtered_lmm_ob = lmm.LMM(filtered_pheno,filtered_kinship_matrix,X0=filtered_covariate_matrix) - if refit: - filtered_lmm_ob.fit(X=xs) - else: - #try: - filtered_lmm_ob.fit() - #except: pdb.set_trace() - ts,ps,beta,betaVar = Ls.association(xs,returnBeta=True) - else: - if x.var() == 0: - p_values.append(np.nan) - t_stats.append(np.nan) - continue - if refit: - lmm_ob.fit(X=x) - ts, ps, beta, betaVar = lmm_ob.association(x) - p_values.append(ps) t_stats.append(ts) return p_values, t_stats +def human_association(snp, + n, + keep, + lmm_ob, + pheno_vector, + covariate_matrix, + kinship_matrix, + refit): + + x = snp[keep].reshape((n,1)) + #x[[1,50,100,200,3000],:] = np.nan + v = np.isnan(x).reshape((-1,)) + + # Check SNPs for missing values + if v.sum(): + keeps = True - v + xs = x[keeps,:] + # If no variation at this snp or all genotypes missing + if keeps.sum() <= 1 or xs.var() <= 1e-6: + return np.nan, np.nan + #p_values.append(np.nan) + #t_stats.append(np.nan) + #continue + + # Its ok to center the genotype - I used options.normalizeGenotype to + # force the removal of missing genotypes as opposed to replacing them with MAF. + + #if not options.normalizeGenotype: + # xs = (xs - xs.mean()) / np.sqrt(xs.var()) + + filtered_pheno = pheno_vector[keeps] + filtered_covariate_matrix = covariate_matrix[keeps,:] + filtered_kinship_matrix = kinship_matrix[keeps,:][:,keeps] + filtered_lmm_ob = lmm.LMM(filtered_pheno,filtered_kinship_matrix,X0=filtered_covariate_matrix) + if refit: + filtered_lmm_ob.fit(X=xs) + else: + #try: + filtered_lmm_ob.fit() + #except: pdb.set_trace() + ts,ps,beta,betaVar = Ls.association(xs,returnBeta=True) + else: + if x.var() == 0: + return np.nan, np.nan + #p_values.append(np.nan) + #t_stats.append(np.nan) + #continue + if refit: + lmm_ob.fit(X=x) + ts, ps, beta, betaVar = lmm_ob.association(x) + return ps, ts + + def run(pheno_vector, genotype_matrix, restricted_max_likelihood=True, |