### R code from vignette source 'vignettes/QuasR/inst/doc/QuasR-Overview.Rnw' ################################################### ### code chunk number 1: options ################################################### options(width=65) options('useFancyQuotes' = FALSE) ################################################### ### code chunk number 2: cite ################################################### citation("QuasR") ################################################### ### code chunk number 3: install (eval = FALSE) ################################################### ## source("http://www.bioconductor.org/biocLite.R") ## biocLite("QuasR") ################################################### ### code chunk number 4: loadLibraries ################################################### library(QuasR) library(BSgenome) library(Rsamtools) library(rtracklayer) library(GenomicFeatures) library(Gviz) ################################################### ### code chunk number 5: help1 (eval = FALSE) ################################################### ## help.start() ################################################### ### code chunk number 6: loadQuasRLibrary (eval = FALSE) ################################################### ## library(QuasR) ################################################### ### code chunk number 7: help2 (eval = FALSE) ################################################### ## ?preprocessReads ################################################### ### code chunk number 8: help3 (eval = FALSE) ################################################### ## help("preprocessReads") ################################################### ### code chunk number 9: assign (eval = FALSE) ################################################### ## x <- 2 ################################################### ### code chunk number 10: ls (eval = FALSE) ################################################### ## ls() ################################################### ### code chunk number 11: printObject (eval = FALSE) ################################################### ## x ################################################### ### code chunk number 12: SampleSession1 ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 13: SampleSession2 ################################################### sampleFile <- "extdata/samples_chip_single.txt" genomeFile <- "extdata/hg19sub.fa" proj <- qAlign(sampleFile, genomeFile) proj ################################################### ### code chunk number 14: SampleSession3 ################################################### qQCReport(proj, "extdata/qc_report.pdf") ################################################### ### code chunk number 15: SampleSession4 ################################################### library(rtracklayer) library(GenomicFeatures) annotFile <- "extdata/hg19sub_annotation.gtf" txStart <- import.gff(annotFile, format="gtf", asRangedData=FALSE, feature.type="start_codon") promReg <- promoters(txStart, upstream=500, downstream=500) names(promReg) <- mcols(promReg)$transcript_name promCounts <- qCount(proj, query=promReg) promCounts ################################################### ### code chunk number 16: sampleFile (eval = FALSE) ################################################### ## sampleFile1 <- system.file(package="QuasR", "extdata", ## "samples_chip_single.txt") ## sampleFile2 <- system.file(package="QuasR", "extdata", ## "samples_rna_paired.txt") ################################################### ### code chunk number 17: auxiliaryFile ################################################### auxFile <- system.file(package="QuasR", "extdata", "auxiliaries.txt") ################################################### ### code chunk number 18: selectGenomeBSgenome ################################################### available.genomes() genomeName <- "BSgenome.Hsapiens.UCSC.hg19" ################################################### ### code chunk number 19: selectGenomeFile (eval = FALSE) ################################################### ## genomeFile <- system.file(package="QuasR", "extdata", "hg19sub.fa") ################################################### ### code chunk number 20: preprocessReadsSingle ################################################### td <- tempdir() infiles <- system.file(package="QuasR", "extdata", c("rna_1_1.fq.bz2","rna_2_1.fq.bz2")) outfiles <- file.path(td, basename(infiles)) res <- preprocessReads(filename = infiles, outputFilename = outfiles, truncateEndBases = 3, Lpattern = "AAAAAAAAAA", minLength = 14, nBases = 2) res unlink(outfiles) ################################################### ### code chunk number 21: preprocessReadsPaired ################################################### td <- tempdir() infiles1 <- system.file(package="QuasR", "extdata", "rna_1_1.fq.bz2") infiles2 <- system.file(package="QuasR", "extdata", "rna_1_2.fq.bz2") outfiles1 <- file.path(td, basename(infiles1)) outfiles2 <- file.path(td, basename(infiles2)) res <- preprocessReads(filename=infiles1, filenameMate=infiles2, outputFilename=outfiles1, outputFilenameMate=outfiles2, nBases=0) res unlink(c(outfiles1,outfiles2)) ################################################### ### code chunk number 22: ChIP_copyExtdata ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 23: ChIP_qAlign ################################################### sampleFile <- "extdata/samples_chip_single.txt" auxFile <- "extdata/auxiliaries.txt" genomeFile <- "extdata/hg19sub.fa" proj1 <- qAlign(sampleFile, genome=genomeFile, auxiliaryFile=auxFile) proj1 ################################################### ### code chunk number 24: ChIP_bamfiles1 ################################################### list.files("extdata", pattern=".bam$") ################################################### ### code chunk number 25: ChIP_bamfiles2 ################################################### list.files("extdata", pattern="^chip_1_1_")[1:3] ################################################### ### code chunk number 26: ChIP_qcplot1 ################################################### qcdat1 <- qQCReport(proj1, pdfFilename="extdata/qc_report.pdf") ################################################### ### code chunk number 27: ChIP_qcplots2 ################################################### qQCReport(proj1, pdfFilename="extdata/qc_report.pdf") ################################################### ### code chunk number 28: ChIP_alignmentStats ################################################### alignmentStats(proj1) ################################################### ### code chunk number 29: ChIP_qExportWig ################################################### qExportWig(proj1, binsize=100L, scaling=TRUE, collapseBySample=TRUE) ################################################### ### code chunk number 30: ChIP_GenomicFeatures ################################################### library(GenomicFeatures) annotFile <- "extdata/hg19sub_annotation.gtf" chrLen <- scanFaIndex(genomeFile) chrominfo <- data.frame(chrom=as.character(seqnames(chrLen)), length=width(chrLen), is_circular=rep(FALSE, length(chrLen))) txdb <- makeTranscriptDbFromGFF(file=annotFile, format="gtf", exonRankAttributeName="exon_number", gffGeneIdAttributeName="gene_name", chrominfo=chrominfo, dataSource="Ensembl", species="Homo sapiens") promReg <- promoters(txdb, upstream=1000, downstream=500, columns=c("gene_id","tx_id")) gnId <- sapply(mcols(promReg)$gene_id, paste, collapse=",") promRegSel <- promReg[ match(unique(gnId), gnId) ] names(promRegSel) <- unique(gnId) head(promRegSel) ################################################### ### code chunk number 31: ChIP_qCount ################################################### cnt <- qCount(proj1, promRegSel) cnt ################################################### ### code chunk number 32: ChIP_visualize ################################################### gr1 <- import("Sample1.wig.gz", asRangedData=FALSE) gr2 <- import("Sample2.wig.gz", asRangedData=FALSE) library(Gviz) axisTrack <- GenomeAxisTrack() dTrack1 <- DataTrack(range=gr1, name="Sample 1", type="h") dTrack2 <- DataTrack(range=gr2, name="Sample 2", type="h") txTrack <- GeneRegionTrack(txdb, name="Transcripts", showId=TRUE) plotTracks(list(axisTrack, dTrack1, dTrack2, txTrack), chromosome="chr3", extend.left=1000) ################################################### ### code chunk number 33: ChIP_rtracklayer ################################################### library(rtracklayer) annotationFile <- "extdata/hg19sub_annotation.gtf" tssRegions <- import.gff(annotationFile, format="gtf", asRangedData=FALSE, feature.type="start_codon", colnames=c("strand", "gene_name")) tssRegions <- tssRegions[!duplicated(tssRegions)] names(tssRegions) <- rep("TSS", length(tssRegions)) head(tssRegions) ################################################### ### code chunk number 34: ChIP_qProfile ################################################### prS <- qProfile(proj1, tssRegions, upstream=3000, downstream=3000, orientation="same") prO <- qProfile(proj1, tssRegions, upstream=3000, downstream=3000, orientation="opposite") lapply(prS, "[", , 1:10) ################################################### ### code chunk number 35: ChIP__visualizeProfile ################################################### prCombS <- do.call("+", prS[-1]) /prS[[1]] prCombO <- do.call("+", prO[-1]) /prO[[1]] plot(as.numeric(colnames(prCombS)), filter(prCombS[1,], rep(1/100,100)), type='l', xlab="Position relative to TSS", ylab="Mean no. of alignments") lines(as.numeric(colnames(prCombO)), filter(prCombO[1,], rep(1/100,100)), type='l', col="red") legend(title="strand", legend=c("same as query","opposite of query"), x="topleft", col=c("black","red"), lwd=1.5, bty="n", title.adj=0.1) ################################################### ### code chunk number 36: ChIP_BSgenomeProject (eval = FALSE) ################################################### ## file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ## ## sampleFile <- "extdata/samples_chip_single.txt" ## auxFile <- "extdata/auxiliaries.txt" ## ## available.genomes() # list available genomes ## genomeName <- "BSgenome.Hsapiens.UCSC.hg19" ## ## proj1 <- qAlign(sampleFile, genome=genomeName, auxiliaryFile=auxFile) ## proj1 ################################################### ### code chunk number 37: RNA_qAlign ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) sampleFile <- "extdata/samples_rna_paired.txt" genomeFile <- "extdata/hg19sub.fa" proj2 <- qAlign(sampleFile, genome=genomeFile, splicedAlignment=TRUE) proj2 ################################################### ### code chunk number 38: RNA_alignmentStats ################################################### proj2unspl <- qAlign(sampleFile, genome=genomeFile, splicedAlignment=FALSE) alignmentStats(proj2) alignmentStats(proj2unspl) ################################################### ### code chunk number 39: RNA_qCount ################################################### geneLevels <- qCount(proj2, txdb, reportLevel="gene") exonLevels <- qCount(proj2, txdb, reportLevel="exon") head(geneLevels) head(exonLevels) ################################################### ### code chunk number 40: RNA_RPMK ################################################### geneRPKM <- t(t(geneLevels[,-1] /geneLevels[,1] *1000) /colSums(geneLevels[,-1]) *1e6) geneRPKM ################################################### ### code chunk number 41: RNA_junction ################################################### exonJunctions <- qCount(proj2, NULL, reportLevel="junction") exonJunctions ################################################### ### code chunk number 42: RNA_junction2 ################################################### knownIntrons <- unlist(intronsByTranscript(txdb)) isKnown <- overlapsAny(exonJunctions, knownIntrons, type="equal") table(isKnown) tapply(rowSums(as.matrix(mcols(exonJunctions))), isKnown, summary) ################################################### ### code chunk number 43: RNA_qcplot1 ################################################### qcdat2 <- qQCReport(proj2unspl, pdfFilename="qc_report.pdf") ################################################### ### code chunk number 44: Bis_qAlign ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) sampleFile <- "extdata/samples_bis_single.txt" genomeFile <- "extdata/hg19sub.fa" proj3 <- qAlign(sampleFile, genomeFile, bisulfite="dir") proj3 ################################################### ### code chunk number 45: Bis_qMeth1 ################################################### meth <- qMeth(proj3, mode="CpGcomb", collapseBySample=TRUE) meth ################################################### ### code chunk number 46: Bis_qMeth2 ################################################### chrs <- readDNAStringSet(genomeFile) sum(vcountPattern("CG",chrs)) length(qMeth(proj3)) length(qMeth(proj3, keepZero=FALSE)) ################################################### ### code chunk number 47: Bis_visualize ################################################### percMeth <- mcols(meth)[,2] *100 /mcols(meth)[,1] summary(percMeth) axisTrack <- GenomeAxisTrack() dTrack1 <- DataTrack(range=gr1, name="H3K4me3", type="h") dTrack2 <- DataTrack(range=meth, data=percMeth, name="Methylation", type="p") txTrack <- GeneRegionTrack(txdb, name="Transcripts", showId=TRUE) plotTracks(list(axisTrack, dTrack1, dTrack2, txTrack), chromosome="chr3", extend.left=1000) ################################################### ### code chunk number 48: Bis_query ################################################### qMeth(proj3, query=GRanges("chr1",IRanges(start=31633,width=2)), collapseBySample=TRUE) qMeth(proj3, query=promRegSel, collapseByQueryRegion=TRUE, collapseBySample=TRUE) ################################################### ### code chunk number 49: Alelle_Extdata ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 50: Allele_qAlign ################################################### sampleFile <- "extdata/samples_chip_single.txt" genomeFile <- "extdata/hg19sub.fa" snpFile <- "extdata/hg19sub_snp.txt" proj1SNP <- qAlign(sampleFile, genome=genomeFile, snpFile=snpFile) proj1SNP ################################################### ### code chunk number 51: Allele_qCount ################################################### head(qCount(proj1, promRegSel)) head(qCount(proj1SNP, promRegSel)) ################################################### ### code chunk number 52: Allele_Bis ################################################### sampleFile <- "extdata/samples_bis_single.txt" genomeFile <- "extdata/hg19sub.fa" proj3SNP <- qAlign(sampleFile, genomeFile, snpFile=snpFile, bisulfite="dir") head(qMeth(proj3SNP, mode="CpGcomb", collapseBySample=TRUE)) ################################################### ### code chunk number 53: qcplotsFig1 ################################################### QuasR:::plotQualByCycle(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 54: qcplotsFig2 ################################################### QuasR:::plotNuclByCycle(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 55: qcplotsFig3 ################################################### QuasR:::plotDuplicated(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 56: qcplotsFig4 ################################################### QuasR:::plotMappings(qcdat1$raw$mapdata, a4layout=FALSE) ################################################### ### code chunk number 57: qcplotsFig5 ################################################### QuasR:::plotUniqueness(qcdat1$raw$unique, a4layout=FALSE) ################################################### ### code chunk number 58: qcplotsFig6 ################################################### QuasR:::plotErrorsByCycle(qcdat1$raw$mm, lmat=rbind(1:2)) ################################################### ### code chunk number 59: qcplotsFig7 ################################################### QuasR:::plotMismatchTypes(qcdat1$raw$mm, lmat=rbind(1:2)) ################################################### ### code chunk number 60: qcplotsFig8 ################################################### QuasR:::plotFragmentDistribution(qcdat2$raw$frag, lmat=rbind(1:2)) ################################################### ### code chunk number 61: alignmentStats ################################################### # using bam files alignmentStats(alignments(proj1)$genome$FileName) alignmentStats(unlist(alignments(proj1)$aux)) # using a qProject object alignmentStats(proj1) ################################################### ### code chunk number 62: sessionInfo ################################################### sessionInfo() ################################################### ### code chunk number 63: cleanUp ################################################### unlink("extdata", recursive=TRUE, force=TRUE)