Correction of gene model annotations improves isoform abundance estimates: the example of ketohexokinase (Khk)

Khk isoforms expression is tissue-specific

In order to assess tissue-specific expression patterns of Khk in mouse, we downloaded RNA-seq data generated from 14 different mouse tissues³³. The availability of 4 biological replicates (2 males and 2 females) per tissue enabled us to conduct a differential gene expression analysis using standard, gene-count based analytical workflows. As previously described in gene specific studies³⁵, we identified a strong tissue specific expression of Khk (adjusted p-value of 0 from DESeq2), with significantly higher expression levels in the liver, small intestine and kidney when compared to other tissues (Figure 1B Underlying data: Table 1). Interestingly, gender did not impact the overall Khk expression levels. These results were all confirmed using the gene level estimates based on Salmon quantifications (Extended Data Figure 1A and B). In particular, patterns of Khk expression across tissues were highly concordant between count-based and Salmon estimates. We also sought to evaluate changes in relative exon usage for this gene using JunctionSeq¹⁴, including both exon and junction counts in our analysis. Therefore, we selected five tissues (liver, spleen, lung, heart, kidney) with known variations in Khk isoform usage^30–32,35 and different levels of overall expression. We clearly identified a preferential inclusion of exon 3 C (Ensembl ID ENSMUSE00000186455) in liver and kidney, as previously described, while exon 3A (Ensembl ID ENSMUSE00001361691) is preferentially retained in heart, spleen and lung (p value < 0.001) (Figure 1C). This trend was also reflected in junctions spanning these exons (Underlying data: Table 2). Variations in 5’ UTR were also observed, most likely resulting from alternative transcription start site usage in liver and kidney. In summary, using traditional count-based methods, we confirmed the tissue-specificity of Khk expression and identified exons preferentially retained in some tissues.

An annotated, unexpressed retained intron biases Khk isoform quantification

As count-based methods might not always reflect the full complexity of splicing patterns and isoform diversity, we sought to quantify relative proportions of Khk isoforms in each tissue to fully capture the complexity of its expression patterns. The GENCODE annotation together with the Salmon¹⁵ quantifier were used to obtain these estimations (Underlying data: Table 3). Interestingly, the quantification results showed that the annotated retained intron (Khk.RI) accounts for more than 15% of expressed isoforms in 8 tissues (Figure 2A). This trend was consistently observed in DRIMSeq²³ proportion estimates (Figure 2A; Underlying data: Table 3) and the raw TPM (transcripts per million) values output by Salmon (Extended data: Figure 2). Nevertheless, examination of the coverage tracks generated from the alignment of the reads to the reference genome showed hardly any detectable expression of the transcript (Figure 2B; Extended data: Figure 3). This observation was also supported by the differential exon usage analysis which clearly revealed very low expression levels for the only Khk.RI-specific exonic region, E012 (Figure 1C; Extended data: Figure 4A). To experimentally confirm the absence of Khk.RI expression, we designed PCR primers amplifying fragments specific to this transcript (see Extended data: Figure 4B; Extended data: Table 1 for a list of primers used in the study) and used RT-qPCR and semiquantitative RT-PCR to measure its expression levels. Khk.RI could hardly be detected using this sensitive method (Extended data: Figure 4C and 4D) and comparison with the expression levels in a Khk-A/C^-/- mouse model showed that it is hardly expressed in heart, kidney, liver, lung and spleen, whereas a product could be amplified using genomic DNA as template (Figure 2C). This experimental validation further demonstrates that the contribution of Khk.RI to the overall Khk expression level was overestimated during the quantification step. Taken together, these results suggest that the presence of non-expressed transcripts in the “raw” gene annotation may result in erroneous detection by a transcript quantification software.

Figure 2. Khk.RI expression levels are overestimated using Salmon quantifications.

A – Estimated relative Khk isoform expression across all tissues in mouse. Proportions were output by DRIMseq using Salmon quantification as an input. B – RNA-seq coverage tracks and sashimi plots highlighting the absence of Khk.RI expression in thymus, spleen and bone marrow samples. Data obtained from all biological replicates were merged prior to plotting. C – Products derived from semiquantitative RT-PCR analysis on cDNAs prepared from total RNA of different mouse organs using the primers specific for Khk.RI. Genomic DNA isolated from the liver was used as positive and RNA isolated from organs of Khk-A/C knockout mice as negative control. (n = 5 and products from two representative mice are shown).

Manual update of the Khk transcript annotations improved quantification results

Since the current genomic annotation did not reflect Khk isoform expression and introduced biases in quantifications, we manually removed the Khk.RI transcript (Ensembl ID ENSMUST00000200978.1) prior to the quantification step. Despite this adjustment, inspection of the junction reads, coverage tracks, and normalised exon and junction counts derived from QoRTs (Figure 1C) revealed discrepancies between quantification estimates and results derived from alignment-based methods (Figure 3A and 3B; Extended Data: Figure 5). These quantification estimates also did not reflect the observations made by previous reports focusing on the characterisation of Khk expression patterns³⁵. An example of such discrepancy may be found examining the heart coverage tracks: while it is quite clear that the KhkC-specific exon 3C (Ensembl ID ENSMUSE00000186455) is hardly captured in comparison with exon 3A (Ensembl ID ENSMUSE00001361691), quantification estimates yielded a score of 39.9% and 16.7% for KhkC and KhkA, respectively (Underlying data: Table 3). Similarly, KhkC estimates were inflated in lung and spleen (14.9% and 18.6% to be compared with the absence of coverage of exon 3C), as were Khk.Skip estimates (32%) in liver (Underlying data: Table 3). Since such quantifications are relying on the reference transcripts provided during the indexing step¹⁵, we reasoned that incorrect transcript models might be the cause of the observed discrepancies and therefore compared coverage tracks, normalised exon counts and isoform models. We identified annotated differences in 3’ end annotations between all isoforms which were not reflected on our coverage tracks (Figure 1C and Figure 3B). As Khk isoforms can be identified unambiguously based on the exclusion patterns of the exons 3A and 3C and regardless of differences in UTRs, we could investigate the impact of these UTR variations on transcript quantifications. We therefore manually updated Khk isoform annotations to provide an identical 3’ end to all isoforms (Underlying data: File 1 and 2) and re-estimated isoform proportions (Figure 3A, second panel). Despite this adjustment, inspection of the proportion estimates still revealed erroneous estimations of isoform expression in particular in the case of liver where KhkA was detected in levels similar to KhkC. Further examination of the results revealed that, while some differences in 5’ end coverage in the dataset were concordant with the current gene annotation, they were not always reflected in the gene model (Figure 1C and Figure 3A, heart, lung and spleen 5’ UTR coverage). Following a similar approach to the one described earlier for 3’ UTRs, we finally manually modified Khk isoforms to provide an identical 5’ and 3’ end to all listed isoforms (Underlying data: File 3). Transcript quantification performed using this updated annotation yielded more concordant results when compared to coverage tracks and in light of previous reports, in particular for tissues such as liver³⁰, small intestine³⁶, heart³², spleen and lung (Figure 3A and 3B; Extended data: Figure 5). Interestingly, a substantial fraction of isoforms detected in kidney were still attributed to the Khk.Skip isoform while both junction count analysis and single gene studies reported a prevalence of KhkC³⁵. Additionally, we computed the relative Khk isoform usage using a new annotation with identical 5’ and 3’ end for all isoforms except Khk.RI which was retained as such in the gene model (Figure 3A). This modification was not sufficient to remove Khk.RI estimation biases, with the retained intron predicted to erroneously account for 20% of the overall gene expression in bone marrow or spleen. We therefore confirmed the impact and importance of both Khk.RI and UTRs annotations on Khk isoform expression estimates.

Figure 3. Khk isoform quantification estimates are improved after manual adjustment of the genomic annotation.

A - Estimated relative Khk isoform expression across all tissues in mouse, using four different genomic annotations: a raw annotation after removal of the Khk.RI transcript, an annotation without Khk.RI and with adjusted 3’ UTR, an annotation without Khk.RI and adjusted 3’ and 5’ UTR, and an annotation with adjusted 3’ and 5’ UTR for all transcripts except Khk.RI. The choice of annotation greatly impacts the quantification results. B – RNA-seq coverage track and sashimi plots illustrating the predominance of different isoforms in spleen, lung, liver, heart, and kidney. C – RT-qPCR analysis of total Khk expression in different mouse organs. Values are expressed as fold-change compared to the expression levels obtained for the heart, which was arbitrarily defined as 1. β-actin was used as the invariant reference gene. Data are mean ± SD (n = 5). D - RT-qPCR analysis of total Khk, KhkA + Khk.A.C, and KhkC + KhkA.C expression in different mouse organs. Values are expressed as fold-change compared to the expression levels obtained for total Khk, which was arbitrarily defined as 1. Data are mean ± SD (n = 5).

To confirm the biological relevance of our newly estimated proportions, we designed primer pairs to specifically target Khk isoforms (see Extended data: Figure 6A and 6B; Extended data: Table 1) and evaluated their relative expression using RT-qPCR. The expression of total Khk was ~60-fold higher in liver and kidney compared to heart, lung, and spleen (Figure 3C), confirming previously described findings³⁵. While we did not manage to achieve a reliable quantitative evaluation of Khk.Skip and KhkA.C expression alone, we accurately measured the expression of (KhkC + KhkA.C) and (KhkA + KhkA.C) in five mouse tissues and normalized it to the total Khk expression (Figure 3C and 3D). We thereby confirmed a strong prevalence of (KhkA + KhkA.C) expression in heart while (KhkC + KhkA.C) accounted for most of the expression measured in kidney and liver (Figure 3D). Therefore, we concluded that KhkA.C levels of expression were much lower than KhkA and KhkC in these three tissues and that the proportion estimates derived from our adjusted genomic annotation reflected the RT-qPCR measurements. The isoform measurements in lung and spleen highlighted low levels of KhkA, KhkC and KhkA.C compared to total Khk expression, strongly suggesting the prevalence of Khk.Skip expression, as observed on coverage tracks (Figure 3B) and in the proportion estimates derived from the updated genomic annotation (Figure 3A).

Finally, we evaluated whether the changes brought to the Khk transcript annotations affected the estimations of overall Khk expression levels. Gene level estimates obtained using quantifications based on 3 of the modified annotations were strongly correlated (Pearson, r > 0.99) with estimates derived from the raw GENCODE annotation (Extended data: Figure 7A). In addition, the high tissue specificity of Khk expression was observed in all cases, with identical expression patterns between annotations (Figure 1A, Extended data: Figure 7B).

Altogether, these findings demonstrate that UTRs and more generally 5’ and 3’ end annotation may greatly influence transcript quantification results. The resulting biases lead to the identification of isoforms that can hardly be detected in biological samples, therefore highlighting the importance of inspecting results for any given gene of interest.

Erroneous genomic annotation of Khk increases discrepancies between sequencing strategies and datasets

We next investigated whether such discrepancies could be observed when using different sequencing library strategies. To do so, we artificially created a single-end dataset by removing one read for each tissue sample. We also created a short-read dataset by trimming the remaining reads to only retain the first 50bp. In both cases, we aligned reads to the reference genome and independently quantified transcript expressions using Salmon as previously described. Regardless of the sequencing strategy considered, we observed a similar overestimation of the Khk.RI fraction as well as discrepancies between proportion estimates and junction counts (Extended data: Figure 8A and 8B). Both differences were corrected using the aforementioned manually curated annotations. We then compared proportion estimates between datasets for each annotation. Quite strikingly, we noted a much better agreement in transcript estimates using our manually modified annotation (Figure 4A–D). While the use of the “raw” annotation only resulted in a 0.56 correlation (Pearson) between estimates from the paired-end and 50bp single-end libraries, the use of an updated annotation resulted in a 0.97 correlation between both platforms. This trend was also observed between paired-end and single-end and between single-end and 50bp single-end respectively (Extended data: Figure 9A-F). To further evaluate the impact of annotations on estimate concordance across datasets, we downloaded RNA-seq data from another study profiling mRNA expression across 13 tissues³⁴. We quantified isoform usage for each tissue and compared those estimates with the ones from the 50bp single-end dataset from Li et al. 2017³³ in order to avoid biases due to differences in read length. In total, 8 tissues could be compared between both studies. Quite strikingly, the agreement between both datasets was not high (Pearson correlation 0.72) when considering fractions derived using the original annotation (Figure 4E). While the removal of Khk.RI without an adjustment of the UTR did further hinder reproducibility between both datasets (Figure 4F), we observed a much stronger consistency of proportion estimates after UTR adjustment (Figure 4H, Pearson correlation 0.73). However, the highest consistency between datasets was reached when using the fully updated annotation (Figure 4G, Pearson correlation 0.91). These results therefore further underscore the importance of appropriate annotations during transcript isoform quantifications. The use of erroneous gene models may further increase discrepancies between sequencing libraries and datasets as exemplified in the case of the Khk gene.

Figure 4. Manual adjustment of the Khk gene model improves the concordance between estimates across sequencing library strategies and datasets.

A – Comparison of relative Khk isoforms expression estimates between the full length paired-end dataset from Li et al. 2017³³ and the short read, single-end dataset. The estimates were generated using the naive Ensembl annotation. B – Comparison of relative Khk isoforms expression estimates between the full length paired-end dataset from Li et al., 2017³³ and the short read, single-end dataset. The estimates were generated using the modified Ensembl annotation where the Khk.RI transcript has been removed. C – Comparison of relative Khk isoforms expression estimates between the full length paired-end dataset from Li et al., 2017³³ and the short read, single-end dataset. The estimates were generated using the modified Ensembl annotation where the Khk.RI transcript has been removed and the 5’ and 3’ UTR of other transcripts were all adjusted. D - Comparison of relative Khk isoforms expression estimates between the full length paired-end dataset from Li et al., 2017³³ and the short read, single-end dataset. The estimates were generated using the modified Ensembl annotation where the 5’ and 3’ UTR of all transcripts except Khk.RI were all adjusted. E – Comparison of relative Khk isoforms expression estimates between the Li et al., 2017³³ dataset and the Söllner et al., 2017³⁴ dataset, using single-end, short (50bp) reads in each case. The estimates were generated using the naive Ensembl annotation. F – Comparison of relative Khk isoforms expression estimates between the Li et al., 2017³³ dataset and the Söllner et al. 2017³⁴ dataset, using single-end, short (50bp) reads in each case. The estimates were generated using the modified Ensembl annotation where the Khk.RI transcript has been removed. G – Comparison of relative Khk isoforms expression estimates between the Li et al., 2017³³ dataset and the Söllner et al. dataset, using single-end, short (50bp) reads in each case. The estimates were generated using the modified Ensembl annotation where the Khk.RI transcript has been removed and the 5’ and 3’ UTR of other transcripts were all adjusted. H - Comparison of relative Khk isoforms expression estimates between the Li et al., 2017³³ dataset and the Söllner et al. dataset, using single-end, short (50bp) reads in each case. The estimates were generated using the modified Ensembl annotation where the 5’ and 3’ UTR of all transcripts except Khk.RI were all adjusted.

Mice

C57BL/6J were obtained from The Jackson Laboratory, while KhkA/C^-/- mice, which are of C57BL/6 background and are lacking both ketohexokinase-A and ketohexokinase-C, were obtained from R. Johnson (University of Colorado) and used as negative control. All mice were housed in a pathogen-free facility at the ETH Phenomics Center (EPIC) under standard conditions (12 h light and 12 h dark cycle) with free access to food and water. 3 female and 2 male C57BL/6J mice and 2 male Khk-A/C^-/- mice were euthanized with CO₂ at the age of 6 weeks and heart, kidney, liver, lung, and spleen were subsequently removed and shock-frozen in liquid nitrogen. The mice did not suffer during the euthanasia with CO₂; the mice were placed into a chamber that contained room air and then CO₂ was gradually introduced with no more than 6 psi to displace at least 20% of the chamber volume per minute. All protocols for animal use and experiments were reviewed and approved by the Veterinary Office of Zurich (Switzerland).

Evaluation of Khk isoform expression in vivo

Total RNA from C57BL/6J and KhkA/C^-/- mice was prepared from frozen tissues with RNeasy Mini Kit (QIAGEN, Hilden, Germany) and treated with DNase I to remove traces of DNA. First-strand complementary DNA (cDNA) was synthesized with random hexamer primers using the High-Capacity cDNA Reverse Transcription Kit (Cat. No. 4368813; Applied Biosystems). Quantitative reverse transcription PCR (RT-qPCR) was performed on a Roche LightCycler 480 in duplicates using 10 ng cDNA and the 2x KAPA SYBR FAST qPCR Master Mix LC480 (Sigma). Thermal cycling was carried out with a 5 min denaturation step at 95 °C, followed by 45 three-step cycles: 10 sec at 95 °C, 10 sec at 60 °C, and 10 sec at 72 °C. Finally, melt curve analysis was carried out to confirm the specific amplification of a target gene and absence of primer dimers. Relative mRNA amount was calculated using the comparative threshold cycle (C_T) method. β -actin was used as the invariant reference gene. The PCR amplification efficiency of RT-qPCR primer sets was determined with serial dilutions of liver cDNAs and was similar for all primer sets. Semiquantitative RT-PCR was performed using Phusion High-Fidelity DNA Polymerase (New England Biolabs) and the following 3-step amplification protocol: 30 sec at 98 °C (denaturation), 30 cycles of 10 sec at 98 °C, 30 sec at 63 °C, and 40 sec at 72 °C, and a final elongation step for 5 min at 72 °C. PCR products were evaluated after gel-electrophoresis. Primer sequences are listed in Extended data: Table 1.

Modification of Khk transcript annotation

A fasta file containing all manually modified Khk transcripts was created. All exonic sequences used to modify the gene model were downloaded from the Ensembl website. A fasta file containing nucleotide sequences of all transcripts from the GENCODE M14 annotation was loaded into R using the Biostrings v 2.46.0 package. All original Khk transcripts were removed from the DNAStringSet object and the updated transcripts were then added. The resulting annotation was written to a fasta file to generate Salmon indexes (see following sections for more details).

RNA-seq data processing and alignment

Fastq files from Li et al., 2017³³ and Söllner et al., 2017³⁴ were downloaded from SRA using the sra-tools software v2.7.0⁴². As they only provided two replicates (instead of 4), we excluded the testis and ovary samples from Li et al., 2017³³. Additionally, due to very low library complexities, we also excluded the pancreatic samples from Söllner et al., 2017³⁴. Reads were aligned to the M14_GRCm38.p5 reference genome using STAR 2.4.2a⁴³ together with the GENCODE M14 annotation (STAR was run with default parameters). The search for novel junctions was allowed during the mapping step. Gene, exon and junction level read counts were generated using the QoRTs software v1.2.42⁴⁴ after excluding reads with multiple alignments (MAPK score less than 255). All workflows were orchestrated using Snakemake v 3.13.3⁴⁵.

Transcript quantifications

Salmon 0.9.1¹⁵ and StringTie 1.3.3b⁴⁰ were used for transcript quantifications. In the case of Salmon, indexes were built from each fasta files using the default quasi-mapping mode and a k-mer of length 31 as recommended in the software documentation. Transcript isoforms were quantified using the default VEBM algorithm. Library types were inferred by the software. Sequence and GC bias corrections were performed during each quantification (--seqBias and –gcBias options) and 100 bootstraps were run to estimate the variance of abundance estimates.

In order to quantify transcripts using StringTie, alignments to the reference genome were first re-generated as described in the “RNA-seq data processing and alignment” section with the additional STAR flag “--outSAMstrandField intronMotif”. Quantifications were then performed with and without the GENCODE annotation used as a guide (-G option in StringTie). All other parameters were set to default values. Newly assembled transcripts were then merged using the transcript merge mode (--merge) with default parameters.

Differential gene expression

Gene count tables were loaded into R (v 3.4.1) as a DESeq2⁴ object to conduct differential gene expression analysis. For the purpose of differential gene expression analysis, we only retained features labelled as “gene” in the GENCODE annotation and of type protein_coding, antisense, sense_intronic, 3prime_overlapping_ncRNA, sense_overlapping or non_coding in order to exclude transcription products requiring specific library preparations to be accurately measured. Genes with very low counts were excluded from downstream analysis: the threshold was set at 50 mapped reads across all samples, corresponding to less than one read per sample on average. Estimated size factors were used to correct for differences in library size. Following the standard DESeq2 workflow⁴⁶, changes in gene expression were modelled using a variable accounting for differences in tissue of origin for each sample. A Likelihood Ratio Test (“LRT” option in DESeq2) was performed to compare this model to a reduced model consisting of only an intercept. Results were extracted with the DESeq2 results function and multiple testing correction was performed using the Benjamini Hochberg procedure⁴⁷.

Differential exon and junction usage

Exon and junction count tables from QoRTs were loaded in R (v 3.4.1) using the JunctionSeq (v 1.8.0) package¹⁴. Changes in exon usage were modelled using the tissue of origin as a main variable. Size factors and dispersions were estimated using default parameters (options “byGenes” and “advanced” respectively, see the JunctionSeq vignette for more details). Dispersion function fits, test for differential usage and estimation of effect sizes were also run using default parameters. Final test results were extracted using the writeCompleteResults function. Feature-level p values were adjusted using the Benjamini Hochberg procedure⁴⁷. Only features with an adjusted p-value below 0.05 were retained.

Estimation of relative isoform proportions

Transcript isoform quantifications from Salmon were loaded into R (v 3.5.1) using the tximport package²¹. Relative isoform proportions and differential transcript usage were then modelled using a Dirichlet-multinomial model as described in Nowicka et al., 2016²³ and implemented in the DRIMSeq package (v 1.6.0). Model precisions were estimated by the dmPrecision function run with default parameters, except for the search grid ranges which were set to -15 and 15. For each sample, feature proportions were then computed using the dmFit function (default parameters).

Data visualization

Genomic annotations and coverage tracks were plotted using the Gviz package (v 1.22.3)⁴⁸. Data from all available biological replicates were pooled together prior to plotting each track. Results from differential exon and junction usage were visualized using JunctionSeq. Other plots were generated with the ggplot2 package.

Underlying data

Fastq files from Li et al., 2017³³ are available from: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA375882

Fastq files from Söllner et al., 2017³⁴ are available from: https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-6081/

Scripts used to generate the analysis presented in this paper are freely and publicly available on Github: https://github.com/chbtchris/Khk_quantifications (Archived scripts: http://doi.org/10.5281/zenodo.2583233⁴⁹).

All underlying data are available at: https://doi.org/10.17605/OSF.IO/NMKFA⁵⁰. Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication). Files available are as follows:

Extended data

All extended data are available at: https://doi.org/10.17605/OSF.IO/NMKFA⁵⁰. Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication). Files available are as follows: