In RNA-Seq, reads are aligned or mapped to a reference genome or transcriptome, and gene-level expression values are then calculated.
However, in highly similar gene families, genes with corresponding pseudogenes, or genes that overlap one another, the same read may correspond to multiple genes.
In Case Study No. 414, “RNA-Seq Is Not a Simple Upgrade from Microarrays in Gene Expression Analysis”, we explained this issue as follows.
In gene families that share highly homologous regions, short-read RNA-Seq reads may correspond to multiple genes. Depending on the analysis method, multi-mapping reads may be excluded, or they may be probabilistically assigned among multiple candidate genes to estimate expression levels.
For genes containing many such reads, the assignment of reads to individual genes depends on the analysis method and annotation used, making it difficult in some cases to accurately estimate gene-specific expression levels or expression differences. In other words, for groups of genes with high sequence similarity, RNA-Seq measurements may not fully reflect the actual differences in expression among individual genes.
To examine how this problem has been investigated, we asked AI about relevant studies and reviewed the papers it suggested.
How are ambiguous reads handled?
The handling of reads that correspond to multiple genes can be broadly classified as follows.
| Method | How the reads are handled | Points to consider |
|---|---|---|
| Exclude them | Reads corresponding to multiple genes are not counted | Genes with few uniquely assignable regions may appear to have lower expression |
| Count them for all candidates | Each read is added to every corresponding gene | The same read is counted multiple times, which can inflate expression estimates |
| Distribute them equally | Each read is divided equally among candidate genes | Reads are divided equally even when the true expression levels differ |
| Distribute them using a model | Reads are assigned probabilistically using unique reads, estimated expression levels, and other information | When little information is available to distinguish genes, the result depends on the model and annotation |
Different processing methods produce different results
In a study describing CoCo, the authors compared quantification results for nested genes obtained using HTSeq-count, featureCounts, RSEM, and other methods.
With standard settings that exclude reads assignable to multiple genes, genes with supporting reads were sometimes not detected or were underestimated. RSEM and featureCounts with modified settings improved detection, but this did not necessarily mean that reads were assigned to each gene in the correct proportions.
A study describing alevin, a quantification method for single-cell RNA-Seq, also showed that discarding multi-mapping reads systematically reduces counts for the genes affected.
Thus, excluding ambiguous reads may lead to a loss of expression signal, whereas distributing them equally may distort the true expression ratio.
Does Model-Based Assignment Produce the Correct Result?
Salmon and RSEM probabilistically assign ambiguous reads among candidate transcripts using information such as uniquely assignable reads, estimated expression levels, and fragment length.
When enough reads are available to distinguish genes, this approach may provide more reasonable estimates than simply excluding ambiguous reads or distributing them equally.
However, if there are almost no sequences that distinguish two genes, the true expression ratio cannot be determined from the RNA-Seq data alone. Even so, analysis software will still output expression values based on its model.
Even when ambiguous reads are distributed using a statistical model, different programs do not necessarily produce the same result. This is because they differ in how reads are matched to candidate transcripts, how expression levels are estimated, how biases are corrected, and how annotations are handled.
A simulation-based comparison of highly similar multigene families in the repeat-rich genome of Trypanosoma cruzi showed that read assignment and estimation error varied substantially depending on the quantification method used, and that accurate quantification became increasingly difficult as sequence similarity approached 100%.
When the RNA-Seq data do not contain enough information to distinguish genes, the reported expression values may differ substantially depending on the program and annotation used. Furthermore, none of these estimates can necessarily be assumed to accurately represent the true expression ratio between the individual genes.
This issue is even more important in single-cell RNA-Seq, where relatively few reads are obtained per cell and only limited regions near the 3′ or 5′ end are often measured.
When individual cells do not contain reads that uniquely distinguish genes, some analysis methods use information about which candidate genes are compatible with each read or UMI to infer the origin of ambiguous reads. In the alevin study, the authors investigated a method that uses reads that are ambiguous between genes rather than discarding them, thereby reducing the underestimation of gene-level counts.
Such inference may improve the estimated counts, but it does not mean that the gene was uniquely identified and directly measured in each individual cell. When highly similar genes are used as cell-type markers, it is therefore important to consider how their expression values were calculated and how ambiguous reads were assigned.
The effect is not limited to expression levels
A study of cross-mappability between genes showed that read misassignment caused by sequence similarity can produce apparent co-expression relationships and false-positive trans-eQTL signals.
The effects of multi-mapping may therefore extend beyond individual gene expression estimates to correlation analysis and network analysis.
How serious is the problem in conventional bulk RNA-Seq?
With 150 bp paired-end bulk RNA-Seq, more sequence information is available than with short single-end reads, making it easier to distinguish the origin of a read. However, ambiguity may still remain for genes with high sequence similarity, even with paired-end sequencing.
A read that corresponds to multiple transcripts does not necessarily create a major problem in gene-level analysis when those transcripts belong to the same gene.
The important case is when the same fragment corresponds to different genes. In particular, paralogs, genes with corresponding pseudogenes, and genes in overlapping regions may remain sensitive to the analysis method even at the gene level.
Could 3′ RNA-Seq or microarrays distinguish these genes more clearly?
Even when genes are highly similar across much of their sequence, the 3′ region may contain sequences unique to one of the genes.
In such cases, 3′ RNA-Seq, which measures sequences near the 3′ end, or a microarray with probes designed against unique regions may distinguish the genes more clearly than conventional RNA-Seq. The 3′ UTR can show greater sequence divergence and has therefore been used in microarray probe design.
However, this will not help if the 3′ regions are also highly similar. Microarrays can also show cross-hybridization when probes bind to similar sequences.
Ideally, genes that are difficult to quantify should be examined not only across the full transcript, but also at the 3′ end and at the actual microarray probe sequences, so that their distinguishability can be compared across measurement methods. However, within the scope of this investigation, we did not identify a study that comprehensively evaluated all of these factors.
Do not accept calculated values uncritically
The studies reviewed here show that, for highly similar genes, estimated expression levels can vary depending on how multi-mapping reads are handled.
Even when analysis software outputs a numerical value, that value does not necessarily accurately represent the true expression level of the individual gene.
When comparing genes with similar sequences, results should therefore be interpreted while considering the quantification method used, the handling of multi-mapping reads, and whether uniquely assignable sequences are available.
When an important candidate gene is identified by RNA-Seq, it can also be useful to inspect the BAM file using IGV or a similar tool and determine which exons or splice junctions the reads correspond to.
For a more detailed examination, STAR or RSEM can be used to generate a transcript-coordinate BAM file, allowing read positions to be viewed on mature mRNA sequences in which the exons have been joined together.
In addition, sequences from regions where reads appear to be concentrated, or where coverage differs substantially between sample groups, can be retrieved and searched using NCBI BLAST to determine whether highly similar regions are also present in other genes or pseudogenes.
For genes whose expression level or differential expression is central to the conclusions of a study, it is important to confirm that the regions covered by the reads contain sequences unique to that gene and, when necessary, to validate the result using qPCR or another method targeting a gene-specific region.
