diff options
Diffstat (limited to 'wqflask')
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/__init__.py | 1 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/benchmark.py | 44 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/chunks.py | 96 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py | 184 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/genotype.py | 51 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/gn2.py | 110 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/gwas.py | 165 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/input.py | 267 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/kinship.py | 168 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/lmm.py | 995 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/lmm2.py | 433 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py | 55 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/phenotype.py | 65 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/plink.py | 107 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 229 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/standalone.py | 110 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/temp_data.py | 25 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py | 122 |
18 files changed, 0 insertions, 3227 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py b/wqflask/wqflask/my_pylmm/pyLMM/__init__.py deleted file mode 100644 index f33c4e74..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/__init__.py +++ /dev/null @@ -1 +0,0 @@ -PYLMM_VERSION="0.51-gn2" diff --git a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py b/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py deleted file mode 100644 index 6c6c9f88..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/benchmark.py +++ /dev/null @@ -1,44 +0,0 @@ -from __future__ import print_function, division, absolute_import - -import collections -import inspect -import time - -class Bench(object): - entries = collections.OrderedDict() - - def __init__(self, name=None): - self.name = name - - def __enter__(self): - if self.name: - print("Starting benchmark: %s" % (self.name)) - else: - print("Starting benchmark at: %s [%i]" % (inspect.stack()[1][3], inspect.stack()[1][2])) - self.start_time = time.time() - - def __exit__(self, type, value, traceback): - if self.name: - name = self.name - else: - name = "That" - - time_taken = time.time() - self.start_time - print(" %s took: %f seconds" % (name, (time_taken))) - - if self.name: - Bench.entries[self.name] = Bench.entries.get(self.name, 0) + time_taken - - - @classmethod - def report(cls): - total_time = sum((time_taken for time_taken in cls.entries.itervalues())) - print("\nTiming report\n") - for name, time_taken in cls.entries.iteritems(): - percent = int(round((time_taken/total_time) * 100)) - print("[{}%] {}: {}".format(percent, name, time_taken)) - print() - - def reset(cls): - """Reset the entries""" - cls.entries = collections.OrderedDict() diff --git a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py b/wqflask/wqflask/my_pylmm/pyLMM/chunks.py deleted file mode 100644 index 9565fb96..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/chunks.py +++ /dev/null @@ -1,96 +0,0 @@ -from __future__ import absolute_import, print_function, division - -import math -import time - - -def divide_into_chunks(the_list, number_chunks): - """Divides a list into approximately number_chunks smaller lists - - >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 3) - [[1, 2, 7], [3, 22, 8], [5, 22, 333]] - >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 4) - [[1, 2, 7], [3, 22, 8], [5, 22, 333]] - >>> divide_into_chunks([1, 2, 7, 3, 22, 8, 5, 22, 333], 5) - [[1, 2], [7, 3], [22, 8], [5, 22], [333]] - >>> - - """ - length = len(the_list) - - if length == 0: - return [[]] - - if length <= number_chunks: - number_chunks = length - - chunksize = int(math.ceil(length / number_chunks)) - - chunks = [] - for counter in range(0, length, chunksize): - chunks.append(the_list[counter:counter+chunksize]) - - return chunks - -def _confirm_chunk(original, result): - all_chunked = [] - for chunk in result: - all_chunked.extend(chunk) - print("length of all chunked:", len(all_chunked)) - assert original == all_chunked, "You didn't chunk right" - - -def _chunk_test(divide_func): - import random - random.seed(7) - - number_exact = 0 - total_amount_off = 0 - - for test in range(1, 1001): - print("\n\ntest:", test) - number_chunks = random.randint(1, 20) - number_elements = random.randint(0, 100) - the_list = list(range(1, number_elements)) - result = divide_func(the_list, number_chunks) - - print("Dividing list of length {} into approximately {} chunks - got {} chunks".format( - len(the_list), number_chunks, len(result))) - print("result:", result) - - _confirm_chunk(the_list, result) - - amount_off = abs(number_chunks - len(result)) - if amount_off == 0: - number_exact += 1 - else: - total_amount_off += amount_off - - - print("\n{} exact out of {} [Total amount off: {}]".format(number_exact, - test, - total_amount_off)) - assert number_exact == 558 - assert total_amount_off == 1580 - return number_exact, total_amount_off - - -def _main(): - info = dict() - #funcs = (("sam", sam_divide_into_chunks), ("zach", zach_divide_into_chunks)) - funcs = (("only one", divide_into_chunks),) - for name, func in funcs: - start = time.time() - number_exact, total_amount_off = _chunk_test(func) - took = time.time() - start - info[name] = dict(number_exact=number_exact, - total_amount_off=total_amount_off, - took=took) - - print("info is:", info) - -if __name__ == '__main__': - _main() - print("\nConfirming doctests...") - import doctest - doctest.testmod() diff --git a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py deleted file mode 100644 index 4312fed0..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/convertlmm.py +++ /dev/null @@ -1,184 +0,0 @@ -# This is a converter for common LMM formats, so as to keep file -# reader complexity outside the main routines. - -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -from __future__ import print_function -from optparse import OptionParser -import sys -import os -import numpy as np -# from lmm import LMM, run_other -# import input -import plink - -usage = """ -python convertlmm.py [--plink] [--prefix out_basename] [--kinship kfile] [--pheno pname] [--geno gname] - - Convert files for runlmm.py processing. Writes to stdout by default. - - try --help for more information - -Examples: - - python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos > test.pheno - - python ./my_pylmm/pyLMM/convertlmm.py --plink --pheno data/test_snps.132k.clean.noX.fake.phenos --geno data/test_snps.132k.clean.noX > test.geno -""" - -# if len(args) == 0: -# print usage -# sys.exit(1) - -option_parser = OptionParser(usage=usage) -option_parser.add_option("--kinship", dest="kinship", - help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program") -option_parser.add_option("--pheno", dest="pheno", - help="Parse a phenotype file (use with --plink only)") -option_parser.add_option("--geno", dest="geno", - help="Parse a genotype file (use with --plink only)") -option_parser.add_option("--plink", dest="plink", action="store_true", default=False, - help="Parse PLINK style") -# option_parser.add_option("--kinship",action="store_false", dest="kinship", default=True, -# help="Parse a kinship file. This is an nxn plain text file and can be computed with the pylmmKinship program.") -option_parser.add_option("--prefix", dest="prefix", - help="Output prefix for output file(s)") -option_parser.add_option("-q", "--quiet", - action="store_false", dest="verbose", default=True, - help="don't print status messages to stdout") -option_parser.add_option("-v", "--verbose", - action="store_true", dest="verbose", default=False, - help="Print extra info") - -(options, args) = option_parser.parse_args() - -writer = None -num_inds = None -snp_names = [] -ind_names = [] - -def msg(s): - sys.stderr.write("INFO: ") - sys.stderr.write(s) - sys.stderr.write("\n") - -def wr(s): - if writer: - writer.write(s) - else: - sys.stdout.write(s) - -def wrln(s): - wr(s) - wr("\n") - - -if options.pheno: - if not options.plink: - raise Exception("Use --plink switch") - # Because plink does not track size we need to read the whole thing first - msg("Converting pheno "+options.pheno) - phenos = [] - count = 0 - count_pheno = None - for line in open(options.pheno,'r'): - count += 1 - list = line.split() - pcount = len(list)-2 - assert(pcount > 0) - if count_pheno == None: - count_pheno = pcount - assert(count_pheno == pcount) - row = [list[0]]+list[2:] - phenos.append(row) - - writer = None - if options.prefix: - writer = open(options.prefix+".pheno","w") - wrln("# Phenotype format version 1.0") - wrln("# Individuals = "+str(count)) - wrln("# Phenotypes = "+str(count_pheno)) - for i in range(count_pheno): - wr("\t"+str(i+1)) - wr("\n") - for i in range(count): - wr("\t".join(phenos[i])) - wr("\n") - num_inds = count - msg(str(count)+" pheno lines written") - -if options.kinship: - is_header = True - count = 0 - msg("Converting kinship "+options.kinship) - writer = None - if options.prefix: - writer = open(options.prefix+".kin","w") - for line in open(options.kinship,'r'): - count += 1 - if is_header: - size = len(line.split()) - wrln("# Kinship format version 1.0") - wrln("# Size="+str(size)) - for i in range(size): - wr("\t"+str(i+1)) - wr("\n") - is_header = False - wr(str(count)) - wr("\t") - wr("\t".join(line.split())) - wr("\n") - num_inds = count - msg(str(count)+" kinship lines written") - -if options.geno: - msg("Converting geno "+options.geno+'.bed') - if not options.plink: - raise Exception("Use --plink switch") - if not num_inds: - raise Exception("Can not figure out the number of individuals, use --pheno or --kinship") - bim_snps = plink.readbim(options.geno+'.bim') - num_snps = len(bim_snps) - writer = None - if options.prefix: - writer = open(options.prefix+".geno","w") - wrln("# Genotype format version 1.0") - wrln("# Individuals = "+str(num_inds)) - wrln("# SNPs = "+str(num_snps)) - wrln("# Encoding = HAB") - for i in range(num_inds): - wr("\t"+str(i+1)) - wr("\n") - - m = [] - def out(i,x): - # wr(str(i)+"\t") - # wr("\t".join(x)) - # wr("\n") - m.append(x) - - snps = plink.readbed(options.geno+'.bed',num_inds, ('A','H','B','-'), out) - - msg("Write transposed genotype matrix") - for g in range(num_snps): - wr(bim_snps[g][1]+"\t") - for i in range(num_inds): - wr(m[g][i]) - wr("\n") - - msg(str(count)+" geno lines written (with "+str(snps)+" snps)") - -msg("Converting done") diff --git a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py b/wqflask/wqflask/my_pylmm/pyLMM/genotype.py deleted file mode 100644 index 49f32e3a..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/genotype.py +++ /dev/null @@ -1,51 +0,0 @@ -# Genotype routines - -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -import numpy as np -from collections import Counter -import operator - -def replace_missing_with_MAF(snp_g): - """ - Replace the missing genotype with the minor allele frequency (MAF) - in the snp row. It is rather slow! - """ - cnt = Counter(snp_g) - tuples = sorted(cnt.items(), key=operator.itemgetter(1)) - l2 = [t for t in tuples if not np.isnan(t[0])] - maf = l2[0][0] - res = np.array([maf if np.isnan(snp) else snp for snp in snp_g]) - return res - -def normalize(ind_g): - """ - Run for every SNP list (for one individual) and return - normalized SNP genotype values with missing data filled in - """ - g = np.copy(ind_g) # copy to avoid side effects - missing = np.isnan(g) - values = g[True - missing] - mean = values.mean() # Global mean value - stddev = np.sqrt(values.var()) # Global stddev - g[missing] = mean # Plug-in mean values for missing data - if stddev == 0: - g = g - mean # Subtract the mean - else: - g = (g - mean) / stddev # Normalize the deviation - return g - diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py b/wqflask/wqflask/my_pylmm/pyLMM/gn2.py deleted file mode 100644 index 821195c8..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py +++ /dev/null @@ -1,110 +0,0 @@ -# Standalone specific methods and callback handler -# -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# Set the log level with -# -# logging.basicConfig(level=logging.DEBUG) - -from __future__ import absolute_import, print_function, division - -import numpy as np -import sys -import logging - -# logger = logging.getLogger(__name__) -logger = logging.getLogger('lmm2') -logging.basicConfig(level=logging.DEBUG) -np.set_printoptions(precision=3,suppress=True) - -progress_location = None -progress_current = None -progress_prev_perc = None - -def progress_default_func(location,count,total): - global progress_current - value = round(count*100.0/total) - progress_current = value - -progress_func = progress_default_func - -def progress_set_func(func): - global progress_func - progress_func = func - -def progress(location, count, total): - global progress_location - global progress_prev_perc - - perc = round(count*100.0/total) - if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5): - progress_func(location, count, total) - logger.info("Progress: %s %d%%" % (location,perc)) - progress_location = location - progress_prev_perc = perc - -def mprint(msg,data): - """ - Array/matrix print function - """ - m = np.array(data) - if m.ndim == 1: - print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:]) - if m.ndim == 2: - print(msg,m.shape,"=\n[", - m[0][0:3]," ... ",m[0][-3:],"\n ", - m[1][0:3]," ... ",m[1][-3:],"\n ...\n ", - m[-2][0:3]," ... ",m[-2][-3:],"\n ", - m[-1][0:3]," ... ",m[-1][-3:],"]") - -def fatal(msg): - logger.critical(msg) - raise Exception(msg) - -def callbacks(): - return dict( - write = sys.stdout.write, - writeln = print, - debug = logger.debug, - info = logger.info, - warning = logger.warning, - error = logger.error, - critical = logger.critical, - fatal = fatal, - progress = progress, - mprint = mprint - ) - -def uses(*funcs): - """ - Some sugar - """ - return [callbacks()[func] for func in funcs] - -# ----- Minor test cases: - -if __name__ == '__main__': - # logging.basicConfig(level=logging.DEBUG) - logging.debug("Test %i" % (1)) - d = callbacks()['debug'] - d("TEST") - wrln = callbacks()['writeln'] - wrln("Hello %i" % 34) - progress = callbacks()['progress'] - progress("I am half way",50,100) - list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15] - mprint("list",list) - matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]] - mprint("matrix",matrix) - ix,dx = uses("info","debug") - ix("ix") - dx("dx") diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py b/wqflask/wqflask/my_pylmm/pyLMM/gwas.py deleted file mode 100644 index 247a8729..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/gwas.py +++ /dev/null @@ -1,165 +0,0 @@ -# pylmm-based GWAS calculation -# -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. -#!/usr/bin/python - -import pdb -import time -import sys -import lmm2 - -import os -import numpy as np -import input -from optmatrix import matrix_initialize -from lmm2 import LMM2 - -import multiprocessing as mp # Multiprocessing is part of the Python stdlib -import Queue - -# ---- A trick to decide on the environment: -try: - from wqflask.my_pylmm.pyLMM import chunks - from gn2 import uses -except ImportError: - has_gn2=False - from standalone import uses - -progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal') - - -def formatResult(id,beta,betaSD,ts,ps): - return "\t".join([str(x) for x in [id,beta,betaSD,ts,ps]]) + "\n" - -def compute_snp(j,n,snp_ids,lmm2,REML,q = None): - result = [] - for snp_id in snp_ids: - snp,id = snp_id - x = snp.reshape((n,1)) # all the SNPs - # if refit: - # L.fit(X=snp,REML=REML) - ts,ps,beta,betaVar = lmm2.association(x,REML=REML,returnBeta=True) - # result.append(formatResult(id,beta,np.sqrt(betaVar).sum(),ts,ps)) - result.append( (ts,ps) ) - if not q: - q = compute_snp.q - q.put([j,result]) - return j - -def f_init(q): - compute_snp.q = q - -def gwas(Y,G,K,restricted_max_likelihood=True,refit=False,verbose=True): - """ - GWAS. The G matrix should be n inds (cols) x m snps (rows) - """ - info("In gwas.gwas") - matrix_initialize() - cpu_num = mp.cpu_count() - numThreads = None # for now use all available threads - kfile2 = False - reml = restricted_max_likelihood - - mprint("G",G) - n = G.shape[1] # inds - inds = n - m = G.shape[0] # snps - snps = m - info("%s SNPs",snps) - assert snps>=inds, "snps should be larger than inds (snps=%d,inds=%d)" % (snps,inds) - - # CREATE LMM object for association - # if not kfile2: L = LMM(Y,K,Kva,Kve,X0,verbose=verbose) - # else: L = LMM_withK2(Y,K,Kva,Kve,X0,verbose=verbose,K2=K2) - - lmm2 = LMM2(Y,K) # ,Kva,Kve,X0,verbose=verbose) - if not refit: - info("Computing fit for null model") - lmm2.fit() # follow GN model in run_other - info("heritability=%0.3f, sigma=%0.3f" % (lmm2.optH,lmm2.optSigma)) - - res = [] - - # Set up the pool - # mp.set_start_method('spawn') - q = mp.Queue() - p = mp.Pool(numThreads, f_init, [q]) - collect = [] - - count = 0 - job = 0 - jobs_running = 0 - jobs_completed = 0 - for snp in G: - snp_id = (snp,'SNPID') - count += 1 - if count % 1000 == 0: - job += 1 - debug("Job %d At SNP %d" % (job,count)) - if numThreads == 1: - debug("Running on 1 THREAD") - compute_snp(job,n,collect,lmm2,reml,q) - collect = [] - j,lst = q.get() - debug("Job "+str(j)+" finished") - jobs_completed += 1 - progress("GWAS2",jobs_completed,snps/1000) - res.append((j,lst)) - else: - p.apply_async(compute_snp,(job,n,collect,lmm2,reml)) - jobs_running += 1 - collect = [] - while jobs_running > cpu_num: - try: - j,lst = q.get_nowait() - debug("Job "+str(j)+" finished") - jobs_completed += 1 - progress("GWAS2",jobs_completed,snps/1000) - res.append((j,lst)) - jobs_running -= 1 - except Queue.Empty: - time.sleep(0.1) - pass - if jobs_running > cpu_num*2: - time.sleep(1.0) - else: - break - - collect.append(snp_id) - - if numThreads==1 or count<1000 or len(collect)>0: - job += 1 - debug("Collect final batch size %i job %i @%i: " % (len(collect), job, count)) - compute_snp(job,n,collect,lmm2,reml,q) - collect = [] - j,lst = q.get() - res.append((j,lst)) - debug("count=%i running=%i collect=%i" % (count,jobs_running,len(collect))) - for job in range(jobs_running): - j,lst = q.get(True,15) # time out - debug("Job "+str(j)+" finished") - jobs_completed += 1 - progress("GWAS2",jobs_completed,snps/1000) - res.append((j,lst)) - - mprint("Before sort",[res1[0] for res1 in res]) - res = sorted(res,key=lambda x: x[0]) - mprint("After sort",[res1[0] for res1 in res]) - info([len(res1[1]) for res1 in res]) - ts = [item[0] for j,res1 in res for item in res1] - ps = [item[1] for j,res1 in res for item in res1] - return ts,ps diff --git a/wqflask/wqflask/my_pylmm/pyLMM/input.py b/wqflask/wqflask/my_pylmm/pyLMM/input.py deleted file mode 100644 index 7063fedf..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/input.py +++ /dev/null @@ -1,267 +0,0 @@ -# pylmm is a python-based linear mixed-model solver with applications to GWAS - -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) - -#The program is free for academic use. Please contact Nick Furlotte -#<nick.furlotte@gmail.com> if you are interested in using the software for -#commercial purposes. - -#The software must not be modified and distributed without prior -#permission of the author. - -#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -#A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR -#CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -#EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -#PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -#PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF -#LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING -#NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -#SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -import os -import sys -import numpy as np -import struct -import pdb - -class plink: - def __init__(self,fbase,kFile=None,phenoFile=None,type='b',normGenotype=True,readKFile=False): - self.fbase = fbase - self.type = type - self.indivs = self.getIndivs(self.fbase,type) - self.kFile = kFile - self.phenos = None - self.normGenotype = normGenotype - self.phenoFile = phenoFile - # Originally I was using the fastLMM style that has indiv IDs embedded. - # NOW I want to use this module to just read SNPs so I'm allowing - # the programmer to turn off the kinship reading. - self.readKFile = readKFile - - if self.kFile: - self.K = self.readKinship(self.kFile) - elif os.path.isfile("%s.kin" % fbase): - self.kFile = "%s.kin" %fbase - if self.readKFile: - self.K = self.readKinship(self.kFile) - else: - self.kFile = None - self.K = None - - self.getPhenos(self.phenoFile) - - self.fhandle = None - self.snpFileHandle = None - - def __del__(self): - if self.fhandle: self.fhandle.close() - if self.snpFileHandle: self.snpFileHandle.close() - - def getSNPIterator(self): - if not self.type == 'b': - sys.stderr.write("Have only implemented this for binary plink files (bed)\n") - return - - # get the number of snps - file = self.fbase + '.bim' - i = 0 - f = open(file,'r') - for line in f: i += 1 - f.close() - self.numSNPs = i - self.have_read = 0 - self.snpFileHandle = open(file,'r') - - self.BytestoRead = self.N / 4 + (self.N % 4 and 1 or 0) - self._formatStr = 'c'*self.BytestoRead - - file = self.fbase + '.bed' - self.fhandle = open(file,'rb') - - magicNumber = self.fhandle.read(2) - order = self.fhandle.read(1) - if not order == '\x01': - sys.stderr.write("This is not in SNP major order - you did not handle this case\n") - raise StopIteration - - return self - - def __iter__(self): - return self.getSNPIterator() - - def next(self): - if self.have_read == self.numSNPs: - raise StopIteration - X = self.fhandle.read(self.BytestoRead) - XX = [bin(ord(x)) for x in struct.unpack(self._formatStr,X)] - self.have_read += 1 - return self.formatBinaryGenotypes(XX,self.normGenotype),self.snpFileHandle.readline().strip().split()[1] - - def formatBinaryGenotypes(self,X,norm=True): - D = { \ - '00': 0.0, \ - '10': 0.5, \ - '11': 1.0, \ - '01': np.nan \ - } - - D_tped = { \ - '00': '1 1', \ - '10': '1 2', \ - '11': '2 2', \ - '01': '0 0' \ - } - - #D = D_tped - - G = [] - for x in X: - if not len(x) == 10: - xx = x[2:] - x = '0b' + '0'*(8 - len(xx)) + xx - a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4]) - L = [D[y] for y in [a,b,c,d]] - G += L - # only take the leading values because whatever is left should be null - G = G[:self.N] - G = np.array(G) - if norm: - G = self.normalizeGenotype(G) - return G - - def normalizeGenotype(self,G): - # print "Before",G - # print G.shape - print "call input.normalizeGenotype" - raise "This should not be used" - x = True - np.isnan(G) - m = G[x].mean() - s = np.sqrt(G[x].var()) - G[np.isnan(G)] = m - if s == 0: G = G - m - else: G = (G - m) / s - # print "After",G - return G - - def getPhenos(self,phenoFile=None): - if not phenoFile: - self.phenoFile = phenoFile = self.fbase+".phenos" - if not os.path.isfile(phenoFile): - sys.stderr.write("Could not find phenotype file: %s\n" % (phenoFile)) - return - f = open(phenoFile,'r') - keys = [] - P = [] - for line in f: - v = line.strip().split() - keys.append((v[0],v[1])) - P.append([(x == 'NA' or x == '-9') and np.nan or float(x) for x in v[2:]]) - f.close() - P = np.array(P) - - # reorder to match self.indivs - D = {} - L = [] - for i in range(len(keys)): - D[keys[i]] = i - for i in range(len(self.indivs)): - if not D.has_key(self.indivs[i]): - continue - L.append(D[self.indivs[i]]) - P = P[L,:] - - self.phenos = P - return P - - def getIndivs(self,base,type='b'): - if type == 't': - famFile = "%s.tfam" % base - else: - famFile = "%s.fam" % base - keys = [] - i = 0 - f = open(famFile,'r') - for line in f: - v = line.strip().split() - famId = v[0] - indivId = v[1] - k = (famId.strip(),indivId.strip()) - keys.append(k) - i += 1 - f.close() - - self.N = len(keys) - sys.stderr.write("Read %d individuals from %s\n" % (self.N, famFile)) - - return keys - - def readKinship(self,kFile): - # Assume the fastLMM style - # This will read in the kinship matrix and then reorder it - # according to self.indivs - additionally throwing out individuals - # that are not in both sets - if self.indivs == None or len(self.indivs) == 0: - sys.stderr.write("Did not read any individuals so can't load kinship\n") - return - - sys.stderr.write("Reading kinship matrix from %s\n" % (kFile) ) - - f = open(kFile,'r') - # read indivs - v = f.readline().strip().split("\t")[1:] - keys = [tuple(y.split()) for y in v] - D = {} - for i in range(len(keys)): D[keys[i]] = i - - # read matrix - K = [] - for line in f: - K.append([float(x) for x in line.strip().split("\t")[1:]]) - f.close() - K = np.array(K) - - # reorder to match self.indivs - L = [] - KK = [] - X = [] - for i in range(len(self.indivs)): - if not D.has_key(self.indivs[i]): - X.append(self.indivs[i]) - else: - KK.append(self.indivs[i]) - L.append(D[self.indivs[i]]) - K = K[L,:][:,L] - self.indivs = KK - self.indivs_removed = X - if len(self.indivs_removed): - sys.stderr.write("Removed %d individuals that did not appear in Kinship\n" % (len(self.indivs_removed))) - return K - - def getCovariates(self,covFile=None): - if not os.path.isfile(covFile): - sys.stderr.write("Could not find covariate file: %s\n" % (phenoFile)) - return - f = open(covFile,'r') - keys = [] - P = [] - for line in f: - v = line.strip().split() - keys.append((v[0],v[1])) - P.append([x == 'NA' and np.nan or float(x) for x in v[2:]]) - f.close() - P = np.array(P) - - # reorder to match self.indivs - D = {} - L = [] - for i in range(len(keys)): - D[keys[i]] = i - for i in range(len(self.indivs)): - if not D.has_key(self.indivs[i]): continue - L.append(D[self.indivs[i]]) - P = P[L,:] - - return P diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py deleted file mode 100644 index 1c157fd8..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py +++ /dev/null @@ -1,168 +0,0 @@ -# pylmm kinship calculation -# -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -# env PYTHONPATH=$pylmm_lib_path:./lib python $pylmm_lib_path/runlmm.py --pheno test.pheno --geno test9000.geno kinship --test - -import sys -import os -import numpy as np -from scipy import linalg -import multiprocessing as mp # Multiprocessing is part of the Python stdlib -import Queue -import time - -from optmatrix import matrix_initialize, matrixMultT - -# ---- A trick to decide on the environment: -try: - from wqflask.my_pylmm.pyLMM import chunks - from gn2 import uses, progress_set_func -except ImportError: - has_gn2=False - import standalone as handlers - from standalone import uses, progress_set_func - -progress,debug,info,mprint = uses('progress','debug','info','mprint') - -def kinship_full(G): - """ - Calculate the Kinship matrix using a full dot multiplication - """ - # mprint("kinship_full G",G) - m = G.shape[0] # snps - n = G.shape[1] # inds - info("%d SNPs",m) - assert m>n, "n should be larger than m (%d snps > %d inds)" % (m,n) - # m = np.dot(G.T,G) - m = matrixMultT(G.T) - m = m/G.shape[0] - # mprint("kinship_full K",m) - return m - -def compute_W(job,G,n,snps,compute_size): - """ - Read 1000 SNPs at a time into matrix and return the result - """ - m = compute_size - W = np.ones((n,m)) * np.nan # W matrix has dimensions individuals x SNPs (initially all NaNs) - for j in range(0,compute_size): - pos = job*m + j # real position - if pos >= snps: - W = W[:,range(0,j)] - break - snp = G[job*compute_size+j] - if snp.var() == 0: - continue - W[:,j] = snp # set row to list of SNPs - return W - -def compute_matrixMult(job,W,q = None): - """ - Compute Kinship(W)*j - - For every set of SNPs matrixMult is used to multiply matrices T(W)*W - """ - res = matrixMultT(W) - if not q: q=compute_matrixMult.q - q.put([job,res]) - return job - -def f_init(q): - compute_matrixMult.q = q - -# Calculate the kinship matrix from G (SNPs as rows!), returns K -# -def kinship(G,computeSize=1000,numThreads=None,useBLAS=False): - - matrix_initialize(useBLAS) - - mprint("G",G) - n = G.shape[1] # inds - inds = n - m = G.shape[0] # snps - snps = m - info("%i SNPs" % (m)) - assert snps>=inds, "snps should be larger than inds (%i snps, %i inds)" % (snps,inds) - - q = mp.Queue() - p = mp.Pool(numThreads, f_init, [q]) - cpu_num = mp.cpu_count() - info("CPU cores: %i" % cpu_num) - iterations = snps/computeSize+1 - - results = [] - K = np.zeros((n,n)) # The Kinship matrix has dimension individuals x individuals - - completed = 0 - for job in range(iterations): - info("Processing job %d first %d SNPs" % (job, ((job+1)*computeSize))) - W = compute_W(job,G,n,snps,computeSize) - if numThreads == 1: - # Single-core - compute_matrixMult(job,W,q) - j,x = q.get() - debug("Job "+str(j)+" finished") - progress("kinship",j,iterations) - K_j = x - K = K + K_j - else: - # Multi-core - results.append(p.apply_async(compute_matrixMult, (job,W))) - # Do we have a result? - while (len(results)-completed>cpu_num*2): - time.sleep(0.1) - try: - j,x = q.get_nowait() - debug("Job "+str(j)+" finished") - K_j = x - K = K + K_j - completed += 1 - progress("kinship",completed,iterations) - except Queue.Empty: - pass - - if numThreads == None or numThreads > 1: - for job in range(len(results)-completed): - j,x = q.get(True,15) - debug("Job "+str(j)+" finished") - K_j = x - K = K + K_j - completed += 1 - progress("kinship",completed,iterations) - - K = K / float(snps) - return K - -def kvakve(K): - """ - Obtain eigendecomposition for K and return Kva,Kve where Kva is cleaned - of small values < 1e-6 (notably smaller than zero) - """ - info("Obtaining eigendecomposition for %dx%d matrix" % (K.shape[0],K.shape[1]) ) - Kva,Kve = linalg.eigh(K) - mprint("Kva",Kva) - mprint("Kve",Kve) - - if sum(Kva < 0): - info("Cleaning %d eigen values (Kva<0)" % (sum(Kva < 0))) - Kva[Kva < 1e-6] = 1e-6 - return Kva,Kve - - - - diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py deleted file mode 100644 index 2a0c7fdc..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py +++ /dev/null @@ -1,995 +0,0 @@ -# pylmm is a python-based linear mixed-model solver with applications to GWAS - -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -from __future__ import absolute_import, print_function, division - -import sys -import time -import uuid - -import numpy as np -from scipy import linalg -from scipy import optimize -from scipy import stats -# import pdb - -# import gzip -# import zlib -import datetime -# import cPickle as pickle -from pprint import pformat as pf - -# Add local dir to PYTHONPATH -import os -cwd = os.path.dirname(__file__) -if sys.path[0] != cwd: - sys.path.insert(1,cwd) - -# pylmm imports -from kinship import kinship, kinship_full, kvakve -import genotype -import phenotype -import gwas -from benchmark import Bench - -# The following imports are for exchanging data with the webserver -import simplejson as json -from redis import Redis -Redis = Redis() -import temp_data - -has_gn2=None - -# sys.stderr.write("INFO: pylmm system path is "+":".join(sys.path)+"\n") -sys.stderr.write("INFO: pylmm file is "+__file__+"\n") - -# ---- A trick to decide on the environment: -try: - sys.stderr.write("INFO: lmm try loading module\n") - import utility.formatting # this is never used, just to check the environment - sys.stderr.write("INFO: This is a genenetwork2 environment\n") - from gn2 import uses, progress_set_func - has_gn2=True -except ImportError: - # Failed to load gn2 - has_gn2=False - import standalone as handlers - from standalone import uses, progress_set_func - sys.stderr.write("WARNING: LMM standalone version missing the Genenetwork2 environment\n") - -progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal') - -#np.seterr('raise') - -#def run_human(pheno_vector, -# covariate_matrix, -# plink_input_file, -# kinship_matrix, -# refit=False, -# loading_progress=None): - -def run_human(pheno_vector, - covariate_matrix, - plink_input_file, - kinship_matrix, - refit=False): - - v = np.isnan(pheno_vector) - keep = True - v - keep = keep.reshape((len(keep),)) - - identifier = str(uuid.uuid4()) - - #print("pheno_vector: ", pf(pheno_vector)) - #print("kinship_matrix: ", pf(kinship_matrix)) - #print("kinship_matrix.shape: ", pf(kinship_matrix.shape)) - - #lmm_vars = pickle.dumps(dict( - # pheno_vector = pheno_vector, - # covariate_matrix = covariate_matrix, - # kinship_matrix = kinship_matrix - #)) - #Redis.hset(identifier, "lmm_vars", lmm_vars) - #Redis.expire(identifier, 60*60) - - if v.sum(): - pheno_vector = pheno_vector[keep] - print("pheno_vector shape is now: ", pf(pheno_vector.shape)) - covariate_matrix = covariate_matrix[keep,:] - print("kinship_matrix shape is: ", pf(kinship_matrix.shape)) - print("keep is: ", pf(keep.shape)) - kinship_matrix = kinship_matrix[keep,:][:,keep] - - print("kinship_matrix:", pf(kinship_matrix)) - - n = kinship_matrix.shape[0] - print("n is:", n) - lmm_ob = LMM(pheno_vector, - kinship_matrix, - covariate_matrix) - lmm_ob.fit() - - - # Buffers for pvalues and t-stats - p_values = [] - t_stats = [] - - #print("input_file: ", plink_input_file) - - with Bench("Opening and loading pickle file"): - with gzip.open(plink_input_file, "rb") as input_file: - data = pickle.load(input_file) - - plink_input = data['plink_input'] - - #plink_input.getSNPIterator() - with Bench("Calculating numSNPs"): - total_snps = data['numSNPs'] - - with Bench("snp iterator loop"): - count = 0 - - with Bench("Create list of inputs"): - inputs = list(plink_input) - - with Bench("Divide into chunks"): - results = chunks.divide_into_chunks(inputs, 64) - - result_store = [] - - key = "plink_inputs" - - # Todo: Delete below line when done testing - Redis.delete(key) - - timestamp = datetime.datetime.utcnow().isoformat() - - # Pickle chunks of input SNPs (from Plink interator) and compress them - #print("Starting adding loop") - for part, result in enumerate(results): - #data = pickle.dumps(result, pickle.HIGHEST_PROTOCOL) - holder = pickle.dumps(dict( - identifier = identifier, - part = part, - timestamp = timestamp, - result = result - ), pickle.HIGHEST_PROTOCOL) - - #print("Adding:", part) - Redis.rpush(key, zlib.compress(holder)) - #print("End adding loop") - #print("***** Added to {} queue *****".format(key)) - for snp, this_id in plink_input: - #with Bench("part before association"): - #if count > 1000: - # break - count += 1 - progress("human",count,total_snps) - - #with Bench("actual association"): - ps, ts = human_association(snp, - n, - keep, - lmm_ob, - pheno_vector, - covariate_matrix, - kinship_matrix, - refit) - - #with Bench("after association"): - p_values.append(ps) - t_stats.append(ts) - - return p_values, t_stats - - -#class HumanAssociation(object): -# def __init__(self): -# - -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,:] - - print("kinship_matrix shape is: ", pf(kinship_matrix.shape)) - print("keeps is: ", pf(keeps.shape)) - 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, -# refit=False, -# temp_data=None): - -def run_other_old(pheno_vector, - genotype_matrix, - restricted_max_likelihood=True, - refit=False): - - """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics - - restricted_max_likelihood -- whether to use restricted max likelihood; True or False - refit -- whether to refit the variance component for each marker - - """ - - print("Running the original LMM engine in run_other (old)") - print("REML=",restricted_max_likelihood," REFIT=",refit) - with Bench("Calculate Kinship"): - kinship_matrix,genotype_matrix = calculate_kinship_new(genotype_matrix) - - print("kinship_matrix: ", pf(kinship_matrix)) - print("kinship_matrix.shape: ", pf(kinship_matrix.shape)) - - # with Bench("Create LMM object"): - # lmm_ob = LMM(pheno_vector, kinship_matrix) - - # with Bench("LMM_ob fitting"): - # lmm_ob.fit() - - print("run_other_old genotype_matrix: ", genotype_matrix.shape) - print(genotype_matrix) - - with Bench("Doing GWAS"): - t_stats, p_values = GWAS(pheno_vector, - genotype_matrix.T, - kinship_matrix, - restricted_max_likelihood=True, - refit=False) - Bench().report() - return p_values, t_stats - -def run_other_new(n,m,pheno_vector, - geno, - restricted_max_likelihood=True, - refit=False): - - """Takes the phenotype vector and genotype matrix and returns a set of p-values and t-statistics - - restricted_max_likelihood -- whether to use restricted max likelihood; True or False - refit -- whether to refit the variance component for each marker - - """ - - print("Running the new LMM2 engine in run_other_new") - print("REML=",restricted_max_likelihood," REFIT=",refit) - - # Adjust phenotypes - n,Y,keep = phenotype.remove_missing_new(n,pheno_vector) - - # if options.maf_normalization: - # G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g ) - # print "MAF replacements: \n",G - # if not options.skip_genotype_normalization: - # G = np.apply_along_axis( genotype.normalize, axis=1, arr=G) - - geno = geno[:,keep] - with Bench("Calculate Kinship"): - K,G = calculate_kinship_new(geno) - - print("kinship_matrix: ", pf(K)) - print("kinship_matrix.shape: ", pf(K.shape)) - - # with Bench("Create LMM object"): - # lmm_ob = lmm2.LMM2(Y,K) - # with Bench("LMM_ob fitting"): - # lmm_ob.fit() - - print("run_other_new genotype_matrix: ", G.shape) - print(G) - - with Bench("Doing GWAS"): - t_stats, p_values = gwas.gwas(Y, - G, - K, - restricted_max_likelihood=True, - refit=False,verbose=True) - Bench().report() - return p_values, t_stats - -# def matrixMult(A,B): -# return np.dot(A,B) - -def matrixMult(A,B): - - # If there is no fblas then we will revert to np.dot() - - try: - linalg.fblas - except AttributeError: - return np.dot(A,B) - - #print("A is:", pf(A.shape)) - #print("B is:", pf(B.shape)) - - # If the matrices are in Fortran order then the computations will be faster - # when using dgemm. Otherwise, the function will copy the matrix and that takes time. - if not A.flags['F_CONTIGUOUS']: - AA = A.T - transA = True - else: - AA = A - transA = False - - if not B.flags['F_CONTIGUOUS']: - BB = B.T - transB = True - else: - BB = B - transB = False - - return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB) - -def calculate_kinship_new(genotype_matrix): - """ - Call the new kinship calculation where genotype_matrix contains - inds (columns) by snps (rows). - """ - assert type(genotype_matrix) is np.ndarray - info("call genotype.normalize") - G = np.apply_along_axis( genotype.normalize, axis=1, arr=genotype_matrix) - mprint("G",genotype_matrix) - info("call calculate_kinship_new") - return kinship(G),G # G gets transposed, we'll turn this into an iterator (FIXME) - -def calculate_kinship_iter(geno): - """ - Call the new kinship calculation where genotype_matrix contains - inds (columns) by snps (rows). - """ - assert type(genotype_matrix) is iter - info("call genotype.normalize") - G = np.apply_along_axis( genotype.normalize, axis=0, arr=genotype_matrix) - info("call calculate_kinship_new") - return kinship(G) - -def calculate_kinship_old(genotype_matrix): - """ - genotype_matrix is an n x m matrix encoding SNP minor alleles. - - This function takes a matrix oF SNPs, imputes missing values with the maf, - normalizes the resulting vectors and returns the RRM matrix. - - """ - info("call calculate_kinship_old") - fatal("THE FUNCTION calculate_kinship_old IS OBSOLETE, use calculate_kinship_new instead - see Genotype Normalization Problem on Pjotr's blog") - n = genotype_matrix.shape[0] - m = genotype_matrix.shape[1] - info("genotype 2D matrix n (inds) is: %d" % (n)) - info("genotype 2D matrix m (snps) is: %d" % (m)) - assert m>n, "n should be larger than m (snps>inds)" - keep = [] - mprint("G (before old normalize)",genotype_matrix) - for counter in range(m): - #print("type of genotype_matrix[:,counter]:", pf(genotype_matrix[:,counter])) - #Checks if any values in column are not numbers - not_number = np.isnan(genotype_matrix[:,counter]) - - #Gets vector of values for column (no values in vector if not all values in col are numbers) - marker_values = genotype_matrix[True - not_number, counter] - #print("marker_values is:", pf(marker_values)) - - #Gets mean of values in vector - values_mean = marker_values.mean() - - genotype_matrix[not_number,counter] = values_mean - vr = genotype_matrix[:,counter].var() - if vr == 0: - continue - keep.append(counter) - genotype_matrix[:,counter] = (genotype_matrix[:,counter] - values_mean) / np.sqrt(vr) - progress('kinship_old normalize genotype',counter,m) - - genotype_matrix = genotype_matrix[:,keep] - mprint("G (after old normalize)",genotype_matrix.T) - kinship_matrix = np.dot(genotype_matrix, genotype_matrix.T) * 1.0/float(m) - return kinship_matrix,genotype_matrix - # return kinship_full(genotype_matrix.T),genotype_matrix - -def GWAS(pheno_vector, - genotype_matrix, - kinship_matrix, - kinship_eigen_vals=None, - kinship_eigen_vectors=None, - covariate_matrix=None, - restricted_max_likelihood=True, - refit=False, - temp_data=None): - """ - Performs a basic GWAS scan using the LMM. This function - uses the LMM module to assess association at each SNP and - does some simple cleanup, such as removing missing individuals - per SNP and re-computing the eigen-decomp - - pheno_vector - n x 1 phenotype vector - genotype_matrix - n x m SNP matrix - kinship_matrix - n x n kinship matrix - kinship_eigen_vals, kinship_eigen_vectors = linalg.eigh(K) - or the eigen vectors and values for K - covariate_matrix - n x q covariate matrix - restricted_max_likelihood - use restricted maximum likelihood - refit - refit the variance component for each SNP - - """ - if kinship_eigen_vals is None: - kinship_eigen_vals = [] - if kinship_eigen_vectors is None: - kinship_eigen_vectors = [] - - n = genotype_matrix.shape[0] - m = genotype_matrix.shape[1] - - if covariate_matrix == None: - covariate_matrix = np.ones((n,1)) - - # Remove missing values in pheno_vector and adjust associated parameters - v = np.isnan(pheno_vector) - if v.sum(): - keep = True - v - print(pheno_vector.shape,pheno_vector) - print(keep.shape,keep) - pheno_vector = pheno_vector[keep] - #genotype_matrix = genotype_matrix[keep,:] - #covariate_matrix = covariate_matrix[keep,:] - #kinship_matrix = kinship_matrix[keep,:][:,keep] - kinship_eigen_vals = [] - kinship_eigen_vectors = [] - - lmm_ob = LMM(pheno_vector, - kinship_matrix, - kinship_eigen_vals, - kinship_eigen_vectors, - covariate_matrix) - if not refit: - lmm_ob.fit() - - p_values = [] - t_statistics = [] - - n = genotype_matrix.shape[0] - m = genotype_matrix.shape[1] - - for counter in range(m): - x = genotype_matrix[:,counter].reshape((n, 1)) - v = np.isnan(x).reshape((-1,)) - if v.sum(): - keep = True - v - xs = x[keep,:] - if xs.var() == 0: - p_values.append(0) - t_statistics.append(np.nan) - continue - - print(genotype_matrix.shape,pheno_vector.shape,keep.shape) - - pheno_vector = pheno_vector[keep] - covariate_matrix = covariate_matrix[keep,:] - kinship_matrix = kinship_matrix[keep,:][:,keep] - lmm_ob_2 = LMM(pheno_vector, - kinship_matrix, - X0=covariate_matrix) - if refit: - lmm_ob_2.fit(X=xs) - else: - lmm_ob_2.fit() - ts, ps, beta, betaVar = lmm_ob_2.association(xs, REML=restricted_max_likelihood) - else: - if x.var() == 0: - p_values.append(0) - t_statistics.append(np.nan) - continue - - if refit: - lmm_ob.fit(X=x) - ts, ps, beta, betaVar = lmm_ob.association(x, REML=restricted_max_likelihood) - - progress("gwas_old",counter,m) - - p_values.append(ps) - t_statistics.append(ts) - - return t_statistics, p_values - - -class LMM: - - """ - This is a simple version of EMMA/fastLMM. - The main purpose of this module is to take a phenotype vector (Y), a set of covariates (X) and a kinship matrix (K) - and to optimize this model by finding the maximum-likelihood estimates for the model parameters. - There are three model parameters: heritability (h), covariate coefficients (beta) and the total - phenotypic variance (sigma). - Heritability as defined here is the proportion of the total variance (sigma) that is attributed to - the kinship matrix. - - For simplicity, we assume that everything being input is a numpy array. - If this is not the case, the module may throw an error as conversion from list to numpy array - is not done consistently. - - """ - def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=True): - - """ - The constructor takes a phenotype vector or array of size n. - It takes a kinship matrix of size n x n. Kva and Kve can be computed as Kva,Kve = linalg.eigh(K) and cached. - If they are not provided, the constructor will calculate them. - X0 is an optional covariate matrix of size n x q, where there are q covariates. - When this parameter is not provided, the constructor will set X0 to an n x 1 matrix of all ones to represent a mean effect. - """ - - if X0 is None: X0 = np.ones(len(Y)).reshape(len(Y),1) - self.verbose = verbose - - #x = Y != -9 - x = True - np.isnan(Y) - #pdb.set_trace() - if not x.sum() == len(Y): - print("Removing %d missing values from Y\n" % ((True - x).sum())) - if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum())) - Y = Y[x] - print("x: ", len(x)) - print("K: ", K.shape) - #K = K[x,:][:,x] - X0 = X0[x,:] - Kva = [] - Kve = [] - self.nonmissing = x - - print("this K is:", K.shape, pf(K)) - - if len(Kva) == 0 or len(Kve) == 0: - # if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) ) - begin = time.time() - # Kva,Kve = linalg.eigh(K) - Kva,Kve = kvakve(K) - end = time.time() - if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin)) - print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve)) - - self.K = K - self.Kva = Kva - self.Kve = Kve - print("self.Kva is: ", self.Kva.shape, pf(self.Kva)) - print("self.Kve is: ", self.Kve.shape, pf(self.Kve)) - self.Y = Y - self.X0 = X0 - self.N = self.K.shape[0] - - # ----> Below moved to kinship.kvakve(K) - # if sum(self.Kva < 1e-6): - # if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0))) - # self.Kva[self.Kva < 1e-6] = 1e-6 - - self.transform() - - def transform(self): - - """ - Computes a transformation on the phenotype vector and the covariate matrix. - The transformation is obtained by left multiplying each parameter by the transpose of the - eigenvector matrix of K (the kinship). - """ - - self.Yt = matrixMult(self.Kve.T, self.Y) - self.X0t = matrixMult(self.Kve.T, self.X0) - self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))]) - self.q = self.X0t.shape[1] - - def getMLSoln(self,h,X): - - """ - Obtains the maximum-likelihood estimates for the covariate coefficients (beta), - the total variance of the trait (sigma) and also passes intermediates that can - be utilized in other functions. The input parameter h is a value between 0 and 1 and represents - the heritability or the proportion of the total variance attributed to genetics. The X is the - covariate matrix. - """ - - S = 1.0/(h*self.Kva + (1.0 - h)) - Xt = X.T*S - XX = matrixMult(Xt,X) - XX_i = linalg.inv(XX) - beta = matrixMult(matrixMult(XX_i,Xt),self.Yt) - Yt = self.Yt - matrixMult(X,beta) - Q = np.dot(Yt.T*S,Yt) - sigma = Q * 1.0 / (float(self.N) - float(X.shape[1])) - return beta,sigma,Q,XX_i,XX - - def LL_brent(self,h,X=None,REML=False): - #brent will not be bounded by the specified bracket. - # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search. - if h < 0: return 1e6 - return -self.LL(h,X,stack=False,REML=REML)[0] - - def LL(self,h,X=None,stack=True,REML=False): - - """ - Computes the log-likelihood for a given heritability (h). If X==None, then the - default X0t will be used. If X is set and stack=True, then X0t will be matrix concatenated with - the input X. If stack is false, then X is used in place of X0t in the LL calculation. - REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True. - """ - - if X is None: - X = self.X0t - elif stack: - self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0] - X = self.X0t_stack - - n = float(self.N) - q = float(X.shape[1]) - beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X) - LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q) - LL = -0.5 * LL - - if REML: - LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(linalg.det(matrixMult(X.T,X))) - np.log(linalg.det(XX)) - LL = LL + 0.5*LL_REML_part - - return LL,beta,sigma,XX_i - - def getMax(self,H, X=None,REML=False): - - """ - Helper functions for .fit(...). - This function takes a set of LLs computed over a grid and finds possible regions - containing a maximum. Within these regions, a Brent search is performed to find the - optimum. - - """ - n = len(self.LLs) - HOpt = [] - for i in range(1,n-2): - if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]: - HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1]))) - if np.isnan(HOpt[-1][0]): - HOpt[-1][0] = [self.LLs[i-1]] - - if len(HOpt) > 1: - if self.verbose: - sys.stderr.write("NOTE: Found multiple optima. Returning first...\n") - return HOpt[0] - elif len(HOpt) == 1: - return HOpt[0] - elif self.LLs[0] > self.LLs[n-1]: - return H[0] - else: - return H[n-1] - - def fit(self,X=None,ngrids=100,REML=True): - - """ - Finds the maximum-likelihood solution for the heritability (h) given the current parameters. - X can be passed and will transformed and concatenated to X0t. Otherwise, X0t is used as - the covariate matrix. - - This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum. - Given this optimum, the function computes the LL and associated ML solutions. - """ - - if X == None: - X = self.X0t - else: - #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)]) - self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0] - X = self.X0t_stack - - H = np.array(range(ngrids)) / float(ngrids) - L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H]) - self.LLs = L - - hmax = self.getMax(H,X,REML) - L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML) - - self.H = H - self.optH = hmax - self.optLL = L - self.optBeta = beta - self.optSigma = sigma - - return hmax,beta,sigma,L - - def association(self,X, h = None, stack=True,REML=True, returnBeta=True): - - """ - Calculates association statitics for the SNPs encoded in the vector X of size n. - If h == None, the optimal h stored in optH is used. - - """ - if stack: - #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)]) - self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0] - X = self.X0t_stack - - if h == None: - h = self.optH - - L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML) - q = len(beta) - ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q) - - if returnBeta: - return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma - return ts,ps - - def tstat(self,beta,var,sigma,q): - - """ - Calculates a t-statistic and associated p-value given the estimate of beta and its standard error. - This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test. - """ - - ts = beta / np.sqrt(var * sigma) - ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q)) - if not len(ts) == 1 or not len(ps) == 1: - print("ts=",ts) - print("ps=",ps) - raise Exception("Something bad happened :(") - return ts.sum(),ps.sum() - - def plotFit(self,color='b-',title=''): - - """ - Simple function to visualize the likelihood space. It takes the LLs - calcualted over a grid and normalizes them by subtracting off the mean and exponentiating. - The resulting "probabilities" are normalized to one and plotted against heritability. - This can be seen as an approximation to the posterior distribuiton of heritability. - - For diagnostic purposes this lets you see if there is one distinct maximum or multiple - and what the variance of the parameter looks like. - """ - import matplotlib.pyplot as pl - - mx = self.LLs.max() - p = np.exp(self.LLs - mx) - p = p/p.sum() - - pl.plot(self.H,p,color) - pl.xlabel("Heritability") - pl.ylabel("Probability of data") - pl.title(title) - -def run_gwas(species,n,m,k,y,geno,cov=None,reml=True,refit=False,inputfn=None,new_code=True): - """ - Invoke pylmm using genotype as a matrix or as a (SNP) iterator. - """ - info("run_gwas") - print('pheno', y) - - if species == "human" : - print('kinship', k ) - ps, ts = run_human(pheno_vector = y, - covariate_matrix = cov, - plink_input_file = inputfn, - kinship_matrix = k, - refit = refit) - else: - print('geno', geno.shape, geno) - - if new_code: - ps, ts = run_other_new(n,m,pheno_vector = y, - geno = geno, - restricted_max_likelihood = reml, - refit = refit) - else: - ps, ts = run_other_old(pheno_vector = y, - genotype_matrix = geno, - restricted_max_likelihood = reml, - refit = refit) - return ps,ts - -def gwas_with_redis(key,species,new_code=True): - """ - Invoke pylmm using Redis as a container. new_code runs the new - version. All the Redis code goes here! - """ - json_params = Redis.get(key) - - params = json.loads(json_params) - - tempdata = temp_data.TempData(params['temp_uuid']) - def update_tempdata(loc,i,total): - """ - This is the single method that updates Redis for percentage complete! - """ - tempdata.store("percent_complete",round(i*100.0/total)) - debug("Updating REDIS percent_complete=%d" % (round(i*100.0/total))) - progress_set_func(update_tempdata) - - def narray(t): - info("Type is "+t) - v = params.get(t) - if v is not None: - # Note input values can be array of string or float - v1 = [x if x != 'NA' else 'nan' for x in v] - v = np.array(v1).astype(np.float) - return v - - def marray(t): - info("Type is "+t) - v = params.get(t) - if v is not None: - m = [] - for r in v: - # Note input values can be array of string or float - r1 = [x if x != 'NA' else 'nan' for x in r] - m.append(np.array(r1).astype(np.float)) - return np.array(m) - return np.array(v) - - def marrayT(t): - m = marray(t) - if m is not None: - return m.T - return m - - # We are transposing before we enter run_gwas - this should happen on the webserver - # side (or when reading data from file) - k = marray('kinship_matrix') - g = marrayT('genotype_matrix') - mprint("geno",g) - y = narray('pheno_vector') - n = len(y) - m = params.get('num_genotypes') - if m is None: - m = g.shape[0] - info("m=%d,n=%d" % (m,n)) - ps,ts = run_gwas(species,n,m,k,y,g,narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params.get('input_file_name'),new_code) - - results_key = "pylmm:results:" + params['temp_uuid'] - - # fatal(results_key) - json_results = json.dumps(dict(p_values = ps, - t_stats = ts)) - - #Pushing json_results into a list where it is the only item because blpop needs a list - Redis.rpush(results_key, json_results) - Redis.expire(results_key, 60*60) - return ps, ts - -def gn2_load_redis(key,species,kinship,pheno,geno,new_code=True): - """ - This function emulates current GN2 behaviour by pre-loading Redis (note the input - genotype is transposed to emulate GN2 (FIXME!) - """ - info("Loading Redis from parsed data") - if kinship == None: - k = None - else: - k = kinship.tolist() - params = dict(pheno_vector = pheno.tolist(), - genotype_matrix = geno.T.tolist(), - num_genotypes = geno.shape[0], - kinship_matrix = k, - covariate_matrix = None, - input_file_name = None, - restricted_max_likelihood = True, - refit = False, - temp_uuid = "testrun_temp_uuid", - - # meta data - timestamp = datetime.datetime.now().isoformat()) - - json_params = json.dumps(params) - Redis.set(key, json_params) - Redis.expire(key, 60*60) - - return gwas_with_redis(key,species,new_code) - -def gn2_load_redis_iter(key,species,kinship,pheno,geno_iterator): - """ - This function emulates GN2 behaviour by pre-loading Redis with - a SNP iterator, for this it sets a key for every genotype (SNP) - """ - print("Loading Redis using a SNP iterator") - for i,genotypes in enumerate(geno_iterator): - gkey = key+'_geno_'+str(i) - Redis.set(gkey, genotypes) - Redis.expire(gkey, 60*60) - - if kinship == None: - k = None - else: - k = kinship.tolist() - params = dict(pheno_vector = pheno.tolist(), - genotype_matrix = "iterator", - num_genotypes = i, - kinship_matrix = k, - covariate_matrix = None, - input_file_name = None, - restricted_max_likelihood = True, - refit = False, - temp_uuid = "testrun_temp_uuid", - - # meta data - timestamp = datetime.datetime.now().isoformat(), - ) - - json_params = json.dumps(params) - Redis.set(key, json_params) - Redis.expire(key, 60*60) - return gwas_with_redis(key,species) - -# This is the main function used by Genenetwork2 (with environment) -# -# Note that this calling route will become OBSOLETE (we should use runlmm.py -# instead) -def gn2_main(): - import argparse - parser = argparse.ArgumentParser(description='Run pyLMM') - parser.add_argument('-k', '--key') - parser.add_argument('-s', '--species') - - opts = parser.parse_args() - - key = opts.key - species = opts.species - - gwas_with_redis(key,species) - - -if __name__ == '__main__': - print("WARNING: Calling pylmm from lmm.py will become OBSOLETE, use runlmm.py instead!") - gn2_main() - diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py deleted file mode 100644 index d871d8d2..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py +++ /dev/null @@ -1,433 +0,0 @@ -# pylmm is a python-based linear mixed-model solver with applications to GWAS - -# Copyright (C) 2013,2014 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -import sys -import time -import numpy as np -from scipy.linalg import eigh, inv, det -import scipy.stats as stats # t-tests -from scipy import optimize -from optmatrix import matrixMult -import kinship - -sys.stderr.write("INFO: pylmm (lmm2) system path is "+":".join(sys.path)+"\n") -sys.stderr.write("INFO: pylmm (lmm2) file is "+__file__+"\n") - -# ---- A trick to decide on the environment: -try: - sys.stderr.write("INFO: lmm2 try loading module\n") - import utility.formatting # this is never used, just to check the environment - sys.stderr.write("INFO: This is a genenetwork2 environment (lmm2)\n") - from gn2 import uses, progress_set_func -except ImportError: - # Failed to load gn2 - has_gn2=False - import standalone as handlers - from standalone import uses, progress_set_func - sys.stderr.write("WARNING: LMM2 standalone version missing the Genenetwork2 environment\n") - -progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal') - - -def calculateKinship(W,center=False): - """ - W is an n x m matrix encoding SNP minor alleles. - - This function takes a matrix oF SNPs, imputes missing values with the maf, - normalizes the resulting vectors and returns the RRM matrix. - """ - n = W.shape[0] - m = W.shape[1] - keep = [] - for i in range(m): - mn = W[True - np.isnan(W[:,i]),i].mean() - W[np.isnan(W[:,i]),i] = mn - vr = W[:,i].var() - if vr == 0: continue - - keep.append(i) - W[:,i] = (W[:,i] - mn) / np.sqrt(vr) - - W = W[:,keep] - K = matrixMult(W,W.T) * 1.0/float(m) - if center: - P = np.diag(np.repeat(1,n)) - 1/float(n) * np.ones((n,n)) - S = np.trace(matrixMult(matrixMult(P,K),P)) - K_n = (n - 1)*K / S - return K_n - return K - -def GWAS(Y, X, K, Kva=[], Kve=[], X0=None, REML=True, refit=False): - """ - - Performs a basic GWAS scan using the LMM. This function - uses the LMM module to assess association at each SNP and - does some simple cleanup, such as removing missing individuals - per SNP and re-computing the eigen-decomp - - Y - n x 1 phenotype vector - X - n x m SNP matrix (genotype matrix) - K - n x n kinship matrix - Kva,Kve = linalg.eigh(K) - or the eigen vectors and values for K - X0 - n x q covariate matrix - REML - use restricted maximum likelihood - refit - refit the variance component for each SNP - - """ - n = X.shape[0] - m = X.shape[1] - prins("Initialize GWAS") - print("genotype matrix n is:", n) - print("genotype matrix m is:", m) - - if X0 is None: - X0 = np.ones((n,1)) - - # Remove missing values in Y and adjust associated parameters - v = np.isnan(Y) - if v.sum(): - keep = True - v - keep = keep.reshape((-1,)) - Y = Y[keep] - X = X[keep,:] - X0 = X0[keep,:] - K = K[keep,:][:,keep] - Kva = [] - Kve = [] - - if len(Y) == 0: - return np.ones(m)*np.nan,np.ones(m)*np.nan - - L = LMM(Y,K,Kva,Kve,X0) - if not refit: L.fit() - - PS = [] - TS = [] - - n = X.shape[0] - m = X.shape[1] - - for i in range(m): - x = X[:,i].reshape((n,1)) - v = np.isnan(x).reshape((-1,)) - if v.sum(): - keep = True - v - xs = x[keep,:] - if xs.var() == 0: - PS.append(np.nan) - TS.append(np.nan) - continue - - Ys = Y[keep] - X0s = X0[keep,:] - Ks = K[keep,:][:,keep] - Ls = LMM(Ys,Ks,X0=X0s) - if refit: - Ls.fit(X=xs) - else: - Ls.fit() - ts,ps = Ls.association(xs,REML=REML) - else: - if x.var() == 0: - PS.append(np.nan) - TS.append(np.nan) - continue - - if refit: - L.fit(X=x) - ts,ps = L.association(x,REML=REML) - - PS.append(ps) - TS.append(ts) - - return TS,PS - -class LMM2: - - """This is a simple version of EMMA/fastLMM. - - The main purpose of this module is to take a phenotype vector (Y), - a set of covariates (X) and a kinship matrix (K) and to optimize - this model by finding the maximum-likelihood estimates for the - model parameters. There are three model parameters: heritability - (h), covariate coefficients (beta) and the total phenotypic - variance (sigma). Heritability as defined here is the proportion - of the total variance (sigma) that is attributed to the kinship - matrix. - - For simplicity, we assume that everything being input is a numpy - array. If this is not the case, the module may throw an error as - conversion from list to numpy array is not done consistently. - - """ - def __init__(self,Y,K,Kva=[],Kve=[],X0=None,verbose=False): - - """The constructor takes a phenotype vector or array Y of size n. It - takes a kinship matrix K of size n x n. Kva and Kve can be - computed as Kva,Kve = linalg.eigh(K) and cached. If they are - not provided, the constructor will calculate them. X0 is an - optional covariate matrix of size n x q, where there are q - covariates. When this parameter is not provided, the - constructor will set X0 to an n x 1 matrix of all ones to - represent a mean effect. - """ - - if X0 is None: - X0 = np.ones(len(Y)).reshape(len(Y),1) - self.verbose = verbose - - x = True - np.isnan(Y) - x = x.reshape(-1,) - if not x.sum() == len(Y): - if self.verbose: sys.stderr.write("Removing %d missing values from Y\n" % ((True - x).sum())) - Y = Y[x] - K = K[x,:][:,x] - X0 = X0[x,:] - Kva = [] - Kve = [] - self.nonmissing = x - - print("this K is:", K.shape, K) - - if len(Kva) == 0 or len(Kve) == 0: - # if self.verbose: sys.stderr.write("Obtaining eigendecomposition for %dx%d matrix\n" % (K.shape[0],K.shape[1]) ) - begin = time.time() - # Kva,Kve = linalg.eigh(K) - Kva,Kve = kinship.kvakve(K) - end = time.time() - if self.verbose: sys.stderr.write("Total time: %0.3f\n" % (end - begin)) - print("sum(Kva),sum(Kve)=",sum(Kva),sum(Kve)) - - self.K = K - self.Kva = Kva - self.Kve = Kve - self.N = self.K.shape[0] - self.Y = Y.reshape((self.N,1)) - self.X0 = X0 - - if sum(self.Kva < 1e-6): - if self.verbose: sys.stderr.write("Cleaning %d eigen values\n" % (sum(self.Kva < 0))) - self.Kva[self.Kva < 1e-6] = 1e-6 - - self.transform() - - def transform(self): - - """ - Computes a transformation on the phenotype vector and the covariate matrix. - The transformation is obtained by left multiplying each parameter by the transpose of the - eigenvector matrix of K (the kinship). - """ - - self.Yt = matrixMult(self.Kve.T, self.Y) - self.X0t = matrixMult(self.Kve.T, self.X0) - self.X0t_stack = np.hstack([self.X0t, np.ones((self.N,1))]) - self.q = self.X0t.shape[1] - - def getMLSoln(self,h,X): - - """ - Obtains the maximum-likelihood estimates for the covariate coefficients (beta), - the total variance of the trait (sigma) and also passes intermediates that can - be utilized in other functions. The input parameter h is a value between 0 and 1 and represents - the heritability or the proportion of the total variance attributed to genetics. The X is the - covariate matrix. - """ - - S = 1.0/(h*self.Kva + (1.0 - h)) - Xt = X.T*S - XX = matrixMult(Xt,X) - XX_i = inv(XX) - beta = matrixMult(matrixMult(XX_i,Xt),self.Yt) - Yt = self.Yt - matrixMult(X,beta) - Q = np.dot(Yt.T*S,Yt) - sigma = Q * 1.0 / (float(self.N) - float(X.shape[1])) - return beta,sigma,Q,XX_i,XX - - def LL_brent(self,h,X=None,REML=False): - #brent will not be bounded by the specified bracket. - # I return a large number if we encounter h < 0 to avoid errors in LL computation during the search. - if h < 0: return 1e6 - return -self.LL(h,X,stack=False,REML=REML)[0] - - def LL(self,h,X=None,stack=True,REML=False): - - """ - Computes the log-likelihood for a given heritability (h). If X==None, then the - default X0t will be used. If X is set and stack=True, then X0t will be matrix concatenated with - the input X. If stack is false, then X is used in place of X0t in the LL calculation. - REML is computed by adding additional terms to the standard LL and can be computed by setting REML=True. - """ - - if X is None: X = self.X0t - elif stack: - self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0] - X = self.X0t_stack - - n = float(self.N) - q = float(X.shape[1]) - beta,sigma,Q,XX_i,XX = self.getMLSoln(h,X) - LL = n*np.log(2*np.pi) + np.log(h*self.Kva + (1.0-h)).sum() + n + n*np.log(1.0/n * Q) - LL = -0.5 * LL - - if REML: - LL_REML_part = q*np.log(2.0*np.pi*sigma) + np.log(det(matrixMult(X.T,X))) - np.log(det(XX)) - LL = LL + 0.5*LL_REML_part - - - LL = LL.sum() - return LL,beta,sigma,XX_i - - def getMax(self,H, X=None,REML=False): - - """ - Helper functions for .fit(...). - This function takes a set of LLs computed over a grid and finds possible regions - containing a maximum. Within these regions, a Brent search is performed to find the - optimum. - - """ - n = len(self.LLs) - HOpt = [] - for i in range(1,n-2): - if self.LLs[i-1] < self.LLs[i] and self.LLs[i] > self.LLs[i+1]: - HOpt.append(optimize.brent(self.LL_brent,args=(X,REML),brack=(H[i-1],H[i+1]))) - if np.isnan(HOpt[-1]): HOpt[-1] = H[i-1] - #if np.isnan(HOpt[-1]): HOpt[-1] = self.LLs[i-1] - #if np.isnan(HOpt[-1][0]): HOpt[-1][0] = [self.LLs[i-1]] - - if len(HOpt) > 1: - if self.verbose: sys.stderr.write("NOTE: Found multiple optima. Returning first...\n") - return HOpt[0] - elif len(HOpt) == 1: return HOpt[0] - elif self.LLs[0] > self.LLs[n-1]: return H[0] - else: return H[n-1] - - - def fit(self,X=None,ngrids=100,REML=True): - - """ - Finds the maximum-likelihood solution for the heritability (h) given the current parameters. - X can be passed and will transformed and concatenated to X0t. Otherwise, X0t is used as - the covariate matrix. - - This function calculates the LLs over a grid and then uses .getMax(...) to find the optimum. - Given this optimum, the function computes the LL and associated ML solutions. - """ - - if X is None: X = self.X0t - else: - #X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)]) - self.X0t_stack[:,(self.q)] = matrixMult(self.Kve.T,X)[:,0] - X = self.X0t_stack - - H = np.array(range(ngrids)) / float(ngrids) - L = np.array([self.LL(h,X,stack=False,REML=REML)[0] for h in H]) - self.LLs = L - - hmax = self.getMax(H,X,REML) - L,beta,sigma,betaSTDERR = self.LL(hmax,X,stack=False,REML=REML) - - self.H = H - self.optH = hmax.sum() - self.optLL = L - self.optBeta = beta - self.optSigma = sigma.sum() - - return hmax,beta,sigma,L - - def association(self,X,h=None,stack=True,REML=True,returnBeta=False): - """ - Calculates association statitics for the SNPs encoded in the vector X of size n. - If h is None, the optimal h stored in optH is used. - - """ - if False: - print "X=",X - print "h=",h - print "q=",self.q - print "self.Kve=",self.Kve - print "X0t_stack=",self.X0t_stack.shape,self.X0t_stack - - if stack: - # X = np.hstack([self.X0t,matrixMult(self.Kve.T, X)]) - m = matrixMult(self.Kve.T,X) - # print "m=",m - m = m[:,0] - self.X0t_stack[:,(self.q)] = m - X = self.X0t_stack - - if h is None: h = self.optH - - L,beta,sigma,betaVAR = self.LL(h,X,stack=False,REML=REML) - q = len(beta) - ts,ps = self.tstat(beta[q-1],betaVAR[q-1,q-1],sigma,q) - - if returnBeta: return ts,ps,beta[q-1].sum(),betaVAR[q-1,q-1].sum()*sigma - return ts,ps - - def tstat(self,beta,var,sigma,q,log=False): - - """ - Calculates a t-statistic and associated p-value given the estimate of beta and its standard error. - This is actually an F-test, but when only one hypothesis is being performed, it reduces to a t-test. - """ - - ts = beta / np.sqrt(var * sigma) - #ps = 2.0*(1.0 - stats.t.cdf(np.abs(ts), self.N-q)) - # sf == survival function - this is more accurate -- could also use logsf if the precision is not good enough - if log: - ps = 2.0 + (stats.t.logsf(np.abs(ts), self.N-q)) - else: - ps = 2.0*(stats.t.sf(np.abs(ts), self.N-q)) - if not len(ts) == 1 or not len(ps) == 1: - raise Exception("Something bad happened :(") - return ts.sum(),ps.sum() - - def plotFit(self,color='b-',title=''): - - """ - Simple function to visualize the likelihood space. It takes the LLs - calcualted over a grid and normalizes them by subtracting off the mean and exponentiating. - The resulting "probabilities" are normalized to one and plotted against heritability. - This can be seen as an approximation to the posterior distribuiton of heritability. - - For diagnostic purposes this lets you see if there is one distinct maximum or multiple - and what the variance of the parameter looks like. - """ - import matplotlib.pyplot as pl - - mx = self.LLs.max() - p = np.exp(self.LLs - mx) - p = p/p.sum() - - pl.plot(self.H,p,color) - pl.xlabel("Heritability") - pl.ylabel("Probability of data") - pl.title(title) - - def meanAndVar(self): - - mx = self.LLs.max() - p = np.exp(self.LLs - mx) - p = p/p.sum() - - mn = (self.H * p).sum() - vx = ((self.H - mn)**2 * p).sum() - - return mn,vx - diff --git a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py b/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py deleted file mode 100644 index 5c71db6a..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/optmatrix.py +++ /dev/null @@ -1,55 +0,0 @@ -import sys -import time -import numpy as np -from numpy.distutils.system_info import get_info; -from scipy import linalg -from scipy import optimize -from scipy import stats - -useNumpy = None -hasBLAS = None - -def matrix_initialize(useBLAS=True): - global useNumpy # module based variable - if useBLAS and useNumpy == None: - print get_info('blas_opt') - try: - linalg.fblas - sys.stderr.write("INFO: using linalg.fblas\n") - useNumpy = False - hasBLAS = True - except AttributeError: - sys.stderr.write("WARNING: linalg.fblas not found, using numpy.dot instead!\n") - useNumpy = True - else: - sys.stderr.write("INFO: using numpy.dot\n") - useNumpy=True - -def matrixMult(A,B): - global useNumpy # module based variable - - if useNumpy: - return np.dot(A,B) - - # If the matrices are in Fortran order then the computations will be faster - # when using dgemm. Otherwise, the function will copy the matrix and that takes time. - if not A.flags['F_CONTIGUOUS']: - AA = A.T - transA = True - else: - AA = A - transA = False - - if not B.flags['F_CONTIGUOUS']: - BB = B.T - transB = True - else: - BB = B - transB = False - - return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB) - -def matrixMultT(M): - # res = np.dot(W,W.T) - # return linalg.fblas.dgemm(alpha=1.,a=M.T,b=M.T,trans_a=True,trans_b=False) - return matrixMult(M,M.T) diff --git a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py b/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py deleted file mode 100644 index 7b652515..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/phenotype.py +++ /dev/null @@ -1,65 +0,0 @@ -# Phenotype routines - -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -import sys -import numpy as np - -# ---- A trick to decide on the environment: -try: - from wqflask.my_pylmm.pyLMM import chunks - from gn2 import uses, progress_set_func -except ImportError: - has_gn2=False - import standalone as handlers - from standalone import uses, progress_set_func - -progress,debug,info,mprint = uses('progress','debug','info','mprint') - -def remove_missing(n,y,g): - """ - Remove missing data from matrices, make sure the genotype data has - individuals as rows - """ - assert(y is not None) - assert y.shape[0] == g.shape[0],"y (n) %d, g (n,m) %s" % (y.shape[0],g.shape) - - y1 = y - g1 = g - v = np.isnan(y) - keep = True - v - if v.sum(): - info("runlmm.py: Cleaning the phenotype vector and genotype matrix by removing %d individuals...\n" % (v.sum())) - y1 = y[keep] - g1 = g[keep,:] - n = y1.shape[0] - return n,y1,g1,keep - -def remove_missing_new(n,y): - """ - Remove missing data. Returns new n,y,keep - """ - assert(y is not None) - y1 = y - v = np.isnan(y) - keep = True - v - if v.sum(): - info("runlmm.py: Cleaning the phenotype vector by removing %d individuals" % (v.sum())) - y1 = y[keep] - n = y1.shape[0] - return n,y1,keep - diff --git a/wqflask/wqflask/my_pylmm/pyLMM/plink.py b/wqflask/wqflask/my_pylmm/pyLMM/plink.py deleted file mode 100644 index 7bd2df91..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/plink.py +++ /dev/null @@ -1,107 +0,0 @@ -# Plink module -# -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# Some of the BED file parsing came from pylmm: -# Copyright (C) 2013 Nicholas A. Furlotte (nick.furlotte@gmail.com) - -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -# According to the PLINK information - -# Parse a textual BIM file and return the contents as a list of tuples -# -# Extended variant information file accompanying a .bed binary genotype table. -# -# A text file with no header line, and one line per variant with the following six fields: -# -# Chromosome code (either an integer, or 'X'/'Y'/'XY'/'MT'; '0' indicates unknown) or name -# Variant identifier -# Position in morgans or centimorgans (safe to use dummy value of '0') -# Base-pair coordinate (normally 1-based, but 0 ok; limited to 231-2) -# Allele 1 (corresponding to clear bits in .bed; usually minor) -# Allele 2 (corresponding to set bits in .bed; usually major) -# -# Allele codes can contain more than one character. Variants with negative bp coordinates are ignored by PLINK. Example -# -# 1 mm37-1-3125499 0 3125499 1 2 -# 1 mm37-1-3125701 0 3125701 1 2 -# 1 mm37-1-3187481 0 3187481 1 2 - -import struct -# import numpy as np - -def readbim(fn): - res = [] - for line in open(fn): - list = line.split() - if len([True for e in list if e == 'nan']) == 0: - res.append( (list[0],list[1],int(list[2]),int(list[3]),int(list[4]),int(list[5])) ) - else: - res.append( (list[0],list[1],list[2],float('nan'),float('nan'),float('nan')) ) - return res - -# .bed (PLINK binary biallelic genotype table) -# -# Primary representation of genotype calls at biallelic variants. Must -# be accompanied by .bim and .fam files. Basically contains num SNP -# blocks containing IND (compressed 4 IND into a byte) -# -# Since it is a biallelic format it supports for every individual -# whether the first allele is homozygous (b00), the second allele is -# homozygous (b11), it is heterozygous (b10) or that it is missing -# (b01). - -# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bed -# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#fam -# http://pngu.mgh.harvard.edu/~purcell/plink2/formats.html#bim - -def readbed(fn,inds,encoding,func=None): - - # For every SNP block fetch the individual genotypes using values - # 0.0 and 1.0 for homozygous and 0.5 for heterozygous alleles - def fetchGenotypes(X): - # D = { \ - # '00': 0.0, \ - # '10': 0.5, \ - # '11': 1.0, \ - # '01': float('nan') \ - # } - - Didx = { '00': 0, '10': 1, '11': 2, '01': 3 } - G = [] - for x in X: - if not len(x) == 10: - xx = x[2:] - x = '0b' + '0'*(8 - len(xx)) + xx - a,b,c,d = (x[8:],x[6:8],x[4:6],x[2:4]) - L = [encoding[Didx[y]] for y in [a,b,c,d]] - G += L - G = G[:inds] - # G = np.array(G) - return G - - bytes = inds / 4 + (inds % 4 and 1 or 0) - format = 'c'*bytes - count = 0 - with open(fn,'rb') as f: - magic = f.read(3) - assert( ":".join("{:02x}".format(ord(c)) for c in magic) == "6c:1b:01") - while True: - count += 1 - X = f.read(bytes) - if not X: - return(count-1) - XX = [bin(ord(x)) for x in struct.unpack(format,X)] - xs = fetchGenotypes(XX) - func(count,xs) diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py deleted file mode 100644 index 6b241cd6..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py +++ /dev/null @@ -1,229 +0,0 @@ -# This is the LMM runner that calls the possible methods using command line -# switches. It acts as a multiplexer where all the invocation complexity -# is kept outside the main LMM routines. -# -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -from optparse import OptionParser -import sys -import tsvreader -import numpy as np - -# Add local dir to PYTHONPATH -import os -cwd = os.path.dirname(__file__) -if sys.path[0] != cwd: - sys.path.insert(1,cwd) - -# pylmm modules -from lmm import gn2_load_redis, gn2_load_redis_iter, calculate_kinship_new, run_gwas -from kinship import kinship, kinship_full -import genotype -import phenotype -from standalone import uses - -progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal') - -usage = """ -python runlmm.py [options] command - - runlmm.py processing multiplexer reads standardised input formats - and calls the different routines (writes to stdout) - - Current commands are: - - parse : only parse input files - redis : use Redis to call into GN2 - kinship : calculate (new) kinship matrix - - try --help for more information -""" - - -parser = OptionParser(usage=usage) -# parser.add_option("-f", "--file", dest="input file", -# help="In", metavar="FILE") -parser.add_option("--kinship",dest="kinship", - help="Kinship file format 1.0") -parser.add_option("--pheno",dest="pheno", - help="Phenotype file format 1.0") -parser.add_option("--geno",dest="geno", - help="Genotype file format 1.0") -parser.add_option("--maf-normalization", - action="store_true", dest="maf_normalization", default=False, - help="Apply MAF genotype normalization") -parser.add_option("--genotype-normalization", - action="store_true", dest="genotype_normalization", default=False, - help="Force genotype normalization") -parser.add_option("--remove-missing-phenotypes", - action="store_true", dest="remove_missing_phenotypes", default=False, - help="Remove missing phenotypes") -parser.add_option("-q", "--quiet", - action="store_false", dest="verbose", default=True, - help="don't print status messages to stdout") -parser.add_option("--blas", action="store_true", default=False, dest="useBLAS", help="Use BLAS instead of numpy matrix multiplication") -parser.add_option("-t", "--threads", - type="int", dest="numThreads", - help="Threads to use") -parser.add_option("--saveKvaKve", - action="store_true", dest="saveKvaKve", default=False, - help="Testing mode") -parser.add_option("--test", - action="store_true", dest="testing", default=False, - help="Testing mode") -parser.add_option("--test-kinship", - action="store_true", dest="test_kinship", default=False, - help="Testing mode for Kinship calculation") - -(options, args) = parser.parse_args() - -if len(args) != 1: - print usage - sys.exit(1) - -cmd = args[0] -print "Command: ",cmd - -k = None -y = None -g = None - -if options.kinship: - k = tsvreader.kinship(options.kinship) - print k.shape - -if options.pheno: - y = tsvreader.pheno(options.pheno) - print y.shape - -if options.geno and cmd != 'iterator': - g = tsvreader.geno(options.geno) - print g.shape - -def check_results(ps,ts): - print np.array(ps) - print len(ps),sum(ps) - p1 = round(ps[0],4) - p2 = round(ps[-1],4) - if options.geno == 'data/small.geno': - info("Validating results for "+options.geno) - assert p1==0.7387, "p1=%f" % p1 - assert p2==0.7387, "p2=%f" % p2 - if options.geno == 'data/small_na.geno': - info("Validating results for "+options.geno) - assert p1==0.062, "p1=%f" % p1 - assert p2==0.062, "p2=%f" % p2 - if options.geno == 'data/test8000.geno': - info("Validating results for "+options.geno) - assert round(sum(ps)) == 4070 - assert len(ps) == 8000 - info("Run completed") - -if y is not None: - n = y.shape[0] - -if cmd == 'run': - if options.remove_missing_phenotypes: - raise Exception('Can not use --remove-missing-phenotypes with LMM2') - n = len(y) - m = g.shape[1] - ps, ts = run_gwas('other',n,m,k,y,g) # <--- pass in geno by SNP - check_results(ps,ts) -elif cmd == 'iterator': - if options.remove_missing_phenotypes: - raise Exception('Can not use --remove-missing-phenotypes with LMM2') - geno_iterator = tsvreader.geno_iter(options.geno) - ps, ts = gn2_load_redis_iter('testrun_iter','other',k,y,geno_iterator) - check_results(ps,ts) -elif cmd == 'redis_new': - # The main difference between redis_new and redis is that missing - # phenotypes are handled by the first - if options.remove_missing_phenotypes: - raise Exception('Can not use --remove-missing-phenotypes with LMM2') - Y = y - G = g - print "Original G",G.shape, "\n", G - # gt = G.T - # G = None - ps, ts = gn2_load_redis('testrun','other',k,Y,G,new_code=True) - check_results(ps,ts) -elif cmd == 'redis': - # Emulating the redis setup of GN2 - G = g - print "Original G",G.shape, "\n", G - if y is not None and options.remove_missing_phenotypes: - gnt = np.array(g).T - n,Y,g,keep = phenotype.remove_missing(n,y,gnt) - G = g.T - print "Removed missing phenotypes",G.shape, "\n", G - else: - Y = y - if options.maf_normalization: - G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g ) - print "MAF replacements: \n",G - if options.genotype_normalization: - G = np.apply_along_axis( genotype.normalize, axis=1, arr=G) - g = None - gnt = None - - # gt = G.T - # G = None - ps, ts = gn2_load_redis('testrun','other',k,Y,G, new_code=False) - check_results(ps,ts) -elif cmd == 'kinship': - G = g - print "Original G",G.shape, "\n", G - if y != None and options.remove_missing_phenotypes: - gnt = np.array(g).T - n,Y,g,keep = phenotype.remove_missing(n,y,g.T) - G = g.T - print "Removed missing phenotypes",G.shape, "\n", G - if options.maf_normalization: - G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g ) - print "MAF replacements: \n",G - if options.genotype_normalization: - G = np.apply_along_axis( genotype.normalize, axis=1, arr=G) - g = None - gnt = None - - if options.test_kinship: - K = kinship_full(np.copy(G)) - print "Genotype",G.shape, "\n", G - print "first Kinship method",K.shape,"\n",K - k1 = round(K[0][0],4) - K2,G = calculate_kinship_new(np.copy(G)) - print "Genotype",G.shape, "\n", G - print "GN2 Kinship method",K2.shape,"\n",K2 - k2 = round(K2[0][0],4) - - print "Genotype",G.shape, "\n", G - K3 = kinship(G) - print "third Kinship method",K3.shape,"\n",K3 - sys.stderr.write(options.geno+"\n") - k3 = round(K3[0][0],4) - if options.geno == 'data/small.geno': - assert k1==0.8333, "k1=%f" % k1 - assert k2==0.9375, "k2=%f" % k2 - assert k3==0.9375, "k3=%f" % k3 - if options.geno == 'data/small_na.geno': - assert k1==0.8333, "k1=%f" % k1 - assert k2==0.7172, "k2=%f" % k2 - assert k3==0.7172, "k3=%f" % k3 - if options.geno == 'data/test8000.geno': - assert k3==1.4352, "k3=%f" % k3 - -else: - fatal("Doing nothing") diff --git a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py b/wqflask/wqflask/my_pylmm/pyLMM/standalone.py deleted file mode 100644 index 40b2021d..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/standalone.py +++ /dev/null @@ -1,110 +0,0 @@ -# Standalone specific methods and callback handler -# -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# Set the log level with -# -# logging.basicConfig(level=logging.DEBUG) - -from __future__ import absolute_import, print_function, division - -import numpy as np -import sys -import logging - -# logger = logging.getLogger(__name__) -logger = logging.getLogger('lmm2') -logging.basicConfig(level=logging.DEBUG) -np.set_printoptions(precision=3,suppress=True) - -progress_location = None -progress_current = None -progress_prev_perc = None - -def progress_default_func(location,count,total): - global progress_current - value = round(count*100.0/total) - progress_current = value - -progress_func = progress_default_func - -def progress_set_func(func): - global progress_func - progress_func = func - -def progress(location, count, total): - global progress_location - global progress_prev_perc - - perc = round(count*100.0/total) - if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5): - progress_func(location, count, total) - logger.info("Progress: %s %d%%" % (location,perc)) - progress_location = location - progress_prev_perc = perc - -def mprint(msg,data): - """ - Array/matrix print function - """ - m = np.array(data) - if m.ndim == 1: - print(msg,m.shape,"=\n",m[0:3]," ... ",m[-3:]) - if m.ndim == 2: - print(msg,m.shape,"=\n[", - m[0][0:3]," ... ",m[0][-3:],"\n ", - m[1][0:3]," ... ",m[1][-3:],"\n ...\n ", - m[-2][0:3]," ... ",m[-2][-3:],"\n ", - m[-1][0:3]," ... ",m[-1][-3:],"]") - -def fatal(msg): - logger.critical(msg) - raise Exception(msg) - -def callbacks(): - return dict( - write = sys.stdout.write, - writeln = print, - debug = logger.debug, - info = logger.info, - warning = logger.warning, - error = logger.error, - critical = logger.critical, - fatal = fatal, - progress = progress, - mprint = mprint - ) - -def uses(*funcs): - """ - Some sugar - """ - return [callbacks()[func] for func in funcs] - -# ----- Minor test cases: - -if __name__ == '__main__': - # logging.basicConfig(level=logging.DEBUG) - logging.debug("Test %i" % (1)) - d = callbacks()['debug'] - d("TEST") - wrln = callbacks()['writeln'] - wrln("Hello %i" % 34) - progress = callbacks()['progress'] - progress("I am half way",50,100) - list = [0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15, - 0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15] - mprint("list",list) - matrix = [[1,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [2,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [3,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [4,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [5,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15], - [6,0.5,0.6096595 , -0.31559815, -0.52793285, 1.16573418e-15]] - mprint("matrix",matrix) - ix,dx = uses("info","debug") - ix("ix") - dx("dx") diff --git a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py b/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py deleted file mode 100644 index 004d45c6..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/temp_data.py +++ /dev/null @@ -1,25 +0,0 @@ -from __future__ import print_function, division, absolute_import -from redis import Redis - -import simplejson as json - -class TempData(object): - - def __init__(self, temp_uuid): - self.temp_uuid = temp_uuid - self.redis = Redis() - self.key = "tempdata:{}".format(self.temp_uuid) - - def store(self, field, value): - self.redis.hset(self.key, field, value) - self.redis.expire(self.key, 60*15) # Expire in 15 minutes - - def get_all(self): - return self.redis.hgetall(self.key) - - -if __name__ == "__main__": - redis = Redis() - for key in redis.keys(): - for field in redis.hkeys(key): - print("{}.{}={}".format(key, field, redis.hget(key, field))) diff --git a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py b/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py deleted file mode 100644 index 66b34ee2..00000000 --- a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py +++ /dev/null @@ -1,122 +0,0 @@ -# Standard file readers -# -# Copyright (C) 2015 Pjotr Prins (pjotr.prins@thebird.nl) -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU Affero General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Affero General Public License for more details. - -# You should have received a copy of the GNU Affero General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -import sys -import os -import numpy as np -import csv - -def kinship(fn): - K1 = [] - print fn - with open(fn,'r') as tsvin: - assert(tsvin.readline().strip() == "# Kinship format version 1.0") - tsvin.readline() - tsvin.readline() - tsv = csv.reader(tsvin, delimiter='\t') - for row in tsv: - ns = np.genfromtxt(row[1:]) - K1.append(ns) # <--- slow - K = np.array(K1) - return K - -def pheno(fn): - Y1 = [] - print fn - with open(fn,'r') as tsvin: - assert(tsvin.readline().strip() == "# Phenotype format version 1.0") - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsv = csv.reader(tsvin, delimiter='\t') - for row in tsv: - ns = np.genfromtxt(row[1:]) - Y1.append(ns) # <--- slow - Y = np.array(Y1) - return Y - -def geno(fn): - G1 = [] - hab_mapper = {'A':0,'H':1,'B':2,'-':3} - pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ] - - print fn - with open(fn,'r') as tsvin: - line = tsvin.readline().strip() - assert line == "# Genotype format version 1.0", line - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsv = csv.reader(tsvin, delimiter='\t') - for row in tsv: - # print(row) - id = row[0] - gs = list(row[1]) - # print id,gs - gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs] - # print id,gs2 - # ns = np.genfromtxt(row[1:]) - G1.append(gs2) # <--- slow - G = np.array(G1) - return G - -def geno(fn): - G1 = [] - for id,values in geno_iter(fn): - G1.append(values) # <--- slow - G = np.array(G1) - return G - -def geno_callback(fn,func): - hab_mapper = {'A':0,'H':1,'B':2,'-':3} - pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ] - - print fn - with open(fn,'r') as tsvin: - assert(tsvin.readline().strip() == "# Genotype format version 1.0") - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsv = csv.reader(tsvin, delimiter='\t') - for row in tsv: - id = row[0] - gs = list(row[1]) - gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs] - func(id,gs2) - -def geno_iter(fn): - """ - Yield a tuple of snpid and values - """ - hab_mapper = {'A':0,'H':1,'B':2,'-':3} - pylmm_mapper = [ 0.0, 0.5, 1.0, float('nan') ] - - print fn - with open(fn,'r') as tsvin: - assert(tsvin.readline().strip() == "# Genotype format version 1.0") - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsvin.readline() - tsv = csv.reader(tsvin, delimiter='\t') - for row in tsv: - id = row[0] - gs = list(row[1]) - gs2 = [pylmm_mapper[hab_mapper[g]] for g in gs] - yield (id,gs2) |