Skip to content

Instantly share code, notes, and snippets.

View CnrLwlss's full-sized avatar

Conor Lawless CnrLwlss

View GitHub Profile
@CnrLwlss
CnrLwlss / pyMC_Hierarchy.py
Last active February 5, 2016 21:57
Hierarchical QFA with pyMC
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
CnrLwlss / SHM.c
Created December 4, 2015 16:34
Truncating priors
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
CnrLwlss / Okazaki.R
Last active November 11, 2015 17:29
# 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
CnrLwlss / OkazakiAnalysis.sh
Last active January 25, 2016 11:28
Shell script for aligning reads from Smith & Whitehouse 2012.
# 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
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
CnrLwlss / MultivariateNormalGauss.R
Created May 15, 2015 16:01
Attempt to simulate Gaussian Random Markov Field
#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
@CnrLwlss
CnrLwlss / Harmonograph.py
Created May 6, 2015 21:57
Simulate an arbitrary number of linked oscillators (harmonograph), drawing output to .png file.
import numpy as np
import math
import matplotlib.pyplot as plt
import webbrowser
def makePlot(xvals,yvals,fname="harmonograph.png",fopen=False,w=50,h=50,width=0.1):
maxs=[np.max(np.abs(xvals)),np.max(np.abs(yvals))]
farthest=max(maxs)
drw=plt.figure(1)
@CnrLwlss
CnrLwlss / Vicky_Ztest.py
Created March 13, 2015 16:30
Z-test in python and alternative...
# Demonstrating significant differences between a
# vector of measurements and a single value
# Using the statsmodels package for doing test
# Using numpy to generate some fake data
from statsmodels.stats import weightstats as stests
import numpy as np
data=np.random.normal(loc=3.4,scale=0.1,size=100)
singleValue=3.3
@CnrLwlss
CnrLwlss / telo.hs
Created December 2, 2014 21:00
Simulating cell population dynamics with Haskell
data Cell = Cell {teloLength :: Double, birthDate :: Double, lifeSpan :: Double, seedCell :: Integer}
instance Eq Cell where (Cell t1 _ _ _) == (Cell t2 _ _ _) = t1 == t2
instance Ord Cell where (Cell t1 _ _ _) `compare` (Cell t2 _ _ _) = t1 `compare` t2
instance Show Cell where show (Cell t b _ _) = show t ++ " " ++ show b
data CellLine x = Senescent | Mother x (CellLine x) (CellLine x) deriving (Show, Read, Eq)
-- Functions modelling cell biology
shortenTelo :: Cell -> Cell
shortenTelo (Cell telo x y oldseed) = Cell (telo - delta) x y newseed
import Control.Parallel (par, pseq)
import Control.DeepSeq
import System.Environment (getArgs)
import Data.Time.Clock (diffUTCTime, getCurrentTime)
data Tree x = Empty | Node x (Tree x) (Tree x) deriving (Show, Read, Eq)
-- Create instance of NFData for Tree data type (defining "normal form")
instance NFData a => NFData (Tree a) where
rnf Empty = ()