### R code from vignette source 'easyRNASeq.Rnw'
###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()
###################################################
### code chunk number 2: Load the library and create the object
###################################################
library(easyRNASeq)
library(RnaSeqTutorial)
count.table <- easyRNASeq(filesDirectory=system.file(
"extdata",
package="RnaSeqTutorial"),
pattern="[A,C,T,G]{6}\\.bam$",
readLength=30L,
organism="Dmelanogaster",
annotationMethod="rda",
annotationFile=system.file(
"data",
"gAnnot.rda",
package="RnaSeqTutorial"),
count="exons"
)
###################################################
### code chunk number 3: Show the results
###################################################
head(count.table)
dim(count.table)
###################################################
### code chunk number 4: Create the object with only 2 files (eval = FALSE)
###################################################
## count.table <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=30L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## count="exons",
## filenames=c("ACACTG.bam", "ACTAGC.bam"))
###################################################
### code chunk number 5: read a gapped bam file (eval = FALSE)
###################################################
## count.table <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=36L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## gapped=TRUE,
## count="exons",
## filenames="gapped.bam")
###################################################
### code chunk number 6: read an export (eval = FALSE)
###################################################
## 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="aln",
## type="SolexaExport",
## count="exons",
## pattern="subset_export",
## filter=compose(
## chromosomeFilter(regex="chr"),
## chastityFilter()))
###################################################
### code chunk number 7: GFF3 (eval = FALSE)
###################################################
## system.file(
## "extdata",
## "Dmel-mRNA-exon-r5.52.gff3",
## package="RnaSeqTutorial")
###################################################
### code chunk number 8: load annotation from biomaRt (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=36L,
## annotationMethod="biomaRt",
## gapped=TRUE,
## count="exons",
## filenames="gapped.bam",
## outputFormat="RNAseq")
###################################################
### code chunk number 9: access the annotation (eval = FALSE)
###################################################
## gAnnot <- genomicAnnotation(rnaSeq)
###################################################
### code chunk number 10: load the annotation from a gff (eval = FALSE)
###################################################
## count.table <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=36L,
## annotationMethod="gff",
## annotationFile=system.file(
## "extdata",
## "Dmel-mRNA-exon-r5.52.gff3",
## package="RnaSeqTutorial"),
## gapped=TRUE,
## count="exons",
## filenames="gapped.bam")
###################################################
### code chunk number 11: RangedData structure
###################################################
library(IRanges)
gAnnot <- RangedData(
IRanges(
start=c(10,30,100),
end=c(21,53,123)),
space=c("chr01","chr01","chr02"),
strand=c("+","+","-"),
transcript=c("trA1","trA2","trB"),
gene=c("gA","gA","gB"),
exon=c("e1","e2","e3"),
universe = "Hs19"
)
gAnnot
###################################################
### code chunk number 12: look at the gAnnot (eval = FALSE)
###################################################
## library(RnaSeqTutorial)
## data(gAnnot)
## gAnnot
###################################################
### code chunk number 13: GRangesList coercion
###################################################
grngs <- as(gAnnot,"GRanges")
grngsList<-split(grngs,seqnames(grngs))
grngsList
###################################################
### code chunk number 14: easyRNASeq output
###################################################
rnaSeq <- easyRNASeq(system.file(
"extdata",
package="RnaSeqTutorial"),
organism="Dmelanogaster",
readLength=30L,
annotationMethod="rda",
annotationFile=system.file(
"data",
"gAnnot.rda",
package="RnaSeqTutorial"),
count="exons",
pattern="[A,C,T,G]{6}\\.bam$",
outputFormat="RNAseq")
###################################################
### code chunk number 15: Transcript counts
###################################################
rnaSeq <- transcriptCounts(rnaSeq)
head(readCounts(rnaSeq,'transcripts'))
###################################################
### code chunk number 16: RPKM normalization (eval = FALSE)
###################################################
## count.table <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=30L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## count="exons",
## filenames=c("ACACTG.bam", "ACTAGC.bam",
## "ATGGCT.bam", "TTGCGA.bam"),
## normalize=TRUE
## )
###################################################
### code chunk number 17: RPKM counts
###################################################
feature.size = width(genomicAnnotation(rnaSeq))
names(feature.size) = genomicAnnotation(rnaSeq)$exon
lib.size=c(
"ACACTG.bam"=56643,
"ACTAGC.bam"=42698,
"ATGGCT.bam"=55414,
"TTGCGA.bam"=60740)
head(RPKM(readCounts(rnaSeq,summarization="exons")$exons,
NULL,
lib.size=lib.size,
feature.size=feature.size))
###################################################
### code chunk number 18: RPKM on RNAseq
###################################################
head(RPKM(rnaSeq,from="transcripts"))
###################################################
### code chunk number 19: define the conditions
###################################################
conditions=c("A","A","B","B")
names(conditions) <- c("ACACTG.bam", "ACTAGC.bam",
"ATGGCT.bam", "TTGCGA.bam")
###################################################
### code chunk number 20: DESeq normalization (eval = FALSE)
###################################################
## countDataSet <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=30L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## count="exons",
## filenames=c("ACACTG.bam", "ACTAGC.bam",
## "ATGGCT.bam", "TTGCGA.bam"),
## normalize=TRUE,
## outputFormat="DESeq",
## conditions=conditions,
## fitType="local"
## )
###################################################
### code chunk number 21: edgeR normalization (eval = FALSE)
###################################################
## dgeList <- easyRNASeq(system.file(
## "extdata",
## package="RnaSeqTutorial"),
## organism="Dmelanogaster",
## readLength=30L,
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## count="exons",
## filenames=c("ACACTG.bam", "ACTAGC.bam",
## "ATGGCT.bam", "TTGCGA.bam"),
## normalize=TRUE,
## outputFormat="edgeR",
## conditions=conditions
## )
###################################################
### code chunk number 22: cleanup
###################################################
file2clean=c(
"ACACTG.fastq",
"ACTAGC.fastq",
"ATGGCT.fastq",
"TTGCGA.fastq")
silent <- sapply(
file2clean,
function(f2c){
if(file.exists(f2c)){file.remove(f2c)}
})
###################################################
### code chunk number 23: Demultiplexing sample
###################################################
alns <- readAligned(
system.file(
"extdata",
package="RnaSeqTutorial"),
pattern="multiplex_export",
filter=compose(
chastityFilter(),
nFilter(2),
chromosomeFilter(regex="chr")),
type="SolexaExport",
withAll=TRUE)
###################################################
### code chunk number 24: Illumina example (eval = FALSE)
###################################################
## alignData(alns)$multiplexIndex
###################################################
### code chunk number 25: barcodes
###################################################
barcodes=c("ACACTG","ACTAGC","ATGGCT","TTGCGA")
###################################################
### code chunk number 26: barcodePlot
###################################################
barcodePlot(alns,
barcodes=barcodes,
type="within",
barcode.length=6,
show.barcode=20,
main="All samples",
xlim=c(0,0.5))
###################################################
### code chunk number 27: demultiplex
###################################################
dem.alns <- demultiplex(alns,
barcodes=barcodes,
edition.dist=2,
barcodes.qty=4,
type="within")
###################################################
### code chunk number 28: check the objects (eval = FALSE)
###################################################
## dem.alns$reads[[1]]
## dem.alns$barcodes[[1]]
###################################################
### code chunk number 29: demPlot
###################################################
par(mfrow=c(2,2))
barcode.frequencies <- lapply(
names(dem.alns$barcodes),
function(barcode,alns){
barcodePlot(
alns$barcodes[[barcode]],
barcodes=barcode,
type="within",barcode.length=6,
show.barcode=20,
main=paste(
"Expected barcode:",
barcode))
},dem.alns)
###################################################
### code chunk number 30: write the file
###################################################
status <- lapply(
names(dem.alns$barcodes),
function(barcode,alns){
writeFastq(
alns$reads[[barcode]],
file=paste(
barcode,
"fastq",sep="."))
},dem.alns)
###################################################
### code chunk number 31: cleanup
###################################################
file2clean=c(
"ACACTG.fastq",
"ACTAGC.fastq",
"ATGGCT.fastq",
"TTGCGA.fastq",
"Rplots.pdf")
silent <- sapply(
file2clean,
function(f2c){
if(file.exists(f2c)){file.remove(f2c)}
})
###################################################
### code chunk number 32: count function example (eval = FALSE)
###################################################
## ## creating a RangedSummarizedExperiment from 4 bam files
## sumExp <- easyRNASeq(filesDirectory=system.file(
## "extdata",
## package="RnaSeqTutorial"),
## pattern="[A,C,T,G]{6}\\.bam$",
## readLength=30L,
## organism="Dmelanogaster",
## annotationMethod="rda",
## annotationFile=system.file(
## "data",
## "gAnnot.rda",
## package="RnaSeqTutorial"),
## count="exons",
## outputFormat="SummarizedExperiment"
## )
## ## the counts
## assays(sumExp)
## ## the sample info
## colData(sumExp)
## ## the 'features' info
## rowRanges(sumExp)
###################################################
### code chunk number 33: load the lib (eval = FALSE)
###################################################
## library(easyRNASeq)
###################################################
### code chunk number 34: synthetic-transcript (eval = FALSE)
###################################################
## ## lib - loaded by easyRNASeq already
## library(genomeIntervals)
##
## ## read in the gff
## gff <- readGff3("dmel_r5-54.gff3")
##
## ## look at the gff
## gff
## nrow(gff[gff$type=="exon",])
## nrow(gff[gff$type=="mRNA",])
## nrow(gff[gff$type=="gene",])
##
## ## map the mRNA ID to the gene ID
## transcriptGeneMapping <- data.frame(
## getGffAttribute(gff[gff$type=="mRNA",],"ID"),
## getGffAttribute(gff[gff$type=="mRNA",],"Parent"))
## head(transcriptGeneMapping)
##
## ## extract the exons IRangesList
## ## ie exons are grouped by gene
## ## and their coordinates are stored as an IRanges object
## sel <- gff$type=="exon"
## rngList<- split(IRanges(start=gff[sel,1],end=gff[sel,2]),
## transcriptGeneMapping[match(
## sapply(strsplit(
## getGffAttribute(gff[sel,],"Parent"),
## ","),"[",1),
## transcriptGeneMapping$ID),"Parent"])
## rngList
##
## ## check what's the gene with the max number of exons
## mostExons <- rev(names(table(elementNROWS(rngList))))[1]
## mostExons
##
## ## work the magic; collapse the genes IRanges
## rngList<- reduce(rngList)
## rngList
##
## ## what's the max now?
## rev(names(table(elementNROWS(rngList))))[1]
##
## ## create the gff
## ## get the gff information; here we simply duplicate the
## ## first exon of every gene by the number of synthetic
## ## exons per gene. The content will be updated afterwards.
## exons <- gff[sel,]
## syntheticGeneModel<- exons[rep(
## match(names(rngList),
## transcriptGeneMapping[
## match(sapply(
## strsplit(getGffAttribute(exons,"Parent"),
## ","),"[",1),
## transcriptGeneMapping$ID),"Parent"]),
## elementNROWS(rngList)),]
##
## ## update the coordinates
## syntheticGeneModel[,1]<- unlist(start(rngList))
## syntheticGeneModel[,2]<- unlist(end(rngList))
##
## ## change the source
## levels(syntheticGeneModel$source)<- "inhouse"
##
## ## get the exon number for the minus strand
## exonNumber<- lapply(elementNROWS(rngList),":",1)
##
## ## reverse them on the plus strand
## sel<- strand(syntheticGeneModel)[cumsum(elementNROWS(rngList))] == "+"
## exonNumber[sel]<- sapply(exonNumber[sel],rev)
##
## ## update the attributes
## syntheticGeneModel$gffAttributes<- paste("ID=",
## rep(names(rngList),elementNROWS(rngList)),
## ":",unlist(exonNumber),";Parent=",
## rep(names(rngList),elementNROWS(rngList)),".0",sep="")
##
## ## write the file
## writeGff3(syntheticGeneModel,file="dmel_synthetic_transcript_r5-54.gff3")
###################################################
### code chunk number 35: easyTrx (eval = FALSE)
###################################################
## sumExp <- easyRNASeq(
## filesDirectory=system.file(
## "extdata",
## package="RnaSeqTutorial"),
## pattern="[A,C,T,G]{6}\\.bam$",
## readLength=30L,
## chr.sel="chr2L",
## organism="Dmelanogaster",
## annotationMethod="gff",
## annotationFile="dmel_synthetic_transcript_r5-54.gff3",
## count="transcripts",
## outputFormat="SummarizedExperiment"
## )
###################################################
### code chunk number 36: Access the overlap flag (eval = FALSE)
###################################################
## genomicAnnotation(rnaSeq)$overlap
###################################################
### code chunk number 37: Access the overlap flag from an SE (eval = FALSE)
###################################################
## rowRanges(sumExp)$overlap
###################################################
### code chunk number 38: bam prep (eval = FALSE)
###################################################
## library(BSgenome.Hsapiens.UCSC.hg19)
## library(GEOquery)
## library(SRAdb)
## library(Rsamtools)
## library(Rsubread)
##
## ## create a temp dir
## dir.create(file.path(getwd(),"tmp1234"))
##
## ## change the working directory
## setwd(file.path(getwd(),"tmp1234"))
##
## ## init SRA
## sqlfile <- getSRAdbFile()
##
## ## init a connection
## sra_con <- dbConnect(SQLite(),sqlfile)
##
## ## list the files
## acc <- c("SRR490224","SRR490225")
## getFASTQinfo( in_acc = acc, srcType = 'ftp' )
##
## ## get the read files
## getSRAfile( in_acc=acc, sra_con, destDir = getwd(),
## fileType = 'fastq', srcType = 'ftp' )
##
## ## close the connection
## dbDisconnect( sra_con )
##
## ## write the human genome sequences
## writeXStringSet(Reduce(append,
## lapply(seqnames(Hsapiens),
## function(nam){
## dss <- DNAStringSet(unmasked(Hsapiens[[nam]]))
## names(dss) <- nam
## dss
## })),
## file="hg19.fa")
##
## ## create the indexes
## dir.create("indexes")
## buildindex(basename=file.path("indexes","hg19"),
## reference="hg19.fa")
##
## ## align the reads
## sapply(dir(pattern="*\\.gz$"),function(fil){
## ## decompress the files
## gunzip(fil)
##
## ## align
## align(index=file.path("indexes","hg19"),
## readfile1=sub("\\.gz$","",fil),
## nsubreads=2,TH1=1,
## output_file=sub("\\.fastq\\.gz$","\\.sam",fil)
## )
##
## ## create bam files
## asBam(file=sub("\\.fastq\\.gz$","\\.sam",fil),
## destination=sub("\\.fastq\\.gz$","",fil),
## indexDestination=TRUE)
## })
##
###################################################
### code chunk number 39: chromosome size (eval = FALSE)
###################################################
## library(BSgenome.Hsapiens.UCSC.hg19)
## chr.sizes=seqlengths(Hsapiens)
###################################################
### code chunk number 40: bam files name (eval = FALSE)
###################################################
## bamfiles=dir(getwd(),pattern="*\\.bam$")
###################################################
### code chunk number 41: Fetching with biomaRt (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(filesDirectory=getwd(),
## organism="Hsapiens",
## readLength=58L,
## annotationMethod="biomaRt",
## count="exons",
## filenames=bamfiles[1],
## outputFormat="RNAseq"
## )
## gAnnot <- genomicAnnotation(rnaSeq)
###################################################
### code chunk number 42: filtering the annot (eval = FALSE)
###################################################
## gAnnot <- gAnnot[space(gAnnot) %in% paste("chr",c(1:22,"X","Y","M"),sep=""),]
## save(gAnnot,file="gAnnot.rda")
###################################################
### code chunk number 43: get a count table (eval = FALSE)
###################################################
## countTable <- easyRNASeq(filesDirectory=getwd(),
## organism="Hsapiens",
## readLength=58L,
## annotationMethod="rda",
## annotationFile="gAnnot.rda",
## count="exons",
## filenames=bamfiles[1]
## )
###################################################
### code chunk number 44: using different annotation (eval = FALSE)
###################################################
## library(GenomicFeatures)
## hg19.tx <- makeTxDbFromUCSC(genome="hg19", tablename="refGene")
##
## gAnnot <- exons(hg19.tx)
##
## colnames(elementMetadata(gAnnot)) <- "exon"
##
## gAnnot <- split(gAnnot,seqnames(gAnnot))
###################################################
### code chunk number 45: sub setting (eval = FALSE)
###################################################
## countTable <- easyRNASeq(filesDirectory=getwd(),
## organism="Hsapiens",
## readLength=58L,
## annotationMethod="env",
## annotationObject=gAnnot,
## count="exons",
## filenames=bamfiles[1],
## chr.sel="chr1"
## )
###################################################
### code chunk number 46: using the count function (eval = FALSE)
###################################################
## sumExp <- easyRNASeq(filesDirectory=getwd(),
## organism="Hsapiens",
## readLength=58L,
## annotationMethod="env",
## annotationObject=gAnnot,
## count="exons",
## filenames=bamfiles[1],
## chr.sel="chr1",
## outputFormat="SummarizedExperiment"
## )
###################################################
### code chunk number 47: DESeq run (eval = FALSE)
###################################################
## conditions <- c("A","B")
## names(conditions) <- bamfiles
##
## countDataSet <- easyRNASeq(filesDirectory=getwd(),
## organism="Hsapiens",
## annotationMethod="env",
## annotationObject=gAnnot,
## count="exons",
## filenames=bamfiles,
## chr.sel="chr1",
## outputFormat="DESeq",
## normalize=TRUE,
## conditions=conditions,
## fitType="local",
## method="blind"
## )
###################################################
### code chunk number 48: read the data (eval = FALSE)
###################################################
## aln <- readAligned("data",type="SolexaExport",pattern="*.txt.gz")
## gc()
## levels(chromosome(aln))
###################################################
### code chunk number 49: fake the data
###################################################
c("chr1.fa","chr10.fa","...","chrY.fa")
###################################################
### code chunk number 50: fetch the annotation (eval = FALSE)
###################################################
## obj <- fetchAnnotation(new('RNAseq',
## organismName="Mmusculus"
## ),
## method="biomaRt")
##
## gAnnot <- genomicAnnotation(obj)
##
## length(grep("NT_",space(gAnnot)))
##
###################################################
### code chunk number 51: fake number 2
###################################################
1181
###################################################
### code chunk number 52: modify the annotation (eval = FALSE)
###################################################
## names(gAnnot) <- paste("chr",names(gAnnot),".fa",sep="")
###################################################
### code chunk number 53: clean (eval = FALSE)
###################################################
## rm(aln,obj)
## gc()
###################################################
### code chunk number 54: chrsize (eval = FALSE)
###################################################
## library(BSgenome.Mmusculus.UCSC.mm9)
## chr.sizes<- seqlengths(Mmusculus)
###################################################
### code chunk number 55: create the chr.map (eval = FALSE)
###################################################
## chr.map <- data.frame(
## from=paste("chr",c(1:19,"X","Y"),".fa",sep=""),
## to=paste("chr",c(1:19,"X","Y"),sep=""))
###################################################
### code chunk number 56: no filter (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(filesDirectory="data",
## organism="custom",
## chr.map=chr.map,
## chr.sizes=chr.sizes,
## filter=compose(
## naPositionFilter(),
## chastityFilter()),
## readLength=50L,
## annotationMethod="env",
## annotationObject=gAnnot,
## format="aln",
## count="genes",
## summarization= "geneModels",
## filenames="1-Feb_ATCACG_L003_R1_001_export.txt.gz",
## outputFormat="RNAseq",
## nbCore=2
## )
###################################################
### code chunk number 57: selected chr (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(filesDirectory="data",
## organism="custom",
## chr.map=chr.map,
## chr.sizes=chr.sizes,
## chr.sel=chr.map$from,
## filter=compose(
## naPositionFilter(),
## chastityFilter()),
## readLength=50L,
## annotationMethod="env",
## annotationObject=gAnnot,
## format="aln",
## count="genes",
## summarization= "geneModels",
## filenames="1-Feb_ATCACG_L003_R1_001_export.txt.gz",
## outputFormat="RNAseq",
## nbCore=2
## )
###################################################
### code chunk number 58: remove annotation (eval = FALSE)
###################################################
## sel <- grep("NT_",names(gAnnot))
## gAnnot <- RangedData(ranges=ranges(gAnnot)[-sel,],values=values(gAnnot)[-sel,])
## colnames(gAnnot) <- gsub("values\\.","",colnames(gAnnot))
###################################################
### code chunk number 59: final call (eval = FALSE)
###################################################
## rnaSeq <- easyRNASeq(filesDirectory="data",
## organism="custom",
## chr.map=chr.map,
## chr.sizes=chr.sizes,
## chr.sel=chr.map$from,
## filter=compose(
## naPositionFilter(),
## chastityFilter()),
## readLength=50L,
## annotationMethod="env",
## annotationObject=gAnnot,
## format="aln",
## count="genes",
## summarization= "geneModels",
## filenames="1-Feb_ATCACG_L003_R1_001_export.txt.gz",
## outputFormat="RNAseq",
## nbCore=2
## )
###################################################
### code chunk number 60: SessionInfo
###################################################
library(easyRNASeq)
library(BSgenome.Dmelanogaster.UCSC.dm3)
library(RnaSeqTutorial)
sessionInfo()
###################################################
### code chunk number 61: some test
###################################################
grngs = GRanges(seqnames=c("chr1", "chr2"),
ranges=IRanges(start=1:2, end=2:3),
strand=c("+","-"))
silent <- strand(grngs)
silent <- reduce(grngs)