The aim of Rmmquant is to quantify the expression of the genes. It is similar to featureCounts (Liao, Smyth, and Shi 2014) and Rsubread, or HTSeq-counts (Anders, Pyl, and Huber 2015) and the countOverlaps of GenomicRanges.
The main difference with other approaches is that Rmmquant explicitely handles multi-mapping reads, in a way described by (Robert and Watson 2015).
Rmmquant is the R port of the C++ mmquant, which has been published previously (Zytnicki 2017).
The easiest method to get the expression of the genes stored in a GTF file, using RNA-Seq data stored in a BAM file is:
dir <- system.file("extdata", package="Rmmquant", mustWork = TRUE)
gtfFile <- file.path(dir, "test.gtf")
bamFile <- file.path(dir, "test.bam")
output <- RmmquantRun(gtfFile, bamFile)In this example, the output is a SummarizedExperiment.
The matrix of counts can be accessed through:
## test
## geneA 1Rmmquant expects at least two kinds of data:
Many aspects of Rmmquant can be controlled with parameters.
All these notions are detailed hereafter.
The Rmmquant package mainly consists in one function, RmmquantRun function, that supports the options described hereafter,
Rmmquant also internally implements an S4 class, RmmquantClass, that stores all the inputs and the outputs of the RmmquantRun function, as well as a validator (that checks the consistency of the inputs).
The annotation file should be in GTF. GFF might work too. The tool only uses the gene/transcript/exon types.
Alternatively, the structure can be given using GenomicRanges, or GenomicRangesList.
When a GenomicRanges is provided, the quantification will be performed on each element of the GenomicRanges.
When a GenomicRangesList is provided, the quantification will be performed on each element too, i.e. on each GenomicRanges. Implicitely, each GenomicRanges of a GenomicRangesList is interpreted as a series of exons, and only the transcript (or the gene) will be quantified.
One or several reads files can be given.
The reads should be given in SAM or BAM format, and preferably be sorted by position, but it is not compulsory. The reads can be single end or paired-end (or a mixture thereof).
You can use the samtools (Li et al. 2009) or the Rsamtools to sort them. Rmmquant uses the NH flag (that provides the number of hits for each read, see the SAM format specification), so be sure that your mapping tool sets it adequately (yes, TopHat2 (Kim et al. 2013) and STAR (Dobin et al. 2013) do it fine). You should also check how your mapping tool handles multi-mapping reads (this can usually be tuned using the appropriate parameters).
The output SummarizedExperiment contains:
The count table can be accessed through the assays method of SummarizedExperiment. This methods returns a list, with only one element: counts.
## test
## geneA 1The columns are the samples (here, only test.bam), and the rows the gene counts (here, only gene_A).
The count table can be used by DESeq2 (Love, Huber, and Anders 2014), for instance, using the DESeqDataSetFromMatrix function.
If the user provided \(n\) reads files, the output will contain \(n\) columns:
| Gene | sample_1 | sample_2 | … |
|---|---|---|---|
| gene_A | … | … | … |
| gene_B | … | … | … |
| gene_B–gene_C | … | … | … |
The row names are the ID of the genes/features. The column names are the sample names. If a read maps several genes (say, gene_B and gene_C), a new feature is added to the matrix, gene_B--gene_C. The reads that can be mapped to these genes will be counted there (but not in the gene_B nor gene_C lines).
The statistics can be accessed using the colData method of SummarizedExperiment:
## DataFrame with 1 row and 5 columns
## n.hits n.uniquely.mapped.reads n.ambiguously.mapped.hits
## <numeric> <numeric> <numeric>
## test 2 1 0
## n.non.uniquely.mapped.hits n.unassigned.hits
## <numeric> <numeric>
## test 0 1This is a DataFrame (defined in S4Vectors), with one column per sample. The content of each column is:
n.hits: the number of hitsn.uniquely.mapped.reads: the number of uniquely mapped readsn.non.uniquely.mapped.hits: the number of non-uniquely mapped hitsn.ambiguously.mapped.hits: the number of hits with several corresponding featuresn.unassigned.hits: the number of hits with no corresponding featureHere, a hit is a possible mapping for a read. When a read maps several times, it means that several hits correspond to a read. A read may have zero, one, or several hits, which correspond to unmapped, uniquely mapped, and non-uniquely mapped reads.
If printGeneName is set to TRUE, the row names of the count table are the gene names, instead of the gene IDs. If two different genes have the same name, the systematic name is added, like: Mat2a (ENSMUSG00000053907).
The gene IDs and gene names should be given in the GTF file after the gene_id and gene_name tags respectively.
The column names of the count matrix should be given in the sampleNames. If not given, the column names are inferred from the file names of the SAM/BAM files.
The library types can be specified using the strands parameter. This parameter can be:
F: the reads are single-end, and the forward strand is sequenced,R: the reads are single-end, and the reverse strand is sequenced,FR: the reads are paired-end, the forward strand is sequenced first, then the reverse strand,RF, FF, FF: other similar cases,U: the reads are single- or pair-ends, and the strand is unknown.The strands parameters expect one value per SAM/BAM file. If only one value is given, it is recycled.
The format is usually inferred from the file name, but you can mention it using the formats option.
This parameters expect one value per SAM/BAM file. If only one value is given, it is recycled.
The way a read \(R\) is mapped to a gene \(A\) depends on the overlap=\(n\) value:
| if \(n\) is | then \(R\) is mapped to \(A\) iff |
|---|---|
| a negative value | \(R\) is included in \(A\) |
| a positive integer | they have at least \(n\) nucleotides in common |
| a float value (0, 1) | \(n\)% of the nucleotides of \(R\) are shared with \(A\) |
We will suppose here that the overlap=1 strategy is used (i.e. a read is attributed to a gene as soon as at least 1 nucleotide overlap). The example can be extended to other strategies as well.
If a read (say, of size 100), maps unambiguously and overlaps with gene A and B, it will be counted as 1 for the new “gene” gene_A--gene_B. However, suppose that only 1 nucleotide overlaps with gene A, whereas 100 nucleotides overlap with gene B (yes, genes A and B overlap). You probably would like to attribute the read to gene B.
The options nOverlapDiff and pcOverlapDiff control this. We compute the number of overlapping nucleotides between a read and the overlapping genes. If a read overlaps “significantly” more with one gene than with all the other genes, they will attribute the read to the former gene only.
The option nOverlapDiff=\(n\) computes the differences of overlapping nucleotides. Let us name \(N_A\) and \(N_B\) the number of overlapping nucleotides with genes A and B respectively. If \(N_A \geq N_B + n\), then the read will be attributed to gene A only.
The option pcOverlapDiff=\(m\) compares the ratio of overlapping nucleotides. If \(N_A / N_B \geq m\%\), then the read will be attributed to gene A only.
If both option nOverlapDiff=\(n\) and pcOverlapDiff=\(m\) are used, then the read will be attributed to gene A only iff both \(N_A \geq N_B + n\) and \(N_A / N_B \geq m\%\).
If the maximum number of reads for a gene is less than countThreshold (a non-negative integer), then the corresponding line is discarded.
Sometimes, there are very few reads that can be mapped unambiguously to a gene A, because it is very similar to gene B.
| Gene | sample_1 |
|---|---|
| gene_A | \(x\) |
| gene_B | \(y\) |
| gene_A–gene_B | \(z\) |
In the previous example, suppose that \(x \ll z\). In this case, you can move all the reads from gene_A to gene_A--gene_B, using the mergeThreshold=\(t\), a float in (0, 1). If \(x < t \cdot y\), then the reads are transferred.
The next examples uses data generated in TBX20BamSubset, where the expression of the genes is compared between the wild type and the TXB20 knock-out mice. The data have been mapped to the mm9 reference, and restricted to chromosome 19.
You can extract the annotation a BioConductor package, such as TxDb.Mmusculus.UCSC.mm9.knownGene.
The default gene IDs are Entrez ID (Maglott et al. 2015). If you prefer Ensembl IDs (Zerbino et al. 2018), you can modify the GenomicRanges accordingly. Notice that Entrez IDs may have zero, one, or more Ensembl counterparts.
DESeq2 can be executed right after Rmmquant.
output <- RmmquantRun(genomicRanges=gr, readsFiles=bamFiles,
sampleNames=sampleNames, sorts=FALSE)
coldata <- data.frame(condition=factor(c(rep("control", 3), rep("KO", 3)),
levels=c("control", "KO")),
row.names=sampleNames)
dds <- DESeqDataSetFromMatrix(countData=assays(output)$counts,
colData =coldata,
design =~ condition)
dds <- DESeq(dds)
res <- lfcShrink(dds, coef=2)
res$padj <- ifelse(is.na(res$padj), 1, res$padj)
res[res$padj < 0.05, ]## log2 fold change (MAP): condition KO vs control
## Wald test p-value: condition KO vs control
## DataFrame with 1 row and 5 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000024997 1083.54 -2.01881 0.338944 3.30589e-10 2.64471e-09While installing the package, if the compiler complains and says
#error This file requires compiler and library support for the ISO C++ 2011 standard.
This support is currently experimental, and must be enabled with the -std=c++11 or -std=gnu++11 compiler options.Add this line
Sys.setenv("PKG_CXXFLAGS"="-std=c++11")before installing the package.
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.3 (2020-10-10)
## os Ubuntu 18.04.5 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate C
## ctype en_US.UTF-8
## tz America/New_York
## date 2021-01-05
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## annotate 1.68.0 2021-01-05 [2] Bioconductor
## AnnotationDbi * 1.52.0 2021-01-05 [2] Bioconductor
## apeglm 1.12.0 2021-01-05 [2] Bioconductor
## askpass 1.1 2019-01-13 [2] CRAN (R 4.0.3)
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.3)
## bbmle 1.0.23.1 2020-02-03 [2] CRAN (R 4.0.3)
## bdsmatrix 1.3-4 2020-01-13 [2] CRAN (R 4.0.3)
## Biobase * 2.50.0 2021-01-05 [2] Bioconductor
## BiocFileCache 1.14.0 2021-01-05 [2] Bioconductor
## BiocGenerics * 0.36.0 2021-01-05 [2] Bioconductor
## BiocManager 1.30.10 2019-11-16 [2] CRAN (R 4.0.3)
## BiocParallel 1.24.1 2021-01-05 [2] Bioconductor
## BiocStyle * 2.18.1 2021-01-05 [2] Bioconductor
## biomaRt 2.46.0 2021-01-05 [2] Bioconductor
## Biostrings * 2.58.0 2021-01-05 [2] Bioconductor
## bit 4.0.4 2020-08-04 [2] CRAN (R 4.0.3)
## bit64 4.0.5 2020-08-30 [2] CRAN (R 4.0.3)
## bitops 1.0-6 2013-08-17 [2] CRAN (R 4.0.3)
## blob 1.2.1 2020-01-20 [2] CRAN (R 4.0.3)
## callr 3.5.1 2020-10-13 [2] CRAN (R 4.0.3)
## cli 2.2.0 2020-11-20 [2] CRAN (R 4.0.3)
## coda 0.19-4 2020-09-30 [2] CRAN (R 4.0.3)
## colorspace 2.0-0 2020-11-11 [2] CRAN (R 4.0.3)
## crayon 1.3.4 2017-09-16 [2] CRAN (R 4.0.3)
## curl 4.3 2019-12-02 [2] CRAN (R 4.0.3)
## DBI 1.1.0 2019-12-15 [2] CRAN (R 4.0.3)
## dbplyr 2.0.0 2020-11-03 [2] CRAN (R 4.0.3)
## DelayedArray 0.16.0 2021-01-05 [2] Bioconductor
## desc 1.2.0 2018-05-01 [2] CRAN (R 4.0.3)
## DESeq2 * 1.30.0 2021-01-05 [2] Bioconductor
## devtools 2.3.2 2020-09-18 [2] CRAN (R 4.0.3)
## digest 0.6.27 2020-10-24 [2] CRAN (R 4.0.3)
## dplyr 1.0.2 2020-08-18 [2] CRAN (R 4.0.3)
## ellipsis 0.3.1 2020-05-15 [2] CRAN (R 4.0.3)
## emdbook 1.3.12 2020-02-19 [2] CRAN (R 4.0.3)
## evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.3)
## fansi 0.4.1 2020-01-08 [2] CRAN (R 4.0.3)
## fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.3)
## genefilter 1.72.0 2021-01-05 [2] Bioconductor
## geneplotter 1.68.0 2021-01-05 [2] Bioconductor
## generics 0.1.0 2020-10-31 [2] CRAN (R 4.0.3)
## GenomeInfoDb * 1.26.2 2021-01-05 [2] Bioconductor
## GenomeInfoDbData 1.2.4 2020-10-12 [2] Bioconductor
## GenomicAlignments 1.26.0 2021-01-05 [2] Bioconductor
## GenomicFeatures * 1.42.1 2021-01-05 [2] Bioconductor
## GenomicRanges * 1.42.0 2021-01-05 [2] Bioconductor
## ggplot2 3.3.3 2020-12-30 [2] CRAN (R 4.0.3)
## glue 1.4.2 2020-08-27 [2] CRAN (R 4.0.3)
## gtable 0.3.0 2019-03-25 [2] CRAN (R 4.0.3)
## hms 0.5.3 2020-01-08 [2] CRAN (R 4.0.3)
## htmltools 0.5.0 2020-06-16 [2] CRAN (R 4.0.3)
## httr 1.4.2 2020-07-20 [2] CRAN (R 4.0.3)
## IRanges * 2.24.1 2021-01-05 [2] Bioconductor
## knitr * 1.30 2020-09-22 [2] CRAN (R 4.0.3)
## lattice 0.20-41 2020-04-02 [2] CRAN (R 4.0.3)
## lifecycle 0.2.0 2020-03-06 [2] CRAN (R 4.0.3)
## locfit 1.5-9.4 2020-03-25 [2] CRAN (R 4.0.3)
## magrittr 2.0.1 2020-11-17 [2] CRAN (R 4.0.3)
## MASS 7.3-53 2020-09-09 [2] CRAN (R 4.0.3)
## Matrix 1.3-0 2020-12-22 [2] CRAN (R 4.0.3)
## MatrixGenerics * 1.2.0 2021-01-05 [2] Bioconductor
## matrixStats * 0.57.0 2020-09-25 [2] CRAN (R 4.0.3)
## memoise 1.1.0 2017-04-21 [2] CRAN (R 4.0.3)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.0.3)
## mvtnorm 1.1-1 2020-06-09 [2] CRAN (R 4.0.3)
## numDeriv 2016.8-1.1 2019-06-06 [2] CRAN (R 4.0.3)
## openssl 1.4.3 2020-09-18 [2] CRAN (R 4.0.3)
## org.Mm.eg.db * 3.12.0 2020-10-15 [2] Bioconductor
## pillar 1.4.7 2020-11-20 [2] CRAN (R 4.0.3)
## pkgbuild 1.2.0 2020-12-15 [2] CRAN (R 4.0.3)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.3)
## pkgload 1.1.0 2020-05-29 [2] CRAN (R 4.0.3)
## plyr 1.8.6 2020-03-03 [2] CRAN (R 4.0.3)
## prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.0.3)
## processx 3.4.5 2020-11-30 [2] CRAN (R 4.0.3)
## progress 1.2.2 2019-05-16 [2] CRAN (R 4.0.3)
## ps 1.5.0 2020-12-05 [2] CRAN (R 4.0.3)
## purrr 0.3.4 2020-04-17 [2] CRAN (R 4.0.3)
## R6 2.5.0 2020-10-28 [2] CRAN (R 4.0.3)
## rappdirs 0.3.1 2016-03-28 [2] CRAN (R 4.0.3)
## RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 4.0.3)
## Rcpp 1.0.5 2020-07-06 [2] CRAN (R 4.0.3)
## RCurl 1.98-1.2 2020-04-18 [2] CRAN (R 4.0.3)
## remotes 2.2.0 2020-07-21 [2] CRAN (R 4.0.3)
## rlang 0.4.10 2020-12-30 [2] CRAN (R 4.0.3)
## rmarkdown * 2.6 2020-12-14 [2] CRAN (R 4.0.3)
## Rmmquant * 1.8.1 2021-01-05 [1] Bioconductor
## rprojroot 2.0.2 2020-11-15 [2] CRAN (R 4.0.3)
## Rsamtools * 2.6.0 2021-01-05 [2] Bioconductor
## RSQLite 2.2.1 2020-09-30 [2] CRAN (R 4.0.3)
## rtracklayer 1.50.0 2021-01-05 [2] Bioconductor
## S4Vectors * 0.28.1 2021-01-05 [2] Bioconductor
## scales 1.1.1 2020-05-11 [2] CRAN (R 4.0.3)
## sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.3)
## stringi 1.5.3 2020-09-09 [2] CRAN (R 4.0.3)
## stringr 1.4.0 2019-02-10 [2] CRAN (R 4.0.3)
## SummarizedExperiment * 1.20.0 2021-01-05 [2] Bioconductor
## survival 3.2-7 2020-09-28 [2] CRAN (R 4.0.3)
## TBX20BamSubset * 1.26.0 2021-01-04 [2] Bioconductor
## testthat 3.0.1 2020-12-17 [2] CRAN (R 4.0.3)
## tibble 3.0.4 2020-10-12 [2] CRAN (R 4.0.3)
## tidyselect 1.1.0 2020-05-11 [2] CRAN (R 4.0.3)
## TxDb.Mmusculus.UCSC.mm9.knownGene * 3.2.2 2020-10-12 [2] Bioconductor
## usethis 2.0.0 2020-12-10 [2] CRAN (R 4.0.3)
## vctrs 0.3.6 2020-12-17 [2] CRAN (R 4.0.3)
## withr 2.3.0 2020-09-22 [2] CRAN (R 4.0.3)
## xfun 0.19 2020-10-30 [2] CRAN (R 4.0.3)
## XML 3.99-0.5 2020-07-23 [2] CRAN (R 4.0.3)
## xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.3)
## xtable 1.8-4 2019-04-21 [2] CRAN (R 4.0.3)
## XVector * 0.30.0 2021-01-05 [2] Bioconductor
## yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.3)
## zlibbioc 1.36.0 2021-01-05 [2] Bioconductor
##
## [1] /tmp/RtmpTloDiV/Rinst28ff4b5f18fe
## [2] /home/biocbuild/bbs-3.12-bioc/R/libraryAnders, Simon, Paul Theodor Pyl, and Wolfgang Huber. 2015. “HTSeq—a Python Framework to Work with High-Throughput Sequencing Data.” Bioinformatics 31 (2): 166–69.
Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal Rna-Seq Aligner.” Bioinformatics 29 (1): 15–21.
Kim, Daehwan, Geo Pertea, Cole Trapnell, Harold Pimentel, Ryan Kelley, and Steven L. Salzberg. 2013. “TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions.” Genome Biology 14 (4): R36.
Li, Heng, Bob Handsaker, Alec Wysoker, Tim Fennell, Jue Ruan, Nils Homer, Gabor Marth, Goncalo Abecasis, Richard Durbin, and 1000 Genome Project Data Processing Subgroup. 2009. “The Sequence Alignment/Map Format and Samtools.” Bioinformatics 25 (16): 2078–9.
Liao, Yang, Gordon K. Smyth, and Wei Shi. 2014. “FeatureCounts: An Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7): 923–30.
Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15 (12): 550.
Maglott, Donna, Jim Ostell, Kim D Pruitt, and Tatiana Tatusova. 2015. “Entrez Gene: gene-centered information at NCBI.” Nucleic Acids Research 33: D54–D58.
Robert, Christelle, and Mick Watson. 2015. “Errors in Rna-Seq Quantification Affect Genes of Relevance to Human Disease.” Genome Biology 16 (1): 177.
Zerbino, Daniel R, Premanand Achuthan, Wasiu Akanni, M Ridwan Amode, Daniel Barrell, Jyothish Bhai, Konstantinos Billis, et al. 2018. “Ensembl 2018.” Nucleic Acids Research 46 (D1): D754–D761.
Zytnicki, Matthias. 2017. “Mmquant: How to Count Multi-Mapping Reads?” BMC Bioinformatics 18 (1): 411.