### R code from vignette source 'QuasR.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(width=80, bibstyle="unsrturl")


###################################################
### 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", 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: <sampleFileSeqToBam (eval = FALSE)
###################################################
## sampleFile1 <- "samples_fastq.txt"
## sampleFile2 <- "samples_bam.txt"
## 
## proj1 <- qAlign(sampleFile1, genomeFile)
## 
## write.table(alignments(proj1)$genome, sampleFile2, sep="\t", row.names=FALSE)
## 
## proj2 <- qAlign(sampleFile2, genomeFile)


###################################################
### code chunk number 19: auxiliaryFile
###################################################
auxFile <- system.file(package="QuasR", "extdata", "auxiliaries.txt")


###################################################
### code chunk number 20: selectGenomeBSgenome
###################################################
available.genomes()
genomeName <- "BSgenome.Hsapiens.UCSC.hg19"


###################################################
### code chunk number 21: selectGenomeFile (eval = FALSE)
###################################################
## genomeFile <- system.file(package="QuasR", "extdata", "hg19sub.fa")


###################################################
### code chunk number 22: 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 23: 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 24: ChIP_copyExtdata
###################################################
file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)


###################################################
### code chunk number 25: 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 26: ChIP_bamfiles1
###################################################
list.files("extdata", pattern=".bam$")


###################################################
### code chunk number 27: ChIP_bamfiles2
###################################################
list.files("extdata", pattern="^chip_1_1_")[1:3]


###################################################
### code chunk number 28: ChIP_qcplot1
###################################################
qcdat1 <- qQCReport(proj1, pdfFilename="extdata/qc_report.pdf")


###################################################
### code chunk number 29: ChIP_qcplots2
###################################################
qQCReport(proj1, pdfFilename="extdata/qc_report.pdf")


###################################################
### code chunk number 30: ChIP_alignmentStats
###################################################
alignmentStats(proj1)


###################################################
### code chunk number 31: ChIP_qExportWig
###################################################
qExportWig(proj1, binsize=100L, scaling=TRUE, collapseBySample=TRUE)


###################################################
### code chunk number 32: 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 <- makeTxDbFromGFF(file=annotFile, format="gtf",
                        chrominfo=chrominfo,
                        dataSource="Ensembl",
                        organism="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 33: ChIP_qCount
###################################################
cnt <- qCount(proj1, promRegSel)
cnt


###################################################
### code chunk number 34: ChIP_visualize
###################################################
gr1 <- import("Sample1.wig.gz")
gr2 <- import("Sample2.wig.gz")

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 35: ChIP_rtracklayer
###################################################
library(rtracklayer)
annotationFile <- "extdata/hg19sub_annotation.gtf"
tssRegions <- import.gff(annotationFile, format="gtf",
                         feature.type="start_codon",
                         colnames="gene_name")
tssRegions <- tssRegions[!duplicated(tssRegions)]
names(tssRegions) <- rep("TSS", length(tssRegions))
head(tssRegions)


###################################################
### code chunk number 36: 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 37: 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 38: 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 39: 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 40: RNA_alignmentStats
###################################################
proj2unspl <- qAlign(sampleFile, genome=genomeFile,
                     splicedAlignment=FALSE)

alignmentStats(proj2)
alignmentStats(proj2unspl)


###################################################
### code chunk number 41: RNA_qCount
###################################################
geneLevels <- qCount(proj2, txdb, reportLevel="gene")
exonLevels <- qCount(proj2, txdb, reportLevel="exon")

head(geneLevels)
head(exonLevels)


###################################################
### code chunk number 42: RNA_RPMK
###################################################
geneRPKM <- t(t(geneLevels[,-1] /geneLevels[,1] *1000)
              /colSums(geneLevels[,-1]) *1e6)
geneRPKM


###################################################
### code chunk number 43: RNA_junction
###################################################
exonJunctions <- qCount(proj2, NULL, reportLevel="junction")
exonJunctions


###################################################
### code chunk number 44: 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 45: RNA_includeSpliced
###################################################
exonBodyLevels <- qCount(proj2, txdb, reportLevel="exon", includeSpliced=FALSE)
summary(exonLevels - exonBodyLevels)


###################################################
### code chunk number 46: RNA_qcplot1
###################################################
qcdat2 <- qQCReport(proj2unspl, pdfFilename="qc_report.pdf")


###################################################
### code chunk number 47: miRNA_extdata
###################################################
file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)


###################################################
### code chunk number 48: 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 49: miRNA_lengthes
###################################################
# get read lengths
library(Biostrings)
oldwd <- setwd(dirname(sampleFile))
lens <- fasta.seqlengths(tab$FileName, nrec=1e5)
lens2 <- fasta.seqlengths(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 50: miRNA_qAlign
###################################################
proj3 <- qAlign(sampleFile2, genomeFile, maxHits=50)
alignmentStats(proj3)


###################################################
### code chunk number 51: miRNA_prepareQuery
###################################################
mirs <- import("extdata/mirbaseXX_qsr.gff3")
names(mirs) <- mirs$Name
preMirs <- mirs[ mirs$type=="miRNA_primary_transcript" ]
matureMirs <- mirs[ mirs$type=="miRNA" ]


###################################################
### code chunk number 52: 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 53: 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 54: 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 55: 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 56: Bis_qMeth1
###################################################
meth <- qMeth(proj4, mode="CpGcomb", collapseBySample=TRUE)
meth


###################################################
### code chunk number 57: Bis_qMeth2
###################################################
chrs <- readDNAStringSet(genomeFile)
sum(vcountPattern("CG",chrs))
length(qMeth(proj4))
length(qMeth(proj4, keepZero=FALSE))


###################################################
### code chunk number 58: 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 59: 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 60: Alelle_Extdata
###################################################
file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)


###################################################
### code chunk number 61: 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 62: Allele_qCount
###################################################
head(qCount(proj1, promRegSel))
head(qCount(proj1SNP, promRegSel))


###################################################
### code chunk number 63: 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 64: qcplotsFig1
###################################################
QuasR:::plotQualByCycle(qcdat1$raw$qa, lmat=rbind(1:2))


###################################################
### code chunk number 65: qcplotsFig2
###################################################
QuasR:::plotNuclByCycle(qcdat1$raw$qa, lmat=rbind(1:2))


###################################################
### code chunk number 66: qcplotsFig3
###################################################
QuasR:::plotDuplicated(qcdat1$raw$qa, lmat=rbind(1:2))


###################################################
### code chunk number 67: qcplotsFig4
###################################################
QuasR:::plotMappings(qcdat1$raw$mapdata, a4layout=FALSE)


###################################################
### code chunk number 68: qcplotsFig5
###################################################
QuasR:::plotUniqueness(qcdat1$raw$unique, a4layout=FALSE)


###################################################
### code chunk number 69: qcplotsFig6
###################################################
QuasR:::plotErrorsByCycle(qcdat1$raw$mm, lmat=rbind(1:2))


###################################################
### code chunk number 70: qcplotsFig7
###################################################
QuasR:::plotMismatchTypes(qcdat1$raw$mm, lmat=rbind(1:2))


###################################################
### code chunk number 71: qcplotsFig8
###################################################
QuasR:::plotFragmentDistribution(qcdat2$raw$frag, lmat=rbind(1:2))


###################################################
### code chunk number 72: alignmentStats
###################################################
# using bam files
alignmentStats(alignments(proj1)$genome$FileName)
alignmentStats(unlist(alignments(proj1)$aux))

# using a qProject object
alignmentStats(proj1)


###################################################
### code chunk number 73: sessionInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 74: cleanUp
###################################################
unlink("extdata", recursive=TRUE, force=TRUE)