###################################################
### chunk number 1: setup
###################################################
#line 17 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
options(width = 70)
jnk <- capture.output(library(locfit))
library(IWB2011)


###################################################
### chunk number 2: sequence-fastqFile
###################################################
#line 44 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(ShortRead)
library(IWB2011)
fastqFile <- system.file("extdata", "SRR002051.reads1-50k.fastq",
                         package="IWB2011")


###################################################
### chunk number 3: sequence-readFastq
###################################################
#line 54 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
fq <- readFastq(fastqFile)
fq
head(sread(fq))
head(quality(fq))


###################################################
### chunk number 4: sequence-matplot
###################################################
#line 65 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
abc <- alphabetByCycle(sread(fq))
matplot(t(abc[1:4,]), type="l", xlab="Cycle", ylab="Count", lwd=2)


###################################################
### chunk number 5: sequence-quality
###################################################
#line 75 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
m <- as(quality(fq), "matrix")
baseQualPerRead <- apply(m, 1, mean)
fq[baseQualPerRead > 20]


###################################################
### chunk number 6: gapped-file
###################################################
#line 87 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(GenomicRanges)
bamFile <- system.file("extdata", "SRR002051.chrI-V.bam",
                       package="IWB2011")
galn <- readGappedAlignments(bamFile)
head(galn)


###################################################
### chunk number 7: gapped-input
###################################################
#line 99 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
tail(sort(table(cigar(galn))))


###################################################
### chunk number 8: annotation
###################################################
#line 109 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
readScript("create-txdb.R")


###################################################
### chunk number 9: annotation-load
###################################################
#line 115 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(GenomicFeatures)
txdbFile <- system.file("extdata", "sacCer2_sgdGene.sqlite",
                        package="IWB2011")
txdb <- loadFeatures(txdbFile)
txdb


###################################################
### chunk number 10: annotation-transcripts
###################################################
#line 131 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
tscript <- transcripts(txdb, columns="gene_id")
head(tscript)
geneIds <- as(values(tscript)$gene_id, "character")
table(table(geneIds))
## following line changes how 'gene_id' is stored, internally
values(tscript)$gene_id <- geneIds


###################################################
### chunk number 11: annotation-unique-tscript
###################################################
#line 146 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
x <- table(geneIds)
ugeneIds <- names(x)[x==1]
idx <- values(tscript)$gene_id %in% ugeneIds
tscript1 <- tscript[idx]
length(tscript1)


###################################################
### chunk number 12: count-strand-blat
###################################################
#line 163 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
strand(tscript1) <- "*"


###################################################
### chunk number 13: count-union
###################################################
#line 173 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(IRanges)
hits <- countOverlaps(galn, tscript1)
table(hits)
galn1 <- galn[hits==1]
counts <- countOverlaps(tscript1, galn1)


###################################################
### chunk number 14: count-sanity-check
###################################################
#line 187 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
table(seqnames(tscript1[counts!=0]))
hist(asinh(counts[counts!=0]))


###################################################
### chunk number 15: differential-representation-explore
###################################################
#line 209 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
data("tscriptCounts")
head(tscriptCounts)
elementMetadata(tscriptCounts)
df <- as(tscriptCounts, "data.frame")
print(splom(asinh(df), pch="."))


###################################################
### chunk number 16: differential-representation-select
###################################################
#line 225 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
idx <- elementMetadata(tscriptCounts)$Replicate
tscriptCounts1 <- tscriptCounts[,idx != "Technical"]
ridx <- rowSums(as(tscriptCounts1, "data.frame")) != 0
tscriptCounts1 <- tscriptCounts1[ridx,]


###################################################
### chunk number 17: differential-representation
###################################################
#line 246 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(edgeR)
df <- as(tscriptCounts1, "data.frame")
grps <- elementMetadata(tscriptCounts1)$Protocol
dge <- DGEList(df, group=grps)
dge <- calcNormFactors(dge)
dge <- estimateCommonDisp(dge)
etest <- exactTest(dge)
ttags <- topTags(etest, n=10, adjust.method="BH")


###################################################
### chunk number 18: differential-representation-hist
###################################################
#line 260 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
hist(etest$table$p.value, 21)


###################################################
### chunk number 19: annotation-biology
###################################################
#line 271 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
geneIds <- rownames(ttags$table)
idx <- values(tscript)$gene_id %in% geneIds
(topScripts <- tscript[idx])


###################################################
### chunk number 20: annotation-BSgenome
###################################################
#line 286 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(BSgenome.Scerevisiae.UCSC.sacCer2)
Scerevisiae
(topSeqs <- getSeq(Scerevisiae, topScripts, as.character=FALSE))
(abc <- alphabetFrequency(topSeqs, baseOnly=TRUE))
rowSums(abc[,c("G", "C")]) / rowSums(abc)


###################################################
### chunk number 21: annotation-org
###################################################
#line 301 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
library(org.Sc.sgd.db)
(namedIds <- Filter(Negate(is.na), mget(geneIds, org.Sc.sgdGENENAME)))


###################################################
### chunk number 22: appendix-differential-representation-create
###################################################
#line 318 "E:/SinglePackageBuilder/job_42_merlot2/IWB2011/inst/doc/Annotation-Lab.Rnw"
readScript("create-tscriptCounts.R")


