withShortRead.R
In Genominator: Analyze, manage and store genomic data

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

###################################################
### code chunk number 1: withShortRead.Rnw:5-6
###################################################
options(width=80)


###################################################
### code chunk number 2: packages
###################################################
require(Genominator)
require(ShortRead)
require(yeastRNASeq)


###################################################
### code chunk number 3: yeastAligned
###################################################
data(yeastAligned)
yeastAligned[[1]]
sapply(yeastAligned, length)


###################################################
### code chunk number 4: filesInYeastRNASeq
###################################################
list.files(file.path(system.file(package = "yeastRNASeq"), "reads"), 
           pattern = "bowtie")


###################################################
### code chunk number 5: chrMap
###################################################
chrMap <- paste("Scchr", sprintf("%02d", 1:16), sep = "")
unique(chromosome((yeastAligned[[1]])))


###################################################
### code chunk number 6: filesArgument
###################################################
files <- list.files(file.path(system.file(package = "yeastRNASeq"), "reads"), 
                    pattern = "bowtie", full.names = TRUE)
names(files) <- sub("_f\\.bowtie\\.gz", "", basename(files))
names(files)


###################################################
### code chunk number 7: importFromAlignedReads
###################################################
eData <- importFromAlignedReads(files, chrMap = chrMap,
                                dbFilename = "my.db",
                                tablename = "raw", type = "Bowtie")
eData
head(eData)


###################################################
### code chunk number 8: biomaRt (eval = FALSE)
###################################################
## require(biomaRt)
## mart <- useMart("ensembl", "scerevisiae_gene_ensembl")
## attributes.gene <- c("ensembl_gene_id", "chromosome_name",  "start_position", 
##                      "end_position", "strand", "gene_biotype")
## attributes.tr <- c("ensembl_gene_id", "ensembl_transcript_id", "ensembl_exon_id", "chromosome_name",  "start_position", 
##                    "end_position", "strand", "gene_biotype", "exon_chrom_start", "exon_chrom_end", "rank")
## ensembl.gene <- getBM(attributes = attributes.gene,  mart = mart)
## ensembl.transcript <- getBM(attributes = attributes.tr,  mart = mart)


###################################################
### code chunk number 9: ensemblAnno
###################################################
data(yeastAnno.sources)
ensembl.gene <- yeastAnno.sources$ensembl.gene
ensembl.transcript <- yeastAnno.sources$ensembl.transcript
head(ensembl.gene, n = 2)
head(ensembl.transcript, n = 2)
dim(ensembl.gene)
dim(ensembl.transcript)
subset(ensembl.gene, ensembl_gene_id == "YPR098C")
subset(ensembl.transcript, ensembl_gene_id == "YPR098C")
length(unique(ensembl.transcript$ensembl_transcript_id))


###################################################
### code chunk number 10: rtracklayer (eval = FALSE)
###################################################
## require(rtracklayer)
## session <- browserSession()
## genome(session) <- "sacCer2"
## ucsc.sgdGene <- getTable(ucscTableQuery(session, "sgdGene"))
## ucsc.ensGene <- getTable(ucscTableQuery(session, "ensGene"))


###################################################
### code chunk number 11: yeastAnno.sources (eval = FALSE)
###################################################
## yeastAnno.sources <- list(ensembl.gene = ensembl.gene,
##                           ensembl.transcript = ensembl.transcript,
##                           ucsc.sgdGene = ucsc.sgdGene,
##                           ucsc.ensGene = ucsc.ensGene)
## save(yeastAnno.sources, file = "yeastAnno.sources.rda")


###################################################
### code chunk number 12: ucscAnno
###################################################
data(yeastAnno.sources)
ucsc.sgdGene <- yeastAnno.sources$ucsc.sgdGene
ucsc.ensGene <- yeastAnno.sources$ucsc.ensGene
head(ucsc.sgdGene, n = 2)
head(ucsc.ensGene, n = 2)
dim(ucsc.sgdGene)
dim(ucsc.ensGene)
subset(ucsc.sgdGene, name == "YPR098C")
subset(ucsc.ensGene, name == "YPR098C")
subset(ucsc.sgdGene, name == "YER102W")
subset(ucsc.ensGene, name == "YER102W")


###################################################
### code chunk number 13: yAnno
###################################################
yAnno <- yeastAnno.sources$ensembl.transcript
yAnno$chr <- match(yAnno$chr, c(as.character(as.roman(1:16)), "MT", "2-micron"))
yAnno$start <- yAnno$start_position
yAnno$end <- yAnno$end_position
rownames(yAnno) <- yAnno$ensembl_exon_id
yAnno.simple <- yAnno[yAnno$chr %in% 1:16, c("chr", "start", "end", "strand")]
head(yAnno.simple, n = 2)
head(yAnno, n = 2)


###################################################
### code chunk number 14: validAnno
###################################################
validAnnotation(yAnno)


###################################################
### code chunk number 15: geneCounts.1
###################################################
geneCounts.1 <- summarizeByAnnotation(eData, yAnno, ignoreStrand = TRUE)
head(geneCounts.1)


###################################################
### code chunk number 16: geneCounts.2
###################################################
geneCounts.2 <- summarizeByAnnotation(eData, yAnno, ignoreStrand = TRUE, 
                                      groupBy = "ensembl_gene_id")
head(geneCounts.2)


###################################################
### code chunk number 17: makeGeneRep
###################################################
yAnno.UI <- makeGeneRepresentation(yAnno, type = "UIgene", gene.id = "ensembl_gene_id",
                                   transcript.id = "ensembl_transcript_id")
head(yAnno.UI)


###################################################
### code chunk number 18: gof (eval = FALSE)
###################################################
## groups <- gsub("_[0-9]_f", "", colnames(geneCounts))
## groups
## plot(regionGoodnessOfFit(geneCounts, groups), chisq = TRUE)


###################################################
### code chunk number 19: computePrimingWeights
###################################################
weightsList <- lapply(yeastAligned, computePrimingWeights, 
                      unbiasedIndex = 20:21, weightsLength = 6L)
sapply(weightsList, summary)


###################################################
### code chunk number 20: addPrimingWeights
###################################################
yeastAligned2 <- mapply(addPrimingWeights, yeastAligned, weightsList)
alignData(yeastAligned2[[1]])
head(alignData(yeastAligned2[[1]])$weights)


###################################################
### code chunk number 21: importWithWeights
###################################################
eData2 <- importFromAlignedReads(yeastAligned2, chrMap = chrMap,
                                 dbFilename = "my.db", tablename = "weights")


###################################################
### code chunk number 22: reweightedCounts
###################################################
reweightedCounts <- summarizeByAnnotation(eData2, yAnno, ignoreStrand = TRUE, 
                                          groupBy = "ensembl_gene_id")
head(reweightedCounts)


###################################################
### code chunk number 23: sessionInfo
###################################################
toLatex(sessionInfo())