rscript.R

# This is example r script for punctaanalysis tool kit.
# If you use batch program, you can use this script.
# To use it, start R and change the working directory to where your pictures and data are stored, 
# and just copy everything here and paste it on R terminal and push return key.
# Once you have run this script, you will get some pdf density graphs,
# and ready to test data sets which has been log2 transformed. ex. wtfixwidth, daffixwidth...
# Then you can test with something like this 
#
# t.test(wtfixwidth, daffixwidth)
# 
# I don't know what kind of test is best. Making a choice is your job.
# 
# Attention: log transformation affect the result of parametric test.
#
# Also, here is an example of density graph
# plot(density(wtwidth), col=2, xlim=c(-2, 2), ylim=c(0,1))
# par(new =T)
# plot(density(dafwidth), xlim=c(-2, 2), ylim=c(0,1))
# 
# I hope this example file helps you to learn R scripting.  20080407 Taizo

######## read data
count <- read.table("count.txt", header = TRUE);
width <- read.table("width.txt", header = TRUE);
gap <- read.table("gap.txt", header = TRUE);
distance <- read.table("distance.txt", header = TRUE);
intensity <- read.table("intensity.txt", header = TRUE);
lineardensity <- read.table("lineardensity.txt", header =TRUE);
fixwidth <- read.table("fixwidth.txt", header =TRUE);
fixvolume <- read.table("fixvolume.txt", header =TRUE);
fixgap <- read.table("fixgap.txt", header =TRUE);
#fixwidthlineardensity <- read.table("fixwidthlineardensity", header =TRUE)
#fixgaplineardensity <- read.table("fixgaplineardensity", header =TRUE)

###########  spadework 
# log2transform. If you need to use raw data, comment out this part
count$count <- log2(count$count)
width$width <- log2(width$width)
gap$gap <- log2(gap$gap)
distance$distance <- log2(distance$distance)
#log2 transform of intensiy case NaN error?
#intensity$intensity <- log2(intensity$intensity)
lineardensity$lineardensity <- log2(lineardensity$lineardensity)
fixwidth$fixwidth <- log2(fixwidth$fixwidth)
fixvolume$fixvolume <- log2(fixvolume$fixvolume)
fixgap$fixgap <- log2(fixgap$fixgap)
#fixwidthlineardensity$fixwidthlineardensity <- log2(fixwidthlineardensity$fixwidthlineardensity)
#fixgaplineardensity$fixgaplineardensity <- log2(fixgaplineardensity$fixgaplineardensity)

#list up data types what you want to analyse
analysisdata<- c("width", "fixwidth","gap","fixvolume", "intensity", "lineardensity", "distance", "fixgap")
datanum <- length(analysisdata)

#Detect genotypes from files automatically
genotype.list <- levels(count$genotype);
numberofgenotype <- length(genotype.list);

# dataframe operation. New vector data are made for each genotype
# follwing makes genotypewidh, genotypegaps etc. (ex: wtwidth, e5ewidth, fc16witdh,....wtvolume, e5volme.)
# "eval(parse(text = ....)))" run the text as command.

for (j in 1:datanum)
{
for (i in 1:numberofgenotype)
{
# This is doing like this. wtwidth <- width[width$genotype == "wt",]
eval(parse(text = paste(genotype.list[i], analysisdata[j], "<-", analysisdata[j], "[", analysisdata[j], "$genotype == genotype.list[", i, "], ]", sep ="")))
# This is doing like this. wtwidth <- wtwidth$width
eval(parse(text = paste(genotype.list[i], analysisdata[j], "<-", genotype.list[i], analysisdata[j], "$",  analysisdata[j], sep = "")))
}
}

######## density graphs
library(KernSmooth)
#list up data types what you want to analyse with density graph here
analysisdata<- c("width", "fixwidth","gap","fixvolume", "intensity", "distance", "fixgap")
datanum <- length(analysisdata)
# genotypes that you want to analyse. everything or can be designate
##################### You may list up genotypes what you want to analyse here.
#genotype.list <- c("wt", "e719","e5", "fc16", "e1598e5", "e1598fc16");
numberofgenotype <- length(genotype.list);


for (j in 1:datanum)
{
xmin <- NULL
xmax <- NULL
ymin <- NULL
ymax <- NULL
for (i in 1:numberofgenotype)
{
# This is doing something like this "wt <- bkde(wtwidth)"
eval(parse(text = paste(genotype.list[i], "<- bkde(", genotype.list[i], analysisdata[j], ")", sep="")))
# appending min and max values to use as axis depiction later
eval(parse(text = paste("xmin <- append(xmin, min(", genotype.list[i], "$x))", sep="")))
eval(parse(text = paste("xmax <- append(xmax, max(", genotype.list[i], "$x))", sep="")))
eval(parse(text = paste("ymin <- append(ymin, min(", genotype.list[i], "$y))", sep="")))
eval(parse(text = paste("ymax <- append(ymax, max(", genotype.list[i], "$y))", sep="")))
}

filename <- paste("density",analysisdata[j], ".pdf", sep="")
pdf(file=filename)
# prepare axis of graph
figtitle <- paste("Density graph of ", analysisdata[j], sep="")
plot(0, main=figtitle, type = 'n', xlab = analysisdata[j], ylab = "kerneldensity", xlim = c(min(xmin), max(xmax)), ylim = c(min(ymin), max(ymax)))

# drow lines 
for (i in 1:numberofgenotype)
{
eval(parse(text = paste("lines(", genotype.list[i], ",type = 'l', col = ", i, ")", sep="")))
}

# write legend
par(usr = c(0, 1, 0, 1))
legend(0.7, 1, genotype.list, col = 1:numberofgenotype, lty = 1)
dev.off()
}
detach("package:KernSmooth")