### R code from vignette source 'RnaSeqTutorial.Rnw'
###################################################
### code chunk number 1: Lood the tutorial library
###################################################
library(RnaSeqTutorial)
###################################################
### code chunk number 2: Read the data
###################################################
library(ShortRead)
aln<-readAligned(
system.file("extdata",package="RnaSeqTutorial"),
pattern="subset_export",type="SolexaExport")
show(aln)
###################################################
### code chunk number 3: The AlignedRead object (eval = FALSE)
###################################################
## chromosome(aln)
## levels(chromosome(aln))
## position(aln)[1:100]
## width(aln)[1:100]
## strand(aln)[1:100]
## sread(aln)
## quality(aln)
###################################################
### code chunk number 4: The chastity field (eval = FALSE)
###################################################
## alignData(aln)$filtering[1:100]
###################################################
### code chunk number 5: Indexing a BAM file with ShortRead
###################################################
library(Rsamtools)
file.copy(
system.file("extdata",
"subset.bam",
package="RnaSeqTutorial"),
getwd())
indexFile <- indexBam("subset.bam")
basename(indexFile)
###################################################
### code chunk number 6: Loading a file using Rsamtools
###################################################
aln2b <- scanBam(
"subset.bam",
index="subset.bam"
)
names(aln2b[[1]])
###################################################
### code chunk number 7: Reading a Tophat BAM file
###################################################
library(GenomicAlignments)
aln3 <- readGAlignments(
system.file("extdata",
"gapped.bam",
package="RnaSeqTutorial"))
###################################################
### code chunk number 8: missing ids
###################################################
head(id(aln))
###################################################
### code chunk number 9: with ids
###################################################
aln4<-readAligned(
system.file("extdata",package="RnaSeqTutorial"),
pattern="subset_export",type="SolexaExport",
withId=TRUE)
head(id(aln4))
###################################################
### code chunk number 10: Filtering the reads (eval = FALSE)
###################################################
## library(easyRNASeq)
## nFilt <- nFilter(2)
## chrFilt <- chromosomeFilter(regex="chr")
## cFilt <- chastityFilter()
## filt <- compose(nFilt,chrFilt,cFilt)
## aln <- aln[filt(aln)]
###################################################
### code chunk number 11: Really filter but do not show
###################################################
nFilt <- nFilter(2)
chrFilt <- chromosomeFilter(regex="chr")
filt <- compose(nFilt,chrFilt)
aln <- aln[filt(aln)]
aln <- aln[alignData(aln)$filtering=="Y"]
###################################################
### code chunk number 12: the filtered aln (eval = FALSE)
###################################################
## show(aln)
###################################################
### code chunk number 13: cleanup
###################################################
rm(aln2b,aln3,aln4)
silent<-gc(verbose=FALSE)
###################################################
### code chunk number 14: Getting the chromosome names and sizes
###################################################
library(BSgenome.Dmelanogaster.UCSC.dm3)
chrSizes <-seqlengths(Dmelanogaster)
chrSizes
###################################################
### code chunk number 15: Exon Transcript Annotation
###################################################
library(biomaRt)
ensembl <- useMart("ensembl",
dataset="dmelanogaster_gene_ensembl")
exon.annotation<-getBM(
c("ensembl_gene_id",
"strand",
"ensembl_transcript_id",
"chromosome_name",
"ensembl_exon_id",
"exon_chrom_start",
"exon_chrom_end"),
mart=ensembl,
filters="chromosome_name",
values=c("2L","2R","3L","3R","4","X"))
###################################################
### code chunk number 16: Exon Transcript object transformation
###################################################
exon.annotation$chromosome <- paste(
"chr",
exon.annotation$chromosome_name,
sep="")
exon.range <- RangedData(
IRanges(
start=exon.annotation$exon_chrom_start,
end=exon.annotation$exon_chrom_end),
space=exon.annotation$chromosome,
strand=exon.annotation$strand,
transcript=exon.annotation$ensembl_transcript_id,
gene=exon.annotation$ensembl_gene_id,
exon=exon.annotation$ensembl_exon_id,
universe = "Dm3"
)
###################################################
### code chunk number 17: Check the objects (eval = FALSE)
###################################################
## show(exon.range)
###################################################
### code chunk number 18: Read in the gff
###################################################
library(genomeIntervals)
gInterval<-readGff3(system.file("extdata",
"annot.gff",
package="RnaSeqTutorial"),
quiet=TRUE)
###################################################
### code chunk number 19: post process the gff
###################################################
exon.range2 <- RangedData(
IRanges(
start=gInterval[,1],
end=gInterval[,2]),
space=gInterval$seq_name,
strand=gInterval$strand,
transcript=as.vector(
getGffAttribute(gInterval,"Parent")),
gene=as.vector(
getGffAttribute(gInterval,"Name")),
exon=as.vector(
getGffAttribute(gInterval,"ID")),
universe = "Dm3"
)
###################################################
### code chunk number 20: importing with rtracklayer (eval = FALSE)
###################################################
## library(rtracklayer)
## exon.range3<-import.gff3(
## system.file("extdata",
## "annot.gff",
## package="RnaSeqTutorial")
## )
###################################################
### code chunk number 21: importing with GenomicFeatures (eval = FALSE)
###################################################
## library(GenomicFeatures)
###################################################
### code chunk number 22: GenomicFeature taste (eval = FALSE)
###################################################
## dm3.tx <- makeTxDbFromUCSC(
## genome="dm3",
## tablename="refGene")
## exon.range4 <- exons(dm3.tx)
## exon.range4
###################################################
### code chunk number 23: clean
###################################################
rm(exon.range,gInterval,ensembl,exon.annotation)
silent<-gc(verbose=FALSE)
###################################################
### code chunk number 24: Genomic coverage
###################################################
cover <- coverage(aln,width=chrSizes)
###################################################
### code chunk number 25: Have a look at the coverage (eval = FALSE)
###################################################
## show(cover)
## show(cover$chr2R)
###################################################
### code chunk number 26: Additional commands for cov
###################################################
runLength(cover$chr4)[1:3]
runValue(cover$chr4)[1:3]
r.start <- runLength(cover$chr4)[1]+1
r.end <- sum(runLength(cover$chr4)[1:2])
as.integer(cover$chr4)[r.start:r.end]
as.integer(window(cover$chr4,r.start,r.end))
###################################################
### code chunk number 27: Exon Transcript coverage
###################################################
exon.coverage<-aggregate(
cover[match(names(exon.range2),names(cover))],
ranges(exon.range2),
sum)
exon.coverage <- ceiling(unlist(exon.coverage)/unique(width(aln)))
names(exon.coverage) <- exon.range2$exon
###################################################
### code chunk number 28: Ex. tr. cov.p.2 (eval = FALSE)
###################################################
## show(exon.coverage)
###################################################
### code chunk number 29: faster (eval = FALSE)
###################################################
## viewSums(
## Views(
## cover[match(names(exon.range2),names(cover))],
## ranges(exon.range2)))
###################################################
### code chunk number 30: clean up
###################################################
rm(aln,chrSizes)
silent<-gc(verbose=FALSE)
###################################################
### code chunk number 31: Load the library and create the object (eval = FALSE)
###################################################
## ## load the library
## library("easyRNASeq")
## library(BSgenome.Dmelanogaster.UCSC.dm3)
##
## count.table <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## chr.sizes=seqlengths(Dmelanogaster),
## readLength=36L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## format="bam",
## count="exons",
## pattern="[A,C,T,G]{6}\\.bam$")
###################################################
### code chunk number 32: Show (eval = FALSE)
###################################################
## head(count.table)
## dim(count.table)
###################################################
### code chunk number 33: open the vignette (eval = FALSE)
###################################################
## vignette("easyRNASeq")
###################################################
### code chunk number 34: RPKM counts (eval = FALSE)
###################################################
## feature.size = width(exon.range2)
## names(feature.size) = exon.range2$exon
## feature.size <- feature.size[!duplicated(names(feature.size))]
## lib.size=c("ACACTG.bam"=56643,
## "ACTAGC.bam"=42698,
## "ATGGCT.bam"=55414,
## "TTGCGA.bam"=60740)
## head(RPKM(count.table,NULL,
## lib.size=lib.size,
## feature.size=feature.size))
###################################################
### code chunk number 35: Chromosome 4 wig export (eval = FALSE)
###################################################
## library(rtracklayer)
## window.size <- 51
## rngs <- breakInChunks(length(cover[["chr4"]]),window.size)
## vals <- viewSums(Views(cover[["chr4"]],rngs))
## #width(rngs)[width(rngs) != width(rngs)[1]] <- width(rngs)[1]
## silent <- export(
## RangedData(rngs,score=vals,universe="Dmelanogaster",space="chr4"),
## con="chr4.wig"
## )
###################################################
### code chunk number 36: normalized exon bed export (eval = FALSE)
###################################################
## exon.RPKM <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## chr.sizes=seqlengths(Dmelanogaster),
## readLength=36L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## format="aln",
## count="exons",
## normalize=TRUE,
## pattern="subset_export",
## type="SolexaExport",
## filter=compose(
## chastityFilter(),
## nFilter(2),
## chromosomeFilter(regex="chr")))
##
## exons <- exon.range2
## exons <- exons[!duplicated(exons$exon),]
## exons$score <- exon.RPKM[,1]
## exons$name <- rownames(exon.RPKM)
## exons <- exons[exons$score>0,]
## export(exons,con="exons.bed")
###################################################
### code chunk number 37: final.cleanup
###################################################
file2clean=c(
"chr4.wig",
"Rplots.pdf",
"subset.bam",
"subset.bam.bai")
silent <- sapply(
file2clean,
function(f2c){
if(file.exists(f2c)){file.remove(f2c)}
})
###################################################
### code chunk number 38: SessionInfo
###################################################
##library(rtracklayer)
library(BSgenome.Dmelanogaster.UCSC.dm3)
library(easyRNASeq)
sessionInfo()
###################################################
### code chunk number 39: solution 3 (eval = FALSE)
###################################################
## table(ngap(aln3))
## table(cigarOpTable(cigar(aln3))[,"N"])
###################################################
### code chunk number 40: solution 4 (eval = FALSE)
###################################################
## exon.grange <- GRanges(
## IRanges(
## start=exon.annotation$exon_chrom_start,
## end=exon.annotation$exon_chrom_end),
## seqnames=Rle(exon.annotation$chromosome),
## strand=Rle(exon.annotation$strand),
## transcript=exon.annotation$ensembl_transcript_id,
## gene=exon.annotation$ensembl_gene_id,
## exon=exon.annotation$ensembl_exon_id,
## seqlengths = chrSizes[
## match(
## unique(exon.annotation$chromosome),
## names(chrSizes))]
## )
###################################################
### code chunk number 41: solution 5 (eval = FALSE)
###################################################
## levels(gInterval$strand) <- c("-","+")
## exon.grange2 <- GRanges(
## IRanges(
## start=gInterval[,1],
## end=gInterval[,2]),
## seqnames=Rle(gInterval$seq_name),
## strand=Rle(gInterval$strand),
## transcript=as.vector(getGffAttribute(
## gInterval,"Parent")),
## gene=as.vector(getGffAttribute(
## gInterval,"Name")),
## exon=as.vector(getGffAttribute(
## gInterval,"ID")),
## seqlengths = chrSizes[
## match(
## unique(gInterval$seq_name),
## names(chrSizes))]
## )
###################################################
### code chunk number 42: solution 5b (eval = FALSE)
###################################################
## exon.grange2 <- as(gInterval,"GRanges")
###################################################
### code chunk number 43: solution 5c (eval = FALSE)
###################################################
## as(as(gInterval,"GRanges"),"RangedData")
###################################################
### code chunk number 44: solution 6 (eval = FALSE)
###################################################
## library(rtracklayer)
## load(system.file("data",
## "gAnnot.rda",
## package="RnaSeqTutorial"))
## export.gff3(gAnnot,con="annot.gff")
## export.gff(gAnnot,con="annot.gff",version="3")
## export(gAnnot,con="annot.gff",version="3")
###################################################
### code chunk number 45: solution7 (eval = FALSE)
###################################################
## as.integer(cover$chr4)
###################################################
### code chunk number 46: solution8 (eval = FALSE)
###################################################
## runLength(cover$chr4)
## runValue(cover$chr4)
###################################################
### code chunk number 47: solution9 better safe than sorry (eval = FALSE)
###################################################
## names(exon.range2)
## names(cover)
## match(names(exon.range2),names(cover))
###################################################
### code chunk number 48: solution10 (eval = FALSE)
###################################################
## exon.coverage <- unlist(exon.coverage)
## names(exon.coverage) <- exon.range2$exon
## head(
## sort(
## exon.coverage[!duplicated(names(exon.coverage))],
## decreasing=TRUE))
###################################################
### code chunk number 49: solution11 (eval = FALSE)
###################################################
## sel <- chromosome(aln) != "chrM"
## aln <- aln[sel]
## exon.counts <- countOverlaps(
## exon.range2,
## split(IRanges(
## start=position(aln),
## width=width(aln)),
## chromosome(aln))
## )
###################################################
### code chunk number 50: solution11.2 (eval = FALSE)
###################################################
##
## plot(
## unlist(exon.coverage),
## unlist(exon.counts),
## log="xy",
## main="countOverlap vs. aggregate",
## xlab="aggregate",
## ylab="CountOverlap",
## pch="+",col=6)
## abline(0,1,lty=2,col="orange")
##
## table(unlist(exon.coverage) - unlist(exon.counts))
###################################################
### code chunk number 51: solution12 (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## chr.sizes=as.list(seqlengths(Dmelanogaster)),
## readLength=36L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## format="bam",
## count="exons",
## pattern="bam$",
## outputFormat="RNAseq")
## show(rnaSeq)
###################################################
### code chunk number 52: Transcript counts (eval = FALSE)
###################################################
## rnaSeq <- transcriptCounts(rnaSeq)
## head(readCounts(rnaSeq,'transcripts'))
###################################################
### code chunk number 53: solution13 (eval = FALSE)
###################################################
## rnaSeq<-geneCounts(rnaSeq,summarization='geneModels')
###################################################
### code chunk number 54: solution13B (eval = FALSE)
###################################################
## rnaSeq2 <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## chr.sizes=as.list(seqlengths(Dmelanogaster)),
## readLength=36L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## format="bam",
## count="genes",
## summarization="geneModels",
## pattern="bam$",
## outputFormat="RNAseq")
## head(readCounts(rnaSeq2,'genes','geneModels'))
###################################################
### code chunk number 55: solution14 (eval = FALSE)
###################################################
## RPKM(rnaSeq,from="transcripts")
## RPKM(rnaSeq2,from="geneModels")
###################################################
### code chunk number 56: solution17 (eval = FALSE)
###################################################
## alignData(alns)$multiplexIndex
###################################################
### code chunk number 57: solution17b (eval = FALSE)
###################################################
## varLabels(alignData(alns))