ericminikel

#!/broad/software/free/Linux/redhat_5_x86_64/pkgs/r_3.0.2/bin/Rscript

# Eric Vallabh Minikel
# CureFFI.org
# 2014-01-14
# example of how to use optparse in R scripts

# usage: ./exampleRScript1.r -a thisisa -b hiagain
#        ./exampleRScript1.r --avar thisisa --bvar hiagain

ericminikel

# Eric Minikel
# 2014-01-12
# explore relationship between inoculated dose and incubation time
# this is all the code from this post: 
# http://www.cureffi.org/2014/01/12/prion-kinetic-models-relationship-inoculum-titer-incubation-time/

logdose = runif(n=100,min=0,max=10) # simulate 100 different experiments with diff titers
days = 10^( (logdose-26.66)/(-12.99) ) + 40 + rnorm(n=100,m=0,s=2) # apply Prusiner's model
plot(logdose,days,xlab="log10(LD50 units)",ylab="days to illness",pch=19,
     main="incubation time interval bioassay\nPrusiner 1982b - Sc237 in hamsters",

ericminikel

primes = []
x = 2
nthprime = 1000 # set to 1000 to get the 1000th prime number
while len(primes) < nthprime:
    isPrime = True
    for item in primes:
        if item > x**.5:
            break
        if x % item == 0:
            isPrime = False

ericminikel

require(rcdk) # chemical informatics functionality in R
require(gap)  # for qq plots later
options(stringsAsFactors=FALSE)
setwd('c:/sci/037cinf/analysis/sar/')

# plot molecules in R plot window instead of separate Java window
rcdkplot = function(molecule,width=500,height=500,marg=0,main='') {
    par(mar=c(marg,marg,marg,marg)) # set margins to zero since this isn't a real plot
    temp1 = view.image.2d(molecule,width,height) # get Java representation into an image matrix. set number of pixels you want horiz and vertical
    plot(NA,NA,xlim=c(1,10),ylim=c(1,10),xaxt='n',yaxt='n',xlab='',ylab='',main=main) # create an empty plot

ericminikel

# simulate PrPC downregulation
kdeg = log(2)/1.5 # 36 half life for PrPSc
kconv = 1.7 #  a conversion rate that works out well for wt mice
PrPSc0 = 10^1.5 # initial titer as of 4 dpi
symptomatic_threshold = 10^9 # titer when mice become asymptomatic
downreg_threshold = sqrt(PrPSc0 * symptomatic_threshold) # geometric mean for 1/2way through disease course

plot(NA,NA,xlim=range(times),ylim=c(1,10^10),log="y", # empty plot over 300 days and 10 log10s
     xlab='days post-inoculation', ylab='prion titer', main=paste("kconversion = ",kconv,sep='')) 

ericminikel

########### what percent of possible degradation is achieved by 1 day, 5 days and 6 days?

maxint = .Machine$integer.max # longest possible time in future to model

# transcription inhibition
r = l0 * 0.50 # 
l = l0 * 1.00 # 
a = m0 * 1.00
m = m0 * 1.00 
b = v0 * 1.00

ericminikel

# Eric Minikel
# CureFFI.org
# 2013-12-03
# Comparison of numerically simulated and analytically derived PrP degradation models

# half lives from literature
mrna_thalf  =  7 # Pfeiffer 1993  http://www.ncbi.nlm.nih.gov/pubmed/8095862
prpc_thalf  =  5 # Borcheldt 1990 http://www.ncbi.nlm.nih.gov/pubmed/1968466/
prpsc_thalf = 30 # Peretz 2001    http://www.ncbi.nlm.nih.gov/pubmed/11507642
division_time = 24 # Ghaemmaghami 2007  http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2084281/

ericminikel

mRNA
====

R(0) = R0
R(infty) = R1

R(t) = R1 + (R0 - R1)exp(-lambda*t)

PrPC
====

ericminikel

# Eric Minikel
# CureFFI.org
# 2013-11-06
# Anti-prion compounds plotted in the space of drugs and chemical libraries

setwd('c:/sci/037cinf/analysis/1')
options(stringsAsFactors=FALSE)
library(stringr)
library(rcdk) # cheminformatics
library(sp) # for point-in-polygon operations

ericminikel

# Eric Minikel
# CureFFI.org
# 2013-07-04
# Bash script to generate FPKMs from Ensembl Human Bodymap 2.0 RNA-seq data

# Announcement of Human BodyMap 2.0 data availability: http://www.ensembl.info/blog/2011/05/24/human-bodymap-2-0-data-from-illumina/
# BAM list: ftp://ftp.ensembl.org/pub/release-70/bam/homo_sapiens/genebuild/

# STEP 1: DOWNLOAD THE BAMS FROM ENSEMBL