### 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))