library("Rmmquant") library("BiocStyle") library("S4Vectors") library("SummarizedExperiment") library("knitr") library("rmarkdown") library("TBX20BamSubset") library("TxDb.Mmusculus.UCSC.mm9.knownGene") library("org.Mm.eg.db") library("DESeq2") opts_chunk$set(message = FALSE, error = FALSE, warning = FALSE, cache = FALSE, fig.width = 5, fig.height = 5)
The aim of Rmmquant is to quantify the expression of the genes.
It is similar to
r Biocpkg("Rsubread"), or
The main difference with other approaches is that Rmmquant explicitely handles multi-mapping reads, in a way described by [@watson].
Rmmquant is the R port of the C++ mmquant, which has been published previously [@mmquant].
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
The matrix of counts can be accessed through:
Rmmquant 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,
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
r Biocpkg("GenomicRanges") is provided, the quantification will be
performed on each element of the
r Biocpkg("GenomicRangesList") is provided, the quantification will be
performed on each element too, i.e. on each
r Biocpkg("GenomicRanges") of a
r Biocpkg("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 [@samtools] or the
r Biocpkg("Rsamtools") to sort them. Rmmquant uses the
(that provides the number of hits for each read, see the
SAM format specification),
so be sure that your mapping tool sets it adequately
do it fine).
You should also check how your mapping tool handles multi-mapping reads
(this can usually be tuned using the appropriate parameters).
r Biocpkg("SummarizedExperiment") contains:
The count table can be accessed through the
assays method of
This methods returns a list, with only one element:
The columns are the samples (here, only
and the rows the gene counts (here, only
The count table can be used by
r Biocpkg("DESeq2") [@deseq2], for instance, using the
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_C), a new feature is
added to the matrix,
The reads that can be mapped to these genes will be counted there
(but not in the
The statistics can be accessed using the
colData method of
This is a
DataFrame (defined in
with one column per sample.
The content of each column is:
n.hits: the number of hits
n.uniquely.mapped.reads: the number of uniquely mapped reads
n.non.uniquely.mapped.hits: the number of non-uniquely mapped hits
n.ambiguously.mapped.hits: the number of hits with several corresponding features
n.unassigned.hits: the number of hits with no corresponding feature
Here, 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.
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,
The gene IDs and gene names should be given in the GTF file after the
gene_name tags respectively.
The column names of the count matrix should be given in the
If not given, the column names are inferred from the file names of the SAM/BAM
The library types can be specified using the
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,
FF: other similar cases,
U: the reads are single- or pair-ends, and the strand is unknown.
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
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
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"
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.
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.
nOverlapDiff=$n$ computes the differences of overlapping
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.
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
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
(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
mergeThreshold=$t$, a float in (0, 1). If $x < t \cdot y$,
then the reads are transferred.
The next examples uses data generated in
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.
bamFiles <- getBamFileList() sampleNames <- names(bamFiles)
You can extract the annotation a BioConductor package, such as
gr <- genes(TxDb.Mmusculus.UCSC.mm9.knownGene, filter=list(tx_chrom="chr19"))
The default gene IDs are Entrez ID [@entrez].
If you prefer Ensembl IDs [@ensembl],
you can modify the
r Biocpkg("GenomicRanges") accordingly.
Notice that Entrez IDs may have zero, one, or more Ensembl counterparts.
ensemblIds <- sapply(as.list(org.Mm.egENSEMBL[mappedkeys(org.Mm.egENSEMBL)]) [mcols(gr)$gene_id], `[[`, 1) gr <- gr[! is.na(names(ensemblIds)), ] names(gr) <- unlist(ensemblIds)
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, ]
While 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
before installing the package.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.