# data[].attributes.attachment_link.readable_text
Research

Open access

Published:

Show authors

BMC Bioinformatics 21, Article number: 348 (2020) Cite this article

5095 Accesses

1 Citations

15 Altmetric

Metrics

Abstract

Background

Bioinformatics has pervaded all fields of biology and has become an indispensable tool for almost all research projects. Although teaching bioinformatics has been incorporated in all traditional life science curricula, practical hands-on experiences in tight combination with wet-lab experiments are needed to motivate students.

Results

We present a tutorial that starts from a practical problem: finding novel enzymes from marine environments. First, we introduce the idea of metagenomics, a recent approach that extends biotechnology to non-culturable microbes. We presuppose that a probe for the screening of metagenomic cosmid library is needed. The students start from the chemical structure of the substrate that should be acted on by the novel enzyme and end with the sequence of the probe. To attain their goal, they discover databases such as BRENDA and programs such as BLAST and Clustal Omega.

Students' answers to a satisfaction questionnaire show that a multistep tutorial integrated into a research wet-lab project is preferable to conventional lectures illustrating bioinformatics tools.

Conclusion

Experimental biologists can better operate basic bioinformatics if a problem-solving approach is chosen.

Background

At present any biologists should capitalize on the resources, data and programs, that are available online to make their experimental plans more efficient and cost-effective [ 1, 2]. For this reason, it is desirable to train students using a problem-solving approach that integrates in silico work into a multidisciplinary experimental project(a few examples [ 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. We provide a tutorial that was administered to students of the courses of marine and environmental biology.

Biodegradation of environmental pollutants by marine prokaryotic enzymes provides the frame into which the bioinformatics tutorial is inserted. In particular, the experimental project proposed to the students, aims at finding an enzyme that is active on a scaffold commonly found among pollutants and synthetic compounds, a so-called "privileged scaffold" [ 14]. In this exercise, the scaffold taken into consideration is indole, an N-heterocyclic aromatic pollutant released in the aquatic environment through the industrial wastewater [ 15].

The students will look for an enzyme with broad specificity that is able to degrade indole and more in general aromatic compounds [ 16]. As it is widely known, only a small fraction of environmental microbes grow under conventional laboratory conditions [ 17, 18, 19]. For this reason, several authors suggested that metagenomes might be a big reservoir of novel enzymes for applications in biocatalysis, biofuels, and bioremediation (for reviews on this subject [ 20, 21, 22, 23]).

The project consists of two parts: a laboratory part, that we won't discuss in detail here, to construct a DNA clones library from the metagenomic marine sample, and a bioinformatics part that is the object of this paper.

We presuppose that DNA has been extracted from prokaryotes present in a seawater sample and that a metagenomic library has been constructed in cosmids. In principle, the screening could be carried out testing the enzymatic activity of interest directly. However, assays on plates can be cumbersome and we propose an alternative strategy that takes advantage of the simplicity of DNA colony hybridization. Only after the clone has been identified by a DNA labelled probe and isolated, the activity will be confirmed in the recombinant E.coli extracts.

The tutorial focuses on the design of a suitable probe for the screening. It combines information deriving from a highly annotated enzyme database such as BRENDA [ 24] with data on uncharacterized open reading frames (ORF) deriving from large scale metagenomic sequencing projects.

BRENDA is freely available for academic users and educational purposes [ 25]. Enzymes are classified according to the catalyzed reactions and well-characterized proteins from different organisms can be found in each class. BRENDA [ 24] can be searched with a structure-based query as well as with a text-based query. It provides a substructure search algorithm that is very useful to use a scaffold as a query. This choice is very convenient because chemical structures identify molecules uniquely whereas names are not unique. For instance, indole is also known as 1-Benzazole,1H-Benzo [b]pyrrole,1H-Indol,1H-Indole, nomenclature can be even more confusing for other molecules. BRENDA can be searched for several specific information at the same time with "Advanced options". The tutorial shows how a bacterial enzyme that is active on a scaffold commonly found among pollutants and synthetic compounds can be found. We imagine that the sequence cannot be used as such for colony hybridization of the cosmid library under stringent conditions. Hence we propose to use the enzyme found in BRENDA [ 24] as a query sequence to look for homologous proteins among uncharacterized ORF from marine metagenomes. This is possible because BLASTp [ 26] consents to search protein sequences from large environmental sequencing projects such as the Malaspina expedition [ 27], Global Ocean Sampling (GOS) campaign [ 28] and Tara Oceans expedition [ 29].

The tutorial was administered to 23 students with no previous knowledge of bioinformatics who filled a satisfaction test at the end of the exercise.

Methods

Aims

This tutorial is intended for laboratory biologists with no previous knowledge of bioinformatics. We envisaged a multistep bioinformatics protocol that is integrated into a project of marine biotechnology. The protocol focuses on the use of four bioinformatics web applications (BLASTp [ 26], Clustal Omega, Cons, Reverse Translate [ 30]) and two biochemical databases (UniProt [ 31] and BRENDA [ 24]). The learning goals are summarized in Table 1.

Experiment overview

We presuppose that the students are engaged in a project for the identification of a novel bacterial enzyme from a seawater metagenomic sample. For this purpose, we assume that a library has been constructed in cosmids and must be screened. In this tutorial, we present a bioinformatics protocol to design a probe to isolate the clone of interest (Fig. 1).

In the first step, the students will search BRENDA [ 24], the enzyme database [ 25], to find a prokaryotic enzyme that is able to use a given class of pollutants as substrates. BRENDA [ 24] consents looking among enzyme ligands by chemical similarity and by substructures. In the tutorial, indole is chosen as an example because it represents a privileged scaffold, i.e. its chemically active structure is common to many natural and synthetic compounds with the ability to bind different targets [ 14]. Unfortunately, it is not (yet) possible to carry out an advanced search drawing a substructure and for this reason, to get an enzyme class that is active on a given substrate, identified by its chemical structure, AND is expressed in bacteria, it is necessary to follow a two-phase, apparently redundant, protocol. In the first phase, the structure is used as the input to get the exact name (i.e. the one used by BRENDA that does not necessarily coincides with the IUPAC name) of the substrate of interest. In the second phase, the exact name and the class of organisms are used with advanced options. Several classes of enzymes are active on molecules that resemble indole. Naphthalene 1,2-dioxygenase has broad specificity and is involved in the degradation of many aromatic compounds [ 32]. By clicking on the E.C. number the students will get much information that is subdivided into different sections. The most interesting ones for the case under study are those concerning enzyme-ligands interactions where it is possible to learn which are all the possible substrates of this class on enzymes, those concerning organism related information and those related to enzyme structures, where the links to UniProt [ https://www.uniprot.org/] [ 31] are found. The students will not find any Naphthalene 1,2-dioxygenases from a prokaryotic marine organism. At present (March 2019) the only well-characterized prokaryotic enzyme is from Pseudomonas putida [ 33].

In the second step, the sequence of Naphthalene 1,2-dioxygenases from Pseudomonas putida will be obtained from UniProt [ https://www.uniprot.org/] [ 31]. Pseudomonas putida is evolutionarily distant from the marine prokaryotes that are present in the sample of marine water used to construct the library. The students will look for homologous sequences from marine prokaryotes. They will take advantage of a large number of uncharacterized coding sequences obtained by massive genomic and metagenomic sequence projects. To carry out such an analysis, they will use BLASTp [ https://blast.ncbi.nlm.nih.gov/Blast.cgi? PROGRAM = blastp&PAGE_TYPE = BlastSearch&LINK_LOC = blasthome] using Naphthalene 1,2-dioxygenases from Pseudomonas putida as the query and limiting the search to the metagenomic proteins deposited in databases (env_nr).

In the third step, a multiple sequence alignment will be carried out using Clustal Omega [ https://www.ebi.ac.uk/Tools/msa/clustalo/] [ 34]. The students will identify conserved regions among the sequences homologous Naphthalene 1,2-dioxygenases from Pseudomonas putida and will derive a consensus using Cons [ http://www.bioinformatics.nl/cgi-bin/emboss/cons] [ 35, 36, 37].

In the fourth step, the retrotranslation of the aminoacidic consensus sequence will be carried out using [ https://www.bioinformatics.org/sms2/rev_trans.html] [ 30]. The students will learn that retrotranslation does not provide a unique DNA sequence unless the codon most frequently used for each amino acid in prokaryotes is chosen.

The output of this tutorial is a sequence of the probe for the screening by colony hybridization.

Requirements

All software used is free with a user-friendly interface available. The only requirement is a computer with an Internet connection.

The used software and databases are:

BRENDA, https://www.brenda-enzymes.org/index.php

BLASTP, https://blast.ncbi.nlm.nih.gov/Blast.cgi

UniProt, http://www.uniprot.org/

Clustal Omega, https://www.ebi.ac.uk/Tools/msa/clustalo/

Cons, http://www.bioinformatics.nl/ emboss-explorer/

Reverse translate, http://www.bioinformatics.org/sms2/rev_trans.html

Detailed protocol

Step 1. BRENDA: searching the enzyme database to find a microbial enzyme able to degrade indole

Access BRENDA (The Comprehensive Enzyme Information System) at https://www.brenda-enzymes.org/index.php.

Draw the molecule of interest, a bicyclic structure, consisting of a six-membered benzene ring fused to a five-membered pyrrole ring, clicking on "Ligand Structure Search" in the BRENDA homepage (Fig. 2) to obtain the exact name of the compound. A "substructure search" with a maximal search time of 120 s restricted to "Substrates" must be selected before running the search Fig. 3).

You obtain the exact names of several molecules containing the structure you drew, you choose "indole". Go back to BRENDA homepage and run an advanced search (Fig. 4) filling in the kingdom (bacteria) and type (substrate) boxes and using the exact name of the molecule (indole) (Fig. 5).

BRENDA shows a list of all the enzymes in the database that use indole as a substrate. The enzyme Naphthalene 1,2-dioxygenase from Pseudomonas has been chosen from the results for this training. It is present in different bacterial species and has a large substrate specificity. Please click on the EC number (Fig. 6) to obtain all the biochemical data about the enzyme.

In the section "enzyme structures", look for "AA sequences" where the links to UniProt [ https://www.uniprot.org/] [ 31] are found (Fig. 7). Pseudomonas putida has been chosen as an organism; click on the UniProt code to open the link to the proteins sequences database.

On the UniProt sheet choose Sequence and download the amino acid sequence in FASTA format (Fig. 8).

Step 2. Protein BLAST (BLASTp) to carry out a local alignment to compare an amino acid sequence to the metagenomic proteins deposited into a database

Open a new Internet page with BLAST at https://blast.ncbi.nlm.nih.gov/Blast.cgi and click on "Protein BLAST" to align amino acid sequences.

Paste the sequence downloaded from UniProt into the "Enter Query Sequence" field, choosing "Metagenomic proteins" as database. In the section "Algorithm parameters", select 500 as "Max target sequences" and run BLAST (Fig. 9).

Select and download all "marine metagenome" sequences in FASTA format (Fig. 10).

Step 3. Clustal omega to perform multiple sequence alignment, helpful to predict relations and similarity among sequences

Open a new Internet page with Clustal Omega tool at https://www.ebi.ac.uk/Tools/msa/clustalo/ to perform a multiple alignment.

Paste your sequences in the blank using Pearson/FASTA as output format, leave the other parameters as default and submit the work (Fig. 11).

Step 4. Cons: an EMBOSS explorer tool to create a consensus sequence from a multiple alignment

Connect to EMBOSS Explorer [ 35, 36, 37] at http://www.bioinformatics.nl/emboss-explorer/ and on the left side of the web sheet click on "cons". Copy the obtained multiple alignment in FASTA format and paste it in the blank (Fig. 12). An example of consensus sequence output is shown in Fig. 13.

Copy a region (selecting about 30 amino acid residues to obtain a ~ 100 bp DNA probe avoiding regions containing too many "Xs").

Step 5: translate an amino acid sequence into the most probable nucleotide sequence

Search for the web tool Reverse Translate ( http://www.bioinformatics.org/sms2/rev_trans.html) to retrotranslate an amino acid sequence into the most probable nucleotide sequence using the default codon usage from E.coli.

Paste the copied sequence in the blank and submit the work (Fig. 14). In Fig. 15, an example of the obtained nucleotide probe is shown.

Results

The tutorial was administered to 23 graduate students who had a bachelor's degree in biology or natural science with a basic knowledge of biochemistry and molecular biology and no previous knowledge of bioinformatics. This tutorial can be administered to undergraduate students too if they have solid bases of biochemistry and molecular biology. Before administering the tutorial, two lectures were given to introduce the programs and the databases reported in Table 1. All the students without the help of the supervisor concluded the tutorial successfully within one hour. After the practice, students have undergone a satisfaction questionnaire of 7 questions (Additional file 1) to evaluate the general interest and the usefulness of the multistep bioinformatics protocol. All students were satisfied although a minority admitted that they would not be able to apply the same tools to another biological project (30%). Questions indicating the student satisfaction degree are gathered in Fig. 16 and in the supplementary material. Through a computer-based approach, students have been able to search into protein and enzyme databases, performing local and multiple sequence alignments obtaining a consensus sequence and retrofitranslate an amino acid sequence to obtain a DNA sequence.

We are aware of the fact that other types of questionnaires could be administered, for instance, students could be asked to explain the steps that they have performed and why.

An even more convincing proof of the usefulness of the tutorial would have been letting the students prepare the cosmid library and, after that, dividing them in two groups. One group would have been asked to screen the library with an enzymatic assay on plates, the second group would have followed the tutorial, designed the probe for colony hybridization and confirmed the activity in the recombinant E.coli extracts.

Discussion

Table 1 summarizes the learning goals of the tutorial. Beyond the immediate ones, i.e. becoming acquainted with very popular bioinformatics tools (Table I) and, in particular, of BRENDA [ 24], there are far-reaching educational aims. In our opinion the students should recognize the major opportunity offered by big data produced by metagenomics projects and the possibility of deriving protein properties by homology.

The tutorial emphasizes the opportunity offered by metagenomic next-generation sequencing projects. They are like mines we can dig to find what is useful for our research. It is highly probable that proteins homologous to the ones we are interested in are present in the uncharacterized big data stored in databanks.

The tutorial highlights the concept that homologous enzymes share similar functions and that homology can be found comparing sequences. It capitalizes on the principle 'Pairwise alignments whisper while multiple alignments shout out loud' (Arthur Lesk). In fact, the students will identify the conserved regions that are a better fit to design a probe carrying out a multiple alignment of proteins homologous to Naphthalene 1,2-dioxygenase from Pseudomonas putida.

Consulting BRENDA [ 24] permits to retrieve original research papers concerning specific classes of enzymes. This will be precious to set up enzymatic assays and check that the cosmid clone isolated by colony hybridization indeed expresses Naphthalene 1,2-dioxygenase.

Conclusions

Computer-based learning is an excellent method to introduce undergraduate and graduate students with biological and biotechnological background to bioinformatics. We have shown that using bioinformatics tools as steps of a research project is more useful than presenting the same tools separately in stand-alone tutorials and that having a clear experimental objective, possibly related to a "trendy" topic such as green economy [ 38] raises the students' interest.

Availability of data and materials

The programs used during the current practice are available at:

BENDA, https://www.brenda-enzymes.org/index.php Accessed 03/09/2019

BLASTP, https://blast.ncbi.nlm.nih.gov/Blast.cgi Accessed 03/09/2019

UniProt, http://www.uniprot.org/ Accessed 03/09/2019

Clustal Omega, https://www.ebi.ac.uk/Tools/msa/clustalo/ Accessed 03/09/2019

Cons, http://www.bioinformatics.nl/ emboss-explorer/ Accessed 03/09/2019

Reverse translate, http://www.bioinformatics.org/sms2/rev_trans.html Accessed 03/09/2019

References

Download references

Acknowledgements

About this supplement

This article has been published as part of BMC Bioinformatics Volume 21 Supplement 10, 2020: Proceedings from the 13th Bioinformatics and Computational Biology International Conference - BBCC2018. The full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume-21-supplement-10 .

Funding

Publication costs are funded MIUR PRIN 2015 2015JHLY35 to M.V.C. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations

1. Dipartimento di Scienze e Tecnologie Ambientali, Biologiche e Farmaceutiche, Università degli Studi della Campania "Luigi Vanvitelli", 81100, Caserta, Italy

Ludovica Liguori & Mariateresa Allocca

2. Istituto di Chimica Biomolecolare -CNR, 80078, Pozzuoli, Italy

Ludovica Liguori & Mariateresa Allocca

3. Dipartimento di Biologia, Università Federico II, 80126, Naples, Italy

Maria Monticelli, Maria Vittoria Cubellis & Bruno Hay Mele

4. Biology and Evolution of Marine Organisms Department, Stazione Zoologica Anton Dohrn, Villa Comunale, 80121, Naples, Italy

Maria Vittoria Cubellis

5. Integrative Marine Ecology Department, Stazione Zoologica Anton Dohrn, Villa Comunale, 80121, Naples, Italy

Bruno Hay Mele

Contributions

L.L performed the experiments, gathered the questionnaires, analyzed the data. M. M and M.A. performed the experiments and prepared the Figs. B.H.M conceived the experiments. M.V.C. and L.L. wrote the paper. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Maria Vittoria Cubellis.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver ( http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Reprints and permissions

About this article

Cite this article

Liguori, L., Monticelli, M., Allocca, M. et al. Bioinformatics tools for marine biotechnology: a practical tutorial with a metagenomic approach. BMC Bioinformatics 21 (Suppl 10), 348 (2020). https://doi.org/10.1186/s12859-020-03560-z
Download citation
Keywords

# data[].attributes.attachment_link.title
Bioinformatics tools for marine biotechnology: a practical tutorial with a metagenomic approach - BMC Bioinformatics

# data[].attributes.attachment_link.description
Background Bioinformatics has pervaded all fields of biology and has become an indispensable tool for almost all research projects. Although teaching bioinformatics has been incorporated in all traditional life science curricula, practical hands-on experiences in tight combination with wet-lab experiments are needed to motivate students.

# data[].attributes.attachment_link.readable_text
Methodology article

Open access

Published:

BMC Bioinformatics 12, Article number: 228 (2011) Cite this article

11k Accesses

23 Citations

Metrics

Abstract

Background

In genetic transcription research, gene expression is typically reported in a test sample relative to a reference sample. Laboratory assays that measure gene expression levels, from Q-RT-PCR to microarrays to RNA-Seq experiments, will compare two samples to the same genetic sequence of interest. Standard practice is to use the log ₂-ratio as the measure of relative expression. There are drawbacks to using this measurement, including unstable ratios when the denominator is small. This paper suggests an alternative estimate based on a proportion that is just as simple to calculate, just as intuitive, with the added benefit of greater numerical stability.

Results

Analysis of two groups of mice measured with 16 cDNA microarrays found similar results between the previously used methods and our proposed methods. In a study of liver and kidney samples measured with RNA-Seq, we found that proportion statistics could detect additional differentially expressed genes usually classified as missing by ratio statistics. Additionally, simulations demonstrated that one of our proposed proportion-based test statistics was robust to deviations from distributional assumptions where all other methods examined were not.

Conclusions

To measure relative expression between two samples, the proportion estimates that we propose yield equivalent results to the log ₂-ratio under most circumstances and better results than the log ₂-ratio when expression values are close to zero.

Background

Several different bioinformatics technologies exist to quantify gene expression. Regardless of technological platform, laboratory assays of gene expression first extract mRNA from a test sample and a control sample. These samples may be labeled with a tag or dye and hybridized to amplified cloned sequences that represent a gene of interest. The amount of mRNA in each sample is usually measured by examining the amount of dye remaining after hybridization. Researchers use Q-RT-PCR to measure expression when there are only one or a few genes of interest. Several lab protocols from various companies exist to quantify gene expression such as RT-PCR assays using intercalating dyes like SYBR Green, the TaqMan Gene Expression Assays, LightCycler, and QuantiGene [ 1- 3]. When genome-wide levels of expression are of interest, microarrays can measure expression for thousands of genes of interest. Microarray platforms employ either cDNA clones [ 4, 5] or n-mer oligonucleotide probes for many genes at once [ 6].

More recently, sequence-based technologies provide more efficient and accurate expression measurements on a genome-wide scale. Evolving from early techniques such as Serial Analysis of Gene Expression (SAGE) to modern techniques such as Massively Parallel Signature Sequencing (MPSS) and RNA Sequencing (RNA-Seq), these approaches now rival microarray-based gene expression analysis for efficiency, cost, and accuracy [ 7]. Sequence-based techniques are also more flexible, allowing for gene expression measurements on a genome-wide level from any organism with a published genome sequence [ 8]. Sequencing employs systems such as the 454 or Illumina platform with the latter demonstrating greater depth and coverage [ 9]. To illustrate the central motive of this paper, Figure 1 demonstrates a two-color competitive hybridization assay of the kind used in TaqMan assays and cDNA microarrays. Other methods involve single-dye hybridization systems or intercalating dyes that bind to double-stranded DNA (dsDNA) product. The statistical models proposed below can be generalized to any scenario where gene expression is measured comparatively in a test sample and a reference sample.

The competitive hybridization process for a two-color system: The number of PCR products equals the number of possible hybridizations. A proportion of the sequences will bind with matching red labeled strands and the remainder bind with the matching green labeled strands. Some sequences will not match (marked with X's) and should not hybridize.

Researchers commonly use the log ₂-ratio to measure relative mRNA expression between two samples. The estimate is as follows. Let R ij represent a summary expression value for gene j in the reference sample i where i = 1,..., n and j = 1,..., K. Let G ij represent a summary expression value for gene j in the test sample i. The value n is the number of paired samples or experiments and K is the number of genes studied. To summarize relative expression between two samples, the log ₂-ratio is

or other similar variants on the theme. The log ₂-ratio is commonly interpreted as the average "log-fold-change" in gene expression between the reference sample and the test sample. Its estimate will be denoted by . If r j = 1, then the ratio between the two samples is 2 ¹ = 2, meaning that the expression of gene j in the test sample is two-fold that of the reference sample on average. If r j = 2, then the ratio between the two samples is 2 ² = 4, meaning that on average the expression in the test sample is four-fold that of the reference sample. Other values of r j are interpreted similarly.

While the interpretation of the log ₂-ratio is appealing, the statistic has an important drawback. When expression in the reference sample is low, is numerically unstable because the denominators R ij are small. As R ij approaches zero, r j increases drastically, approaching infinity. When R ij = 0, then r j is undefined. Thus, when reference sample expression is low, we get extreme estimates or missing values for r j . This phenomenon is especially common when measuring gene expression in simple organisms. In bacteria, for example, transcription may be binary; either on or off. The log ₂-ratio is least reliable for these systems. This problem persists in human genomics research for certain experimental conditions and genes of interest.

This article proposes a new estimate to compare mRNA expression in two samples. This estimate is the proportion of mRNA in the test sample p j for each gene. The proportion takes the amount of mRNA in the test sample and compares it to the total amount of expressed mRNA represented by the sum of the test and reference samples. One formula for estimating the proportion is

The proportion is well-defined for all values of R ij and G ij . For example, when R ij = 0, then G ij/( G ij + R ij ) = 1. We can interpret the number as follows: mRNA expression is observed only in the test sample and not in the reference sample. Similarly, if G ij = 0, then G ij/( G ij + R ij ) = 0 and this means that mRNA expression is observed in the reference sample only. Figure 2 demonstrates the relationship between the log ₂-ratio and the proportion estimates, which follows a logistic function. The relationship is roughly linear near the center point but non-linear at the extreme values. A detailed description for the estimate of the proportion, , and an alternative derived from a maximum likelihood estimate, is in the Results section.

The relationship between the log2-ratio r j and the proportion p j.

More generally, p j can be interpreted as the proportion of mRNA from gene j expressed in the test sample. As p j deviates from 0.5, then there is differential expression between the test and reference samples. As p j approaches one, then gene j is up-regulated in the test sample. As p j approaches zero, then gene j is down-regulated in the test sample. The proportion statistic p j can also be transformed into a percentage: p j × 100% for reporting. For example, if p j = 0.75 then we can say that 75% of the mRNA expressed in the experiment comes from the test sample. The proportion estimate can easily be used to test for differential expression between groups. Under the null hypothesis of no gene expression, p j = 0.5. The alternative hypothesis is differential expression, p j ≠ 0.5. The log ₂-ratio estimate requires a different hypothesis test. Under the null hypothesis, r j = 0 and under the alternative, r j ≠ 0.

Using a proportion p j to describe relative expression for gene j instead of the log ₂-ratio r j maintains the ability to interpret differential expression and test for differences. The added benefit of the proportion is the ability to preserve all data points, even for experiments with very low expression values. Typically when values of R ij are very small, researchers eliminate the j th probe of the i th experiment from their analysis. Eliminating missing data results in a loss of information and potential bias and loss of power. The proportion estimate does not require the removal of extreme, but legitimate, data points.

The Results section provides details that describe the estimation of statistics for p j . The section also provides several test statistics for hypothesis tests of p j . Estimation and testing are developed in the frequentist context but the Bayesian context can also be used, as described in the Appendix. The Results section compares the testing scenarios in simulations and two datasets. The first dataset consists of expression data from a cDNA microarray platform and the second dataset uses RNA-Seq. Both the log ₂-ratio and proportion statistics achieve roughly equivalent results under usual conditions, but one of the proportion statistics performs better across a variety of distributional assumptions. Proportion statistics also detect differentially expressed genes that would typically be classified as missing data.

Results

Parameter Estimates and Hypothesis Testing

We propose a new strategy for the comparison of expression values that is tied to the underpinnings of the hybridization process and its natural interpretation using a binomial distribution. Figure 1 illustrates the hybridization process in a way that justifies the use of a binomial distribution. The description is specific to a two-color hybridization platform. The same concept extends to any system where both test samples and reference samples are assayed.

For each gene sequence, suppose that researchers amplify sequences resulting in M ij clones, where j is the gene probe index and i is the sample number, in order to co-hybridize the extracted mRNA sequences from the reference and test samples. Usually M ij is in the millions, but the exact value will be unknown. For each probe, suppose it hybridizes to a test target with probability p j and to a reference target with probability 1 - p j . This reflects the proportion of available test sequences versus references sequences. We assume that each probe must hybridize to mRNA extracted from either the test or reference sample. Then, the number of hybridizing test target sequences Y ij follows a binomial distribution with size M ij and probability p j . We wish to estimate p j to calculate the proportion of hybridized test target sequences. The maximum likelihood estimate for p j is . In this scenario, Y ij = G ij and M ij = G ij + R ij where R ij represents the expression value for gene j in the reference sample i and G ij represents an expression value for gene j in the test sample i when there are i = 1,..., n paired experiments and j = 1,..., K genes. Therefore, to summarize n experiments the estimated proportion for each gene j is

To test for differential expression we set up a decision with the null hypothesis H ₀: p = 0.5 versus the alternative hypothesis H ₁: p ≠ 0.5. The test derived for this binomial distribution has a test statistic

The test statistic z j is compared to a quantile from the normal distribution z 1-α/2. If the type I error is α = 0.05, then z 0.975 = 1.96. If |z j| > 1.96, then gene j is declared differentially expressed between test and reference samples. The z 1-α/2 quantile is replaced by a quantile when the variance estimate uses instead of p j . This test of binomial proportions, however, is not robust to deviations from the binomial distribution. Indeed, we do not believe that expression data will always follow a binomial distribution, but we include this derivation to motivate the choice of this statistic. Instead, we recommend an alternative test statistic that can be used whether distributional assumptions are met or not. The alternative test statistic simply uses a normal approximation to the binomial distribution and calculates a sample variance estimate. Then the test statistic for differential expression is

where we estimate the proportion and the sample variance in the usual way, . If | t j | > t 1-α/2,n-1, then gene j is differentially expressed between test and reference samples. This test is valid for sufficiently large sample sizes (see Table 1).

Calculating corresponding confidence intervals for each of the test statistics above is straightforward. Previous research suggests adjusting confidence intervals for binomial proportions. The most popular adjustment of these intervals uses the Agresti-Coull procedure [ 10, 11]. We recommend this procedure to estimate confidence intervals for both of the proportion estimates above.

The proposed statistics are evaluated within a frequentist framework. A Bayesian framework is provided in the Appendix.

Simulation Results

We ran a series of simulations to compare the inference behavior of proportion based statistics, and , to log-ratio based statistics, and . The proportion statistics and are introduced in equations 2 and 3 above and their test statistics are given in equations 4 and 5. The ratio-based statistics that have been used in the literature previously are described in equation 1 ( ) and equation 6 in the Methods section ( ). In preliminary simulation exercises, we found that the performance of some test statistics was heavily dependent on the distribution used to generate the expression data. Thus, we generated expression data under four different distributions. The simulation results in Table 1 present a subset of the sample sizes and fold changes examined. More extensive tables are in Additional file 1.

The performance of the estimators under four different distributions is summarized in Table 2. The table was created after examination of the empirical type I error and power for each statistic in each simulation (Additional file 1). The proportion performs at or above the others with the exception of for the Poisson, although still performs adequately there (cp. Tables 1 and 2). The statistics and the DESeq analysis exhibit unacceptable performance under one or more of the distributional assumptions. Statistics , and limma/empirical Bayes (EBA) and edgeR analysis are always good or acceptable, depending on the distributional assumptions. The edgeR and DESeq analyses have type I errors less than 0.05 in many instances. On the average, the and EBA tests have moderately better type I error and power than , (Additional file 1). Although might be expected to outperform the other methods under the binomial assumption, detection under this assumption is easy and all methods performed equally well. In conclusion, the statistic and limma/EBA have the best inference in our simulations overall. The empirical Bayes approach results in better power than on the average under the Poisson and normal distributions.

Table 2 Comparison of estimators from the simulations.

Analysis of Gene Expression in Mice with apoAI Knockout

To examine the performance of our method on cDNA microarray data, we analyzed the expression values reported in Ge et al (2003) [ 12]. Since the apoAI experiment was a control-treatment experiment that used a third sample as a reference, this data exhibits how the methods of this paper can be extended to the case of a difference of two proportions. When testing for the difference between control and treatment, the p-values from and were very similar in magnitude. This was true for both raw p-values and p-values adjusted for multiple-testing. The order of the p-values was also similar, but not identical (see Figure 3). When using the limma/EBA method, the p-values from and were again similar in magnitude, although the order varied more after the 7th probe (Table 3). The top 8 most differentially expressed probes from the original analysis differed from those selected by using t-statistics in the 8th probe, although the top 9 probes for both sets are the same (Table 3). In the original analysis, the top 8 probes corresponded to four distinct genes, and were confirmed by real time quantitative PCR [ 13].

Scatterplot of p-values for log 2-ratio and proportion for the mouse data. The general ordering of the genes is similar, although not identical, using the two methods.

When using , there were 158 (2.5%) unanalyzable probes because one or more of the samples had both G ij = 0 and R ij = 0, which made undefined. The statistic was defined for all probes because G ij and R ij were never zero for all samples of a specific probe. For this data, none of the 158 unanalyzed probes were in the top eight when using , although if they were a potential discovery they would have been missed using . To avoid this problem, one may add an arbitrary constant to all probes before taking the log-ratio. If merely raw p-values were selected at α = 0.05, then would have selected 850, but there would have been 9 more significant p-values if an arbitrary 0.05 were added to the data to avoid zero denominators when using log-ratios. By comparison, would have selected 871 probes.

Therefore, is able to give comparable results to for this cDNA microarray experiment, with the slight advantage that it provided information for 158 more probes in the study, without an arbitrary constant.

Analysis of Differential Expression in Human Kidney and Liver Cells

To examine the performance of our methods on RNA-Seq data, we analyzed the expression values reported in Marioni et al (2008) [ 9]. This data compared the expression of human kidney and liver cells sampled from the same person. Concentrations of 3 pM of cDNA were sequenced using the Illumina platform in five lanes. The original paper analyzed the expression of 32,000 sequences and reported that 11,493 of the sequences were differentially expressed with q-values less than 0.001 (FDR < 0.1%) [ 14]. Supplemental Table 3 from Marioni et al (2008) provides the results of 17,708 sequences analyzed with both RNA-Seq technology and Affymetrix microarrays. They reported that 8,113 of Affyymetrix probe sets were differentially expressed with q-values less than 0.001.

In order to compare the methods in the original paper to those we are proposing, we used a type I error rate of α = 0.05/32000 for all tests. In this way, the threshold can be universally applied to all genes and methods while controlling the genomewide error rate.

Table 4 shows the total number of significant genes detected for all methods as well as the overlap between each. The least powerful method to detect differential expression was the test of while the most powerful method was the test of . Of our proposed methods, the statistic gave conclusions that overlapped most with the original methods in the paper, the likelihood ratio test (LRT) based on the maximum likelihood estimate [ 9]. In fact, found all of the differentially expressed genes found by the LRT. Tests conducted using the edgeR package gave the second largest overlap. Comparing the proposed methods with the Affymetrix data reported in the original paper [ 9], the most overlap in calls was with the tests, followed by the LRT using . If we look at the most recent methods developed for RNA-Seq data, the results of the DESeq package overlap most with , followed by the edgeR package, the LRT and then .

Table 4 Significant genes detected from the dataset in Marioni et al (2008).

Log ₂-ratio estimates produced much missing data in the RNA-Seq data analysis, with 8,947 sequences eliminated from analysis. Of these missing sequences, the proposed methods detected 4,171 ( ) or 2,406 ( ) significant calls (see Table 5). This means that using a test based on the log ₂-ratio may miss many possibly important differentially expressed genes. When using a LRT as suggested in Marioni et al (2008), 3,979 sequences were omitted from analysis because of missing data. The edgeR and DESeq packages did not generate any missing data. Of the 3,979 missing values generated by the LRT, our proposed methods detected 3,064 ( ) or 2,208 ( ) additional differentially expressed genes while the edgeR and DESeq packages detected an additional 235 and 197 genes. Since the true differential expression in this data is unknown, the differences between these methods are intriguing, but it is not clear whether one method is more accurate than another in this analysis. Overall, these findings suggest that test statistics based on a proportion statistic do not result in missing data, and more importantly, can detect possibly important differentially expressed genes that the log ₂-ratio based methods would miss.

Table 5 A summary of the missing values for each of the tests and the number of significant genes detected by other methods within those missing values.

Discussion

Although log ₂-ratios are widely used to compare two groups of expression data, there are limitations to using these statistics. The largest drawback to ratio statistics is that they are unstable as the denominator gets closer to zero. In addition, frequentist methods for constructing a corresponding variance and formally testing hypotheses of differential expression are unsatisfying and more complicated than typical scenarios [ 15, 16]. Due to these drawbacks, we proposed an alternative to testing for differential expression with all of the advantages of a log ₂-ratio statistic and none of its disadvantages.

We examined the proposed alternative, a proportion statistic, in four sets of simulations and two different sets of expression data. In simulations, the statistic , plus a constant and limma/EBA were robust to changes in distributional assumptions and the others were not. For the case of the Poisson distribution with rate parameter λ = 3, the statistic was underpowered, but otherwise and performed similarly well in simulations. The simulations suggest the use of in differential expression analyses because it uniformly preserved type I error and had competitive power. Note, however, that is not uniformly most powerful, and statistics derived from specific distributions can beat it when the distributional assumptions hold. The performance of the empirical Bayes analysis was competitive with in simulations but not in the analysis of the RNA-Seq dataset. Future research of interest may extend the test within an empirical Bayes framework, akin to what already exists for the log ₂-ratio. This may even further extend the clearly demonstrated feasibility of to detect differentially expressed genes.

Additionally, while the popular statistic performs sufficiently well under many simulation conditions, it suffers from problems with missing data in real data analysis problems when expression values are low. The addition of constants 0.05 and 0.5 appreciably improve simulation results and make the performance of nearly as good as (see Additional file 1). Nevertheless, this ad hoc procedure can be avoided using . In both the analysis of a cDNA microarray set and an RNA-Seq dataset, the log ₂-ratio based statistics led to missing values. Of these genes with missing ratio values, the proportion statistic was able to detect instances of statistically significant differential expression. We therefore recommend for general use over the other statistics discussed.

Methods

This section describes the data generation process in our simulations and the data collection in the two datasets analyzed in this paper.

Simulations

The proposed test statistics were evaluated under four different distributions. Though sophisticated simulations can be used to mimic expression data, the simulations below use simple scenarios so as to examine the performance of test statistics under basic distributions and to compare the eight different methods clearly and meaningfully. The first set of simulated intensity values were sampled from an exponential distribution that mimics the values from a 16-bit TIFF image of a cDNA microarray with respect to center and spread. The reference sample was taken from an Exp(1/4000) and the test sample was taken from a c × Exp(1/4000) where c was the fold-change value, c = 1,2,3,4,5. The four statistics, , and , were calculated for each value of c and sample sizes n = 3, 5, 10, 15, 20, 25, 30, 40, 50. Additionally we evaluated the ratio statistic after shifting values for an arbitrarily small constant set at either 0.05 and 0.5. For further comparison, a standard implementation of the limma/empirical Bayes method of Smyth (2004) was perfomed [ 17]. For simulations of count data values, we evaluated methods that account for overdispersion in the tests of differential expression using the edgeR and DESeq packages in Bioconductor [ 18, 19]. The implementation in both packages fixes a constant library size for each sample so that normalization is not executed. Sample sizes larger than n = 50 give simulation results similar to those for sample sizes of 50. In order to compare results, the p-value for an independent t-test was computed, with a null hypothesis of no difference between the two sample means. The null hypothesis was rejected if the p-value was below α = 0.05 and the proportion of rejections out of 1000 simulations was recorded (Table 1 and Additional file 1). The null hypothesis of no differential expression is equivalent to a fold change of one, c = 1. When the fold change is greater than one, we are calculating the power to detect differential expression. In this way, type I error and power were compared across the different methods. The results would be equivalent when using reciprocal fold changes instead. The simulations for an exponential distribution were repeated for an Exp(1/400) distribution, to study the effects of changing the scale.

A second set of simulated sampled intensity values from a Binomial( M = 10000, p = 0.5) distribution were obtained. The choice of this distribution was motivated by the derivation behind the maximum likelihood estimate of the proportion . The size was chosen to mimic the values from a 16-bit TIFF image of a cDNA microarray with respect to center. Analogous simulations to the exponential above were conducted with respect to statistics, sample sizes, and fold changes. For the binomial distribution, fold-changes of 2, 3, 4, and 5 correspond to binomial probabilities of 2/3, 3/4, 4/5, and 5/6 respectively. The simulations were repeated for a Binomial( M = 100, p = 0.5) distribution, to study the effects of change in the size parameter. A third set of simulated sampled intensity values from a Poisson(λ = 3) distribution were obtained. This distribution is motivated by the derivation behind the likelihood ratio test used in Marioni et al (2008) [ 9]. The parameter λ = 3 was chosen to mimic the number of categories arising from a smooth histogram of values from the RNA-Seq data. The simulations were repeated for a Poisson(λ = 30) distribution, to study the effects of change in the rate parameter. Analogous simulations to the exponential above were conducted with respect to statistics, sample sizes, and fold changes.

A fourth set of simulated sampled intensity values from a Normal( μ = 5, σ = 1) distribution were obtained. This distribution was included since many analyses assume expression data to be normally distributed. The center was chosen to mimic values from cDNA data with mean 5,000 and standard deviation 1,000, scaled to Normal( μ = 5, σ = 1). Analogous simulations to the exponential above were conducted with respect to statistics, sample sizes, and fold changes. For the normal distribution, fold-changes of 2, 3, 4, and 5 correspond to test samples of Normal( c × μ, σ = 1). The simulations were repeated for a Normal( μ = 10, σ = 2) distribution, to study the impact of changing the parameters.

All simulations were conducted using R http://www.r-project.org and the code is available in Additional file 2.

Gene Expression Data from Mice using cDNA Microarrays

We examined our proposed approach in a well-known and often cited set of cDNA microarrays. We chose this set because many research groups have evaluated their methods on this data and consequently the differential expression behavior in this data are better understood. The Apo AI experiment used cDNA microarrays to measure gene expression in the livers of 8 inbred control mice versus 8 mice with the Apo AI gene "knocked out" [ 20]. For each microarray, the reference sample was created from the pooled cDNA of the eight control mice. The goal of the experiment was to detect differential expression in the liver between control mice and the genetic knockout strain [ 12]. Since the Apo AI gene plays a role in HDL metabolism, differentially expressed genes are likely associated with lipid metabolism. The data can be obtained as an Rdata object from http://www.bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/Data/apoai.zip on the Bioconductor website. Welch two-sample t-statistics for each of the 6,384 probes were calculated,

where X trt and X cont were either our proportion estimators, and or the usual log ₂-ratio estimators, and . Since the variability of the cDNA data resembles the exponential distribution, the assumptions for methods and do not hold and therefore they were not used. To account for multiple testing, the original analysis used the maxT step-down procedure based on the t-statistics and found eight significantly differentially expressed probe sequences [ 12]. In order to explore the performance of alternative methods with both of the test statistics, the limma/EBA method of Smyth (2004) was computed [ 17]. Although this method was developed for log ₂-ratio values, we used the same programs on the proportion values as well.

Gene Expression Data from Human Kidney and Liver Cells using RNA-Seq

In order to examine the performance of our new approach on a sequence-based technology, we analyzed a set of RNA-Seq data discussed in Marioni et al (2008) [ 9]. This set of data compared the expression of 32,000 sequences in human kidney and liver cells extracted from the same person. The expression was also measured using Affymetrix U133 oligonucleotide arrays. Data was obtained from Supplemental Table 2 in the original manuscript. To compare our methods with those reported in the Supplemental Table 3 of their manuscript, we extracted the same five lanes of Illumina sequencing data corresponding to 3 pM concentrations of cDNA. We calculated both of the proportion tests outlined in the Results section, the ratio-based test provided in the Background section, and compared them to the methods from the original paper and more recent methods that account for overdispersion [ 18, 19].

The methods to test for differential expression from RNA-Seq data in the original paper used a likelihood ratio test (LRT) for inference [ 9]. Their test assumes that the expression data follows a Poisson distribution where the rate of expression λ is equivalent in kidney (K) and liver (L) cells under the null hypothesis. For gene j, the likelihood ratio test compares H ₀ : λ Kj = λ Lj versus the alternative hypothesis that expression rates differ H ₁ : λ Kj ≠ λ Lj. The likelihood ratio test is

for gene j. The maximum likelihood estimate for the alternative hypothesis from the above LRT is denoted by . The original paper also tested for differential expression on the Affymetrix platform for the same tissue samples. The methods employed were an empirical Bayes analysis with a false discovery rate of 0.1% [ 17]. More recent developments that account for overdispersion in the tests of differential expression were implemented using the edgeR and DESeq packages in Bioconductor [ 18, 19].

Appendix: Bayesian Estimation and Inference

In order to compare our proposed methods to previously suggested test statistics in the data analysis sections, we evaluated the proportion statistics within a frequentist testing framework. It is also possible to conceive the model in a Bayesian framework. Given the binomial assumption presented in the Results section, a Bayesian analysis can be conducted. Let the beta distribution be denoted by β( a, b), with density

Where Γ ( a) is the gamma function with parameter a. We denote the Bayesian estimator of p j by . Using a Beta prior for p j with parameters a and b, the posterior distribution of , is β ( y + a, M - y + b), where and with density

[ 21]. To compare the performance of the Bayesian with frequentist statistics , and , credible intervals and confidence intervals can be constructed and coverage can be examined in simulations. For data where the difference of two proportions is required, the posterior distribution derived in [ 22] can be used.

References

Download references

Acknowledgements

The authors would like to thank Suzanne Grindle in the Department of Microbiology at the University of Minnesota for her feedback, as well as that of the anonymous reviewers, which has resulted in improvements to the manuscript. Support for this research was provided by the Institute for Pure and Applied Mathematics at UCLA.

Author information

Authors and Affiliations

1. Division of Biostatistics, School of Public Health, University of Minnesota, Minneapolis, MN, 55455, USA

Tracy L Bergemann

2. Cardiac Rhythm Disease Management, Medtronic, Mounds View, MN, 55112, USA

Tracy L Bergemann

3. Department of Mathematics and Computer Science, Biola University, La Mirada, CA, 90639, USA

Jason Wilson

Corresponding authors

Correspondence to Tracy L Bergemann or Jason Wilson.

Additional information

Authors' contributions

TLB and JW contributed equally to the research. Both authors read and approved the final manuscript.

Electronic supplementary material

Authors' original submitted files for images

Below are the links to the authors' original submitted files for images.

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Reprints and permissions

About this article

Cite this article

Bergemann, T.L., Wilson, J. Proportion statistics to detect differentially expressed genes: a comparison with log-ratio statistics. BMC Bioinformatics 12, 228 (2011). https://doi.org/10.1186/1471-2105-12-228
Download citation
Keywords

# data[].attributes.attachment_link.title
Proportion statistics to detect differentially expressed genes: a comparison with log-ratio statistics - BMC Bioinformatics

# data[].attributes.attachment_link.description
Background In genetic transcription research, gene expression is typically reported in a test sample relative to a reference sample. Laboratory assays that measure gene expression levels, from Q-RT-PCR to microarrays to RNA-Seq experiments, will compare two samples to the same genetic sequence of interest.

# data[].attributes.attachment_link.title
Bienvenidos al mundo de las estadisticas | Walter Sosa Escudero | TEDxRiodelaPlata

# data[].attributes.attachment_link.title
Our strategy to achieve and document reproducible computing - BMC Bioinformatics

# data[].attributes.attachment_link.description
BMC Bioinformatics volume 14, Article number: A19 (2013) Cite this article

# data[].attributes.attachment_link.readable_text
Proceedings

Open access

Published:

BMC Bioinformatics 13, Article number: S5 (2012) Cite this article

3704 Accesses

8 Citations

Metrics

Abstract

Background

The T box riboswitch controls bacterial transcription by structurally responding to tRNA aminoacylation charging ratios. Knowledge of the thermodynamic stability difference between two competing structural elements within the riboswitch, the terminator and the antiterminator, is critical for effective T box-targeted drug discovery.

Methods

The ΔG of aminoacyl tRNA synthetase (aaRS) T box riboswitch terminators and antiterminators was predicted using DINAMelt and the resulting ΔΔG (ΔG Terminator - ΔG Antiterminator) values were compared.

Results

Average ΔΔG values did not differ significantly between the bacterial species analyzed, but there were significant differences based on the type of aaRS.

Conclusions

The data indicate that, of the bacteria studied, there is little potential for drug targeting based on overall bacteria-specific thermodynamic differences of the T box antiterminator vs. terminator stability, but that aaRS-specific thermodynamic differences could possibly be exploited for designing drug specificity.

Background

The T box riboswitch is an important regulatory mechanism found in Gram-positive bacteria including many pathogens [ 1- 3]. The riboswitch responds to high levels of uncharged (non-aminoacylated) tRNA to control the transcription of cognate genes (e.g., aminoacyl tRNA synthetase, aaRS, genes) [ 3]. Cognate, uncharged tRNA binds the 5'-untranslated region of T box mRNA during transcription and, when present in sufficient quantities, results in antitermination (Figure 1a) [ 3]. The tRNA anticodon binds a specifier sequence in Stem 1, thus providing the cognate specificity, while the uncharged tRNA acceptor end nucleotides bind the first four bases in a seven nucleotide bulge of the highly conserved antiterminator structural element [ 4]. Aminoacylated-tRNA also binds the specifier sequence, but does not bind the antiterminator element [ 5]. In the absence of uncharged tRNA bound to the antiterminator element, a more thermodynamically stable stem-loop structure forms (the terminator) followed by factor-independent transcription termination a few nucleotides later [ 3]. The antiterminator and the terminator are mutually exclusive structural elements due to sharing common nucleotides (Figure 1b). Ligands that target and disrupt the T box riboswitch function could be potential antibacterial agents in light of the critical genes regulated by the T box mechanism [ 1, 2].

T box transcription antitermination mechanism. a) Uncharged, cognate tRNA binds the 5' -untranslated region of the nascent mRNA and stabilizes the antiterminator. b) In the absence of uncharged, cognate tRNA, the more stable terminator forms and transcription is terminated before the translation start site.

We have been investigating the structure-function relationship of the T box antiterminator element and the key recognition features necessary for ligands to specifically bind the antiterminator and disrupt its function. There are few detailed medicinal chemistry studies of ligands targeting RNA [ 6] and this project has required an extensive multidisciplinary approach. The solution structure of antiterminator model RNA AM1A [ 7] was determined using molecular modelling with NMR-derived constraints [ 8]. The structure indicated that the seven-nucleotide bulge of the antiterminator was not pre-ordered for tRNA binding, but rather, that binding of the tRNA acceptor end must require a certain extent of tertiary-structure capture and/or an induced fit. Fluorescence life-time studies confirmed that modest antiterminator structural reorganization occurs upon tRNA binding in a magnesium-dependent manner [ 9]. While some RNAs have specific divalent metal ion binding sites, for the antiterminator RNA, the Mg ²⁺ binds via a diffuse, outer-sphere interaction [ 10]. In vitro selection studies of both the antiterminator [ 11] and tRNA [ 12] indicate that there are likely no direct interactions between the tRNA and the antiterminator other than the known base pairing between the acceptor end nucleotides and the first four nucleotides at the 5'-end of the seven-nucleotide bulge. Given this antiterminator structure-function information, ligands could potentially disrupt tRNA binding simply by competing with the base pairing between the tRNA acceptor end and the antiterminator bulge nucleotides.

Aminoglycosides bind AM1A in a structure-specific manner, most likely via electrostatics [ 13, 14]. In contrast, two novel classes of heterocyclic compounds, triazoles and oxazolidinones, bind AM1A with structure-specificity and high affinity, but without reliance on electrostatics [ 15- 18] and can alter T box transcription antitermination [ 16]. A computational, quantitative structure activity relationship analysis has shown that hydrophobic interactions play a significant role in the binding of these compounds to AM1A [ 19]. Furthermore, the small molecule ligands disrupt tRNA binding to the antiterminator in a structure-specific manner [ 20].

From a drug discovery perspective, a key factor to determine is the range of ligand-induced stabilization that can be accommodated without overly stabilizing the antiterminator element and precluding terminator formation. The goal of our T box drug discovery project is to determine the key ligand features that lead to specific antiterminator binding, but that do not result in excessive stabilization of the antiterminator secondary structure. These ligands could then potentially compete with tRNA for binding to the antiterminator, but still allow terminator formation such that transcription of a T box gene critical for bacterial survival would be terminated (Figure 2). The predicted thermodynamic stability (ΔG) of the terminator and antiterminator structural elements have been compared for the B. subtilis tyrS T box [ 4]. However, there has been no systematic comparison of predicted thermodynamic stability differences for a larger set of T box genes. Using a structural bioinformatics approach, we have analyzed the differences in predicted free energy (ΔΔG) between antiterminators and terminators in a set of aaRS T box genes in order to predict an upper limit of ligand-induced stabilization that can potentially be accommodated.

T box drug discovery goal. As soon as the antiterminator is transcribed (a), a small molecule ligand competes with tRNA for binding the antiterminator (b), but does not overly stabilize the antiterminator such that the terminator can form and transcription is terminated (c).

Methods

The thermodynamic stability of T box antiterminator and terminator structural elements was calculated using the DINAMelt server [ 21]. The DINAMelt server computes the secondary structure and free energy of the folded RNA using a secondary structure folding algorithm and empirically-derived nearest neighbour coefficients [ 21]. The folding algorithm predicts the minimum energy RNA secondary structure using the available thermodynamic data for base pairing, base stacking and destabilizing energies [ 22, 23]. The sequences analyzed were predicted to be involved in aaRS T box antiterminator and terminator structural elements from Bacillus cereus (BC), Bacillus subtilis (BS), Clostridium botulinum (CB), Clostridium difficile (CDF), Clostridium perfringens (CPE), Staphylococcus aureus Mu50 (SA), Streptococcus agalactiae (SAG), and, Streptococcus pyogenes (SPY) [ 2]. The ΔG for each aaRS T box antiterminator and terminator sequence was determined using the Quickfold (RNA 3.0) option on the DINAMelt server and the predicted thermodynamic stability difference calculated from ΔΔG = ΔG Terminator - ΔG Antiterminator. The % suboptimal setting was adjusted as necessary to obtain the lowest free energy antiterminator fold that had the consensus secondary structure of the core seven-nucleotide bulge containing the 5' -UGGN-3' nucleotides that are complementary to the tRNA acceptor end nucleotides [ 4].

Results and discussion

The free energy values for the antiterminator and terminator structural elements of the T box genes analyzed were calculated using the DINAMelt webserver [ 21]. The predicted ΔG values are listed in Additional File 1 and the relative comparisons of terminator vs.antiterminator stability (ΔΔG) are summarized in Figure 3. The overall average ΔΔG for all aaRS studied was -12.8 kcal/mol. The average ΔΔG values did not differ significantly between bacteria when comparing between the pathogenic bacteria nor between pathogenic vs. the non-pathogenic bacteria studied (BS) (Figure 3a). In contrast, there were significant differences in ΔΔG averages between specific aaRS with alanyl aaRS having the smallest average ΔG (-7.8 ± 3.5 kcal/mol) and glycyl aaRS having the largest (-20.1 ± 4.6 kcal/mol) (Figure 3b). Based on these results, the glycyl aaRS, on average, may be best able to accommodate ligand-induced stabilization of the antiterminator and still allow formation of the terminator. An important factor to consider, however, is that the free energy calculations are based on empirically-derived parameters for known RNA structural motifs [ 21]. Structural motifs, especially in loops and bulges, that have not been previously characterized might contribute to the stability of the RNA elements and not be accounted for in the DINAMelt ΔG calculations. Since the loop of the antiterminator is not highly conserved [ 4], most likely there is no structural motif within it that might contribute to the antiterminator stability, however, the possibility cannot be excluded. Further investigation of experimentally-derived free energy values of individual antiterminators and terminators is needed.

Average calculated ΔΔG. Average calculated aaRS ΔΔG values grouped by a) bacteria and b) aaRS. Standard deviation indicated by bar and number of aaRS sequences averaged noted over the individual columns.

Conclusions

The free energy of T box riboswitch antiterminator and terminator elements was predicted and compared for a series of aaRS T box genes. The observed aaRS-specific stability differences between these key riboswitch structural elements could potentially be targeted to effect ligand-specificity in future drug discovery efforts.

Abbreviations

(aaRS):

Aminoacyl tRNA synthetase

(BC):

Bacillus cereus

(BS):

Bacillus subtilis

(CB):

Clostridium botulinum

(CDF):

Clostridium difficile

(CPE):

Clostridium perfringens

Mu50 (SA):

Staphylococcus aureus

(SAG):

Streptococcus agalactiae and, Streptococcus pyogenes (SPY).

References

Download references

Acknowledgements

We wish to thank the National Institutes of Health (GM073188) for support of this work and the office of the VP for Research, Ohio University and the DAAD for support of FJ on an Internationale Studien- und Ausbildungspartnerschaften (ISAP) grant.

This article has been published as part of BMC Bioinformatics Volume 13 Supplement 2, 2012: Proceedings from the Great Lakes Bioinformatics Conference 2011. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/13/S2

Author information

Authors and Affiliations

1. Department of Chemistry, Leipzig University, Leipzig, Germany

Franziska Jentzsch

2. Department of Chemistry & Biochemistry, Ohio University, Athens, OH, USA

Franziska Jentzsch & Jennifer V Hines

Corresponding author

Correspondence to Jennifer V Hines.

Additional information

Competing interests

JH is a co-inventor on U.S. Patent number 7,005,441 for T box binding assays and ligands.

Authors' contributions

FJ carried out the sequence folding and free energy predictions. JH conceived of the study, designed and coordinated the study and prepared the manuscript.

Electronic supplementary material

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Reprints and permissions

About this article

Cite this article

Jentzsch, F., Hines, J.V. Interfacing medicinal chemistry with structural bioinformatics: implications for T box riboswitch RNA drug discovery. BMC Bioinformatics 13 (Suppl 2), S5 (2012). https://doi.org/10.1186/1471-2105-13-S2-S5
Download citation
Keywords

# data[].attributes.attachment_link.title
Interfacing medicinal chemistry with structural bioinformatics: implications for T box riboswitch RNA drug discovery - BMC Bioinformatics

# data[].attributes.attachment_link.description
Background The T box riboswitch controls bacterial transcription by structurally responding to tRNA aminoacylation charging ratios. Knowledge of the thermodynamic stability difference between two competing structural elements within the riboswitch, the terminator and the antiterminator, is critical for effective T box-targeted drug discovery.

# data[].attributes.attachment_link.title
COMO APRENDER COMPUTACION [ GRATIS ] 😱😱 DESDE CERO 2024  | CLASE #01

# data[].attributes.attachment_link.readable_text
The functional properties of a protein are strongly dependent on its structural conformation. The primary question addressed in this work is how to determine structural and therefore functional similarity from 3D protein structures. The approach we take relies on protein structure alignment, which elucidates functional protein relationships that are not depicted by the sequence.

Most current structural alignment tools are based on geometric properties of protein structures, and do not incorporate the effect of physical properties responsible for protein shape and similarity. The effect of physical features, such as hydrogen bonding, bond and torsion angles, hydropathy, and solvent accessibility, are not well understood and hence not used while developing structure alignment tools. In this work, we exploit physical properties in protein structural alignment algorithms. In particular, we incorporate dominant physical features in one of the exact protein structure alignment tools, CMOS [ 1], which is based on the contact map overlap maximization (MAX-CMO) formulation of protein structure alignment. CMOS is a state-of-the-art tool, providing fast structural alignments that are in excellent agreement with the SCOP fold family classification. An implementation of CMOS is accessible at http://eudoxus.cheme.cmu.edu/cmos/cmos.html.

# data[].attributes.attachment_link.title
Exploiting physical properties in protein structure alignment - BMC Bioinformatics

# data[].attributes.attachment_link.description
The functional properties of a protein are strongly dependent on its structural conformation. The primary question addressed in this work is how to determine structural and therefore functional similarity from 3D protein structures. The approach we take relies on protein structure alignment, which elucidates functional protein relationships that are not depicted by the sequence.

# data[].attributes.attachment_link.title
La física cuántica explicada en menos de 5 minutos | BBC Mundo

# data[].attributes.attachment_link.title
Pero... ¿PARA QUÉ SIRVEN LAS MATEMÁTICAS?