metaseqr: The main metaseqr pipeline

In metaseqR: An R package for the analysis and result reporting of RNA-Seq data by combining multiple statistical algorithms.

Description

This function is the main metaseqr workhorse and implements the main metaseqr workflow which performs data read, filtering, normalization and statistical selection, creates diagnostic plots and exports the results and a report if requested. The metaseqr function is responsible for assembling all the steps of the metaseqr pipeline which i) reads the input gene or exon read count table ii) performs prelimininary filtering of data by removing chrM and other non-essential information for a typical differential gene expression analysis as well as a preliminary expression filtering based on the exon counts, if an exon read count file is provided. iii) performs data normalization with one of currently widely used algorithms, including EDASeq (Risso et al., 2011), DESeq (Anders and Huber, 2010), edgeR (Robinson et al., 2010), NOISeq (Tarazona et al., 2012) or no normalization iv) performs a second stage of filtering based on the normalized gene expression according to several gene filters v) performs statistical testing with one or more of currently widely used algorithms, including DESeq (Anders and Huber, 2010), edgeR (Robinson et al., 2010), NOISeq (Tarazona et al., 2012), limma (Smyth et al., 2005) for RNA-Seq data, baySeq (Hardcastle et al., 2012) vi) in the case of multiple statistical testing algorithms, performs meta-analysis using one of five available methods (see the meta.p argument) vii) exports the resulting differentially expressed gene list in text tab-delimited format viii) creates a set of diagnostic plots either available in the aforementioned packages or metaseqr specific ones and ix) creates a comprehensive HTML report which summarizes the run information, the results and the diagnostic plots. Certain diagnostic plots (e.g. the volcano plot) can be interactive with the use of the external Highcharts (http://www.highcharts.com) JavaScript library for interactive graphs. Although the inputs to the metaseqr workflow are many, in practice, setting only very few of them and accepting the defaults as the rest can result in quite comprehensible results for mainstream organisms like mouse, human, fly and rat.

Usage

    metaseqr(counts, sample.list, exclude.list = NULL,
        file.type = c("auto", "sam", "bam", "bed"),
        path = NULL, contrast = NULL, libsize.list = NULL,
        id.col = 4, gc.col = NA, name.col = NA, bt.col = NA,
        annotation = c("download", "embedded"), gene.file = NULL,
        org = c("hg18", "hg19", "hg38", "mm9", "mm10", "rn5", "rn6",
            "dm3", "dm6", "danrer7", "pantro4", "susscr3", "tair10", 
            "equcab2","custom"),
        refdb = c("ensembl", "ucsc", "refseq"),
        trans.level = c("gene", "transcript", "exon"),
        count.type = c("gene", "exon","utr"),
        utr.flank = 500,
        exon.filters = list(min.active.exons = list(exons.per.gene = 5, 
                min.exons = 2, frac = 1/5)),
        gene.filters = list(length = list(length = 500), 
                avg.reads = list(average.per.bp = 100, quantile = 0.25), 
                expression = list(median = TRUE, mean = FALSE, quantile = NA, 
                    known = NA, custom = NA), 
                biotype = get.defaults("biotype.filter", org[1]),
                presence = list(frac = 0.25, min.count = 10,
					per.condition = FALSE)),
        when.apply.filter = c("postnorm", "prenorm"),
        normalization = c("deseq", "edaseq", "edger", "noiseq", "nbpseq", 
                "each", "none"),
        norm.args = NULL,
        statistics = c("deseq", "edger", "noiseq", "bayseq", "limma", 
                "nbpseq"),
        stat.args = NULL,
        adjust.method = sort(c(p.adjust.methods, "qvalue")),
        meta.p = if (length(statistics) > 1) c("simes", "bonferroni", "fisher", 
                "dperm.min", "dperm.max", "dperm.weight", "fperm", "whitlock", 
                "minp", "maxp", "weight", "pandora", "none") else "none",
        weight = rep(1/length(statistics), length(statistics)),
        nperm = 10000, reprod=TRUE, pcut = NA, log.offset = 1,
        preset = NULL,
        qc.plots = c("mds", "biodetection", "countsbio", "saturation", 
                "readnoise", "filtered", "correl", "pairwise", "boxplot", 
                "gcbias", "lengthbias", "meandiff", "meanvar", "rnacomp", 
                "deheatmap", "volcano", "biodist"),
        fig.format = c("png", "jpg", "tiff", "bmp", "pdf", "ps"),
        out.list = FALSE, export.where = NA,
        export.what = c("annotation", "p.value", "adj.p.value", 
                "meta.p.value", "adj.meta.p.value", "fold.change", 
                "stats", "counts","flags"),
        export.scale = c("natural", "log2", "log10", "vst", "rpgm"),
        export.values = c("raw", "normalized"),
        export.stats = c("mean", "median", "sd", "mad", "cv", 
                "rcv"),
        export.counts.table = FALSE,
        restrict.cores = 0.6, report = TRUE, report.top = 0.1,
        report.template = "default", save.gene.model = TRUE,
        verbose = TRUE, run.log = TRUE, progress.fun = NULL, ...)

Arguments

counts	a text tab-delimited file containing gene or exon counts in one of the following formats: i) the first column contains unique gene or exon identifiers and the rest of the columns contain the read counts for each sample. Thus the first cell of each row is a gene or exon accession and the rest are integers representing the counts for that accession. In that case, the annotation parameter should strictly be "download" or an external file in proper format. ii) The first n columns should contain gene or exon annotation elements like chromosomal locations, gene accessions, exon accessions, GC content etc. In that case, the annotation parameter can also be "embedded". The ideal embedded annotation contains 8 columns, chromosome, gene or exon start, gene or exon end, gene or exon accession, GC-content (fraction or percentage), strand, HUGO gene symbol and gene biotype (e.g. "protein_coding" or "ncRNA"). When the annotation parameter is "embedded", certain of these features are mandatory (co-ordinates and accessions). If they are not present, the pipeline will not run. If additional elements are not present (e.g. GC content or biotypes), certain features of metaseqr will not be available. For example, EDASeq normalization will not be performed based on a GC content covariate but based on gene length which is not what the authors of EDASeq suggest. If biotypes are not present, a lot of diagnostic plots will not be available. If the HUGO gene symbols are missing, the final annotation will contain only gene accessions and thus be less comprehensible. Generally, it's best to set the annotation parameter to "download" to ensure the most comprehensible results. Counts can be a data frame satisfying the above conditions. It is a data frame by default when read2count is used. counts can also be an .RData file (output of save function which contains static input elements (list containing the gene model (exon counts for each gene constructed by the construct.gene.model function, gene and exon annotation to avoid re-downloading and/or gene counts depending on count.type). This kind of input facilitates the re-analysis of the same experiment, using different filtering, normalization and statistical algorithms. Finally, counts can be a list representing the gene model (exon counts for each gene) constructed by the construct.gene.model function (provided for backwards compatibility). This .RData file can be generated by setting save.gene.model=TRUE when performing data analysis for the first time.
sample.list	a list containing condition names and the samples under each condition. It should have the format sample.list <- list(ConditionA=c("Sample_A1", "Sample_A2", "Sample_A3"), ConditionB=c("Sample_B1", "Sample_B2"), ConditionC=c("Sample_C1", "Sample_C2")). The names of the samples in list members MUST match the column names containing the read counts in the counts file. If they do not match, the pipeline will either crash or at best, ignore several of your samples. Alternative, sample.list can be a small tab-delimited file structured as follows: the first line of the external tab delimited file should contain column names (names are not important). The first column MUST contain UNIQUE sample names and the second column MUST contain the biological condition where each of the samples in the first column should belong to. In this case, the function make.sample.list is used. If the counts argument is missing, the sample.list argument MUST be a targets text tab-delimited file which contains the sample names, the BAM/BED file names and the biological conditions/groups for each sample/file. The file should be text tab-delimited and structured as follows: the first line of the external tab delimited file should contain column names (names are not important). The first column MUST contain UNIQUE sample names. The second column MUST contain the raw BAM/BED files WITH their full path. Alternatively, the path argument should be provided (see below). The third column MUST contain the biological condition where each of the samples in the first column should belong to.
exclude.list	a list of samples to exclude, in the same (list) format as sample.list above.
path	an optional path where all the BED/BAM files are placed, to be prepended to the BAM/BED file names in the targets file. If not given and if the files in the second column of the targets file do not contain a path to a directory, the current directory is assumed to be the BAM/BED file container.
file.type	the type of raw input files. It can be "auto" for auto-guessing, "bed" for BED files, "sam" for SAM files or "bam" for BAM files.
contrast	a character vector of contrasts to be tested in the statistical testing step(s) of the metaseqr pipeline. Each element of contrast should STRICTLY have the format "ConditionA_vs_ConditionB_vs_...". A valid example based on the sample.list above is contrast <- c("ConditionA_vs_ConditionB", "ConditionA_vs_ConditionC", "ConditionA_vs_ConditionB_vs_ConditionC"). The first element of pairwise contrasts (e.g. "ConditionA" above) MUST be the control condition or any reference that ConditionB is checked against. metaseqr uses this convention to properly calculate fold changes. If it's NULL, a contrast between the first two members of the sample.list will be auto-generated.
libsize.list	an optional named list where names represent samples (MUST be the same as the samples in sample.list) and members are the library sizes (the sequencing depth) for each sample. For example libsize.list <- list(Sample_A1=32456913, Sample_A2=4346818).
id.col	an integer denoting the column number in the file (or data frame) provided with the counts argument, where the unique gene or exon accessions are. Default to 4 which is the standard feature name column in a BED file.
gc.col	an integer denoting the column number in the file (or data frame) provided with the counts argument, where each gene's GC content is given. If not provided, GC content normalization provided by EDASeq will not be available.
name.col	an integer denoting the column number in the file (or data frame) provided with the counts argument, where the HUGO gene symbols are given. If not provided, it will not be available when reporting results. In addition, the "known" gene filter will not be available.
bt.col	an integer denoting the column number in the file (or data frame) provided with the counts argument, where the gene biotypes are given. If not provided, the "biodetection", "countsbio", "saturation", "filtered" and "biodist" plots will not be available.
annotation	instructs metaseqr where to find the annotation for the given counts file. It can be one of i) "download" (default) for automatic downloading of the annotation for the organism specified by the org parameter (using biomaRt), ii) "embedded" if the annotation elements are embedded in the read counts file or iv) a file specified by the user which should be as similar as possible to the "download" case, in terms of column structure.
gene.file	an external gene annotation file required when annotation="embedded", count.type="exon" or count.type="utr" and org="custom". See the result of get.annotation on the format of the external gene file.
org	the supported organisms by metaseqr. These can be, for human genomes "hg18", "hg19" or "hg38", for mouse genomes "mm9", "mm10", for rat genomes "rn5" or "rn6", for drosophila genome "dm3" or "dm6", for zebrafish genome "danrer7", for chimpanzee genome "pantro4", for pig genome "susscr3", for Arabidopsis thaliana genome "tair10" and for Equus caballus genome "equcab2". Finally, "custom" will instruct metaseqR to completely ignore the org argument and depend solely on annotation file provided by the user.
refdb	the reference annotation repository from which to retrieve annotation elements to use with metaseqr. It can be one of "ensembl" (default), "ucsc" or "refseq".
trans.level	perform differential expression analysis at which transcriptional unit, can be one of "gene" (default), "transcript" for reporting differential expression at the transcript level or "exon" for exon level.
count.type	the type of reads inside the counts file. It can be one of "gene", "exon" or "utr" for quant seq (lexogen) protocol. This is a very important and mandatory parameter as it defines the course of the workflow.
utr.flank	the number of base pairs to flank the 3' UTR of transcripts when analyzing Quant-Seq data.
exon.filters	a named list whose names are the names of the supported exon filters and its members the filter parameters. See section "Exon filters" below for details.
gene.filters	a named list whose names are the names of the supported gene filters and its members the filter parameters. See section "Gene filters" below for details.
when.apply.filter	a character string determining when to apply the exon and/or gene filters, relative to normalization. It can be "prenorm" to apply apply the filters and exclude genes from further processing before normalization, or "postnorm" to apply the filters after normalization (default). In the case of when.apply.filter="prenorm", a first normalization round is applied to a copy of the gene counts matrix in order to derive the proper normalized values that will constitute the several expression-based filtering cutoffs.
normalization	the normalization algorithm to be applied on the count data. It can be one of "edaseq" for EDASeq normalization, "deseq" for the normalization algorithm (individual options specified by the norm.args argument) in the DESq package (default), "edger" for the normalization algorithms present in the edgeR package (specified by the norm.args argument), "noiseq" for the normalization algorithms present in the NOISeq package (specified by the norm.args argument), "nbpseq" for the normalization algorithms present in the NBPSeq package (specified by the norm.args argument) or "none" to not normalize the data (highly unrecommended). It can also be "each" where in this case, the normalization applied will be specific to each statistical test used (i.e. the normalization method bundled with each package and used in its examples and documentation). The last choice is for future use!
norm.args	a named list whose names are the names of the normalization algorithm parameters and its members parameter values. See section "Normalization parameters" below for details. Leave NULL for the defaults of normalization. If normalization="each", it must be a named list of lists, where each sub-list contains normalization parameters specific to each statistical test to be used. The last choice is for future use!
statistics	one or more statistical analyses to be performed by the metaseqr pipeline.It can be one or more of "deseq" (default) to conduct statistical test(s) implemented in the DESeq package, "edger" to conduct statistical test(s) implemented in the edgeR package, "limma" to conduct the RNA-Seq version of statistical test(s) implemented in the limma package, "noiseq" to conduct statistical test(s) implemented in the NOISeq package, "bayseq" to conduct statistical test(s) implemented in the baySeq package and "nbpseq" to conduct statistical test(s) implemented in the NBPSeq package. In any case individual algorithm parameters are controlled by the contents of the stat.args list.
stat.args	a named list whose names are the names of the statistical algorithms used in the pipeline. Each member is another named list whose names are the algorithm parameters and its members are the parameter values. See section "Statistics parameters" below for details. Leave NULL for the defaults of statistics.
adjust.method	the multiple testing p-value adjustment method. It can be one of p.adjust.methods or "qvalue" from the qvalue Bioconductor package. Defaults to "BH" for Benjamini-Hochberg correction.
meta.p	the meta-analysis method to combine p-values from multiple statistical tests . It can be one of "simes" (default), "bonferroni", "minp", "maxp", "weight", "pandora", "dperm.min", "dperm.max", "dperm.weight", "fisher", "fperm", "whitlock" or "none". For the "fisher" and "fperm" methods, see the documentation of the R package MADAM. For the "whitlock" method, see the documentation of the survcomp Bioconductor package. With the "maxp" option, the final p-value is the maximum p-value out of those returned by each statistical test. This is equivalent to an "intersection" of the results derived from each algorithm so as to have a final list with the common genes returned by all statistical tests. Similarly, when meta.p="minp", is equivalent to a "union" of the results derived from each algorithm so as to have a final list with all the genes returned by all statistical tests. The latter can be used as a very lose statistical threshold to aggregate results from all methods regardless of their False Positive Rate. With the "simes" option, the method proposed by Simes (Simes, R. J., 1986) is used. With the "dperm.min", "dperm.max", "dperm.weight" options, a permutation procedure is initialed, where nperm permutations are performed across the samples of the normalized counts matrix, producing nperm permuted instances of the initital dataset. Then, all the chosen statistical tests are re-executed for each permutation. The final p-value is the number of times that the p-value of the permuted datasets is smaller than the original dataset. The p-value of the original dataset is created based on the choice of one of dperm.min, dperm.max or dperm.weight options. In case of dperm.min, the intial p-value vector is consists of the minimum p-value resulted from the applied statistical tests for each gene. The maximum p-value is used with the dperm.max option. With the dperm.weight option, the weight weighting vector for each statistical test is used to weight each p-value according to the power of statistical tests (some might work better for a specific dataset). Be careful as the permutation procedure usually requires a lot of time. However, it should be the most accurate. This method will NOT work when there are no replicated samples across biological conditions. In that case, use meta.p="simes" instead. Finally, there are the "minp", "maxp" and "weight" options which correspond to the latter three methods but without permutations. Generally, permutations would be accurate to use when the experiment includes >5 samples per condition (or even better 7-10) which is rather rare in RNA-Seq experiments. Finally, "pandora" is the same as "weight" and is added to be in accordance with the metaseqR paper.
weight	a vector of weights with the same length as the statistics vector containing a weight for each statistical test. It should sum to 1. Use with caution with the dperm.weight parameter! Theoretical background is not yet solid and only experience shows improved results!
nperm	the number of permutations performed to derive the meta p-value when meta.p="fperm" or meta.p="dperm". It defaults to 10000.
reprod	create reproducible permutations when meta.p="dperm.min", meta.p="dperm.max" or meta.p="dperm.weight". Ideally one would want to create the same set of indices for a given dataset so as to create reproducible p-values. If reprod=TRUE, a fixed seed is used by meta.perm for all the datasets analyzed with metaseqr. If reprod=FALSE, then the p-values will not be reproducible, although statistical significance is not expected to change for a large number of resambling. Finally, reprod can be a numeric vector of seeds with the same length as nperm so that the user can supply his/her own seeds.
pcut	a p-value cutoff for exporting differentially genes, default is to export all the non-filtered genes.
log.offset	an offset to be added to values during logarithmic transformations in order to avoid Infinity (default is 1).
preset	an analysis strictness preset. preset can be one of "all.basic", "all.normal", "all.full", "medium.basic", "medium.normal", "medium.full", "strict.basic", "strict.normal" or "strict.full", each of which control the strictness of the analysis and the amount of data to be exported. For an explanation of the presets, see the section "Presets" below.
qc.plots	a set of diagnostic plots to show/create. It can be one or more of "mds", "biodetection", "rnacomp", "countsbio", "saturation", "readnoise", "filtered", "boxplot", "gcbias", "lengthbias", "meandiff", "meanvar", "deheatmap", "volcano", "biodist", "venn". The "mds" stands for Mutlti-Dimensional Scaling and it creates a PCA-like plot but using the MDS dimensionality reduction instead. It has been succesfully used for NGS data (e.g. see the package htSeqTools) and it shows how well samples from the same condition cluster together. For "biodetection", "countsbio", "saturation", "rnacomp", "readnoise", "biodist" see the vignette of NOISeq package. The "saturation" case has been rewritten in order to display more samples in a more simple way. See the help page of diagplot.noiseq.saturation. In addition, the "readnoise" plots represent an older version or the RNA composition plot included in older versions of NOISeq. For "gcbias", "lengthbias", "meandiff", "meanvar" see the vignette of EDASeq package. "lenghtbias" is similar to "gcbias" but using the gene length instead of the GC content as covariate. The "boxplot" option draws boxplots of log2 transformed gene counts. The "filtered" option draws a 4-panel figure with the filtered genes per chromosome and per biotype, as absolute numbers and as fractions of the genome. See also the help page of diagplot.filtered. The "deheatmap" option performs hierarchical clustering and draws a heatmap of differentially expressed genes. In the context of diagnostic plots, it's useful to see if samples from the same groups cluster together after statistical testing. The "volcano" option draws a volcano plot for each contrast and if a report is requested, an interactive volcano plot is presented in the HTML report. The "venn" option will draw an up to 5-way Venn diagram depicting the common and specific to each statistical algorithm genes and for each contrast, when meta-analysis is performed. The "correl" option creates two correlation graphs: the first one is a correlation heatmap (a correlation matrix which depicts all the pairwise correlations between each pair of samples in the counts matrix is drawn as a clustered heatmap) and the second one is a correlogram plot, which summarizes the correlation matrix in the form of ellipses (for an explanation please see the vignette/documentation of the R package corrplot. Set qc.plots=NULL if you don't want any diagnostic plots created.
fig.format	the format of the output diagnostic plots. It can be one or more of "png", "jpg", "tiff", "bmp", "pdf", "ps". The native format "x11" (for direct display) is not provided as an option as it may not render the proper display of some diagnostic plots in some devices.
out.list	a logical controlling whether to export a list with the results in the running environment.
export.where	an output directory for the project results (report, lists, diagnostic plots etc.)
export.what	the content of the final lists. It can be one or more of "annotation", to bind the annoation elements for each gene, "p.value", to bind the p-values of each method, "adj.p.value", to bind the multiple testing adjusted p-values, "meta.p.value", to bind the combined p-value from the meta-analysis, "adj.meta.p.value", to bind the corrected combined p-value from the meta-analysis, "fold.change", to bind the fold changes of each requested contrast, "stats", to bind several statistics calclulated on raw and normalized counts (see the export.stats argument), "counts", to bind the raw and normalized counts for each sample.
export.scale	export values from one or more transformations applied to the data. It can be one or more of "natural", "log2", "log10", "vst" (Variance Stabilizing Transormation, see the documentation of DESeq package) and "rpgm" which is ratio of mapped reads per gene model (either the gene length or the sum of exon lengths, depending on count.type argument). Note that this is not RPKM as reads are already normalized for library size using one of the supported normalization methods. Also, "rpgm" might be misleading when normalization is other than "deseq".
export.values	It can be one or more of "raw" to export raw values (counts etc.) and "normalized" to export normalized counts.
export.stats	calculate and export several statistics on raw and normalized counts, condition-wise. It can be one or more of "mean", "median", "sd", "mad", "cv" for the Coefficient of Variation, "rcv" for a robust version of CV where the median and the MAD are used instead of the mean and the standard deviation.
export.counts.table	exports also the calculated read counts table when input is read from bam files and exports also the normalized count table in all cases. Defaults to FALSE.
restrict.cores	in case of parallel execution of several subfunctions, the fraction of the available cores to use. In some cases if all available cores are used (restrict.cores=1 and the system does not have sufficient RAM, the pipeline running machine might significantly slow down.
report	a logical value controlling whether to produce a summary report or not. Defaults to TRUE.
report.top	a fraction of top statistically significant genes to append to the HTML report. This helps in keeping the size of the report as small as possible, as appending the total gene list might create a huge HTML file. Users can always retrieve the whole gene lists from the report links. Defaults to 0.1 (top 10 genes). Set to NA or NULL to append all the statistically significant genes to the HTML report.
report.template	an HTML template to use for the report. Do not change this unless you know what you are doing.
save.gene.model	in case of exon analysis, a list with exon counts for each gene will be saved to the file export.where/data/gene_model.RData. This file can be used as input to metaseqR for exon count based analysis, in order to avoid the time consuming step of assembling the counts for each gene from its exons
verbose	print informative messages during execution? Defaults to TRUE.
run.log	write a log file of the metaseqr run using package log4r. Defaults to TRUE. The filename will be auto-generated.
progress.fun	a function which updates a Progress object from shiny. This function must accept a detail argument. See http://shiny.rstudio.com/articles/progress.html
...	further arguments that may be passed to plotting functions, related to par.

Value

If out.list is TRUE, a named list whose length is the same as the number of requested contrasts. Each list member is named according to the corresponding contrast and contains a data frame of differentially expressed genes for that contrast. The contents of the data frame are defined by the export.what, export.scale, export.stats, export.values parameters. If report is TRUE, the output list contains two main elements. The first is described above (the analysis results) and the second contains the same results but in HTML formatted tables.

Exon filters

The exon filters are a set of filters which are applied after the gene models are assembled from the read counts of individual exons and before the gene expression is summarized from the exons belonging to each gene. These filters can be applied when the input read counts file contains exon reads. It is not applicable when the input file already contains gene counts. Such filters can be for example "accept genes where all the exons contain more than x reads" or "accept genes where there is read presence in at least m/n exons, n being the total exons of the gene". Such filters are NOT meant for detecting differential splicing as also the whole metaseqr pipeline, thus they should not be used in that context. The exon.filters argument is a named list of filters, where the names are the filter names and the members are the filter parameters (named lists with parameter name, parameter value). See the usage of the metaseqr function for an example of how these lists are structured. The supported exon filters in the current version are: i) min.active.exons which implements a filter for demanding m out of n exons of a gene to have a certain read presence with parameters exons.per.gene, min.exons and frac. The filter is described as follows: if a gene has up to exons.per.gene exons, then read presence is required in at least min.exons of them, else read presence is required in a frac fraction of the total exons. With the default values, the filter instructs that if a gene has up to 5 exons, read presence is required in at least 2, else in at least 20 exons, in order to be accepted. More filters will be implemented in future versions and users are encouraged to propose exon filter ideas to the author by mail. See metaseqr usage for the defaults. Set exon.filters=NULL to not apply any exon filtering.

Gene filters

The gene filters are a set of filters applied to gene expression as this is manifested through the read presence on each gene and are preferably applied after normalization. These filters can be applied both when the input file or data frame contains exon read counts and gene read counts. Such filter can be for example "accept all genes above a certain count threshold" or "accept all genes with expression above the median of the normalized counts distribution" or "accept all with length above a certain threshold in kb" or "exclude the 'pseudogene' biotype from further analysis". The supported gene filters in the current version, which have the same structure as the exon filters (named list of lists with filter names, parameter names and parameter arguments) are: i) length which implements a length filter where genes are accepted for further analysis if they are above length (its parameter) kb. ii) avg.reads which implements a filter where a gene is accepted for further analysis if it has more average reads than the quantile of the average count distribution per average.per.bp base pairs. In summary, the reads of each gene are averaged per average.per.bp based on each gene's length (in case of exons, input the "gene's length" is the sum of the lengths of exons) and the quantile quantile of the average counts distribution is calculated for each sample. Genes passing the filter should have an average read count larger than the maximum of the vector of the quantiles calculated above. iii) expression which implements a filter based on the overall expression of a gene. The parameters of this filter are: median, where genes below the median of the overall count distribution are not accepted for further analysis (this filter has been used to distinguish between "expressed" and "not expressed" genes in several cases, e.g. (Mokry et al., NAR, 2011) with a logical as value, mean which is the same as median but using the mean, quantile which is the same as the previous two but using a specific quantile of the total counts distribution, known, where in this case, a set of known not-expressed genes in the system under investigation are used to estimate an expression cutoff. This can be quite useful, as the genes are filtered based on a "true biological" cutoff instead of a statistical cutoff. The value of this filter is a character vector of HUGO gene symbols (MUST be contained in the annotation, thus it's better to use annotation="download") whose counts are used to build a "null" expression distribution. The 90th quantile of this distribution is then the expression cutoff. This filter can be combined with any other filter. Be careful with gene names as they are case sensitive and must match exactly ("Pten" is different from "PTEN"!). iv) biotype where in this case, genes with a certain biotype (MUST be contained in the annotation, thus it's better to use annotation="download") are excluded from the analysis. This filter is a named list of logical, where names are the biotypes in each genome and values are TRUE or FALSE. If the biotype should be excluded, the value should be TRUE else FALSE. See the result of get.defaults("biotype.filter","hg19") for an example. Finally, in future versions there will be support for user-defined filters in the form of a function. v) presence where in this case, a gene is further considered for statistical testing if frac (x100 for a percentage value) have more than min.count reads across all samples (per.condition=FALSE) or across the samples of each condition (per.condition=TRUE).

Normalization parameters

The normalization parameters are passed again as a named list where the names of the members are the normalization parameter names and the values are the normalization parameter values. You should check the documentation of the packages EDASeq, DESeq, edgeR, NOISeq and NBPSeq for the parameter names and parameter values. There are a few exceptions in parameter names: in case of normalization="edaseq" the only parameter names are within.which and between.which, controlling the withing lane/sample and between lanes/samples normalization algorithm. In the case of normalization="nbpseq", there is one additional parameter called main.method which can take the calues "nbpseq" or "nbsmyth". These values correspond to the two different workflows available in the NBPSeq package. Please, consult the NBPSeq package documentation for further details. For the rest of the algorithms, the parameter names are the same as the names used in the respective packages. For examples, please use the get.defaults function.

Statistics parameters

The statistics parameters as passed to statistical algorithms in metaseqr, exactly with the same way as the normalization parametes above. In this case, there is one more layer in list nesting. Thus, stat.args is a named list whose names are the names the algorithms used (see the statistics parameter). Each member is another named list,with parameters to be used for each statistical algorithm. Again, the names of the member lists are parameter names and the values of the member lists are parameter values. You should check the documentations of DESeq, edgeR, NOISeq, baySeq, limma and NBPSeq for these parameters. There are a few exceptions in parameter names: In case of statistics="edger", apart from the rest of the edgeR statistical testing arguments, there is the argument main.method which can be either "classic" or "glm", again defining whether the binomial test or GLMs will be used for statistical testing. For examples, please use the get.defaults function. When statistics="nbpseq", apart from the rest arguments of the NBPSeq functions estimate.disp and estimate.dispersion, there is the argument main.method which can be "nbpseq" or "nbsmyth". This argument determines the parameters to be used by the estimate.dispersion function or by the estimate.disp function to estimate RNA-Seq count dispersions. The difference between the two is that they constitute different starting points for the two workflows in the package NBPSeq. The first worklfow (with main.method="nbpseq" and the estimate.dispersion function is NBPSeq package specific, while the second (with main.method="nbsmyth" and the estimate.disp function is similar to the workflow of the edgeR package. For additional information regarding the statistical testing in NBPSeq, please consult the documentation of the NBPSeq package. Additinally, please note that there is currently a problem with the NBPSeq package and the workflow that is specific to the NBPSeq package. The problem has to do with function exporting as there are certain functions which are not recognized from the package internally. For this reason and until it is fixed, only the Smyth workflow will be available with the NBPSeq package (thus stat.args$main.method="nbpseq" will not be available)!

Presets

The analysis presets are a set of keywords (only one can be used) that predefine some of the parameters of the metaseqr pipeline. For the time being they are quite simple and they control i) the strictness of filtering and statistical thresholding with three basic levels ("all", "medium", "strict") and ii) the data columns that are exported, again in three basic ways ("basic", "normal", "full") controlling the amount of data to be exported. These keywords can be combined with a dot in the middle (e.g. "all.basic" to define an analysis preset. When using analysis presets, the following argumentsof metaseqr are overriden: exon.filters, gene.filters, pcut, export.what, export.scale, export.values, exon.stats. If you want to explicitly control the above arguments, the preset argument should be set to NULL (default). Following is a synopsis of the different presets and the values of the arguments they moderate:

• "all.basic": use all genes (do not filter) and export all genes and basic annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=NULL

  • gene.filters=NULL

  • pcut=1

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean")

• "all.normal": use all genes (do not filter) and export all genes and normal annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=NULL

  • gene.filters=NULL

  • pcut=1

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean","sd","cv")

  In this case, the above described arguments become:

  • exon.filters=NULL

  • gene.filters=NULL

  • pcut=1

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2","log10","vst")

  • export.values=c("raw","normalized")

  • export.stats=c("mean","median","sd","mad","cv","rcv")

• "medium.basic": apply a medium set of filters and and export statistically significant genes and basic annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=5,min.exons=2,frac=1/5))

  • gene.filters=list(length=list(length=500), avg.reads=list(average.per.bp=100,quantile=0.25), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.05

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean")

• "medium.normal": apply a medium set of filters and and export statistically significant genes and normal annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=5,min.exons=2,frac=1/5))

  • gene.filters=list(length=list(length=500), avg.reads=list(average.per.bp=100,quantile=0.25), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.05

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean","sd","cv")

  and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=5,min.exons=2,frac=1/5))

  • gene.filters=list(length=list(length=500), avg.reads=list(average.per.bp=100,quantile=0.25), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.05

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2","log10","vst")

  • export.values=c("raw","normalized")

  • export.stats=c("mean","median","sd","mad","cv","rcv")

• "strict.basic": apply a strict set of filters and and export statistically significant genes and basic annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=4,min.exons=2,frac=1/4))

  • gene.filters=list(length=list(length=750), avg.reads=list(average.per.bp=100,quantile=0.5), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.01

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean")

• "strict.normal": apply a strict set of filters and and export statistically significant genes and normal annotation and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=4,min.exons=2,frac=1/4))

  • gene.filters=list(length=list(length=750), avg.reads=list(average.per.bp=100,quantile=0.5), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.01

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2")

  • export.values=c("normalized")

  • export.stats=c("mean","sd","cv")

  and statistics elements. In this case, the above described arguments become:

  • exon.filters=list(min.active.exons=list(exons.per.gene=4,min.exons=2,frac=1/4))

  • gene.filters=list(length=list(length=750), avg.reads=list(average.per.bp=100,quantile=0.5), expression=list(median=TRUE,mean=FALSE,quantile=NA,known=NA,custom=NA), biotype=get.defaults("biotype.filter",org[1]))

  • pcut=0.01

  • export.what=c("annotation","p.value","adj.p.value","meta.p.value", "adj.meta.p.value","fold.change","stats","counts")

  • export.scale=c("natural","log2","log10","vst")

  • export.values=c("raw","normalized")

  • export.stats=c("mean","median","sd","mad","cv","rcv")

Note

Please note that currently only gene and exon annotation from Ensembl (http://www.ensembl.org), UCSC and RefSeq are supported. Thus, the unique gene or exon ids in the counts files should correspond to valid Ensembl, UCSC or RefSeq gene or exon accessions for the organism of interest. If you are not sure about the source of your counts file or do not know how to produce it, it's better to start from the original BAM/BED files (metaseqr will use the read2count function to create a counts file). Keep in mind that in the case of BED files, the performance will be significantly lower and the overall running time significanlty higher as the R functions which are used to read BED files to proper structures (GenomicRanges) and calculate the counts are quite slow. An alternative way is maybe the easyRNASeq package (Delhomme et al, 2012). The read2count function does not use this package but rather makes use of standard Bioconductor functions to handle NGS data. If you wish to work outside R, you can work with other popular read counters such as the HTSeq read counter (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html). Please also note that in the current version, the members of the gene.filters and exon.filters lists are not checked for validity so be careful to supply with correct names otherwise the pipeline will crash or at the best case scenario, will ignore the filters. Also note that when you are supplying metaseqr wtih an exon counts table, gene annotation is always downloaded so please be sure to have a working internet connection. In addition to the above, if you have a multiple core system, be very careful on how you are using the restrict.cores argument and generally how many cores you are using with scripts purely written in R. The analysis with exon read data can very easily cause memory problems, so unless you have more than 64Gb of RAM available, consider setting restrict.cores to something like 0.2 when working with exon data. Finally, if you do not wish to download the same annotation again and again when performing multiple analyses, it is best to use the get.annotation function to download and store the resulting data frames in local files and then use these files with the annotation option.

Please note that the meta-analysis feature provided by metaseqr does not satisfy the strict definition of "meta-analysis", which is the combination of multiple similar datasets under the same statistical methodology. Instead it is the use of mulitple statistical tests applied to the same data. For the Simes method, please consult also "Simes, R. J. (1986). "An improved Bonferroni procedure for multiple tests of significance". Biometrika 73 (3): 751–754."

Author(s)

Panagiotis Moulos

Examples

# An example pipeline with exon counts
data("hg19.exon.data",package="metaseqR")
metaseqr(
 counts=hg19.exon.counts,
 sample.list=list(normal="normal",paracancerous="paracancerous",cancerous="cancerous"),
 contrast=c("normal_vs_paracancerous","normal_vs_cancerous",
        "normal_vs_paracancerous_vs_cancerous"),
 libsize.list=libsize.list.hg19,
 id.col=4,
 annotation="download",
 org="hg19",
 count.type="exon",
 normalization="edaseq",
 statistics="deseq",
 pcut=0.05,
 qc.plots=c("mds", "biodetection", "countsbio", "saturation", "rnacomp", 
        "boxplot", "gcbias", "lengthbias", "meandiff", "readnoise","meanvar", 
        "readnoise", "deheatmap", "volcano", "biodist", "filtered"),
 fig.format=c("png","pdf"),
 export.what=c("annotation","p.value","adj.p.value","fold.change","stats",
        "counts"),
 export.scale=c("natural","log2","log10","vst"),
 export.values=c("raw","normalized"),
 export.stats=c("mean","median","sd","mad","cv","rcv"),
 restrict.cores=0.8,
 gene.filters=list(
     length=list(
         length=500
     ),
     avg.reads=list(
         average.per.bp=100,
         quantile=0.25
     ),
     expression=list(
         median=TRUE,
         mean=FALSE
     ),
     biotype=get.defaults("biotype.filter","hg18")
 )
)

# An example pipeline with gene counts
data("mm9.gene.data",package="metaseqR")
result <- metaseqr(
 counts=mm9.gene.counts,
 sample.list=list(e14.5=c("e14.5_1","e14.5_2"), adult_8_weeks=c("a8w_1","a8w_2")),
 contrast=c("e14.5_vs_adult_8_weeks"),
 libsize.list=libsize.list.mm9,
 annotation="download",
 org="mm9",
 count.type="gene",
 normalization="edger",
 statistics=c("deseq","edger","noiseq"),
 meta.p="fisher",
 pcut=0.05,
 fig.format=c("png","pdf"),
 export.what=c("annotation","p.value","meta.p.value","adj.meta.p.value",
        "fold.change"),
 export.scale=c("natural","log2"),
 export.values="normalized",
 export.stats=c("mean","sd","cv"),
 export.where=getwd(),
 restrict.cores=0.8,
 gene.filters=list(
     length=list(
         length=500
     ),
     avg.reads=list(
            average.per.bp=100,
            quantile=0.25
     ),
     expression=list(
            median=TRUE,
            mean=FALSE,
            quantile=NA,
            known=NA,
            custom=NA
     ),
     biotype=get.defaults("biotype.filter","mm9")
 ),
 out.list=TRUE
)
head(result$data[["e14.5_vs_adult_8_weeks"]])

metaseqR documentation built on Nov. 8, 2020, 5:57 p.m.