compute_read_counts: Compute read counts

In recount3: Explore and download data from the recount3 project

Description

As described in the recount workflow, the counts provided by the recount2 project are base-pair counts. You can scale them using transform_counts() or compute the read counts using the area under coverage information (AUC).

Usage

compute_read_counts(
  rse,
  round = TRUE,
  avg_mapped_read_length = "recount_qc.star.average_mapped_length"
)

Arguments

rse	A RangedSummarizedExperiment-class created by create_rse().
round	A logical(1) specifying whether to round the transformed counts or not.
avg_mapped_read_length	A character(1) specifying the metdata column name that contains the average fragment length after aligning. This is typically twice the average read length for paired-end reads.

Details

This function is similar to recount::read_counts(use_paired_end = TRUE, round = TRUE) but more general and with a different name to avoid NAMESPACE conflicts. Note that the default value of round is different than in recount::read_counts(). This was done to match the default value of round in transform_counts().

Value

A matrix() with the read counts. By default this function uses the average read length to the QC annotation.

References

Collado-Torres L, Nellore A and Jaffe AE. recount workflow: Accessing over 70,000 human RNA-seq samples with Bioconductor version 1; referees: 1 approved, 2 approved with reservations. F1000Research 2017, 6:1558 doi: 10.12688/f1000research.12223.1.

See Also

Other count transformation functions: compute_scale_factors(), is_paired_end(), transform_counts()

Examples

## Create a RSE object at the gene level
rse_gene_SRP009615 <- create_rse_manual("SRP009615")
colSums(compute_read_counts(rse_gene_SRP009615)) / 1e6

## Create a RSE object at the gene level
rse_gene_DRP000499 <- create_rse_manual("DRP000499")
colSums(compute_read_counts(rse_gene_DRP000499)) / 1e6

## You can compare the read counts against those from recount::read_counts()
## from the recount2 project which used a different RNA-seq aligner
## If needed, install recount, the R/Bioconductor package for recount2:
# BiocManager::install("recount")
recount2_readsums <- colSums(assay(recount::read_counts(
    recount::rse_gene_SRP009615
), "counts")) / 1e6
recount3_readsums <- colSums(compute_read_counts(rse_gene_SRP009615)) / 1e6
recount_readsums <- data.frame(
    recount2 = recount2_readsums[order(names(recount2_readsums))],
    recount3 = recount3_readsums[order(names(recount3_readsums))]
)
plot(recount2 ~ recount3, data = recount_readsums)
abline(a = 0, b = 1, col = "purple", lwd = 2, lty = 2)

## Repeat for DRP000499, a paired-end study
recount::download_study("DRP000499", outdir = tempdir())
load(file.path(tempdir(), "rse_gene.Rdata"), verbose = TRUE)

recount2_readsums <- colSums(assay(recount::read_counts(
    rse_gene
), "counts")) / 1e6
recount3_readsums <- colSums(compute_read_counts(rse_gene_DRP000499)) / 1e6
recount_readsums <- data.frame(
    recount2 = recount2_readsums[order(names(recount2_readsums))],
    recount3 = recount3_readsums[order(names(recount3_readsums))]
)
plot(recount2 ~ recount3, data = recount_readsums)
abline(a = 0, b = 1, col = "purple", lwd = 2, lty = 2)

recount3 documentation built on Feb. 13, 2021, 2 a.m.