2.1 Citing QuasR

If you use QuasR (Gaidatzis et al. 2015) in your work, you can cite it as follows:

citation("QuasR")

## 
## Please use the QuasR reference below to cite the software itself. If
## you were using qAlign with Rbowtie as aligner, it can be cited as
## Langmead et al. (2009) (unspliced alignments) or Au et al. (2010)
## (spliced alignments). If you were using qAlign with Rhisat2 as
## aligner, it can be cited as Kim et al. (2015).
## 
##   Gaidatzis D, Lerch A, Hahne F, Stadler MB. QuasR: Quantification
##   and annotation of short reads in R. Bioinformatics 31(7):1130-1132
##   (2015).
## 
##   Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and
##   memory-efficient alignment of short DNA sequences to the human
##   genome. Genome Biology 10(3):R25 (2009).
## 
##   Au KF, Jiang H, Lin L, Xing Y, Wong WH. Detection of splice
##   junctions from paired-end RNA-seq data by SpliceMap. Nucleic Acids
##   Research, 38(14):4570-8 (2010).
## 
##   Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with
##   low memory requirements. Nat Methods, 12(4):357-60 (2015).
## 
## This free open-source software implements academic research by the
## authors and co-workers. If you use it, please support the project by
## citing the appropriate journal articles.
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.

2.2 Installation

QuasR is a package for the R computing environment and it is assumed that you have already installed R. See the R project at (http://www.r-project.org). To install the latest version of QuasR, you will need to be using the latest version of R. QuasR is part of the Bioconductor project at (http://www.bioconductor.org). To get QuasR together with its dependencies you can use

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("QuasR")

Bioconductor works on a 6-monthly official release cycle. As with other Bioconductor packages, there are always two versions of QuasR. Most users will use the current official release version, which will be installed by BiocManager::install if you are using the current release version of R. There is also a development version of QuasR that includes new features due for the next official release. The development version will be installed if you are using the development version of Bioconductor (see version = "devel" in BiocManager). The official release version always has an even second number (for example 1.20.1), whereas the developmental version has an odd second number (for example 1.21.4).

2.3 Loading of QuasR and other required packages

In order to run the code examples in this vignette, the QuasR package and a few additional packages need to be loaded:

suppressPackageStartupMessages({
    library(QuasR)
    library(BSgenome)
    library(Rsamtools)
    library(rtracklayer)
    library(GenomicFeatures)
    library(Gviz)
})

2.4 How to get help

Most questions about QuasR will hopefully be answered by the documentation or references. If you’ve run into a question which isn’t addressed by the documentation, or you’ve found a conflict between the documentation and software itself, then there is an active support community which can offer help.

The authors of the package (maintainer: Michael Stadler michael.stadler@fmi.ch) always appreciate receiving reports of bugs in the package functions or in the documentation. The same goes for well-considered suggestions for improvements.

Any other questions or problems concerning QuasR should be posted to the Bioconductor support site (https://support.bioconductor.org). Users posting to the support site for the first time should read the helpful posting guide at (https://support.bioconductor.org/info/faq/). Note that each function in QuasR has it’s own help page, as described in the section 3.1. Posting etiquette requires that you read the relevant help page carefully before posting a problem to the site.

3.1 A brief introduction to R

If you already use R and know about its command line interface, just skip this section and continue with section 3.2.

The structure of this vignette and in particular this section is based on the excellent user guide of the limma package, which we would like to hereby acknowledge. R is a program for statistical computing. It is a command-driven language meaning that you have to type commands into it rather than pointing and clicking using a mouse. In this guide it will be assumed that you have successfully downloaded and installed R from (http://www.r-project.org) as well as QuasR (see section 2.2). A good way to get started is to type

help.start()

at the R prompt or, if you’re using R for Windows, to follow the drop-down menu items Help \(\succ\) Html help. Following the links Packages \(\succ\) QuasR from the html help page will lead you to the contents page of help topics for functions in QuasR.

Before you can use any QuasR commands you have to load the package by typing

library(QuasR)

at the R prompt. You can get help on any function in any loaded package by typing ? and the function name at the R prompt, for example

?preprocessReads

or equivalently

help("preprocessReads")

for detailed help on the preprocessReads function. The individual function help pages are especially important for listing all the arguments which a function will accept and what values the arguments can take.

A key to understanding R is to appreciate that anything that you create in R is an object. Objects might include data sets, variables, functions, anything at all. For example

x <- 2

will create a variable x and will assign it the value 2. At any stage of your R session you can type

ls()

to get a list of all the objects you have created. You can see the contents of any object by typing the name of the object at the prompt. The following command will print out the contents of x:

x

We hope that you can use QuasR without having to spend a lot of time learning about the R language itself but a little knowledge in this direction will be very helpful, especially when you want to do something not explicitly provided for in QuasR or in the other Bioconductor packages. For more details about the R language see An Introduction to R which is available from the online help. For more background on using R for statistical analysis see (Dalgaard 2002).

3.2 Sample QuasR session

This is a quick overview of what an analysis could look like for users preferring to jump right into an analysis. The example uses data that is provided with the QuasR package, which is first copied to the current working directory, into a subfolder called "extdata":

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

The sequence files to be analyzed are listed in sampleFile (see section 5.1 for details). The sequence reads will be aligned using bowtie (Langmead et al. 2009) (from the Rbowtie package (Hahne, Lerch, and Stadler 2012)) to a small reference genome (consisting of three short segments from the hg19 human genome assembly, available in full for example in the BSgenome.Hsapiens.UCSC.hg19 package). Make sure that you have sufficient disk space, both in your R temporary directory (tempdir()) as well as to store the resulting alignments (see section 7.2).

sampleFile <- "extdata/samples_chip_single.txt"
genomeFile <- "extdata/hg19sub.fa"

proj <- qAlign(sampleFile, genomeFile)

## Creating .fai file for: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub.fa

## alignment files missing - need to:

##     create alignment index for the genome

##     create 2 genomic alignment(s)

## Creating an Rbowtie index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 2 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a3fdd8ac0.txt
## Genomic alignments have been created successfully

proj

## Project: qProject
##  Options   : maxHits         : 1
##              paired          : no
##              splicedAlignment: FALSE
##              bisulfite       : no
##              snpFile         : none
##              geneAnnotation  : none
##  Aligner   : Rbowtie v1.24.0 (parameters: -m 1 --best --strata)
##  Genome    : /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vi.../hg19sub.fa (file)
## 
##  Reads     : 2 files, 2 samples (fastq format):
##    1. chip_1_1.fq.bz2  Sample1 (phred33)
##    2. chip_2_1.fq.bz2  Sample2 (phred33)
## 
##  Genome alignments: directory: same as reads
##    1. chip_1_1_d0a6130e4bf.bam
##    2. chip_2_1_d0a723b4f4d.bam
## 
##  Aux. alignments: none

The proj object keeps track of all the information of a sequencing experiment, for example where sequence and alignment files are stored, and what aligner and reference genome was used to generate the alignments.

Now that the alignments have been generated, further analyses can be performed. A quality control report is saved to the "extdata/qc_report.pdf" file using the qQCReport function.

qQCReport(proj, "extdata/qc_report.pdf")

## collecting quality control data

## creating QC plots

The number of alignments per promoter region is quantified using qCount. Genomic coordinates for promoter regions are imported from a gtf file (annotFile) into the GRanges-object with the name promReg:

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)

## counting alignments...done

promCounts

##                width Sample1 Sample2
## TNFRSF18-003    1000      20       4
## TNFRSF18-002    1000      20       4
## TNFRSF18-001    1000      20       4
## TNFRSF4-001     1000       5       2
## SDF4-007        1000       8       2
## SDF4-001        1000       8       2
## SDF4-002        1000       8       2
## SDF4-201        1000       8       2
## B3GALT6-001     1000      25     274
## RPS7-001        1000     121     731
## RPS7-008        1000     121     731
## RPS7-009        1000     121     731
## RPS7-005        1000     121     731
## C3orf10-201     1000     176     496
## C3orf10-001     1000     176     496
## AC034193.1-201  1000       5       2
## VHL-001         1000      61     336
## VHL-002         1000      61     336
## VHL-201         1000      61     336

4.1 File storage locations

Apart from qExportWig and qQCReport, which generate wig files and pdf reports, qAlign is the only function in QuasR that stores files on the disk (see section 7.2 for details). All files generated by qAlign are listed here by type, together with their default location and how locations can be changed.

• Temporary files (default: tempdir()): Temporary files include reference genomes in fasta format, decompressed input sequence files, and temporary alignments in text format, and can require a large amount of disk space. By default, these files will be written to the temporary directory of the R process (as reported by tempdir()). If using clObj for parallel processing, this may be the tempdir() from the cluster node(s). An alternative location can be set using the TMPDIR environment variable (see ?tempdir).
  
• Alignment files (bam format) (default: same directory as the input sequence files): Alignments against reference genome and auxiliary targets are stored in bam format in the same directory that also contains the input sequence file (listed in sampleFile). Please note that if the input sequence file corresponds to a symbolic link, QuasR will follow the link and use the directory of the original file instead. An alternative directory can be specified with the alignmentsDir argument from qAlign, which will store all bam files in that directory even if the input sequence files are located in different directories.
  
• Alignment index files (default: depends on genome and snpFile arguments): Many alignment tools including bowtie require an index of the reference sequence to perform alignments. If necessary, qAlign will build this index automatically and store it in a default location that depends on the genome argument:
  • BSgenome: If genome is the name of a BSgenome package (such as "BSgenome.Hsapiens.UCSC.hg19"), the index will be stored as a new R package in the default library path (as reported by .libPaths()[1], see ?install.packages for details). The name of this index package will be the name of the original BSgenome package with a suffix for the index type, for example "BSgenome.Hsapiens.UCSC.hg19.Rbowtie".
    
  • fasta: If genome refers to a reference genome file in fasta format, the index will be stored in a subdirectory at the same location. Similarly, the indices for files listed in auxiliaryFile are store at the location of these files. For example, the Rbowtie index for the genome at "./genome/mm9.fa" is stored in "./genome/mm9.fa.Rbowtie".
    
  • Allele-specific analysis: A special case is the allele-specific analysis, where reference and alternative alleles listed in snpFile (e.g. "./mySNPs.tab") are injected into the genome (e.g. "BSgenome.Mmusculus.UCSC.mm9") to create two variant genomes to be indexed. These indices are saved at the location of the snpFile in a directory named after snpFile, genome and the index type (e.g. "./mySNPs.tab.BSgenome.Mmusculus.UCSC.mm9.A.fa.Rbowtie").

5.1 Create a sample file

The sample file is a tab-delimited text file with two or three columns. The first row contains the column names: For a single read experiment, these are ‘FileName’ and ‘SampleName’; for a paired-end experiment, these are ‘FileName1’, ‘FileName2’ and ‘SampleName’. If the first row does not contain the correctly spelled column names, QuasR will not accept the samples file. Subsequent rows contain the input sequence files.

Here are examples of such sample files for a single read experiment:

FileName    SampleName
chip_1_1.fq.bz2 Sample1
chip_2_1.fq.bz2 Sample2

and for a paired-end experiment:

FileName1   FileName2   SampleName
rna_1_1.fq.bz2  rna_1_2.fq.bz2  Sample1
rna_2_1.fq.bz2  rna_2_2.fq.bz2  Sample2

These example files are also contained in the QuasR package and may be used as templates. The path of the files can be determined using:

sampleFile1 <- system.file(package="QuasR", "extdata",
                           "samples_chip_single.txt")
sampleFile2 <- system.file(package="QuasR", "extdata",
                           "samples_rna_paired.txt")

The columns FileName for single-read, or FileName1 and FileName2 for paired-end experiments contain paths and names to files containing the sequence data. The paths can be absolute or relative to the location of the sample file. This allows combining files from different directories in a single analysis. For each input sequence file, qAlign will create one alignment file and by default store it in the same directory as the sequence file. Already existing alignment files with identical parameters will not be re-created, so that it is easy to reuse the same sequence files in multiple projects without unnecessarily copying sequence files or recreating alignments.

The SampleName column contains sample names for each sequence file. The same name can be used on several lines to indicate multiple sequence files that belong to the same sample (qCount can use this information to automatically combine counts for one sample from multiple files).

Three file formats are supported for input files (but cannot be mixed within a single sample file):

• fasta files have names that end with ‘.fa’, ‘.fna’ or ‘.fasta’. They contain only sequences (and no base qualities) and will thus by default be aligned on the basis of mismatches (the best alignment is the one with fewest mismatches).
  
• fastq files have names that end with ‘.fq’ or ‘.fastq’. They contain sequences and corresponding base qualities and will be aligned by default using these qualities. The encoding scheme of base qualities is automatically detected for each individual fastq file.
  
• bam files have names that end with ‘.bam’. They can be used if the sequence reads have already been aligned outside of QuasR, and QuasR will only be used for downstream analysis based on the alignments contained in the bam files. This makes it possible to use alignment tools that are not available within QuasR, but making use of this option comes with a risk and should only be used by experienced users. For example, it cannot be guaranteed any more that certain assumptions made by qCount are fulfilled by the external aligner (see below). When using external bam files, we recommend to use files which contain only one alignment per read. This may also include multi-hit reads, for which one of the alignments is randomly selected. This allows QuasR to count the total number of reads by counting the total number of alignments. Furthermore, if the bam files also contain the unmapped reads, QuasR will be able to calculate the fraction of mapped reads. For bisulfite samples we require ungapped alignments stored in unpaired or paired ff orientation (even if the input reads are fr). For allele-specific bam files, QuasR requires an additional tag for each alignment called XV of type A (printable character) with the possible values R (Reference), U (Unknown) and A (Alternative).

fasta and fastq files can be compressed with gzip, bzip2 or xz (file extensions ‘.gz’, ‘.bz2’ or ‘xz’, respectively) and will be automatically decompressed when necessary.

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)

5.2 Create an auxiliary file (optional)

By default QuasR aligns reads only to the reference genome. However, it may be interesting to align non-matching reads to further targets, for example to identify contamination from vectors or a different species, or in order to quantify spike-in material not contained in the reference genome. In QuasR, such supplementary reference files are called auxiliary references and can be specified to qAlign using the auxiliaryFile argument (see section 7.2 for details). The format of the auxiliary file is similar to the one of the sample file described in section 5.1: It contains two columns with column names ‘FileName’ and ‘AuxName’ in the first row. Additional rows contain names and files of one or several auxiliary references in fasta format.

An example auxiliary file looks like this:

FileName    AuxName
NC_001422.1.fa  phiX174

and is available from your QuasR installation at

auxFile <- system.file(package="QuasR", "extdata", "auxiliaries.txt")

5.3 Select the reference genome

Sequence reads are primarily aligned against the reference genome. If necessary, QuasR will create an aligner index for the genome. The reference genome can be provided in one of two different formats:

• a string, referring to the name of a BSgenome package:

  available.genomes()

  ##  [1] "BSgenome.Alyrata.JGI.v1"                  
  ##  [2] "BSgenome.Amellifera.BeeBase.assembly4"    
  ##  [3] "BSgenome.Amellifera.UCSC.apiMel2"         
  ##  [4] "BSgenome.Amellifera.UCSC.apiMel2.masked"  
  ##  [5] "BSgenome.Aofficinalis.NCBI.V1"            
  ##  [6] "BSgenome.Athaliana.TAIR.04232008"         
  ##  [7] "BSgenome.Athaliana.TAIR.TAIR9"            
  ##  [8] "BSgenome.Btaurus.UCSC.bosTau3"            
  ##  [9] "BSgenome.Btaurus.UCSC.bosTau3.masked"     
  ## [10] "BSgenome.Btaurus.UCSC.bosTau4"            
  ## [11] "BSgenome.Btaurus.UCSC.bosTau4.masked"     
  ## [12] "BSgenome.Btaurus.UCSC.bosTau6"            
  ## [13] "BSgenome.Btaurus.UCSC.bosTau6.masked"     
  ## [14] "BSgenome.Btaurus.UCSC.bosTau8"            
  ## [15] "BSgenome.Carietinum.NCBI.v1"              
  ## [16] "BSgenome.Celegans.UCSC.ce10"              
  ## [17] "BSgenome.Celegans.UCSC.ce11"              
  ## [18] "BSgenome.Celegans.UCSC.ce2"               
  ## [19] "BSgenome.Celegans.UCSC.ce6"               
  ## [20] "BSgenome.Cfamiliaris.UCSC.canFam2"        
  ## [21] "BSgenome.Cfamiliaris.UCSC.canFam2.masked" 
  ## [22] "BSgenome.Cfamiliaris.UCSC.canFam3"        
  ## [23] "BSgenome.Cfamiliaris.UCSC.canFam3.masked" 
  ## [24] "BSgenome.Cjacchus.UCSC.calJac3"           
  ## [25] "BSgenome.Dmelanogaster.UCSC.dm2"          
  ## [26] "BSgenome.Dmelanogaster.UCSC.dm2.masked"   
  ## [27] "BSgenome.Dmelanogaster.UCSC.dm3"          
  ## [28] "BSgenome.Dmelanogaster.UCSC.dm3.masked"   
  ## [29] "BSgenome.Dmelanogaster.UCSC.dm6"          
  ## [30] "BSgenome.Drerio.UCSC.danRer10"            
  ## [31] "BSgenome.Drerio.UCSC.danRer5"             
  ## [32] "BSgenome.Drerio.UCSC.danRer5.masked"      
  ## [33] "BSgenome.Drerio.UCSC.danRer6"             
  ## [34] "BSgenome.Drerio.UCSC.danRer6.masked"      
  ## [35] "BSgenome.Drerio.UCSC.danRer7"             
  ## [36] "BSgenome.Drerio.UCSC.danRer7.masked"      
  ## [37] "BSgenome.Ecoli.NCBI.20080805"             
  ## [38] "BSgenome.Gaculeatus.UCSC.gasAcu1"         
  ## [39] "BSgenome.Gaculeatus.UCSC.gasAcu1.masked"  
  ## [40] "BSgenome.Ggallus.UCSC.galGal3"            
  ## [41] "BSgenome.Ggallus.UCSC.galGal3.masked"     
  ## [42] "BSgenome.Ggallus.UCSC.galGal4"            
  ## [43] "BSgenome.Ggallus.UCSC.galGal4.masked"     
  ## [44] "BSgenome.Ggallus.UCSC.galGal5"            
  ## [45] "BSgenome.Hsapiens.1000genomes.hs37d5"     
  ## [46] "BSgenome.Hsapiens.NCBI.GRCh38"            
  ## [47] "BSgenome.Hsapiens.UCSC.hg17"              
  ## [48] "BSgenome.Hsapiens.UCSC.hg17.masked"       
  ## [49] "BSgenome.Hsapiens.UCSC.hg18"              
  ## [50] "BSgenome.Hsapiens.UCSC.hg18.masked"       
  ## [51] "BSgenome.Hsapiens.UCSC.hg19"              
  ## [52] "BSgenome.Hsapiens.UCSC.hg19.masked"       
  ## [53] "BSgenome.Hsapiens.UCSC.hg38"              
  ## [54] "BSgenome.Hsapiens.UCSC.hg38.masked"       
  ## [55] "BSgenome.Mdomestica.UCSC.monDom5"         
  ## [56] "BSgenome.Mfascicularis.NCBI.5.0"          
  ## [57] "BSgenome.Mfuro.UCSC.musFur1"              
  ## [58] "BSgenome.Mmulatta.UCSC.rheMac2"           
  ## [59] "BSgenome.Mmulatta.UCSC.rheMac2.masked"    
  ## [60] "BSgenome.Mmulatta.UCSC.rheMac3"           
  ## [61] "BSgenome.Mmulatta.UCSC.rheMac3.masked"    
  ## [62] "BSgenome.Mmulatta.UCSC.rheMac8"           
  ## [63] "BSgenome.Mmusculus.UCSC.mm10"             
  ## [64] "BSgenome.Mmusculus.UCSC.mm10.masked"      
  ## [65] "BSgenome.Mmusculus.UCSC.mm8"              
  ## [66] "BSgenome.Mmusculus.UCSC.mm8.masked"       
  ## [67] "BSgenome.Mmusculus.UCSC.mm9"              
  ## [68] "BSgenome.Mmusculus.UCSC.mm9.masked"       
  ## [69] "BSgenome.Osativa.MSU.MSU7"                
  ## [70] "BSgenome.Ptroglodytes.UCSC.panTro2"       
  ## [71] "BSgenome.Ptroglodytes.UCSC.panTro2.masked"
  ## [72] "BSgenome.Ptroglodytes.UCSC.panTro3"       
  ## [73] "BSgenome.Ptroglodytes.UCSC.panTro3.masked"
  ## [74] "BSgenome.Ptroglodytes.UCSC.panTro5"       
  ## [75] "BSgenome.Ptroglodytes.UCSC.panTro6"       
  ## [76] "BSgenome.Rnorvegicus.UCSC.rn4"            
  ## [77] "BSgenome.Rnorvegicus.UCSC.rn4.masked"     
  ## [78] "BSgenome.Rnorvegicus.UCSC.rn5"            
  ## [79] "BSgenome.Rnorvegicus.UCSC.rn5.masked"     
  ## [80] "BSgenome.Rnorvegicus.UCSC.rn6"            
  ## [81] "BSgenome.Scerevisiae.UCSC.sacCer1"        
  ## [82] "BSgenome.Scerevisiae.UCSC.sacCer2"        
  ## [83] "BSgenome.Scerevisiae.UCSC.sacCer3"        
  ## [84] "BSgenome.Sscrofa.UCSC.susScr11"           
  ## [85] "BSgenome.Sscrofa.UCSC.susScr3"            
  ## [86] "BSgenome.Sscrofa.UCSC.susScr3.masked"     
  ## [87] "BSgenome.Tgondii.ToxoDB.7.0"              
  ## [88] "BSgenome.Tguttata.UCSC.taeGut1"           
  ## [89] "BSgenome.Tguttata.UCSC.taeGut1.masked"    
  ## [90] "BSgenome.Tguttata.UCSC.taeGut2"           
  ## [91] "BSgenome.Vvinifera.URGI.IGGP12Xv0"        
  ## [92] "BSgenome.Vvinifera.URGI.IGGP12Xv2"        
  ## [93] "BSgenome.Vvinifera.URGI.IGGP8X"

  genomeName <- "BSgenome.Hsapiens.UCSC.hg19"

In this example, the BSgenome package "BSgenome.Hsapiens.UCSC.hg19" refers to an unmasked genome; alignment index and alignments will be performed on the full unmasked genome sequence (recommended). If using a masked genome (e.g. "BSgenome.Hsapiens.UCSC.hg19.masked"), masked regions will be replaced with "N" bases, and this hard-masked version of the genome will be used for creating the alignment index and further alignments.

• a file name, referring to a sequence file containing one or several reference sequences (e.g. chromosomes) in fasta format:

  genomeFile <- system.file(package="QuasR", "extdata", "hg19sub.fa")

5.4 Sequence data pre-processing

The preprocessReads function can be used to prepare the input sequence files prior to alignment. The function takes one or several sequence files (or pairs of files for a paired-end experiment) in fasta or fastq format as input and produces the same number of output files with the processed reads.

In the following example, we truncate the reads by removing the three bases from the 3’-end (the right side), remove the adapter sequence AAAAAAAAAA from the 5’-end (the left side) and filter out reads that, after truncation and adapter removal, are shorter than 14 bases or contain more than 2 N bases:

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)

##   filtering /tmp/RtmpicCqbg/Rinst766929a7604f/QuasR/extdata/rna_1_1.fq.bz2

##   filtering /tmp/RtmpicCqbg/Rinst766929a7604f/QuasR/extdata/rna_2_1.fq.bz2

res

##                  rna_1_1.fq.bz2 rna_2_1.fq.bz2
## totalSequences             3002           3000
## matchTo5pAdapter            466            463
## matchTo3pAdapter              0              0
## tooShort                    107             91
## tooManyN                      0              0
## lowComplexity                 0              0
## totalPassed                2895           2909

unlink(outfiles)

preprocessReads returns a matrix with a summary of the pre-processing. The matrix contains one column per (pair of) input sequence files, and contains the total number of reads (totalSequences), the number of reads that matched to the five prime or three prime adapters (matchTo5pAdapter and matchTo3pAdapter), the number of reads that were too short (tooShort), contained too many non-base characters (tooManyN) or were of low sequence complexity (lowComplexity, deactivated by default). Finally, the number of reads that passed the filtering steps is reported in the last row (totalPassed).

In the example below we process paired-end reads, removing all pairs with one or several N bases. Even if only one sequence in a pair fulfills the filtering criteria, both reads in the pair are removed, thereby preserving the matching order of the sequences in the two files:

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)

##   filtering /tmp/RtmpicCqbg/Rinst766929a7604f/QuasR/ex.../rna_1_1.fq.bz2 and
##     /tmp/RtmpicCqbg/Rinst766929a7604f/QuasR/ext.../rna_1_2.fq.bz2

res

##                  rna_1_1.fq.bz2:rna_1_2.fq.bz2
## totalSequences                            3002
## matchTo5pAdapter                            NA
## matchTo3pAdapter                            NA
## tooShort                                     0
## tooManyN                                     3
## lowComplexity                                0
## totalPassed                               2999

More details on the preprocessReads function can be found in the function documentation (see ?preprocessReads) or in the section 7.1.

6.1 ChIP-seq: Measuring protein-DNA binding and chromatin modifications

Here we show an exemplary single-end ChIP-seq workflow using a small number of reads from a histone 3 lysine 4 trimethyl (H3K4me3) ChIP-seq experiment. This histone modification is known to locate to genomic regions with a high density of CpG dinucleotides (so called CpG islands); about 60% of mammalian genes have such a CpG island close to their transcript start site. All necessary files are included in the QuasR package, and we start the example workflow by copying those files into the current working directly, into a subfolder called "extdata":

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

sampleFile <- "extdata/samples_chip_single.txt"
auxFile <- "extdata/auxiliaries.txt"
genomeFile <- "extdata/hg19sub.fa"

proj1 <- qAlign(sampleFile, genome=genomeFile, auxiliaryFile=auxFile)

## alignment files missing - need to:

##     create 2 auxiliary alignment(s)

## Creating an Rbowtie index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/NC_001422.1.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:

## nodeNames
## malbec2 
##       1

## Performing auxiliary alignments for 2 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a73484c6d.txt
## Auxiliary alignments have been created successfully

proj1

## Project: qProject
##  Options   : maxHits         : 1
##              paired          : no
##              splicedAlignment: FALSE
##              bisulfite       : no
##              snpFile         : none
##              geneAnnotation  : none
##  Aligner   : Rbowtie v1.24.0 (parameters: -m 1 --best --strata)
##  Genome    : /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vi.../hg19sub.fa (file)
## 
##  Reads     : 2 files, 2 samples (fastq format):
##    1. chip_1_1.fq.bz2  Sample1 (phred33)
##    2. chip_2_1.fq.bz2  Sample2 (phred33)
## 
##  Genome alignments: directory: same as reads
##    1. chip_1_1_d0a6130e4bf.bam
##    2. chip_2_1_d0a723b4f4d.bam
## 
##  Aux. alignments: 1 file, directory: same as reads
##    a. /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vig.../NC_001422.1.fa  phiX174
##      1. chip_1_1_d0a74af2d5.bam
##      2. chip_2_1_d0acc99101.bam

list.files("extdata", pattern=".bam$")

## [1] "chip_1_1_d0a6130e4bf.bam"      "chip_1_1_d0a74af2d5.bam"      
## [3] "chip_2_1_d0a723b4f4d.bam"      "chip_2_1_d0acc99101.bam"      
## [5] "phiX_paired_withSecondary.bam"

list.files("extdata", pattern="^chip_1_1_")[1:3]

## [1] "chip_1_1_d0a6130e4bf.bam"     "chip_1_1_d0a6130e4bf.bam.bai"
## [3] "chip_1_1_d0a6130e4bf.bam.txt"

## collecting quality control data

## creating QC plots

qQCReport(proj1, pdfFilename="extdata/qc_report.pdf")

## collecting quality control data

## creating QC plots

alignmentStats(proj1)

##                 seqlength mapped unmapped
## Sample1:genome      95000   2339      258
## Sample2:genome      95000   3609      505
## Sample1:phiX174      5386    251        7
## Sample2:phiX174      5386    493       12

qExportWig(proj1, binsize=100L, scaling=TRUE, collapseBySample=TRUE)

## collecting mapping statistics for scaling...done
## start creating wig files...
##   Sample1.wig.gz (Sample1)
##   Sample2.wig.gz (Sample2)
## done

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

## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OK

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)

## GRanges object with 6 ranges and 2 metadata columns:
##                   seqnames      ranges strand |         gene_id     tx_id
##                      <Rle>   <IRanges>  <Rle> | <CharacterList> <integer>
##   ENSG00000176022     chr1 31629-33128      + | ENSG00000176022         1
##   ENSG00000186891     chr1   6452-7951      - | ENSG00000186891         2
##   ENSG00000186827     chr1 14013-15512      - | ENSG00000186827         6
##   ENSG00000078808     chr1 31882-33381      - | ENSG00000078808         9
##   ENSG00000171863     chr2   1795-3294      + | ENSG00000171863        17
##   ENSG00000252531     chr2   7160-8659      + | ENSG00000252531        26
##   -------
##   seqinfo: 3 sequences from an unspecified genome

cnt <- qCount(proj1, promRegSel)

## counting alignments...done

cnt

##                 width Sample1 Sample2
## ENSG00000176022  1500     157     701
## ENSG00000186891  1500      22       5
## ENSG00000186827  1500      10       3
## ENSG00000078808  1500      73     558
## ENSG00000171863  1500      94     339
## ENSG00000252531  1500      59       9
## ENSG00000247886  1500     172     971
## ENSG00000254999  1500     137     389
## ENSG00000238642  1500       8       3
## ENSG00000134086  1500       9      18
## ENSG00000238345  1500      13      25
## ENSG00000134075  1500       7       3

gr1 <- import("Sample1.wig.gz")

## Warning in asMethod(object): NAs introduced by coercion

gr2 <- import("Sample2.wig.gz")

## Warning in asMethod(object): NAs introduced by coercion

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)

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)

## GRanges object with 6 ranges and 1 metadata column:
##       seqnames      ranges strand |   gene_name
##          <Rle>   <IRanges>  <Rle> | <character>
##   TSS     chr1   6949-6951      - |    TNFRSF18
##   TSS     chr1 14505-14507      - |     TNFRSF4
##   TSS     chr1 29171-29173      - |        SDF4
##   TSS     chr1 32659-32661      + |     B3GALT6
##   TSS     chr2   3200-3202      + |        RPS7
##   TSS     chr3   2386-2388      + |     C3orf10
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

prS <- qProfile(proj1, tssRegions, upstream=3000, downstream=3000, 
                orientation="same")

## profiling alignments...done

prO <- qProfile(proj1, tssRegions, upstream=3000, downstream=3000, 
                orientation="opposite")

## profiling alignments...done

lapply(prS, "[", , 1:10)

## $coverage
## -3000 -2999 -2998 -2997 -2996 -2995 -2994 -2993 -2992 -2991 
##     8     8     8     8     8     8     8     8     8     8 
## 
## $Sample1
## -3000 -2999 -2998 -2997 -2996 -2995 -2994 -2993 -2992 -2991 
##     1     0     0     0     0     0     0     0     0     0 
## 
## $Sample2
## -3000 -2999 -2998 -2997 -2996 -2995 -2994 -2993 -2992 -2991 
##     0     0     0     2     0     0     1     1     1     0

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)

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

6.2 RNA-seq: Gene expression profiling

In QuasR, an analysis workflow for an RNA-seq dataset is very similar to the one described above for a ChIP-seq experiment. The major difference is that here reads are aligned using qAlign(..., splicedAlignment=TRUE, aligner="Rhisat2"), which will cause qAlign to align reads with the HISAT2 aligner (Kim, Langmead, and Salzberg 2015) (via the Rhisat2 package), rather than with bowtie (Langmead et al. 2009). Before the Rhisat2 package was available (introduced in Bioconductor 3.9), qAlign(... splicedAlignment=TRUE) aligned reads using SpliceMap (Au et al. 2010), which is not recommended now but still possible in order to reproduce old results. Spliced paired-end alignments are also supported; the splicedAlignment argument can be freely combined with the paired argument. In addition, HISAT2 also allows the specification of known splice sites, which can help in the read alignment. This is done by specifying the argument geneAnnotation in qAlign(), to either a .gtf file or a sqlite database generated by exporting a TxDb object.

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

sampleFile <- "extdata/samples_rna_paired.txt"
genomeFile <- "extdata/hg19sub.fa"

proj2 <- qAlign(sampleFile, genome=genomeFile, splicedAlignment=TRUE, aligner="Rhisat2")

## alignment files missing - need to:

##     create alignment index for the genome

##     create 2 genomic alignment(s)

## Creating an Rhisat2 index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 2 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a798faa58.txt
## Genomic alignments have been created successfully

proj2

## Project: qProject
##  Options   : maxHits         : 1
##              paired          : fr
##              splicedAlignment: TRUE
##              bisulfite       : no
##              snpFile         : none
##              geneAnnotation  : none
##  Aligner   : Rhisat2 v1.0.1 (parameters: -k 2)
##  Genome    : /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vi.../hg19sub.fa (file)
## 
##  Reads     : 2 pairs of files, 2 samples (fastq format):
##    1. rna_1_1.fq.bz2  rna_1_2.fq.bz2  Sample1 (phred33)
##    2. rna_2_1.fq.bz2  rna_2_2.fq.bz2  Sample2 (phred33)
## 
##  Genome alignments: directory: same as reads
##    1. rna_1_1_d0a32af3a9f.bam
##    2. rna_2_1_d0ab1f2c50.bam 
## 
##  Aux. alignments: none

proj2unspl <- qAlign(sampleFile, genome=genomeFile,
                     splicedAlignment=FALSE)

## alignment files missing - need to:

##     create 2 genomic alignment(s)

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 2 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a7d49de30.txt
## Genomic alignments have been created successfully

alignmentStats(proj2)

##                seqlength mapped unmapped
## Sample1:genome     95000   5961        6
## Sample2:genome     95000   5914        2

alignmentStats(proj2unspl)

##                seqlength mapped unmapped
## Sample1:genome     95000   2258     3746
## Sample2:genome     95000   2652     3348

geneLevels <- qCount(proj2, txdb, reportLevel="gene")

## extracting gene regions from TxDb...done
## counting alignments...done
## collapsing counts by query name...done

exonLevels <- qCount(proj2, txdb, reportLevel="exon")

## extracting exon regions from TxDb...done
## counting alignments...done

head(geneLevels)

##                 width Sample1 Sample2
## ENSG00000078808  4697     710    1083
## ENSG00000134075   589    1173    1303
## ENSG00000134086  4213     279     295
## ENSG00000171863  5583    2924    2224
## ENSG00000176022  2793      62     344
## ENSG00000186827  1721      37       8

head(exonLevels)

##    width Sample1 Sample2
## 1   2793      62     344
## 10   187       3       0
## 11   307       3       0
## 12   300      11       2
## 13   493      19       2
## 14   129       7       0

geneRPKM <- t(t(geneLevels[,-1] /geneLevels[,1] *1000)
              /colSums(geneLevels[,-1]) *1e6)
geneRPKM

##                 Sample1 Sample2
## ENSG00000078808   21350   31786
## ENSG00000134075  281287  304966
## ENSG00000134086    9354    9653
## ENSG00000171863   73974   54915
## ENSG00000176022    3135   16979
## ENSG00000186827    3037     641
## ENSG00000186891    2681     201
## ENSG00000238345       0       0
## ENSG00000238642       0       0
## ENSG00000247886       0       0
## ENSG00000252531    6066    1691
## ENSG00000254999  213296  222826

exonJunctions <- qCount(proj2, NULL, reportLevel="junction")

## counting junctions...done

exonJunctions

## GRanges object with 46 ranges and 2 metadata columns:
##        seqnames      ranges strand |   Sample1   Sample2
##           <Rle>   <IRanges>  <Rle> | <numeric> <numeric>
##    [1]     chr1 12213-12321      + |         3         0
##    [2]     chr1 13085-13371      - |         1         0
##    [3]     chr1 18069-18837      + |         9        16
##    [4]     chr1 18069-18837      - |         7         4
##    [5]     chr1 18185-18837      - |         2         0
##    ...      ...         ...    ... .       ...       ...
##   [42]     chr1 14166-14362      + |         0         1
##   [43]     chr1 19308-23623      - |         0         2
##   [44]     chr1 29327-32271      + |         0         2
##   [45]     chr1 29327-32271      - |         0         1
##   [46]     chr3   2504-5589      - |         0         3
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

knownIntrons <- unlist(intronsByTranscript(txdb))
isKnown <- overlapsAny(exonJunctions, knownIntrons, type="equal")
table(isKnown)

## isKnown
## FALSE  TRUE 
##    25    21

tapply(rowSums(as.matrix(mcols(exonJunctions))),
       isKnown, summary)

## $`FALSE`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       3       7      47      31     340 
## 
## $`TRUE`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     2.0    16.0    50.6    91.0   210.0

exonBodyLevels <- qCount(proj2, txdb, reportLevel="exon", includeSpliced=FALSE)

## extracting exon regions from TxDb...done
## counting alignments...done

summary(exonLevels - exonBodyLevels)

##      width      Sample1       Sample2   
##  Min.   :0   Min.   :  0   Min.   :  0  
##  1st Qu.:0   1st Qu.:  0   1st Qu.:  0  
##  Median :0   Median :  3   Median :  1  
##  Mean   :0   Mean   : 42   Mean   : 35  
##  3rd Qu.:0   3rd Qu.: 48   3rd Qu.: 50  
##  Max.   :0   Max.   :818   Max.   :647

## collecting quality control data

## creating QC plots

6.3 smRNA-seq: small RNA and miRNA expression profiling

Expression profiling of miRNAs differs only slightly from the profiling of mRNAs. There are a few details that need special care, which are outlined in this section.

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

# prepare sample file with processed reads filenames
sampleFile <- file.path("extdata", "samples_mirna.txt")
sampleFile

## [1] "extdata/samples_mirna.txt"

sampleFile2 <- sub(".txt", "_processed.txt", sampleFile)
sampleFile2

## [1] "extdata/samples_mirna_processed.txt"

tab <- read.delim(sampleFile, header=TRUE, as.is=TRUE)
tab

##     FileName SampleName
## 1 mirna_1.fa     miRNAs

tab2 <- tab
tab2$FileName <- sub(".fa", "_processed.fa", tab$FileName)
write.table(tab2, sampleFile2, sep="\t", quote=FALSE, row.names=FALSE)
tab2

##               FileName SampleName
## 1 mirna_1_processed.fa     miRNAs

# remove adapters
oldwd <- setwd(dirname(sampleFile))
res <- preprocessReads(tab$FileName,
                       tab2$FileName,
                       Rpattern="TGGAATTCTCGGGTGCCAAGG")

##   filtering mirna_1.fa

res

##                  mirna_1.fa
## totalSequences         1000
## matchTo5pAdapter          0
## matchTo3pAdapter       1000
## tooShort                  0
## tooManyN                  0
## lowComplexity             0
## totalPassed            1000

setwd(oldwd)

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

proj3 <- qAlign(sampleFile2, genomeFile, maxHits=50)

## alignment files missing - need to:

##     create 1 genomic alignment(s)

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 1 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a154bf20b.txt
## Genomic alignments have been created successfully

alignmentStats(proj3)

##               seqlength mapped unmapped
## miRNAs:genome     95000   1000        0

mirs <- import("extdata/mirbaseXX_qsr.gff3")
names(mirs) <- mirs$Name
preMirs <- mirs[ mirs$type=="miRNA_primary_transcript" ]
matureMirs <- mirs[ mirs$type=="miRNA" ]

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(end(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)

matureMirsExtended <- resize(matureMirs, width=1L, fix="start") + 3L

# quantify mature miRNAs
cnt <- qCount(proj3, matureMirsExtended, orientation="same")

## counting alignments...done

cnt

##                 width miRNAs
## qsr-miR-9876-5p     7     13
## qsr-miR-9876-3p     7    984

# quantify pre-miRNAs
cnt <- qCount(proj3, preMirs, orientation="same")

## counting alignments...done

cnt

##              width miRNAs
## qsr-mir-9876    75   1000

6.4 Bis-seq: Measuring DNA methylation

Sequencing of bisulfite-converted genomic DNA allows detection of methylated cytosines, which in mammalian genomes typically occur in the context of CpG dinucleotides. The treatment of DNA with bisulfite induces deamination of non-methylated cytosines, converting them to uracils. Sequencing and aligning of such bisulfite-converted DNA results in C-to-T mismatches. Both alignment of converted reads, as well as the interpretation of the alignments for calculation of methylation levels require specific approaches and are supported in QuasR by qAlign (bisulfite argument, section 7.2) and qMeth (section 7.9), respectively.

We start the analysis by copying the example data files into the current working directly, into a subfolder called "extdata". Then, bisulfite-specific alignment is selected in qAlign by setting bisulfite to "dir" for a directional experiment, or to "undir" for an undirectional Bis-seq experiment:

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

sampleFile <- "extdata/samples_bis_single.txt"
genomeFile <- "extdata/hg19sub.fa"

proj4 <- qAlign(sampleFile, genomeFile, bisulfite="dir")

## alignment files missing - need to:

##     create alignment index for the genome

##     create 1 genomic alignment(s)

## Creating an RbowtieCtoT index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 1 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a4d852f47.txt
## Genomic alignments have been created successfully

proj4

## Project: qProject
##  Options   : maxHits         : 1
##              paired          : no
##              splicedAlignment: FALSE
##              bisulfite       : dir
##              snpFile         : none
##              geneAnnotation  : none
##  Aligner   : Rbowtie v1.24.0 (parameters: -k 2 --best --strata -v 2)
##  Genome    : /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vi.../hg19sub.fa (file)
## 
##  Reads     : 1 file, 1 sample (fasta format):
##    1. bis_1_1.fa.bz2  Sample1
## 
##  Genome alignments: directory: same as reads
##    1. bis_1_1_d0a2c42c2f7.bam
## 
##  Aux. alignments: none

The resulting alignments are not different from those of non-Bis-seq experiments, apart from the fact that they may contain many C-to-T (or A-to-G) mismatches that are not counted as mismatches when aligning the reads. The number of alignments in specific genomic regions could be quantified using qCount as with ChIP-seq or RNA-seq experiments. For quantification of methylation the qMeth function is used:

meth <- qMeth(proj4, mode="CpGcomb", collapseBySample=TRUE)
meth

## GRanges object with 3110 ranges and 2 metadata columns:
##          seqnames      ranges strand | Sample1_T Sample1_M
##             <Rle>   <IRanges>  <Rle> | <integer> <integer>
##      [1]     chr1       19-20      * |         1         1
##      [2]     chr1       21-22      * |         1         1
##      [3]     chr1       54-55      * |         3         1
##      [4]     chr1       57-58      * |         3         0
##      [5]     chr1       80-81      * |         6         5
##      ...      ...         ...    ... .       ...       ...
##   [3106]     chr3 44957-44958      * |         8         7
##   [3107]     chr3 44977-44978      * |         5         3
##   [3108]     chr3 44981-44982      * |         4         3
##   [3109]     chr3 44989-44990      * |         1         1
##   [3110]     chr3 44993-44994      * |         1         1
##   -------
##   seqinfo: 3 sequences from an unspecified genome

By default, qMeth quantifies methylation for all cytosines in CpG contexts, combining the data from plus and minus strands (mode="CpGcomb"). The results are returned as a GRanges object with coordinates of each CpG, and two metadata columns for each input sequence file in the qProject object. These two columns contain the total number of aligned reads that overlap a given CpG (C-to-C matches or C-to-T mismatches, suffix _T in the column name), and the number of read alignments that had a C-to-C match at that position (methylated events, suffix _M).

Independent of the number of alignments, the returned object will list all CpGs in the target genome including the ones that have zero coverage, unless you set keepZero=FALSE:

chrs <- readDNAStringSet(genomeFile)
sum(vcountPattern("CG",chrs))

## [1] 3110

length(qMeth(proj4))

## [1] 3110

length(qMeth(proj4, keepZero=FALSE))

## [1] 2929

The fraction methylation can easily be obtained as the ratio between _M and _T columns:

percMeth <- mcols(meth)[,2] *100 /mcols(meth)[,1]
summary(percMeth)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0    75.0    90.9    75.4   100.0   100.0     181

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)

If qMeth is called without a query argument, it will by default return methylation states for each C or CpG in the genome. Using a query argument it is possible to restrict the analysis to specific genomic regions, and if using in addition collapseByQueryRegion=TRUE, the single base methylation states will further be combined for all C’s that are contained in the same query region:

qMeth(proj4, query=GRanges("chr1",IRanges(start=31633,width=2)),
      collapseBySample=TRUE)

## GRanges object with 1 range and 2 metadata columns:
##           seqnames      ranges strand | Sample1_T Sample1_M
##              <Rle>   <IRanges>  <Rle> | <integer> <integer>
##   CpGcomb     chr1 31633-31634      * |        10         2
##   -------
##   seqinfo: 3 sequences from an unspecified genome

qMeth(proj4, query=promRegSel, collapseByQueryRegion=TRUE,
      collapseBySample=TRUE)

## GRanges object with 12 ranges and 2 metadata columns:
##        seqnames      ranges strand | Sample1_T Sample1_M
##           <Rle>   <IRanges>  <Rle> | <numeric> <numeric>
##    [1]     chr1 31629-33128      + |       426        74
##    [2]     chr1   6452-7951      - |       388       244
##    [3]     chr1 14013-15512      - |       627       560
##    [4]     chr1 31882-33381      - |       522       232
##    [5]     chr2   1795-3294      + |       997       539
##    ...      ...         ...    ... .       ...       ...
##    [8]     chr3   1276-2775      + |       715       253
##    [9]     chr3 19069-20568      + |       253       204
##   [10]     chr3 26692-28191      + |       934       818
##   [11]     chr3 26834-28333      + |       934       777
##   [12]     chr3 13102-14601      - |       307       287
##   -------
##   seqinfo: 3 sequences from an unspecified genome

Finally, qMeth allows the retrieval of methylation states for individual molecules (per alignment). This is done by using a query containing a single genomic region (typically small, such as a PCR amplicon) and setting reportLevel="alignment". In that case, the return value of qMeth will be a list (over samples) of lists (with four elements giving the identities of alignment, C nucleotide, strand and the methylation state). See the documentation of qMeth for more details.

6.5 Allele-specific analysis

All experiment types supported by QuasR (ChIP-seq, RNA-seq and Bis-seq; only alignments to the genome, but not to auxiliaries) can also be analyzed in an allele-specific manner. For this, a file containing genomic location and the two alleles of known sequence polymorphisms has to be provided to the snpFile argument of qAlign. The file is in tab-delimited text format without a header and contains four columns with chromosome name, position, reference allele and alternative allele.

Below is an example of a SNP file, also available from system.file(package="QuasR", "extdata", "hg19sub_snp.txt"):

chr1    3199    C   T
chr1    3277    C   T
chr1    4162    C   T
chr1    4195    C   T
...

For a given locus, either reference or alternative allele may but does not have to be identical to the sequence of the reference genome (genomeFile in this case). To avoid an alignment bias, all reads are aligned separately to each of the two new genomes, which QuasR generates by injecting the SNPs listed in snpFile into the reference genome. Finally, the two alignment files are combined, only retaining the best alignment for each read. While this procedure takes more than twice as long as aligning against a single genome, it has the advantage to correctly align reads even in regions of high SNP density and has been shown to produce more accurate results.

While combining alignments, each read is classified into one of three groups (stored in the bam files under the XV tag):

• R: the read aligned better to the reference genome
  
• U: the read aligned equally well to both genomes (unknown origin)
  
• A: the read aligned better to the alternative genome

Using these alignment classifications, the qCount and qMeth functions will produce three counts instead of a single count; one for each class. The column names will be suffixed by _R, _U and _A.

The examples below use data provided with the QuasR package, which is first copied to the current working directory, into a subfolder called "extdata":

file.copy(system.file(package="QuasR", "extdata"), ".", recursive=TRUE)

## [1] TRUE

The example below aligns the same reads that were also used in the ChIP-seq workflow (section 6.1), but this time using a snpFile:

sampleFile <- "extdata/samples_chip_single.txt"
genomeFile <- "extdata/hg19sub.fa"
snpFile <- "extdata/hg19sub_snp.txt"
proj1SNP <- qAlign(sampleFile, genome=genomeFile, snpFile=snpFile)

## alignment files missing - need to:

##     create alignment index for the genome

##     create 2 genomic alignment(s)

## Reading and processing the SNP file: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt

## Creating the genome fasta file containing the SNPs: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.R.fa

## Creating the genome fasta file containing the SNPs: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.A.fa

## Creating a .fai file for the snp genome: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.R.fa

## Creating a .fai file for the snp genome: /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.A.fa

## Creating an Rbowtie index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.R.fa

## Finished creating index

## Creating an Rbowtie index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.A.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 2 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a3a628d66.txt
## Genomic alignments have been created successfully

proj1SNP

## Project: qProject
##  Options   : maxHits         : 1
##              paired          : no
##              splicedAlignment: FALSE
##              bisulfite       : no
##              snpFile         : /tmp/RtmpicCqbg/Rbuild7669.../hg19sub_snp.txt
##              geneAnnotation  : none
##  Aligner   : Rbowtie v1.24.0 (parameters: -k 2 --best --strata -v 2)
##  Genome    : /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vi.../hg19sub.fa (file)
## 
##  Reads     : 2 files, 2 samples (fastq format):
##    1. chip_1_1.fq.bz2  Sample1 (phred33)
##    2. chip_2_1.fq.bz2  Sample2 (phred33)
## 
##  Genome alignments: directory: same as reads
##    1. chip_1_1_d0a656e8066.bam
##    2. chip_2_1_d0a5b1c40d2.bam
## 
##  Aux. alignments: none

Instead of one count per promoter region and sample, qCount now returns three (promRegSel was generated in the ChIP-seq example workflow):

head(qCount(proj1, promRegSel))

## counting alignments...done

##                 width Sample1 Sample2
## ENSG00000176022  1500     157     701
## ENSG00000186891  1500      22       5
## ENSG00000186827  1500      10       3
## ENSG00000078808  1500      73     558
## ENSG00000171863  1500      94     339
## ENSG00000252531  1500      59       9

head(qCount(proj1SNP, promRegSel))

## counting alignments...done

##                 width Sample1_R Sample1_U Sample1_A Sample2_R Sample2_U
## ENSG00000176022  1500         0       133         0         0       559
## ENSG00000186891  1500         4        16         0         0         5
## ENSG00000186827  1500         2         8         0         0         2
## ENSG00000078808  1500         0        59         0         0       432
## ENSG00000171863  1500         4        78         0         8       263
## ENSG00000252531  1500         3        50         2         0         6
##                 Sample2_A
## ENSG00000176022         0
## ENSG00000186891         0
## ENSG00000186827         0
## ENSG00000078808         0
## ENSG00000171863         0
## ENSG00000252531         0

The example below illustrates use of a snpFile for Bis-seq experiments. Similarly as for qCount, the count types are labeled by R, U and A. These labels are added to the total and methylated column suffixes _T and _M, resulting in a total of six instead of two counts per feature and sample:

sampleFile <- "extdata/samples_bis_single.txt"
genomeFile <- "extdata/hg19sub.fa"
proj4SNP <- qAlign(sampleFile, genomeFile,
                   snpFile=snpFile, bisulfite="dir")

## alignment files missing - need to:

##     create alignment index for the genome

##     create 1 genomic alignment(s)

## Creating an RbowtieCtoT index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.R.fa

## Finished creating index

## Creating an RbowtieCtoT index for /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/extdata/hg19sub_snp.txt.hg19sub.fa.A.fa

## Finished creating index

## Testing the compute nodes...OK
## Loading QuasR on the compute nodes...OK
## Available cores:
## malbec2: 1
## Performing genomic alignments for 1 samples. See progress in the log file:
## /tmp/RtmpicCqbg/Rbuild7669357602ca/QuasR/vignettes/QuasR_log_d0a4631da66.txt
## Genomic alignments have been created successfully

head(qMeth(proj4SNP, mode="CpGcomb", collapseBySample=TRUE))

## GRanges object with 6 ranges and 6 metadata columns:
##       seqnames    ranges strand | Sample1_TR Sample1_MR Sample1_TU
##          <Rle> <IRanges>  <Rle> |  <integer>  <integer>  <integer>
##   [1]     chr1     19-20      * |          0          0          1
##   [2]     chr1     21-22      * |          0          0          1
##   [3]     chr1     54-55      * |          0          0          3
##   [4]     chr1     57-58      * |          0          0          3
##   [5]     chr1     80-81      * |          0          0          6
##   [6]     chr1   103-104      * |          0          0          6
##       Sample1_MU Sample1_TA Sample1_MA
##        <integer>  <integer>  <integer>
##   [1]          1          0          0
##   [2]          1          0          0
##   [3]          1          0          0
##   [4]          0          0          0
##   [5]          5          0          0
##   [6]          5          0          0
##   -------
##   seqinfo: 3 sequences from an unspecified genome

7.1 preprocessReads

The preprocessReads function can be used to prepare the input sequences before alignment to the reference genome, for example to filter out low quality reads unlikely to produce informative alignments. When working with paired-end experiments, the paired reads are expected to be contained in identical order in two separate files. For this reason, both reads of a pair are filtered out if any of the two reads fulfills the filtering criteria. The following types of filtering tasks can be performed (in the order as listed):

1. Truncate reads: remove nucleotides from the start and/or end of each read.
   
2. Trim adapters: remove nucleotides at the beginning and/or end of each read that match to a defined (adapter) sequence. The adapter trimming is done by calling trimLRPatterns from the Biostrings package (Pages et al., n.d.).
   
3. Filter out low quality reads: Filter out reads that fulfill any of the filtering criteria (contain more than nBases N bases, are shorter than minLength or have a dinucleotide complexity of less than complexity-times the average complexity of the human genome sequence).

The dinucleotide complexity is calculated in bits as Shannon entropy using the following formula \(-\sum_i f_i \cdot \log_2 f_i\), where \(f_i\) is the frequency of dinucleotide \(i\) (\(i=1, 2, ..., 16\)).

7.2 qAlign

qAlign is the function that generates alignment files in bam format, for all input sequence files listed in sampleFile (see section 5.1), against the reference genome (genome argument), and for reads that do not match to the reference genome, against one or several auxiliary target sequences (auxiliaryFile, see section 5.2).

The reference genome can be provided either as a fasta sequence file or as a BSgenome package name (see section 5.3). If a BSgenome package is not found in the installed packages but available from Bioconductor, it will be automatically downloaded.

The alignment program is set by aligner, and parameters by maxHits, paired, splicedAlignment and alignmentParameter. Currently, aligner can only be set to "Rbowtie", which is a wrapper for bowtie (Langmead et al. 2009) and SpliceMap (Au et al. 2010), or "Rhisat2", which is a wrapper for HISAT2 (Kim, Langmead, and Salzberg 2015). When aligner="Rbowtie", SpliceMap will be used if splicedAlignment=TRUE (not recommended anymore except for reproducing older analyses). However, it is recommended to create spliced alignment using splicedAlignment=TRUE, aligner="Rhisat2", which will use the HISAT2 aligner and typically leads to more sensistive alignments and shorter alignment times compared to SpliceMap. The alignment strategy is furthermore affected by the parameters snpFile (alignments to variant genomes containing sequence polymorphisms) and bisulfite (alignment of bisulfite-converted reads). Finally, clObj can be used to enable parallelized alignment, sorting and conversion to bam format.

For each input sequence file listed in sampleFile, one bam file with alignments to the reference genome will be generated, and an additional one for each auxiliary sequence file listed in auxiliaryFile. By default, these bam files are stored at the same location as the sequence files, unless a different location is specified under alignmentsDir. If compatible alignment files are found at this location, they will not be regenerated, which allows re-use of the same sequencing samples in multiple analysis projects by listing them in several project-specific sampleFiles.

The alignment process produces temporary files, such as decompressed input sequence files or raw alignment files before conversion to bam format, which can be several times the size of the input sequence files. These temporary files are stored in the directory specified by cacheDir, which defaults to the R process temporary directory returned by tempdir(). The location of tempdir() can be set using environment variables (see ?tempdir).

qAlign returns a qProject object that contains all file names and paths, as well as all alignment parameters necessary for further analysis (see section 7.3 for methods to access the information contained in a qProject object).

7.3 qProject class

The qProject objects are returned by qAlign and contain all information about a sequencing experiment needed for further analysis. It is the main argument passed to the functions that start with a q letter, such as qCount, qQCReport and qExportWig. Some information inside of a qProject object can be accessed by specific methods (in the examples below, x is a qProject object):

• length(x) gets the number of input files.
  
• genome(x) gets the reference genome as a character(1). The type of genome is stored as an attribute in attr(genome(x),"genomeFormat"): "BSgenome" indicates that genome(x) refers to the name of a BSgenome package, "file" indicates that it contains the path and file name of a genome in fasta format.
  
• auxiliaries(x) gets a data.frame with auxiliary target sequences. The data.frame has one row per auxiliary target file, and two columns “FileName” and “AuxName”.
  
• alignments(x) gets a list with two elements "genome" and "aux". "genome" contains a data.frame with length(x) rows and two columns "FileName" (containing the path to bam files with genomic alignments) and "SampleName". "aux" contains a data.frame with one row per auxiliary target file (with auxiliary names as row names), and length(x) columns (one per input sequence file).
  
• x[i] returns a qProject object instance with i input files, where i can be an NA-free logical, numeric, or character vector.

7.4 qQCReport

The qQCReport function samples a random subset of sequences and alignments from each sample or input file and generates a series of diagnostic plots for estimating data quality. The available plots vary depending on the types of available input (fasta, fastq, bam files or qProject object; paired-end or single-end). The plots below show the currently available plots as produced by the ChIP-seq example in section 6.1 (except for the fragment size distributions which are based on an unspliced alignment of paired-end RNA seq reads):

• Quality score boxplot shows the distribution of base quality values as a box plot for each position in the input sequence. The background color (green, orange or red) indicates ranges of high, intermediate and low qualities.

• Nucleotide frequency plot shows the frequency of A, C, G, T and N bases by position in the read.

• Duplication level plot shows for each sample the fraction of reads observed at different duplication levels (e.g. once, two-times, three-times, etc.). In addition, the most frequent sequences are listed.

• Mapping statistics shows fractions of reads that were (un)mappable to the reference genome.

• Library complexity shows fractions of unique read(-pair) alignment positions. Please note that this measure is not independent from the total number of reads in a library, and is best compared between libraries of similar sizes.

• Mismatch frequency shows the frequency and position (relative to the read sequence) of mismatches in the alignments against the reference genome.

• Mismatch types shows the frequency of read bases that caused mismatches in the alignments to the reference genome, separately for each genome base.

• Fragment size shows the distribution of fragment sizes inferred from aligned read pairs.
  

7.5 alignmentStats

alignmentStats is comparable to the idxstats function from Samtools; it returns the size of the target sequence, as well as the number of mapped and unmapped reads that are contained in an indexed bam file. The function works for arguments of type qProject, as well as on a string with one or several bam file names. There is however a small difference in the two that is illustrated in the following example, which uses the qProject object from the ChIP-seq workflow:

# using bam files
alignmentStats(alignments(proj1)$genome$FileName)

##                          seqlength mapped unmapped
## chip_1_1_d0a6130e4bf.bam     95000   2339      258
## chip_2_1_d0a723b4f4d.bam     95000   3609      505

alignmentStats(unlist(alignments(proj1)$aux))

##                         seqlength mapped unmapped
## chip_1_1_d0a74af2d5.bam      5386    251        0
## chip_2_1_d0acc99101.bam      5386    493        0

# using a qProject object
alignmentStats(proj1)

##                 seqlength mapped unmapped
## Sample1:genome      95000   2339      258
## Sample2:genome      95000   3609      505
## Sample1:phiX174      5386    251        7
## Sample2:phiX174      5386    493       12

If calling alignmentStats on the bam files directly as in the first two expressions of the above example, the returned numbers correspond exactly to what you would obtain by the idxstats function from Samtools, only that the latter would report them separately for each target sequence, while alignmentStats sums them for each bam file. These numbers correctly state that there are zero unmapped reads in the auxiliary bam files. However, if calling alignmentStats on a qProject object, it will report 7 and 12 unmapped reads in the auxiliary bam files. This is because alignmentStats is aware that unmapped reads are removed from auxiliary bam files by QuasR, but can be calculated from the total number of reads to be aligned to the auxiliary target, which equals the number of unmapped reads in the corresponding genomic bam file.

7.6 qExportWig

qExportWig creates fixed-step wig files (see (http://genome.ucsc.edu/goldenPath/help/wiggle.html) for format definition) from the genomic alignments contained in a qProject object. The combine argument controls if several input files are combined into a single multi-track wig file, or if they are exported as individual wig files. Alignments of single read experiments can be shifted towards there 3’-end using shift; paired-end alignments are automatically shifted by half the insert size. The resolution of the created wig file is defines by the binsize argument, and if scaling=TRUE, multiple alignment files in the qProject object are scaled by their total number of reads.

7.7 qCount

qCount is the workhorse for counting alignments that overlap query regions. Usage and details on parameters can be obtained from the qCount function documentation. Two aspects that are of special importance are also discussed here:

7.8 qProfile

The qProfile function differs from qCount in that it returns alignments counts relative to their position in the query region. Except for upstream and downstream, the arguments of qProfile and qCount are the same. This section will describe these two additional arguments; more details on the other arguments are available in section 7.7 and from the qProfile function documentation.

The regions to be profiled are anchored by the biological start position, which are aligned at position zero in the return value. The biological start position is defined as start(query) for regions on the plus strand and end(query) for regions on the minus strand. The anchor positions are extended to the left and right sides by the number of bases indicated in the upstream and downstream arguments.

• upstream indicates the number of bases upstream of the anchor position, which is on the left side of the anchor point for regions on the plus strand and on the right side for regions on the minus strand.
  
• downstream indicates the number of bases downstream of the anchor position, which is on the left side of the anchor point for regions on the plus strand and on the left side for regions on the minus strand.

Be aware that query regions with a * strand are handled the same way as regions on the plus strand.

7.9 qMeth

qMeth is used exclusively for Bis-seq experiments. In contrast to qCount, which counts the number of read alignments per query region, qMeth quantifies the number of C and T bases per cytosine in query regions, in order to determine methylation status.

qMeth can be run in one of four modes, controlled by the mode argument:

• CpGcomb: Only C’s in CpG context are considered. It is assumed that methylation status of the CpG base-pair on both strands is identical. Therefore, the total and methylated counts obtained for the C at position \(i\) and the C on the opposite strand at position \(i+1\) are summed.
  
• CpG: As with CpGcomb, only C’s in CpG context are quantified. However, counts from opposite strand are not summed, resulting in separate output values for C’s on both strands.
  
• allC: All C’s contained in query regions are quantified, keeping C’s from either strand separate. While this mode allows quantification of non-CpG methylation, it should be used with care, as the large result could use up available memory. In that case, a possible work-around is to divide the region of interest (e.g. the genome) into several regions (e.g. chromosomes) and call qMeth separately for each region.
  
• var: In this mode, only alignments on the opposite strand from C’s are analysed in order to collect evidence for sequence polymorphisms. Methylated C’s are hot-spots for C-to-T transitions, which in a Bis-seq experiment cannot be discriminated from completely unmethylated C’s. The information is however contained in alignments to the G from the opposite strand: Reads containing a G are consistent with a non-mutated C, and reads with an A support the presence of a sequence polymorphism. qMeth(..., mode="var") returns counts for total and matching bases for all C’s on both strands. A low fraction of matching bases is an indication of a mutation and can be used as a basis to identify mutated positions in the studied genome relative to the reference genome. Such positions should not be included in the quantification of methylation.

When using qMeth in a allele-specific quantification (see also section 6.5, cytosines (or CpGs) that overlap a sequence polymorphism will not be quantified.