### R code from vignette source 'Rngsapps.Rnw'

###################################################
### code chunk number 1: Rngsapps.Rnw:130-135
###################################################
targets <- data.frame(Samples=c("AP3 domain, flower stage 4", "AP3 domain, flower stage 4", 
                      "Translatome, flower stage 4", "Translatome, flower stage 4"),
                      Factor=c("AP3", "AP3", "TRL", "TRL"),
                      Fastq=c("SRR064154", "SRR064155", "SRR064166", "SRR064167"))
targets


###################################################
### code chunk number 2: Rngsapps.Rnw:143-152
###################################################
library(rtracklayer); library(GenomicRanges); library(Rsamtools)
gff <- import.gff("./data/TAIR10_GFF3_trunc.gff", asRangedData=FALSE)
seqlengths(gff) <- end(ranges(gff[which(elementMetadata(gff)[,"type"]=="chromosome"),]))
subgene_index <- which(elementMetadata(gff)[,"type"] == "gene")
gffsub <- gff[subgene_index,] # Returns only gene ranges
strand(gffsub) <- "*" # For strand insensitive analysis
gffsub[1:4,1:2]
ids <- elementMetadata(gffsub)[, "group"]
gffsub <- split(gffsub) # Coerce to GRangesList


###################################################
### code chunk number 3: Rngsapps.Rnw:160-175
###################################################
samples <- as.character(targets$Fastq)
samplespath <- paste("./data/", samples, sep="")
countDF <- data.frame(row.names=ids)
for(i in samplespath) {
        # aligns <- readBamGappedAlignments(i) # Substitute next two lines with this one.
        aligns <- read.table(i, header=TRUE)
        aligns <- GRanges(seqnames=aligns$seqnames, IRanges(aligns$start, aligns$end), strand=aligns$strand)
        counts <- countOverlaps(gffsub, aligns)
        countDF <- cbind(countDF, counts)
}
colnames(countDF) <- samples
rownames(countDF) <- gsub(".*=", "", rownames(countDF))
countDF[1:4,]
write.table(countDF, "./data/countDF", quote=FALSE, sep="\t", col.names = NA)
countDF <- read.table("./data/countDF")


###################################################
### code chunk number 4: Rngsapps.Rnw:183-192
###################################################
returnRPKM <- function(counts, gffsub) {
        geneLengthsInKB <- sum(width(gffsub))/1000 # Number of bases per exonRanges element in kbp
        millionsMapped <- sum(counts)/1e+06 # Factor for converting to million of mapped reads.
        rpm <- counts/millionsMapped # RPK: reads per kilobase of exon model.
        rpkm <- rpm/geneLengthsInKB # RPKM: reads per kilobase of exon model per million mapped reads.
        return(rpkm)
}
countDFrpkm <- apply(countDF, 2, function(x) returnRPKM(counts=x, gffsub=gffsub))
countDFrpkm[1:4,]


###################################################
### code chunk number 5: Rngsapps.Rnw:200-202
###################################################
d <- cor(countDF, method="pearson")
plot(hclust(dist(1-d))) # Sample tree


###################################################
### code chunk number 6: Rngsapps.Rnw:210-213
###################################################
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/colAg.R")
countDFrpkm_mean <- colAg(myMA=countDFrpkm, group=c(1,1,2,2), myfct=mean)
countDFrpkm_mean[1:4,]


###################################################
### code chunk number 7: Rngsapps.Rnw:216-222
###################################################
countDFrpkm_mean <- cbind(countDFrpkm_mean, log2ratio=log2(countDFrpkm_mean[,1]/countDFrpkm_mean[,2]))
countDFrpkm_mean <- countDFrpkm_mean[is.finite(countDFrpkm_mean[,3]), ]
degs2fold <- countDFrpkm_mean[countDFrpkm_mean[,3] >= 1 | countDFrpkm_mean[,3] <= -1,]
degs2fold[1:4,]
write.table(degs2fold, "./data/degs2fold", quote=FALSE, sep="\t", col.names = NA)
degs2fold <- read.table("./data/degs2fold")


###################################################
### code chunk number 8: Rngsapps.Rnw:230-242
###################################################
library(DESeq)
countDF <- read.table("./data/countDF")
conds <- targets$Factor
cds <- newCountDataSet(countDF, conds) # Creates object of class CountDataSet derived from eSet class
counts(cds)[1:4, ] # CountDataSet has similar accessor methods as eSet class.
cds <- estimateSizeFactors(cds) # Estimates library size factors from count data. Alternatively, one can provide here the true library sizes with sizeFactors(cds) <- c(..., ...)
cds <- estimateDispersions(cds) # Estimates the variance within replicates
res <- nbinomTest(cds, "TRL", "AP3") # Calls DEGs with nbinomTest
res <- na.omit(res)
res2fold <- res[res$log2FoldChange >= 1 | res$log2FoldChange <= -1,]
res2foldpadj <- res2fold[res2fold$padj <= 0.01, ]
res2foldpadj[1:4,1:8]


###################################################
### code chunk number 9: Rngsapps.Rnw:250-260
###################################################
library(edgeR)
countDF <- read.table("./data/countDF")
y <- DGEList(counts=countDF, group=conds) # Constructs DGEList object
y <- estimateCommonDisp(y) # Estimates common dispersion
y <- estimateTagwiseDisp(y) # Estimates tagwise dispersion
et <- exactTest(y, pair=c("TRL", "AP3")) # Computes exact test for the negative binomial distribution.
topTags(et, n=4)
edge <- as.data.frame(topTags(et, n=50000)) 
edge2fold <- edge[edge$logFC >= 1 | edge$logFC <= -1,]
edge2foldpadj <- edge2fold[edge2fold$adj.P.Val <= 0.01, ]


###################################################
### code chunk number 10: Rngsapps.Rnw:267-274
###################################################
bothDF <- merge(res, countDFrpkm_mean, by.x=1, by.y=0, all=TRUE); bothDF <- na.omit(bothDF)
cor(bothDF[,"log2FoldChange"], bothDF[,"log2ratio"], method="spearman")
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/overLapper.R")
setlist <- list(edgeR=rownames(edge2foldpadj), DESeq=as.character(res2foldpadj[,1]), RPKM=rownames(degs2fold))
OLlist <- overLapper(setlist=setlist, sep="_", type="vennsets")
counts <- sapply(OLlist$Venn_List, length)
vennPlot(counts=counts)


###################################################
### code chunk number 11: Rngsapps.Rnw:282-291
###################################################
library(GOstats); library(GO.db); library(ath1121501.db)
geneUniverse <- rownames(countDF) 
geneSample <- res2foldpadj[,1]
params <- new("GOHyperGParams", geneIds = geneSample, universeGeneIds = geneUniverse, 
        annotation="ath1121501", ontology = "MF", pvalueCutoff = 0.5,
	conditional = FALSE, testDirection = "over")
hgOver <- hyperGTest(params)
summary(hgOver)[1:4,]
htmlReport(hgOver, file = "data/MyhyperGresult.html")


###################################################
### code chunk number 12: Rngsapps.Rnw:302-304
###################################################
cov <- coverage(aligns)
cov[1:2]


###################################################
### code chunk number 13: Rngsapps.Rnw:314-319
###################################################
aligns <- read.table(samplespath[3], header=TRUE)
aligns <- GRanges(seqnames=aligns$seqnames, IRanges(aligns$start, aligns$end), strand=aligns$strand)
cov <- coverage(aligns)
islands <- slice(cov, lower = 50)
islands[[1]]


###################################################
### code chunk number 14: Rngsapps.Rnw:326-335
###################################################
plotCov <- function(mycov=cov, mychr=1, mypos=c(1,1000), mymain="Coverage", ...) {
    op <- par(mar=c(8,3,6,1))
    plot(as.numeric(mycov[[mychr]][mypos[1]:mypos[2]]), type="l",
    lwd=1, col="blue", ylab="", main=mymain, xlab="", xaxt="n", ...)
    axis(1, las = 2, at=seq(1,mypos[2]-mypos[1], length.out= 10),
    labels=as.integer(seq(mypos[1], mypos[2], length.out= 10)))
    par(op)
    }
plotCov(mycov=cov, mychr="1", mypos=c(1,length(cov[["1"]])))


###################################################
### code chunk number 15: Rngsapps.Rnw:343-352
###################################################
chip_signal_list <- sapply(samples, list)
for(i in seq(along=samplespath)) {
	aligns <- read.table(samplespath[i], header=TRUE)
	aligns <- GRanges(seqnames=aligns$seqnames, IRanges(aligns$start, aligns$end), strand=aligns$strand)
	chip_signal_list[[i]] <- aligns
}
chip_signal_list[["SRR064154"]][1:4,]
chip_signal_list <- sapply(names(chip_signal_list), function(x) resize(chip_signal_list[[x]], width = 200))
for(i in names(chip_signal_list)) start(chip_signal_list[[i]])[start(chip_signal_list[[i]]) < 0 ] <- 1


###################################################
### code chunk number 16: Rngsapps.Rnw:360-366
###################################################
Nreads <- sapply(names(chip_signal_list), function(x)  length(chip_signal_list[[x]]))
normfactor <- 10^6/Nreads
chip_signal_list <- sapply(names(chip_signal_list), function(x) coverage(chip_signal_list[[x]]) * normfactor[x])
chip_signal_list[["SRR064154"]][1:2,]
chip_peak_list <- sapply(names(chip_signal_list), function(x) slice(chip_signal_list[[x]], lower=5))
chip_peak_list[[1]][[1]][1:3]


###################################################
### code chunk number 17: Rngsapps.Rnw:374-392
###################################################
library(BayesPeak)
sig <- read.table(samplespath[1], header=TRUE)
bgr <- read.table(samplespath[3], header=TRUE)
sig <- RangedData(space=sig[,1], IRanges(start=sig[,2], end=sig[,3]), strand=sig[,5])
bgr <- RangedData(space=bgr[,1], IRanges(start=bgr[,2], end=bgr[,3]), strand=bgr[,5])

raw.output <- bayespeak(treatment=sig, control=bgr)   
# unreliable.jobs <- log(raw.output$QC$lambda1) < 1.5 # Removal of false positives due to overfitting.
# bpeaks <- as.data.frame(summarise.peaks(raw.output, method = "lowerbound", exclude.jobs = unreliable.jobs)) 
bpeaks <- as.data.frame(summarise.peaks(raw.output, method = "lowerbound"))

bpeaksIR <- IRanges(start=bpeaks[, 2], end=bpeaks[, 3])
bpeaksIR <- lapply(unique(bpeaks[,1]), function(x) IRanges(start=bpeaks[bpeaks[,1]==x,2], end=bpeaks[bpeaks[,1]==x,3]))
names(bpeaksIR) <- unique(names(sig))
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/chipseqFct.R") # Imports the rangeCoverage function.
sigcovDF <- rangeCoverage(summaryFct=viewMaxs, myname="sig_", peaksIR=bpeaksIR, sig=sig, readextend=1, normfactor=10^6/length(start(sig))) 
bgrcovDF <- rangeCoverage(summaryFct=viewMaxs, myname="bgr_", peaksIR=bpeaksIR, sig=bgr, readextend=1, normfactor=10^6/length(start(bgr))) 
bpeaksDF <- cbind(bpeaks, sigcovDF[,-1], bgrcovDF[,-1]); bpeaksDF[1:4,] 


###################################################
### code chunk number 18: Rngsapps.Rnw:399-400
###################################################
plotCov(mycov=chip_signal_list[[3]], mychr="1", mypos=c(2000,12000), ylim=c(0,140))


###################################################
### code chunk number 19: Rngsapps.Rnw:408-412
###################################################
simple_peak <- as.data.frame(as(chip_peak_list[[1]], "IRangesList"))
   # simple_peak <- as.data.frame(chip_peak_list[[1]])
commonpeaks <- subsetByOverlaps(as(bpeaks, "RangedData"), as(simple_peak, "RangedData"), minoverlap=100)
bpeaksDF[bpeaksDF$start %in% start(commonpeaks),][1:4,]


###################################################
### code chunk number 20: Rngsapps.Rnw:419-427
###################################################
library(ChIPpeakAnno)
annoRD <- unlist(gffsub)
names(annoRD) <- gsub(".*=", "", elementMetadata(annoRD)[, "group"])
annoRD <- as(annoRD, "RangedData"); 
peaksRD <- RangedData(space=bpeaksDF$space, IRanges(bpeaksDF$start, bpeaksDF$end))
annotatedPeak <- annotatePeakInBatch(peaksRD, AnnotationData = annoRD)
as.data.frame(annotatedPeak)[1:4,1:11]
bpeaksDF[1:4,]


###################################################
### code chunk number 21: Rngsapps.Rnw:434-435
###################################################
sessionInfo()