### R code from vignette source 'vignettes/QuasR/inst/doc/QuasR-Overview.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex(width=80) ################################################### ### code chunk number 2: options ################################################### #options(width=65) options('useFancyQuotes' = FALSE, continue=" ", digits=3) ################################################### ### code chunk number 3: cite ################################################### citation("QuasR") ################################################### ### code chunk number 4: install (eval = FALSE) ################################################### ## source("http://www.bioconductor.org/biocLite.R") ## biocLite("QuasR") ################################################### ### code chunk number 5: loadLibraries ################################################### library(QuasR) library(BSgenome) library(Rsamtools) library(rtracklayer) library(GenomicFeatures) library(Gviz) ################################################### ### code chunk number 6: help1 (eval = FALSE) ################################################### ## help.start() ################################################### ### code chunk number 7: loadQuasRLibrary (eval = FALSE) ################################################### ## library(QuasR) ################################################### ### code chunk number 8: help2 (eval = FALSE) ################################################### ## ?preprocessReads ################################################### ### code chunk number 9: help3 (eval = FALSE) ################################################### ## help("preprocessReads") ################################################### ### code chunk number 10: assign (eval = FALSE) ################################################### ## x <- 2 ################################################### ### code chunk number 11: ls (eval = FALSE) ################################################### ## ls() ################################################### ### code chunk number 12: printObject (eval = FALSE) ################################################### ## x ################################################### ### code chunk number 13: SampleSession1 ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 14: SampleSession2 ################################################### sampleFile <- "extdata/samples_chip_single.txt" genomeFile <- "extdata/hg19sub.fa" proj <- qAlign(sampleFile, genomeFile) proj ################################################### ### code chunk number 15: SampleSession3 ################################################### qQCReport(proj, "extdata/qc_report.pdf") ################################################### ### code chunk number 16: 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 17: 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 18: auxiliaryFile ################################################### auxFile <- system.file(package="QuasR", "extdata", "auxiliaries.txt") ################################################### ### code chunk number 19: selectGenomeBSgenome ################################################### available.genomes() genomeName <- "BSgenome.Hsapiens.UCSC.hg19" ################################################### ### code chunk number 20: selectGenomeFile (eval = FALSE) ################################################### ## genomeFile <- system.file(package="QuasR", "extdata", "hg19sub.fa") ################################################### ### code chunk number 21: 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 22: 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 23: ChIP_copyExtdata ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 24: 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 25: ChIP_bamfiles1 ################################################### list.files("extdata", pattern=".bam$") ################################################### ### code chunk number 26: ChIP_bamfiles2 ################################################### list.files("extdata", pattern="^chip_1_1_")[1:3] ################################################### ### code chunk number 27: ChIP_qcplot1 ################################################### qcdat1 <- qQCReport(proj1, pdfFilename="extdata/qc_report.pdf") ################################################### ### code chunk number 28: ChIP_qcplots2 ################################################### qQCReport(proj1, pdfFilename="extdata/qc_report.pdf") ################################################### ### code chunk number 29: ChIP_alignmentStats ################################################### alignmentStats(proj1) ################################################### ### code chunk number 30: ChIP_qExportWig ################################################### qExportWig(proj1, binsize=100L, scaling=TRUE, collapseBySample=TRUE) ################################################### ### code chunk number 31: 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 32: ChIP_qCount ################################################### cnt <- qCount(proj1, promRegSel) cnt ################################################### ### code chunk number 33: 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 34: 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 35: 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 36: 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 37: 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 38: 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 39: RNA_alignmentStats ################################################### proj2unspl <- qAlign(sampleFile, genome=genomeFile, splicedAlignment=FALSE) alignmentStats(proj2) alignmentStats(proj2unspl) ################################################### ### code chunk number 40: RNA_qCount ################################################### geneLevels <- qCount(proj2, txdb, reportLevel="gene") exonLevels <- qCount(proj2, txdb, reportLevel="exon") head(geneLevels) head(exonLevels) ################################################### ### code chunk number 41: RNA_RPMK ################################################### geneRPKM <- t(t(geneLevels[,-1] /geneLevels[,1] *1000) /colSums(geneLevels[,-1]) *1e6) geneRPKM ################################################### ### code chunk number 42: RNA_junction ################################################### exonJunctions <- qCount(proj2, NULL, reportLevel="junction") exonJunctions ################################################### ### code chunk number 43: 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 44: RNA_includeSpliced ################################################### exonBodyLevels <- qCount(proj2, txdb, reportLevel="exon", includeSpliced=FALSE) summary(exonLevels - exonBodyLevels) ################################################### ### code chunk number 45: RNA_qcplot1 ################################################### qcdat2 <- qQCReport(proj2unspl, pdfFilename="qc_report.pdf") ################################################### ### code chunk number 46: miRNA_extdata ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 47: miRNA_preprocessReads ################################################### # prepare sample file with processed reads filenames sampleFile <- file.path("extdata", "samples_mirna.txt") sampleFile sampleFile2 <- sub(".txt", "_processed.txt", sampleFile) sampleFile2 tab <- read.delim(sampleFile, header=TRUE, as.is=TRUE) tab tab2 <- tab tab2$FileName <- sub(".fa", "_processed.fa", tab$FileName) write.table(tab2, sampleFile2, sep="\t", quote=FALSE, row.names=FALSE) tab2 # remove adapters oldwd <- setwd(dirname(sampleFile)) res <- preprocessReads(tab$FileName, tab2$FileName, Rpattern="TGGAATTCTCGGGTGCCAAGG") res setwd(oldwd) ################################################### ### code chunk number 48: miRNA_lengthes ################################################### # get read lengthes library(Biostrings) oldwd <- setwd(dirname(sampleFile)) lens <- fasta.info(tab$FileName, nrec=1e5) lens2 <- fasta.info(tab2$FileName, nrec=1e5) setwd(oldwd) # plot length distribution lensTab <- rbind(raw=tabulate(lens,50), processed=tabulate(lens2,50)) colnames(lensTab) <- 1:50 barplot(lensTab/rowSums(lensTab)*100, xlab="Read length (nt)", ylab="Percent of reads") legend(x="topleft", bty="n", fill=gray.colors(2), legend=rownames(lensTab)) ################################################### ### code chunk number 49: miRNA_qAlign ################################################### proj3 <- qAlign(sampleFile2, genomeFile, maxHits=50) alignmentStats(proj3) ################################################### ### code chunk number 50: miRNA_prepareQuery ################################################### mirs <- import("extdata/mirbaseXX_qsr.gff3", asRangedData=FALSE) names(mirs) <- mirs$Name preMirs <- mirs[ mirs$type=="miRNA_primary_transcript" ] matureMirs <- mirs[ mirs$type=="miRNA" ] ################################################### ### code chunk number 51: miRNA_coverage ################################################### library(Rsamtools) alns <- scanBam(alignments(proj3)$genome$FileName, param=ScanBamParam(what=scanBamWhat(), which=preMirs[1]))[[1]] alnsIR <- IRanges(start=alns$pos - start(preMirs), width=alns$qwidth) mp <- barplot(as.vector(coverage(alnsIR)), names.arg=1:max(max(alnsIR)), xlab="Relative position in pre-miRNA", ylab="Alignment coverage") rect(xleft=mp[start(matureMirs)-start(preMirs)+1,1], ybottom=-par('cxy')[2], xright=mp[end(matureMirs)-start(preMirs)+1,1], ytop=0, col="#CCAA0088", border=NA, xpd=NA) ################################################### ### code chunk number 52: miRNA_extendQuery ################################################### matureMirsExtended <- matureMirs start(matureMirsExtended) <- ifelse(as.logical(strand(matureMirs)=="+"), start(matureMirs) - 3, end(matureMirs) - 3) end(matureMirsExtended) <- ifelse(as.logical(strand(matureMirs)=="+"), start(matureMirs) + 3, end(matureMirs) + 3) ################################################### ### code chunk number 53: miRNA_quantify ################################################### # quantify mature miRNAs cnt <- qCount(proj3, matureMirsExtended, orientation="same") cnt # quantify pre-miRNAs cnt <- qCount(proj3, preMirs, orientation="same") cnt ################################################### ### code chunk number 54: Bis_qAlign ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) sampleFile <- "extdata/samples_bis_single.txt" genomeFile <- "extdata/hg19sub.fa" proj4 <- qAlign(sampleFile, genomeFile, bisulfite="dir") proj4 ################################################### ### code chunk number 55: Bis_qMeth1 ################################################### meth <- qMeth(proj4, mode="CpGcomb", collapseBySample=TRUE) meth ################################################### ### code chunk number 56: Bis_qMeth2 ################################################### chrs <- readDNAStringSet(genomeFile) sum(vcountPattern("CG",chrs)) length(qMeth(proj4)) length(qMeth(proj4, keepZero=FALSE)) ################################################### ### code chunk number 57: 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 58: Bis_query ################################################### qMeth(proj4, query=GRanges("chr1",IRanges(start=31633,width=2)), collapseBySample=TRUE) qMeth(proj4, query=promRegSel, collapseByQueryRegion=TRUE, collapseBySample=TRUE) ################################################### ### code chunk number 59: Alelle_Extdata ################################################### file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE) ################################################### ### code chunk number 60: 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 61: Allele_qCount ################################################### head(qCount(proj1, promRegSel)) head(qCount(proj1SNP, promRegSel)) ################################################### ### code chunk number 62: Allele_Bis ################################################### sampleFile <- "extdata/samples_bis_single.txt" genomeFile <- "extdata/hg19sub.fa" proj4SNP <- qAlign(sampleFile, genomeFile, snpFile=snpFile, bisulfite="dir") head(qMeth(proj4SNP, mode="CpGcomb", collapseBySample=TRUE)) ################################################### ### code chunk number 63: qcplotsFig1 ################################################### QuasR:::plotQualByCycle(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 64: qcplotsFig2 ################################################### QuasR:::plotNuclByCycle(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 65: qcplotsFig3 ################################################### QuasR:::plotDuplicated(qcdat1$raw$qa, lmat=rbind(1:2)) ################################################### ### code chunk number 66: qcplotsFig4 ################################################### QuasR:::plotMappings(qcdat1$raw$mapdata, a4layout=FALSE) ################################################### ### code chunk number 67: qcplotsFig5 ################################################### QuasR:::plotUniqueness(qcdat1$raw$unique, a4layout=FALSE) ################################################### ### code chunk number 68: qcplotsFig6 ################################################### QuasR:::plotErrorsByCycle(qcdat1$raw$mm, lmat=rbind(1:2)) ################################################### ### code chunk number 69: qcplotsFig7 ################################################### QuasR:::plotMismatchTypes(qcdat1$raw$mm, lmat=rbind(1:2)) ################################################### ### code chunk number 70: qcplotsFig8 ################################################### QuasR:::plotFragmentDistribution(qcdat2$raw$frag, lmat=rbind(1:2)) ################################################### ### code chunk number 71: alignmentStats ################################################### # using bam files alignmentStats(alignments(proj1)$genome$FileName) alignmentStats(unlist(alignments(proj1)$aux)) # using a qProject object alignmentStats(proj1) ################################################### ### code chunk number 72: sessionInfo ################################################### toLatex(sessionInfo()) ################################################### ### code chunk number 73: cleanUp ################################################### unlink("extdata", recursive=TRUE, force=TRUE)