What does gene-level summarization mean?

Microarray data provides probe-level expression measurements and RNA-seq data provides exon-level or transcript-level (i.e. different isoforms of the same gene) expression measurements. However, current functional annotations are mainly assigned at the gene or protein level. Therefore, when multiple probes or transcripts are mapped to the same gene, they need to be summarized into a single value for that gene. At the Gene Annotation step, users can choose to use the averages or medians of multiple probe intensities (microarray), or sums of counts from multiple transcripts (RNA-seq) to perform gene-level summarization.

What is the rationale behind filtering gene expression data?

The purpose of filtering is to increase the statistical power of differential expression analysis be removing any genes that are less likely to be informative. Please refer to the paper Independent filtering increases detection power for high-throughput experiments for detailed discussion and benchmark tests

Low variance filter: genes whose expression values do not change across different samples, and thus have very low variance. Genes are ranked by their variance from low to high, and you can exclude a certain percentile of genes with the lowest variance by adjusting the "Variance filter" slider. The above referenced study has suggested that up to 50% genes can be removed based on their variance with improved results

Low abundance filter: genes with very low abundance are not measured relaibly and amy not be biologically important. You can exclude genes below a certain threshold by adjusting the "Low abundance" slider. The above referenced study has suggested 10% genes can be removed based on their abundance with improved results

How should I choose a suitable normalization procedure?

Normalizing the data accounts for systematic technical sources of variation so that biologically-driven changes in gene expression can be better detected between samples. A gene expression normalization method should be chosen unless the data has already been normalized, in which case the user should select "None".

All of the normalization methods available on EcoToxXplorer are well-established and have been used in many previous studies. They are based on slightly different assumptions about the underlaying distributions, but should produce relatively similar results. If you are concerned about significant differences between normalization methods, you can try out more than one and visualize the results using the provided plots.

Tip: if you are not sure whether the data is already log transformed or not, you can easily figure this out by visualizing the data (i.e. boxplot). For microarray data, log transformed data values are usually less than 16. For RNA-seq data with 1 million reads, log2(1,000,000) is less than 20. Therefore if all data values are all below 20, it is reasonable to assume that the data has already been log transformed.

How to detect and deal with potential outliers?

Potential outlier samples can be identified from PCA plots. The potential outlier will distinguish itself as the one located far away from the major clusters formed by the remaining samples. To deal with outliers, the first thing is to check if the sample was measured properly. In many cases, outliers are the result of operational errors during the analytical process. If those values cannot be corrected, the sample should be removed from the input data and the analysis re-started.

Which statistical method to use for differential analysis?

Limma is a popular method for differential analysis that was first developed for microarray differential analysis. It addresses the problem of low sample sizes typical to whole-transcriptome studies by using the whole expression profile to make more stable estimates of gene expression variance. EdgeR and DESeq2 were both developed to analyze RNAseq data. All three methods are well-established and should give similar results. Please note:

• EdgeR and DESeq2 are only designed for RNAseq data and will be disabled for microarray data.
• Due to high computational resources required, DESeq2 will be disabled when your dataset contain over 50 samples.

My data contains multiple metadata, how should I choose a proper method for differential analysis?

In differential expression analysis, you should first determine whether any of the metadata encode blocking factors, then decide on how to classify individual samples into groups, and finally decide which groups of samples should be compared to each other using statistical tests. Let's assume that none of your metadata are blocking factors (more on that later) and try to understand how selecting primary and secondary factors creates different groups of samples. Consider the "Estrogen" example data, generated in a study that measured gene expression at multiple time points in breast cancer cells in which the estrogen receptor (ER) was either present or absent. Here, the metadata are "ER" and "TIME". As the figure below shows, selecting "ER" as the primary factor divides the data into two groups because "ER" has two different levels ('present' and 'absent'). Selecting "TIME" as the secondary factor results in four groups because the two primary groups are split based on the two time points. If there were three time points, each primary group would be split into three groups, resulting in six groups overall.

The defined groups can now be compared to find genes that are differentially expressed between them (more details on this in later sections). In some experimental designs, we aren't interested in finding the genes that are differentially expressed between the groups defined by the secondary factor because it is a blocking factor. Examples of blocking factors are subject IDs when multiple samples were taken from the same subject (e.g. paired samples, multiple tissue types), or batches of samples that were measured at different times or in different locations. If you indicate that your secondary factor is a blocking factor, EcoToxXplorer will conduct comparisons within the groups that it defines, which typically improves the accuracy of the overall result.

I received the error message "no residual degrees of freedom", what should I do?

This means you do not have enough samples to perform the analysis you specified, usually when combining two metadata in an independent two-factor analysis (no blocking factors). In this case, the total number of groups will be the product of the number of levels in each metadata factor (i.e. if the primary metadata contains 3 levels, and the secondary metadata contains 4, the total number of groups will be 3 * 4 = 12). We recommend a minimum of 3 samples per group, therefore at least 36 samples are required in order to perform a 3 x 4 two-factor analysis.

In this case, you should focus on a single primary metadata and leave the seconday metadata as "Not available", and perform differential analysis with regard to individual metadata. You can then choose the other metadata as the primary metadata and perform the analysis again. If there are no or very few significant genes identified, it is most likely that incorporating the secondary metadata into the analysis will not affect the result.

What are the differences between pair-wise comparisons and time-series comparisons?

A pair-wise comparison tests for genes that are differentially expressed between any pair of groups. For example, take three groups A, B, and C. The "all pairwise" comparison will contrast A-B, A-C, and B-C. A time-series comparison will only contrast consecutive pairs of groups, so in our example only A-B and B-C. Time-series are commonly used when gene expression was measured at multiple time points, or after treatments with varying concentrations/durations.

What is a nested comparison?

A nested comparison allows you to determine which genes respond differently to a treatment condition, respective to some other metadata. For example, consider the experimental design described in the section on multiple metadata where cells with and without an ER were measured at 10hrs and at 48 hrs. To find the genes that respond differently over time in the ER vs. noER cells, you would perform a nested comparison. First, compare ER10-ER48 to find the genes that are differentially expressed in cells with an ER (ERgenes). Next, compare noER10-noER48 to find the genes differentially expressed in cells with no ER (noERgenes). Finally, to find the genes that respond differently over time in ER vs. noER cells, compare ERgenes-noERgenes.

Selecting "Interaction only" will return significant results from only the ERgenes-noERgenes contrast. Otherwise the full model is returned, which is the combination of significant genes from the ER10-ER48, noER10-noER48, and the ERgenes-noERgenes contrasts.

What is the difference between PCA and t-SNE?

PCA (principal component analysis) and t-SNE (t-distributed stochastic neighbor embedding) are both popular dimension reduction techniques. In PCA, each principal component is the linear combination of predictor variables that explains the greatest amount of variability in the outcome variable, after accounting for previously computed principal components. Unlike PCA, t-SNE utilizes random walks to estimate non-linear relationships between predictor variables for each sample. This means that each iteration of t-SNE will generate slightly different results.

How to use and interpret 3D PCA and t-SNE visualization?

The interactive PCA visualization summarizes all the data into the the first three principal components (PCs). Each data point in the Scores Plot represents a sample. Samples that are close together are more similar to each other. The colors of these data points are based on the factor labels. Users can change the colors according to any of the two factor labels.

Each data point in the Loadings Plot represents a feature. When scores and loadings plots are viewed from the identical perspective, the direction of separation on the scores plot can be explained by the corresponding features on the same directions - i.e. features on the two ends of the direction contribute more to the pattern of separation.

How are gene-level benchmark doses (BMDs) calculated?

Gene-level BMDs are calculated using six steps.

1. Multiple statistical models are fit to the gene expression values for each gene
2. The model fits are filtered to remove poor fits using the lack-of-fit p-value
3. From the remaining fits, the best-fit model for each gene is selected, thus different genes can have different best-fit models
4. Benchmark responses are computed for each gene based on the expression of the control samples
5. The BMD is calculated as the dose that corresponds to the benchmark response
6. BMDs are filtered to remove those with low confidence

More details are given for each of these steps in the below FAQs.

Where do the statistical models come from and what are their exact equations?

The statistical models are the same models the the US Environmental Protection Agency and the National Toxicology Program have used in their studies on 'omics dose-response analysis. They are the same models available in BMDExpress, and are the ones recommended in the peer-reviewed report that outlines the National Toxicology Program Approach to Genomic Dose-Response Modeling. Exp2 - Exp4 are four different forms of exponential models. Poly2, Poly3, and Poly4 are polynomial models with degree 2, 3, and 4. More details and the full mathematical forms can be found in the NTP report linked above.

Which statistical models should I select?

The NTP recommendations say to use all of the models other than Poly3 and Poly4. While these higher degree polynomials often give good fits, there are concerns that allowing too many changes of direction may over fit the data. However, of the remaining models (Exp2-5, Linear, Polynomial degree 2, Power, and Hill), only Poly 2 allows for non-monotonic behaviour. Thus, if you expect non-monotonic behaviour, you may wish to include higher order polynomials. The NTP report suggests that future implementations of dose-response software should include the ability to constrain Poly3 and Poly4 models to only change direction once. We plan to include this feature in future updates to EcoToxXplorer.

What is a lack-of-fit p-value and which threshold should I choose?

A statistical model has a "lack-of-fit" if it fails to adequately explain the relationship between the x (dose) and y (response) variables. In a lack-of-fit statistical test, the null hypothesis is that the model fits the data. Thus, in these tests a significant p-value indicates that there is evidence that the model does not adequately explain the data, and so here we check that the p-value is greater than the significance threshold . Significance thresholds commonly range from 0.05 to 0.5 depending on the desired stringency. The NTP recommendations give a threshold of 0.10.

How is the best fit model selected for each gene?

After applying the lack-of-fit p-value threshold, there may be several statistical models that fit the data well. From these remaining models, the one with the lowest AIC (Akaike information criterion) is selected as the best fit. The AIC is a measure of prediction error that penalizes models with more parameters. This means that if there are two models that do an equally good job of explaining the data, the model with fewer parameters will be selected.

The AIC is not displayed in the results table on the curve fitting page, but it is included in the bmd.csv file that can be downloaded from the Analysis Pipeline side panel.

What is the BMR factor?

BMR means "benchmark response", which is the pre-determined response level that is considered "adverse" or "significant". The dose that corresponds to the BMR is defined as the benchmark dose (BMD). In an ideal case, we would know which change in response variable is physiologically significant and potentially toxic, however this is rarely known for individual genes. For 'omics dose-response analysis, the NTP recommended approach defines the BMR for a gene as the mean of the control gene expression values, plus or minus a certain number of standard deviations of the control values. The number of standard deviations is the BMR factor, and increasing this parameter increases the absolute value of the BMR for each gene.

In the figure above, the BMR factor is one since this is the number of standard deviations that was used to calculate the BMR.

How are the BMD_l and BMD_u calculated for each gene?

The BMD_l and BMD_u are the lower and upper limits of the 95% confidence interval of the BMD, computed using the profile liklihood method. The BMD_l is sometimes used instead of the BMD since it is a more conservative estimate that depends on the uncertainty of the model fit to the data.

Why are there some genes with model fits but no BMDs?

There are several quality criteria applied to the BMDs to filter out low-confidence or otherwise undesirable BMD estimates:

• Confidence interval width: a wide confidence interval indicates a low-confidence estimate of the BMD. BMDs are eliminated if the BMD_u/BMD_l is greater than 40.
• High dose extrapolation: BMDs that are higher than the highest dose are sometimes due to unpredictable model behaviour beyond the measured data (especially with higher order polynomials). BMDs are eliminated if they are greater than the highest dose.

Additionally, for some genes, the expression values never exceed the BMR and thus EcoToxXplorer cannot compute a BMD.

What is a transcriptomic BMD (omicBMD) and how is it calculated?

The intial curve fitting and BMD calculation are done at the gene level (geneBMD). However, the main statistic of interest is usually the dose at which the whole-transcriptome is responding to chemical exposure, or the omicBMD. The NTP report defines the omicBMD based on pathway enrichment analysis of the geneBMDs (see the Pathway BMD Analysis FAQ section for more details). However, with many ecological species these pathway-based omicBMDs are unstable due to the low number of annotated gene sets compared to popular mammalian model organisms, and thus it has been proposed to estimate a statistical omicBMD based on the distribution of geneBMDs. There are three different methods for estimating statistical omicBMDs in EcoToxXplorer:

• omicBMD₂₀: the omicBMD is equal to the 20th lowest geneBMD value
• omicBMD_10th: the omicBMD is equal to the 10th percentile of the geneBMD values
• omicBMD_mode: the omicBMD is equal to the mode of the first peak in the distribution of geneBMD values

How are pathway BMDs (pathBMDs) calculated?

Gene set analysis is used to identify significantly overrepresented pathways in the list of geneBMDs. Pathway-level BMDs (pathBMDs) are calculated as the bootstrapped median of the geneBMDs in that pathway. The National Toxicology Program Approach to Genomic Dose-Response Modeling defines the omicBMD as the lowest pathBMD.

How is the gene set analysis done?

EcoToxXplorer uses overrepresentation analysis (ORA) to compute gene set enrichment statistics. ORA is a statistical technique to identify gene sets or pathways that have a significant overlap with a gene list of interest. In EcoToxXplorer, Hypergeometric tests are used to compute the p-values.

Should I change the filters for the enrichment analysis results?

The NTP approach recommends no α cut-off, but to only consider pathways with at least 3 genes or 5% of the total pathway genes.

Which gene set libraries are available for pathway BMD analysis?

EcoToxXplorer supports enrichment analysis with gene sets from the Gene Ontology, PANTHER, KEGG, Reactome, and MSigDB databases. Note - not all gene set libraries are available for all species.

The GO:BP, GO:MF, and GO:CC gene sets include the complete set of Gene Ontology terms (> 45 000) for the biological process, molecular function, and cellular component categories. The PANTHER:BP, PANTHER:MF, and PANTHER:CC are reduced sets of GO terms ("GO slims") that have been manually chosen based on the PANTHER protein classification system. Briefly, the PANTHER project has created > 15 000 phylogenetic trees that encode the evolutionary relationships within protein families. Subsets of GO terms were chosen that best reflect the function gain or loss along the branches of the PANTHER trees for each of the BP, MF, and CC categories. In general, GO slims can simplify the interpretation of enrichment analysis results because they reduce the number of highly similar GO terms.

The KEGG and Reactome gene sets are networks of molecular interactions that represent biological pathways and processes. Reactome pathways are created through a process similar to scientific peer review, where different experts create and review the pathway organization, and all interactions contain references to the primary literature. KEGG pathways are also based on molecular interactions in the primary literature, but are accompanied by an extensive ortholog mapping that allows KEGG pathways to be rapidly extended to additional species based on genome sequence homology.

What do the colours on the pathway heatmaps represent?

The pathway heatmap values are calculated through a series of steps:

1. The fitted model for each gene is evaluated across the range of doses in the uploaded data.
2. The resulting modeled expression values are normalized by subtracting the minimum and then dividing by the maximum, producing a set of numbers between zero and one.
3. Genes with a BMD that is down-regulated compared to the control are identified. Their values are inverted by subtracting 1 and then multiplying by -1, which creates a consistent colour gradient from blue on the left to red on the right. The direction of regulation is noted by the coloured boxes on the left side with red for up-regulated and green for down-regulated.
4. Genes are ordered by their BMD.
5. The heatmap cells where geneBMDs fall are coloured black.

The purpose of these steps are to produce a visually appealing heatmap that clearly shows how the expression of each pathway gene compares to the others.