Intro_Biocom_ND_319_Tutorial9/Tutorial 9 at master · amidlige/Intro_Biocom_ND_319_Tutorial9 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#Import Packages
import pandas
import numpy
from scipy.optimize import minimize
from scipy.stats import norm
import re
import os
from plotnine import *

os.chdir("/Users/madelinebuynak/Desktop")

#load data
data = pandas.read_csv("ponzr1.csv",header=0,sep='\t')

subset=data.loc[data.mutation.isin(['WT','I231N']),:]

#two likelihood functions
def nllike(p,obs):
    B0=p(0)
    B1=p(1)
    sigma=p(2)

def null (p,obs):
    B0=p[0]
    sigma=p[1]

    expected=B0
    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
    return nll

def alter (p,obs):
    B0=p[0]
    B1=p[1]
    sigma=p[2]

expected=B0+B1*obs.x
nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
return nll

initialGuess=numpy.array([1,1,1])
fitNull=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
fitAlter=minimize(nllike,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)


subset=data.loc[data.mutation.sin(['WT','V456D']),:]


def null (p,obs):
    B0=p[0]
    sigma=p[1]

    expected=B0
    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
    return nll

def alter (p,obs):
    B0=p[0]
    B1=p[1]
    sigma=p[2]

    expected=B0+B1*obs.x
    nll=-1*norm(expected,sigma).logpdf(obs.y).sum()
    return nll

initialGuess=numpy.array([1,1,1])
fitNull=minimize(null,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)
fitAlter=minimize(alter,initialGuess,method="Nelder-Mead",options={'disp':True},args=subset)

D=2*(fitNull.fun-fitAlter.fun)
1-ch2.cdf(x=D,df=1)

#M124K: p-value ~ 0.72 (no effect of treatment)

#V456D: p-value ~ 5.6e-6 (effect of treatment)

#I213N: p-value ~ 0.88 (no effect of treatment)