CnrLwlss

# Plot comparing ranks

import string
import random
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

def experiment(genes,vmin=0.0,vmax=1.0):
    '''Simulate an experiment with (average) measurement for a set of genotypes'''

CnrLwlss

import matplotlib as mp
import matplotlib.pyplot as plt
import numpy as np
import string

# Set font for figure according to journal specs
mp.rcParams['font.family'] = "Arial"

# Gene deletion names read from file might not necessarily be in lower case
gnames=["EXO1","RAD9","TMA20"]

CnrLwlss

#http://matplotlib.org/examples/color/colormaps_reference.html
#http://matplotlib.org/users/colormaps.html

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

N = 200
x = np.random.randn(N)
y = np.random.randn(N)

CnrLwlss

# Hierarchical uniform distributions
Nsamps=500000
r_min=1
r_max=10
trunc=FALSE

makePlot=function(){
  xlim=c(-5,15)
  op=par(mfrow=c(1,3))
  hist(r,breaks=100,freq=FALSE,xlim=xlim)

CnrLwlss

def hierarchy_inf(data,par,iter=250000,burn=1000,thin=100):
    priors={}
    x0=mc.Uniform('x0',par.x0_min,par.x0_max)
    tau=mc.Uniform('tau',par.tau_min,par.tau_max)

    priors["x0"]=x0
    priors["tau"]=tau
  
    r=mc.Uniform('r',par.r_min,par.r_max)
    r_delta=mc.Uniform('r_delta',0,(par.r_max-par.r_min)/2)

CnrLwlss

double MCMC_base_truncate_low(double truncate, struct_data *D,struct_para *D_para,struct_priors *D_priors,double *accept,double *h,double para,double (*foo)(struct struct_data *D,struct struct_para *D_para,struct struct_priors *D_priors,double,int,int),int l, int m){
  double logu,logaprob,can;
   //can=rnorm(para,*h); /*can=para+gsl_ran_gaussian(RNG,*h);*/
   // if(can<=(truncate)){
   //   can=para;
   // }

   can=truncate-1.0;
   while(can<=truncate){
	 can=rnorm(para,*h);  

CnrLwlss

# Problem with bedpe file format as generated by bedtools
# https://code.google.com/p/bedtools/issues/detail?id=152
#SRR364781 -> wt_sample
#SRR364782 -> wt_replicate
#SRR364783 -> pol32_sample
#SRR364784 -> pol32_replicate

#source("https://bioconductor.org/biocLite.R")
#biocLite("ShortRead")
#biocLite("org.Sc.sgd.db")

CnrLwlss

# Data from Smith & Whitehouse (2012) http://dx.doi.org/10.1038/nature10895
#SRR364781 -> wt_sample
#SRR364782 -> wt_replicate
#SRR364783 -> pol32_sample
#SRR364784 -> pol32_replicate

roots=(
SRR364781
SRR364782
SRR364783

CnrLwlss

count_DDY835=c(10,13,9,7,0,15,7,0,7,0,0,10,9,5,0,10,12,6,9,8,0,10,10,0,11,0,0,9,0,0,0,0)
count_DDY836=c(15,11,8,17,0,10,7,0,7,0,0,12,8,8,0,9,12,13,10,12,0,13,10,0,8,0,0,10,0,0,0,0)
gtype=c("wt","nmd2D","exo1D","rad24D","cdc13D","nmd2_exo1D","nmd2D_rad24D","nmd2D_cdc13D",
"exo1D_rad24D","exo1D_cdc13D","rad24D_cdc13D","nmd2D_exo1D_rad24D","nmd2D_exo1D_cdc13D","nmd2D_rad24D_cdc13D",
"exo1D_rad24D_cdc13D","nmd2D_exo1D_rad24D_cdc13D","rif1_wt","rif1_nmd2D","rif1_exo1D","rif1_rad24D","rif1_cdc13D",
"rif1_nmd2_exo1D","rif1_nmd2D_rad24D","rif1_nmd2D_cdc13D","rif1_exo1D_rad24D","rif1_exo1D_cdc13D","rif1_rad24D_cdc13D","rif1_nmd2D_exo1D_rad24D",
"rif1_nmd2D_exo1D_cdc13D","rif1_nmd2D_rad24D_cdc13D","rif1_exo1D_rad24D_cdc13D","rif1_nmd2D_exo1D_rad24D_cdc13D")

df=data.frame(gene=gtype,count_DDY835=count_DDY835,count_DDY836=count_DDY836,stringsAsFactors=FALSE)

CnrLwlss

#install.packages(c("mvtnorm","fields"))
library(fields)
library(mvtnorm)

grmf=function(nrow,ncol,sigsq=1,V=1){
	pos=expand.grid(row=1:nrow,col=1:ncol)
	mat=matrix(nrow=nrow,ncol=ncol,0)

	dists=rdist(pos)
	#dists[row(dists)==col(dists)]=1