silencing time stuff

# !!! warning messages can be ignored

d <- read.csv("sil_time.csv",sep=";",dec=".",header=FALSE)
dm <- d[,1]
cc <- d[1:24,2]

dm.shapiro <- shapiro.test(dm)
cc.shapiro <- shapiro.test(cc)

# visual normality test
# qqnorm(dm)
# qqnorm(cc)

dm.conf <- t.test(dm)$conf.int
cc.conf <- t.test(cc)$conf.int

# p-value = 0.3811
# there is greater than 5 % chance that
# population medians are the same
cc.dm.wilcox <- wilcox.test(cc,dm,alternative="two.sided",paired=FALSE)

# box plots with 95 % confidence intervals of means
cairo_pdf(filename="cc_dm_boxplot.pdf",onefile=TRUE,width=8,height=8,pointsize=12)
boxplot(cbind(dm,cc),names=c("Diabetes mellitus type 2","Colorectal carcinoma"),main="Drug target distances")
points(c(dm.conf[1],dm.conf[2],cc.conf[1],cc.conf[2])~c(1,1,2,2),col="red",pch="-",cex=5)
dev.off()

# ####

f <- read.csv("silt.csv",sep=";",dec=".",header=TRUE)

g <- f[,2]

a <- length(g)/(max(g)-min(g))
i <- min(g)
j <- max(g)
e <- 10
hh <- 0
t <- 1
c <- 0

while(TRUE){
    c <- c+1
    j <- j-e
    h <- densty(i,j,g)
    print(paste("iteration: interval [",i,";",j,"] average: ",h,sep=""))
    if(h<hh & t==1){
        j <- j+e
        e <- e*2
        h <- hh
        t <- 2
    }
    if(h>hh & t>1){
        e <- 10
        t <- 1
        hh <- h
    }
    if(h<hh & t>1){
        j <- j+e
        if(t==2){
            e <- e/4
        }else{
            e <- round(e/2)
        }
        if(e==1){
            print(paste("i = ",i,sep=""))
            print(paste("j = ",j,sep=""))
            print(paste("c = ",c,sep=""))
            print(paste("h = ",hh,sep=""))
            print(paste("a = ",a,sep=""))
            high <- list(i,j,h)
            break
        }
        h <- hh
        t <- 3
    }
    if(h>hh & t==1){
        hh <- h
    }
}

# searching for lower limit

a <- length(g)/(max(g)-min(g))
i <- min(g)
j <- max(g)
e <- 10
hh <- 0
t <- 1
c <- 0

while(TRUE){
    c <- c+1
    i <- i+e
    h <- densty(i,j,g)
    print(paste("iteration: interval [",i,";",j,"] average: ",h,sep=""))
    if(h<hh & t==1){
        i <- i-e
        e <- e*2
        h <- hh
        t <- 2
    }
    if(h>hh & t>1){
        e <- 10
        t <- 1
        hh <- h
    }
    if(h<hh & t>1){
        i <- i-e
        if(t==2){
            e <- e/4
        }else{
            e <- round(e/2)
        }
        if(e==1){
            print(paste("i = ",i,sep=""))
            print(paste("j = ",j,sep=""))
            print(paste("c = ",c,sep=""))
            print(paste("h = ",hh,sep=""))
            print(paste("a = ",a,sep=""))
            low <- list(i,j,h)
            break
        }
        h <- hh
        t <- 3
    }
    if(h>hh & t==1){
        hh <- h
    }
}

densty <- function(i,j,d){
    a <- length(d)/(max(d)-min(d))
    ij <- subset(d,d>i & d<j)
    aa <- length(ij)/(max(ij)-min(ij))
    return(aa)
}

cairo_pdf(filename="silt_hist.pdf",onefile=TRUE,width=12,height=8,pointsize=12)
hist(g,breaks=50,xlab="Silencing time",main="Frequency of drug targets\nas a function of silencing time")
abline(v=31,col="red",lwd=2)
abline(v=87,col="red",lwd=2)
dev.off()