R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

## Question 1: Random Scatter
set.seed(296)
data = read.csv("hcmv-263hxkx-1qhtfgz.txt", header=T)
site.original = data.matrix(data)
n.base <- 229354 #number of DNA sequence bases
n.site <- 296 #296 palindrome sites
Title1 = 'Palindrome locations'
library(lattice)
#Create a scatter plot of the locations of the original data
site.original.plot =stripplot(site.origianl, pch=5, cex=0.8,col="black", 
           main=Title1, xlab= 'Base Pairs', jitter.data=T, grid="h") 
site.original.plot

#randomly generate 3 different 296 panlindrome locations from 229354 DNA sequence bases
site.random1 = round(runif(n.site,0,n.base))
site.random2 = round(runif(n.site,0,n.base))
site.random3 = round(runif(n.site,0,n.base))
#plot 3 scatter plots of runif
stripplot(site.random1, pch=5, cex=0.8, jitter.data=T, grid = "h", col="blue", #3 different scatter plots
          main = Title1, xlab = 'Base Pairs') #1-dimensional scatter plot

     
stripplot(site.random2, pch=5, cex=0.8, jitter.data=T, grid = "h", col="red", 
          main = Title1, xlab = 'Base Pairs')

     
stripplot(site.random3, pch=5, cex=0.8, jitter.data=T, grid = "h", col="green", 
          main = Title1, xlab = 'Base Pairs', add=T)

#plot histograms comparing sample 1 and sample 2
hist(site.random1, breaks=85, col=rgb(1,1,0,0.7), main="Histogram of sample 1 and sample 2", xlab="Base Pair",xaxt='n')
hist(site.random2, breaks=85, col=rgb(0,1,1,0.4), main="Histogram of sample 2", xlab="Base Pair", add=T, xaxt='n')

axis(1, at=seq(0, n.base+10000, 10000))
legend("topright", c("sample1", "sample2"), col=c(rgb(1,1,0,0.7), rgb(0,1,1,0.4)), lwd=8)

#plot histograms comparing original data and sample3
hist(site.random3, breaks=85, col=rgb(1,0,0,0.5), main="Histogram of original data and sample 3", xlab="BasePair", xaxt='n')
hist(site.origianl, breaks=85, col=rgb(0,0,1,0.5), main="Histogram of original data", xlab="Base Pair", add=T, xaxt='n')

axis(1, at=seq(0, n.base+10000, 10000))
legend("topright", c("sample3", "data"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), lwd=8)

for(x in c(51,57,66)){
 count.sample = rep(n.site/x, x)
 count.data = as.vector(table(cut(data, breaks=seq(0, n.base, length.out = x + 1), include.lowest=TRUE)))
 
 hist(data, breaks=x, probability=T, col=rgb(1,0,0,0.5))
}

Error in cut.default(data, breaks = seq(0, n.base, length.out = x + 1),  : 
  'x' must be numeric

#chi-squared
data = read.csv("hcmv-263hxkx-1qhtfgz.txt", header=T)
array.data = data.matrix(data)
n.base <- 229354 #number of DNA sequence bases
n.site <- 296 #296 palindrome sites
gene <- seq(1,n.base)
reg.split <- function(n.region, gene, site){ #function to split regions into n regions
  count.int <- table(cut(site, breaks = seq(1, length(gene), length.out=n.region+1), include.lower=TRUE))
  count.vector <- as.vector(count.int)
  count.tab <-table(count.vector)
  return(count.tab)
}
n.region <- 50 # n=50 regions here
reg.split(n.region, gene, site.random) #generate a table for the dataset

count.vector
 1  2  3  4  5  6  7  8  9 10 12 
 1  2  3  6 14  7  4  7  2  2  2 

chisqtable <- function(n.region, site, n.base){ #chi-squared table for the dataset
  n <- length(site) 
  lambda.est <- n/n.region #estimate for lambda
  count.int <- table(cut(site, breaks = seq(1, length(gene), length.out=n.region+1), include.lowest=TRUE)) #divide into n.region number of non-overlapping intervals
  count.vector <- as.vector(count.int)  #get the count levels range
  count.range <- max(count.vector) - min(count.vector) +1
  
  table <- matrix(rep(NA, count.range*3), count.range, 3) #create a table
  for(i in 1:count.range){
    offset <-min(count.vector) - 1
    table[i, 1] <- i + offset #first column is the count level
    table[i, 2] <- sum(count.vector == i + offset) #2nd column is the observed count
    if ((i + offset == min(count.vector) && (min(count.vector)) != 0))
        table[i, 3] <- ppois(i+offset, lambda.est)*n.region
        else if (i + offset == max(count.vector))
          table[i, 3] <- (1 - ppois(i + offset - 1, lambda.est))*n.region
        else 
          table[i, 3] <- (ppois(i + offset, lambda.est) - ppois(i + offset - 1, lambda.est))*n.region
  }
  
  return(table)
}
site.random.table <- chisqtable(n.region, site.random, n.base)
site.random.table

      [,1] [,2]      [,3]
 [1,]    1    1 0.9290793
 [2,]    2    2 2.3526650
 [3,]    3    3 4.6425922
 [4,]    4    6 6.8710365
 [5,]    5   14 8.1353072
 [6,]    6    7 8.0268365
 [7,]    7    4 6.7884103
 [8,]    8    7 5.0234236
 [9,]    9    2 3.3042964
[10,]   10    2 1.9561435
[11,]   11    0 1.0527609
[12,]   12    2 0.9174487

install.packages("gplots")

also installing the dependencies ‘gtools’, ‘gdata’

trying URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.4/gtools_3.8.1.tgz'
Content type 'application/x-gzip' length 262688 bytes (256 KB)
==================================================
downloaded 256 KB

trying URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.4/gdata_2.18.0.tgz'
Content type 'application/x-gzip' length 1141723 bytes (1.1 MB)
==================================================
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.4/gplots_3.0.1.1.tgz'
Content type 'application/x-gzip' length 509475 bytes (497 KB)
==================================================
downloaded 497 KB

The downloaded binary packages are in
    /var/folders/xy/341ywwj15pq82lydhqf89xr80000gn/T//RtmpWrrJrc/downloaded_packages

#Monte Carlo Practice 
set.seed(12)
n.site 

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.