RNA-seq Parent-of-Origin Classification with Machine Learning applied to Alignment Features

doi:10.21203/rs.3.rs-3214264/v1

Download PDF

Research Article

RNA-seq Parent-of-Origin Classification with Machine Learning applied to Alignment Features

https://doi.org/10.21203/rs.3.rs-3214264/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 12 Mar, 2024

Read the published version in BMC Bioinformatics →

You are reading this latest preprint version

Background

Parent-of-origin allele-specific gene expression (ASE) can be detected in interspecies hybrids by virtue of RNA sequence variants between the parental haplotypes. ASE is detectable by differential expression analysis (DEA) applied to the counts of RNA-seq read pairs aligned to parental references, but aligners do not always choose the correct parental reference.

Results

We used public data from four species pairs that are known to hybridize. For each pair, we obtained RNA-seq read pairs from both species and measured our ability to assign each read to its proper species by comparing reads to the transcriptome or genome references. We tested four software packages that assign each read pair to a reference position and found that they often favored the incorrect species reference. To address this problem, we introduce a post process that extracts alignment features and trains a random forest classifier to choose the better alignment. On each dataset tested, our machine-learning post-processor achieved higher accuracy than the aligner by itself at choosing the correct species per RNA-seq read pair.

Conclusions

For the parent-of-origin classification of RNA-seq, machine learning can improve the accuracy of alignment-based methods. This approach should be useful for enhancing ASE detection in interspecies hybrids. We believe this is the first application of machine learning to this problem domain.

Machine Learning

RNA-seq

Allele-specific Expression

Sequence Alignment

RNA sequencing is a ubiquitous tool in molecular biology (1, 2). The technology can detect gene transcription and quantify RNA abundance (2, 3). Paired, short reads, as generated by Illumina sequencing machines, constitute the vast majority of RNA-seq in public databases (1), and several software tools specialize in mapping short reads to references. Mappers such as Salmon (4) and Kallisto (5) use alignment-free algorithms, while mappers such as Bowtie2 (6, 7) and STAR (8, 9) are alignment-based. A comparison of these two approaches found that alignment-free methods are faster, and are equally accurate on simulated data, but fail to achieve the accuracy of alignment-based methods for quantification of human RNA-seq (10). Various short-read aligners have been tested and compared for many applications (11–17). Aligners such as HiSat2 (18, 19) and STAR (8, 9) specialize in mapping RNA-seq to genome sequences, a task that requires virtual splicing of read sequences into exons and generation of alignments that skip over introns in the reference.

Differential expression analysis (DEA) is the statistical and comparative analysis of RNA-seq quantities between environmental conditions, organisms, tissues, or other factors (3, 20, 21). As a special case, DEA is used to detect imbalanced transcription of the maternal and paternal alleles in a diploid genome. The imbalance phenomenon is called allele-specific expression (ASE) (22). Possible epigenetic mechanisms include genomic imprinting by DNA methylation or chromatin modification (23). ASE has been associated with cancer and other human diseases (22). ASE may be of evolutionary importance, as it is seen in interspecies hybrids of animals and plants (24). ASE has been detected in the Neanderthal genes inherited by some modern humans (25). ASE has been documented in seeds of plant hybrids, including hybrids of model organisms (26) as well as important crops (27). MetaImprint (28), the Plant Imprinting Database (29), and ASMdb (30) are examples of databases of imprinted genes.

When using RNA-seq for ASE detection, a critical step is the computational association of each RNA sequence with its putative source allele. At least three basic approaches have been used to accomplish this. (1) A first approach aligns reads to a concatenation of two parental references, trusting the aligner to choose the better target. This was done in a study of crosses involving three species of the flowering plant Arabidopsis, A. thaliana, A. lyrata, and A. halleri (31). In this case, references were available for only two of the three species, but a reference transcriptome for the third species was generated ad hoc by applying the Trinity (32) assembler to RNA-seq from that species. (2) A second approach aligns each read pair to each parental reference separately, then compares the quality of both alignments. This was used in two studies of crosses between ecotypes of Arabidopsis thaliana (26, 33). In these studies, each parent’s transcriptome reference was computed by aligning its RNA-seq reads to a common reference, then customizing that reference with consensus polishing software (34). (3) A third approach aligns reads to a single reference to bin them by gene, and then analyzes the reads for known single-nucleotide polymorphisms (SNPs) between the parental alleles of each gene. This approach has been used in many studies including studies of Arabidopsis thaliana (35, 36) and Mus musculus (37) intra-species crosses. This general approach has also been implemented in software (38, 39).

Machine learning has been applied to similar challenges. Linear regression analysis was recently used for DEA with RNA-seq data (40). Machine learning was recently used to detect novel isoforms in RNA-seq data (41). Machine learning has been used to detect ASE, without RNA-seq data, by searching for methylation signals in genomic DNA (42, 43).

We propose a machine learning approach to binning RNA-seq reads for the purpose of ASE detection in interspecies hybrid crosses. Our approach is inspired by trio binning, a technique that assigns DNA reads to either parental haplotype before computing a haplotype-resolved genome assembly of their cross (44). We train a binary classifier to assign each RNA-seq read pair to its correct parent by analyzing features of its alignments to both parental references. Once trained, the classifier could be deployed on RNA-seq from the hybrid cross to assist ASE detection. Given the alignment of one read pair to both parental references, the classifier would choose the more likely parental alignment, which would indicate the most likely allele and gene of origin. Our results indicate that using a combination of an aligner plus machine learning is more accurate than using an aligner by itself. The resulting boost in parent-of-origin classification accuracy offers potential to boost the quality of ASE detection.

We developed the process illustrated in Fig. 1. The training pipeline (Fig. 1, A) uses RNA-seq read pairs from both parents, as well as a reference genome or transcriptome from each parent. Every read pair is aligned to both parental references. Features are extracted from the alignments by a Python script that parses the aligner output. A machine-learning model is trained to predict the parent-of-origin per read pair based on the alignment features. The prediction pipeline (Fig. 1, B) would apply a similar pipeline to read pairs from the hybrid offspring. In this case, the unknown parent-of-origin per read pair must be predicted. To assess model accuracy in hybrids, we simulated hybrid data by combining real RNA-seq data from both parents in equal quantities. Possible limitations of the simulation are addressed in Discussion.

Fig. 1 Process for classifying RNA-seq by machine learning. Data from parents P1 and P2 are represented by blue and yellow, respectively, while data from their hybrid offspring are green. A) Training uses many RNA-seq read pairs from both parents, P1 and P2. For compactness, only one read pair is shown here, and it is from parent P2. The read pair is aligned to both parent references. Features are extracted from the best alignment to each reference. A binary classifier, such as the random forest shown here, is trained to examine the features for each read pair and assign it to the correct parent. B) The trained classifier can be used to classify RNA-seq read pairs from a hybrid. Classifying a read pair as P2, as shown, implies its RNA is likely transcribed from the P2 haplotype in the hybrid genome. The classified reads can be binned by gene and quantified to detect allele-specific gene expression in the hybrid.

We searched INSDC databases (45) to find pairs of organisms each having (1) a reference genome assembly, and (2) a reference transcriptome assembly, and (3) at least one RNA-seq dataset. We filtered for pairs of organisms known to hybridize and belonging to the same genus. We filtered further for organism pairs having RNA-seq datasets that were comparable by read length and sample collection methods. In some cases, after low map rates were seen, different RNA-seq data was sought out. The search yielded data for: (1) two species of the flowering plant genus Arabidopsis; (2) two species of the cultivated plant genus Brassica; (3) two strains of the mouse species Mus musculus; and (4) two species of the equine genus Equus.

We selected mapping software packages based on the software’s ability to generate SAM/BAM files (46) that include the required and optional fields required for our alignment feature extraction process. For mapping to reference transcripts, we selected the aligners Bowtie2 (6,7) and STAR (8,9), which we configured for mapping to RNA references lacking introns. For mapping to reference genomes, we selected HiSat2 (18,19), which incorporates Bowtie2, and STAR using its default configuration for mapping to DNA containing introns. For an additional comparison, the alignment-free mapping software Salmon was included for transcript mapping, though its outputs did not satisfy the criteria for alignment feature extraction.

For each genus, RNA-seq read pairs from two parental organisms were combined in equal quantities. Reads were mapped to both parental references, using either both transcriptomes or both genomes. For machine-learning approaches, fifty-three features were selected to represent the pairwise sequence alignment generated by one aligner of one RNA-seq read pair to two parental references, where both references were either transcriptomes or genomes. The features are listed in Table 1. A portion of the reads were used to train a random forest model (47). A separate portion of the reads were used to quantify the predictive performance of the trained model. The models’ performance was compared to these other approaches: comparing alignment scores or relying on the mapper directly.

Table 1 Alignment Features used for Machine Learning

Feature type

Extraction

Technical notes

AS: Alignment Score

ED: Edit Distance

MM: Mismatch count

HQMM: HQ mismatch count

GO: Gap Open count

GE: Gap Extend count

INS: Insertion count

HQINS: HQ insertion count

DELS: Deletion count

HQDEL: HQ deletion count

Taken from:

P1 R1,
P1 R2,
P2 R1,
P2 R2,

10 feature types,

40 features total.

High-quality (HQ) means that the base call quality score is the maximal value. The HQ requirement was applied to the one base involved in a mismatch or insertion, and to the two surrounding bases for deletion. INS or DEL refer to an extra or missing base in the read, respectively. GO is the number of separate indels, and GE is the number of bases in indels.

AS diff

ED diff

MM diff

HQMM diff

GO diff

GE diff

INS diff

DELS diff

HQINS diff

HQDEL diff

MAT diff

Subtract

(P1 R1 + P1 R2)

from

(P2 R1 + P2 R2)

Each difference represents the sum over the read pair alignments to parent 2 minus the equivalent sum for parent 1.

MAT is the matched base count. See A) for other feature types.

Span diff

Subtract P1 span

from P2 span.

Span is the length of the read pair alignment along the reference.

Parent choice

Compare P1 to P2

Use -1 or +1 to indicate whether P1 or P2 had the greater alignment score, respectively, or 0 if tied.

After an aligner, such as Bowtie2, aligned an RNA-seq read pair to both parental references, the one best alignment per reference was selected based on flags in the aligner output. Then, 53 features were extracted for machine learning. A) These 10 features, taken directly from the aligner output, describe the alignment of one read to one reference. Four sets of 10 we extracted to represent reads R1 and R2 aligned to references of parents P1 and P2. B, C, D) These features were derived from combinations of those in part A.

Results with genus Arabidopsis

This experiment used public data from A. lyrata and A. halleri, two species of the flowering plant genus Arabidopsis. The transcriptome results are characterized in Table 2. As shown in column A, we achieved 83% accuracy by running Bowtie2 on the concatenation of parental references. As shown in column B, we achieved 81% accuracy by running Bowtie2 twice, once on each parent reference, then choosing the one alignment having the higher alignment score, which is one of the criteria used by aligners internally. As shown in column C, we achieved 96% accuracy by applying the random forest post-process to the alignment generated for column B. Out of these three approaches, the random forest achieved the highest performance by all measures including accuracy, F1-score, and MCC. Columns D, E, and F show the same experiment repeated using the STAR aligner configured for transcript alignments. Again, the highest performance was achieved with random forest. For comparison, column G characterizes the performance of an alignment-free approach.

Table 2 Classification Performance with Reference Transcripts of genus Arabidopsis

Arabidopsis RNA

Bowtie2

Bo_AS

Bo_RF

STAR

St_AS

St_RF

Salmon

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref:

Ties:

82.8%

84.7%

80.9%

81.6%

83.1%

0.656

51.1%

81.3%

71.3%

91.2%

89.0%

79.2%

0.638

40.1%

13.6%

96.0%

92.1%

99.9%

95.8%

0.923

99.7%

46.1%

80.9%

81.6%

80.2%

80.5%

81.1%

0.618

50.5%

80.9%

69.9%

92.0%

89.7%

78.6%

0.634

49.8%

14.7%

92.3%

92.0%

92.5%

92.3%

0.846

98.4%

46.1%

75.5%

76.9%

74.0%

74.8%

75.8%

0.510

51.0%

Performance metrics for parent-of-origin classification in Arabidopsis. In all seven approaches, RNA-seq read pairs were assigned to either of two reference transcriptomes. Whether used with Bowtie2 or STAR, the random forest method demonstrated superior performance. For the sake of directional statistics like sensitivity, species A. lyrata and A. halleri were designated as the negative and positive classes, respectively. A) Parent chosen by the Bowtie2 aligner. B) Parent chosen by comparing Bowtie2 alignment scores. C) Parent chosen by the random forest classifier using Bowtie2 alignment features. D, E, F) Similar to columns A, B, C, but using the STAR aligner, configured to avoid splicing. G) Parent chosen by the alignment-free software, Salmon.

To evaluate alternate machine learning architectures, the Bowtie2 experiment was repeated with three other architectures: a gradient boosting classifier, a support vector machine, and a multi-layer perceptron. The results, described in Table S1, did not exceed those of Table 2. To evaluate whether more trees would help the random forest, the experiment was repeated after increasing the number of trees within the random forest. The results, described in Table S2, did not exceed those of Table 2. Therefore, the default random forest model was used for the remaining experiments.

Whereas Table 2 showed results with the transcriptome references, Table 3 shows the results using the genome references and splice-aware alignment. The highest performance, as measured by accuracy, F1, or MCC, was achieved with the machine-learning post-process (columns C and F). The overall highest accuracy on Arabidopsis data was achieved with HiSat2 plus random forest, an in this case, the accuracy climbed from 67% with the aligner alone to 96% with machine learning.

Table 3 Classification Performance with Reference Genomes of genus Arabidopsis

Arabidopsis DNA

HiSat2

Hi_AS

Hi_RF

STAR

St_AS

St_RF

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref:

Ties:

66.8%

38.3%

95.4%

89.2%

53.6%

0.410

28.7%

82.8%

72.0%

93.5%

91.8%

80.7%

0.671

39.2%

13.2%

96.2%

93.2%

99.1%

99.0%

96.0%

0.925

99.7%

47.1%

67.5%

40.5%

94.5%

88.0%

55.4%

0.415

30.0%

81.4%

69.8%

93.1%

91.0%

79.0%

0.647

38.4%

14.7%

92.6%

93.8%

91.5%

91.7%

92.7%

0.853

98.5%

51.2%

Performance metrics for parent-of-origin classification in Arabidopsis. In all six approaches, RNA-seq read pairs were assigned to either of two reference genomes. Whether used with HiSat2 or STAR, the random forest led to superior accuracy, F1, and MCC. For the sake of directional statistics like sensitivity, species A. lyrata and A. halleri were designated as the negative and positive classes, respectively. A) Parent chosen by the HiSat2 aligner. B) Parent chosen by comparing HiSat2 alignment scores. C) Parent chosen by the random forest classifier using HiSat2 alignment features. D, E, F) Similar to columns A, B, and C but using the STAR aligner, configured for splicing.

Mild class imbalance is seen in Table 2. For example, in column C, sensitivity < specificity and precision > recall. (Recall and sensitivity are the same in binary classification). Also, the “Pos Pref” statistic shows that 46% of reads were assigned to the positive class, though the sample contained 50% from each class. These statistics reveal a bias for choosing the negative class (A. lyrata) and a tendency to mistakenly assign A. halleri reads to the A. lyrata parent. Investigation showed the mistakes were non-random, with a few transcripts attracting large portions of the mistaken alignments. Class bias is seen again in Table 3. The two aligners had a positive preference of 29% and 30%, while the random forest post-processor brought these more in balance with values of 47% and 51%. Thus, the trend in class bias seen with genomes was opposite of that seen with transcriptomes. Possible causes and mitigations for the bias will be addressed in the Discussion.

Results with genus Brassica

This experiment used public data from B. oleracea and B. rapa, two species in genus Brassica which includes many cabbage-like cultivars consumed by humans. The transcriptome and genome results are characterized in Tables 4 and 5 respectively. The aligner-plus-machine-learning combination showed the best performance in every trial, as measured by accuracy, F1, or MCC.

Table 4 Classification Performance with Reference Transcripts of genus Brassica

Brassica RNA

Bowtie2

Bo_AS

Bo_RF

STAR

St_AS

St_RF

Salmon

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref:

Ties:

88.7%

92.3%

85.1%

86.1%

89.1%

0.777

52.0%

91.5%

94.2%

88.9%

89.4%

91.7%

0.832

52.6%

3.2%

93.6%

92.5%

94.6%

94.5%

93.5%

0.872

98.3%

98.2%

49.0%

90.2%

93.4%

87.1%

87.8%

90.5%

0.806

51.8%

92.2%

94.6%

89.9%

90.3%

92.4%

0.846

52.4%

5.0%

93.8%

93.7%

93.9%

93.8%

0.876

98.4%

49.9%

84.6%

88.9%

80.3%

81.9%

85.2%

0.694

52.5%

Performance metrics for transcript-based parent-of-origin classification in Brassica. Whether used with Bowtie2 or STAR, the random forest improved the accuracy, F1, and MCC. For directional statistics, species B. rapa and B. oleracea were considered the negative and positive classes, respectively. A, B, C) Using the Bowtie2 aligner, a parent was chosen by the aligner, or by comparing alignment scores, or by the random forest, respectively. D, E, F) Similar to A, B, and C but using the STAR aligner, configured to avoid splicing. G) Parent chosen by the alignment-free software, Salmon.

Table 5 Classification Performance with Reference Genomes of genus Brassica

Brassica DNA

HiSat2

Hi_AS

Hi_RF

STAR

St_AS

St_RF

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref:

Ties:

91.8%

95.1%

88.6%

89.3%

92.1%

0.839

53.3%

92.9%

95.6%

90.3%

90.8%

93.1%

0.860

52.7%

3.4%

95.1%

94.7%

95.6%

95.5%

95.1%

0.903

99.0%

98.9%

49.5%

93.0%

96.0%

90.1%

90.6%

93.2%

0.862

53.0%

92.3%

94.6%

90.0%

90.4%

92.5%

0.847

52.3%

5.6%

93.9%

93.2%

94.7%

94.6%

93.9%

0.879

98.4%

98.3%

49.2%

Performance metrics for genome-based parent-of-origin classification in Brassica. Whether used with HiSat2 or STAR, the random forest improved the accuracy, F1, and MCC. For directional statistics, species B. rapa and B. oleracea were considered the negative and positive classes, respectively. A, B, C) Parent chosen by the HiSat2 aligner, or by comparing HiSat2 alignment scores, or by the random forest using HiSat2 alignment features, respectively. D, E, F) Similar to columns A, B, and C but using the STAR aligner, configured for splicing.

The overall highest accuracy with Brassica was achieved with HiSat2 plus random forest, and in this case, the accuracy climbed from 89% with the aligner alone to 94% with machine learning. Each mapping software package showed better performance with the Brassica data (Tables 4 and 5), than with the Arabidopsis data (Tables 2 and 3). Consequently, the performance differences between approaches were smaller. Nevertheless, the random forest gains seem smaller in Brassica than in Arabidopsis. A possible explanation is suggested by the smaller numbers of alignment score ties in Brassica (3%-6%) than in Arabidopsis (13%-15%). It may be that the random forest has the most effect when the aligners generate many equally good (or equally bad) alignments between the two parental references.

Mild class imbalance is seen in the Brassica results, though the positive preference was always close to 50%. The alignment-based methods (columns A, B, D, and E) showed the most imbalance between sensitivity and specificity, and these imbalances were reduced by the machine-learning post-process (columns C and F).

Results with genus Mus

This experiment used public data from Mus musculus, the mouse species that serves as a model organism for mammalian genetics. The B6 and D2 strains are inbred strains of the same species, developed and maintained in the laboratory. Parent-of-origin classification accuracy was low, barely exceeding the 50% expected by random guessing. Accuracy did not exceed 57% by any method tested. These results suggest that these within-species hybrids did not present enough genomic variants for parent-of-origin classification. For this reason, the results are shown in supplement, Tables S3 and S4, and the procedure is not recommended for crosses between very similar genotypes. Nevertheless, even on this dataset, the random forest post-processor added value, achieving 56%-57% accuracy, compared to 53%-56% achieved by the aligners directly. As shown in the supplemental tables, the random forest predictions showed positive class bias (toward D2) using transcript references, but negative class bias using genome references.

Results with genus Equus

This experiment used public data from the genus Equus. The equine species of horse, E. caballus, and donkey, E. asinus, can hybridize to yield a mule or hinny, depending on which parent is male or female (48). Results are shown in Tables 6 and 7. The machine-learning method provided the best performance using transcriptomes or genomes. These results indicate that our method is suitable for interspecies hybrids of animals as well as plants.

Table 6 Classification Performance with Reference Transcripts of genus Equus

Equus RNA

Bowtie2

Bo_AS

Bo_RF

STAR

St_AS

St_RF

Salmon

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref

Ties:

71.9%

78.2%

65.6%

69.5%

73.6%

0.441

56.3%

76.6%

74.9%

78.3%

77.6%

76.2%

0.533

48.3%

37.3%

81.7%

91.0%

72.4%

76.7%

83.2%

0.645

92.0%

91.5%

59.3%

71.2%

79.0%

63.3%

68.3%

73.3%

0.428

57.8%

76.6%

75.1%

78.1%

77.4%

76.2%

0.532

48.5%

40.1%

85.2%

79.1%

91.4%

90.2%

84.3%

0.710

93.3%

93.0%

43.9%

63.7%

74.8%

52.5%

61.2%

67.3%

0.281

61.2%

By many measures including accuracy, F1, and MCC, the random forest performance surpassed that of the other methods tested. For directional statistics, donkey and horse were considered the negative and positive classes, respectively. A, B, C) Using the Bowtie2 aligner, a parent was chosen by the aligner, or by comparing alignment scores, or by the random forest, respectively. D, E, F) Similar to A, B, and C but using the STAR aligner, configured to avoid splicing. G) Parent chosen by the alignment-free software, Salmon.

Table 7 Classification Performance with Reference Genomes of genus Equus

Equus DNA

HiSat2

Hi_AS

Hi_RF

STAR

St_AS

St_RF

Accuracy:

Sensitivity:

Specificity:

Precision:

F1-score:

MCC:

AUPRC:

AUROC:

Pos Pref

Ties:

73.1%

72.4%

73.7%

73.3%

72.9%

0.461

49.4%

78.9%

78.2%

79.5%

79.2%

78.7%

0.577

49.4%

33.6%

84.3%

91.0%

77.5%

80.2%

85.3%

0.692

93.9%

93.7%

56.7%

71.6%

71.1%

72.2%

71.9%

71.5%

0.432

49.4%

78.3%

77.6%

79.0%

78.7%

78.1%

0.566

49.3%

36.9%

86.2%

81.5%

90.8%

89.8%

85.5%

0.726

94.1%

94.0%

45.4%

Table 7 Legend: By many measures including accuracy, F1, and MCC, the random forest performance surpassed that of the other methods tested. For directional statistics, donkey and horse were considered the negative and positive classes, respectively. A, B, C) Parent chosen by the HiSat2 aligner, or by comparing HiSat2 alignment scores, or by the random forest using HiSat2 alignment features, respectively. D, E, F) Similar to columns A, B, and C but using the STAR aligner, configured for splicing.

The overall highest accuracy with Equus was achieved with STAR DNA plus random forest, and in this case, the accuracy climbed from 72% with the aligner alone to 86% with machine learning. The random forest post-process introduced mild levels of class bias on the Equus data, though the direction varied. Positive preference increased when the random forest was used with Bowtie2 or HiSat2 but decreased when used with STAR. This indicates that, at least on these data, the bias is a computational artifact.

Results Summary and Generalization

The results on the three interspecies datasets are compared in Fig. 2. For each of the Arabidopsis, Brassica, and Equus datasets, the approach that achieved the highest accuracy is compared to the other approaches tested with the same aligner. On each dataset, the random forest accuracy surpassed that from relying on the alignment scores or the aligner by itself. On Arabidopsis and Brassica, the maximum accuracy was achieved using HiSat2 alignments to the genome reference, plus the random forest. On Equus, the maximum accuracy was achieved using STAR alignments to the genome reference, plus the random forest. Since no one approach worked best on all three datasets, it may be advantageous to experiment with several aligners, as we have done here, when using other datasets.

Fig. 2 Comparison of three models by accuracy. Using the genera Arabidopsis and Brassica, the parent-of-origin read classification accuracy was maximized using HiSat2 alignments to genomes followed by random forest. Using the genus Equus, accuracy was maximized using STAR alignments to genomes followed by random forest. The vertical axis ranges from 50% (guessing) to 100% (perfect).

Unlike many machine learning models, random forest models can be interrogated for indications of which features were most important. Fig. 3 illustrates the top five features reported for the three most performant models in this study. (In Supplement, Table S5 shows the top five features for all twelve models presented in this paper.) The different rankings in this comparison indicate that the feature importance varied with the dataset and the aligner. However, the features in the intersection of these three lists are Parent (a number indicating the parent with higher alignment score) and MAT diff (the P2-P1 difference in matched bases). Also, there was at least one HQ feature in each list. Our HQ features count problems (mismatches, insertions, deletions) that involve positions in reads assigned the maximal quality score by the sequencing instrument software. Notably, the match counts and quality scores are not reflected in the alignment scores generated by the mappers used here. Thus, it seems the random forest models learned to rely on complementary features.

Fig. 3 Alignment Features Ranked by Importance. Top five alignment features, ranked by importance, used by the random forest model. The three columns here correspond to those in Figure 2. The figure indicates that different models relied on different features. Labels like “P2 R1” refer to the alignment of read 1 to parent 2. Labels with “HQ” count only those events that involve a maximal base call quality score in the read. Labels with “diff” include the difference between parent 2 and 1 alignments. Labels with “MM” and “MAT” refer to mismatched or matched base counts, respectively. Labels with “INS” or “DEL” refer to bases inserted or deleted in the read, respectively. The label “PARENT” indicates a feature that was 1 if P2 had the better combined alignment score, or -1 if P1 had the better score, or 0 if they were tied.

In all the experiments shown so far, each model was trained on one set of read pairs and tested on another. However, each test set was carved from the same RNA-seq runs as its cognate training set. In an experiment on actual hybrids, the training RNA-seq from the parents and would come from different sequencing libraries and sequencing runs compared to the test RNA-seq from the hybrid. We wondered how well our models would perform on RNA-seq runs other than those on which it was trained. In other words, could the model generalize.

We performed a test with the model trained on Equus. We downloaded two additional parental RNA-seq runs from the same Equus project. We selected runs for testing that were generated from different individuals, i.e., a different horse and a different donkey, compared t the runs used for training. These secondary runs were aligned to the parental transcriptomes with Bowtie2. The random forest model that was trained on the primary runs was used without retraining to make predictions on the secondary runs. The model achieved 82.5% accuracy, 0.659 MCC, and 91.4% AUROC on the secondary runs. The full results are shown in Table S6. Overall, the model performed similarly whether evaluated with the primary or secondary runs. Thus, the model was able to generalize and maintain accuracy with RNA-seq runs other than those on which it was trained.

We demonstrated process improvement for one part of pipelines that study allele-specific expression (ASE) in hybrids. We examined several procedures for classifying paired, short-read RNA-seq data according to their parent of origin. The best performance was achieved by a novel process that applies machine-learning to features extracted from pairwise sequence alignments. We believe that this is the first application of machine learning to the problem of binning RNA-seq reads by parent-of-origin.

Using public data representing seven species in four genera, we trained classifiers to bin RNA-seq read pairs by their parent of origin. The reason for training on parental reads is that their true parent of origin is known in all cases. For evaluation, we used combinations of parental reads, again because their true parent of origin is known. Thus, the evaluations were conducted on simulated hybrid datasets composed of mixtures of real RNA-seq data.

We tested with one transcript aligner, one genome aligner, and one aligner that worked with transcripts or genomes. For each aligner, we evaluated three configurations: relying on the aligner to choose the parent of origin, or by choosing a parent based on the better alignment score, or by feeding alignment features into a machine-learning post-processor. In all our experiments, the best performance was achieved using the post-processor.

For machine learning, we used random forest models and interrogated them for feature importance. The top five features per model included on based on alignment scores and two that are independent of alignment scores. It appears that the models learned to rely on complementary features. The features that we extracted from alignments included imperfection counts (mismatches, indels), differences in imperfection counts, and the difference in match counts. The models were not given the actual match counts, read lengths, or alignment spans, because these might allow models to distinguish reads by their RNA-seq library or sequencing run. If models had focused on library-specific or run-specific features, they would not generalize to other RNA-seq data. We saw confirmation that our models could generalize when we trained and tested a model on one Equus RNA-seq dataset but also tested the trained model on another RNA-seq dataset. The model achieved similar performance on both datasets.

We relied on simulations for evaluation. We simulated hybrid RNA-seq by combining real RNA-seq from two individuals representing the parents. Since simulation can fail to capture the complexity of real data, we speculate how it might affect our approach. First, real hybrid individuals inherit only one allele per gene per parent, but our simulated hybrid data contained up to four alleles per gene. Therefore, it seems that real data from hybrids would be less complex than what we simulated. So, we expect classifiers trained on simulated hybrids to perform as well on real hybrids. Second, a real hybrid individual could express an allele not expressed in either parent. But our classifiers are trained on alignments not sequences. So, assuming the hybrid’s allele or gene would be represented in the parental references, we expect our training regime would prepare classifiers to handle these cases well. Third, RNA-seq from real hybrids may reflect novel alleles or genes generated, for example, by a mid-gene crossing-over event during meiosis in either parent’s germ cell. These cases should be rare, but we expect the classifier to do well with these, since each section of the new allele is likely to have a better alignment to its true parent than to the other parent.

The fact that machine learning could improve on aligner accuracy should not be taken as criticism of aligners. The aligners are general-purpose tools that have enabled many biological discoveries. It is possible that the aligners could have been tuned for our specific tasks using their command-line parameters for adjusting their alignment score weights or formulas. It is, however, far from clear how to optimize those parameters. Instead, we applied machine learning to the aligner outputs. The training step exploited the large quantities of data provided by millions of reads and their accurate parental labels. It is likely that the machine learning models excelled because they were each trained specifically for one species pair.

In ASE studies, accurate quantification is key. Meta-analyses of published lists of ASE genes in Arabidopsis have noted little overlap between the lists (33,49). Other meta-analyses have called into question published claims of weak imprinting in humans and mice (50). Simulations have shown that ASE discoveries are sensitive to underlying map bias (51) and read trimming (52). Our method appears to increase the portions of RNA-seq reads assigned to the correct parent of a hybrid organism. This improvement could lead to improved sensitivity and specificity and thus higher confidence in ASE detection. However, several of our results indicated map bias. Map bias can result from different levels of completeness or quality between the two references. In our results, the direction of the bias differed by computational approach, hinting that the bias could be a computational artifact. Regardless of its cause, map bias can induce false conclusions about ASE (51). Our machine-learning approach did not eliminate map bias and it exacerbated it on some data. However, our approach would permit map bias detection and mitigation. One mitigation strategy would test several aligners, as we have done here, and choose the one leaving the least bias after machine learning. Another would incorporate prior class weights to a model’s loss function. A heuristic described by King and Zeng (53) is already implemented by the fit function of scikit-learn models. Alternatively, bias could be quantified and compensated for during ASE statistical analysis. During the training step, when we align RNA-seq for which the true parent-of-origin is known, mismatch rates could be measured for that dataset. Statistical compensation could be adjusted to fit the observed bias per gene or per transcript.

To put our approach into practice on a real hybrid organism, an experimenter would need to sequence RNA from one or more hybrids plus both of its parents, then align all the RNA-seq data to both parental references. Either two genome references or two transcriptome references could be used. (For organisms lacking trusted references, a reference transcriptome could be generated by de novo assembly of either parental RNA-seq dataset using e.g., Trinity (32), with the limitation that this ad hoc reference would only reflect genes expressed by the parent.) The experimenter would train and test a classifier on the parental RNA-seq alignments, then apply the trained classifier to the hybrid RNA-seq alignments. The experimenter would infer the allele of origin per hybrid RNA-seq read pair using its alignment to the parental reference predicted by the classifier. The counts per allele per gene could be given to any differential expression analysis pipeline for ASE detection. Implementation of such a pipeline is left for future work. One implementation detail will be the proper assignment of read pairs that align to one parent exclusively; our study considered only the read pairs that aligned to both parental references.

Our approach decides the parent-of-origin before investigating the gene-of-origin. It seems common practice to choose the gene first, based on alignments, and the parental allele second, based on sequence variants. Our results show that the initial alignment often chooses the wrong parent, and that deciding the parent-of-origin first can avoid that pitfall. Furthermore, our approach of deciding the parent first allows for the creation of labeled training sets and the injection of machine learning.

Random forest models were trained on alignment features. When trained and tested on RNA-seq paired-end reads aligned to reference transcriptomes or genomes, these models provided higher accuracy at parent-of-origin classification compared to any aligner by itself. This study establishes that machine learning can play a role in RNA-seq analysis of allele-specific gene expression in hybrids.

All mapping software ran under Linux (Rocky 9.1) on the Saga computing cluster in Norway. Alignment files in SAM/BAM format were manipulated with samtools (46) version 1.16.1. Software compilations used GCC 11.3.0. The machine learning software ran inside Jupyter notebooks on Google CoLab Pro virtual computers with one CPU and 12 GB RAM.

Reference genome files were downloaded in FASTA format from Ensembl (54). Ensembl provides a “primary_assembly” file when chromosome sequences are available. All the assemblies used here, except mouse b6, lacked chromosome assignments. For consistency, the “toplevel” files of scaffolds were used in all cases. Reference transcript files were downloaded from the cDNA directories corresponding to these genomes. The cDNA files contain intron-free transcript references. RNA-seq files were downloaded from NCBI SRA (55). In every case, the SRA normalized file was selected in order to obtain original base call quality values. SRA files were processed with ‘fastq-dump --split-3’ from the SRA-Toolkit version 3.0.3 to create paired FASTQ files. Where two database accessions are given below, the first is a general one given in the publication and the second indicates the data subset used here.

● Arabidopsis. We used the A. halleri reference with accession GCA_900078215 (no publication) and the A. lyrata reference with accession GCA_000004255.1 (56). (Newer references for A. lyrata became available recently (57) but too late for inclusion here.) All the RNA-seq data came from a study of RNA editing across Arabidopsis species (58) in which total leaf RNA samples were subjected to rRNA depletion and 2x100 Illumina sequencing. The A. lyrata and A. halleri RNA-seq datasets have DDJB accessions DRA007657 and DRA007658 and SRA accessions DRR161380 and DRR161381, respectively.

● Brassica. We used the B. rapa reference with accession GCA_000309985.1 (59,60) and the B. oleracea reference with accession GCA_000695525.1 (61). The B. rapa RNA-seq came from a study of heterosis in Chinese cabbage hybrids, Project PRJNA876066 (62). The RNA-seq used here was from one of the inbred, non-hybrid parents, C-1 SRR21735970. The B. oleracea RNA-seq derived from a study of a Chinese kale allotetraploid, Project PRJNA885390 (63). The RNA-seq used here was from the diploid parent, CC SRR21778809. Both RNA-seq datasets reflect 2x150 Illumina sequencing.

● Mus. We used the M. musculus (mouse) C57BL/6J (‘B6’) reference with accession GCA_000001635.9 (64) and the DBA/2J (‘D2’) reference with GCA_001624505.1 (65). All the RNA-seq was 2x100 Illumina from a study of gene expression in the retinas of the two mouse strains (66). We used RNA-seq with accessions SRR8690244 and SRR8690250.

● Equus. We used the E. caballus (horse) reference with accession GCA_002863925.1 (67) and the E. asinus (donkey) reference GCA_016077325.2 (48). Training used 2x151 Illumina RNA-seq runs with accessions SRR23724220 and SRR24443170. Validation used data from different individuals that were parts of the same studies, accessions SRR23724221 and SRR24443174.

Raw reads were end-trimmed to remove adapter, low-quality bases, and N base calls using Trim Galore (68) version 0.6.10 using command line ‘trim_galore --cores 4 --trim-n --paired’ plus read 1 and 2 filenames. The read sets were not subjected to duplicate removal, which can be unhelpful unless universal molecular identifiers are present (69) .

The mappers selected were ones compared in (10). The mappers selected use a range of algorithmic approaches including hash tables of k-mer to position in Salmon (4), the FM-index (70,71) and Burrows-Wheeler transform (72,73) in Bowtie2 (6,7), and a suffix array (74) in STAR (8,9). Some other mappers had to be excluded because their output formats did not support our alignment feature extraction: DART (11), Kallisto (5), and GSNAP (75).

Using Bowtie2 (6,7) version 2.4.5, references were indexed with default parameters. Reads were aligned with command line ‘bowtie2 --no-unal --no-mixed --no-discordant --sensitive --end-to-end --threads 4’ plus options to specify the target sequence filename, the R1 and R2 filenames, and the output filename. We extracted the primary alignment per read using BAM flags. Although Bowtie2 has an option to report multiple alignments per read, the option was not used because the documentation says these alignments are not necessarily the best. HiSat2 (18,19) version 2.2.1 was used to align RNA paired reads to genomic DNA using all the same options as Bowtie2.

Using STAR (8,9) version 2.7.10b, references were indexed with command line ‘STAR --runThreadN 4 --runMode genomeGenerate’, plus options to specify the target directory and read filenames. Reads were aligned with command line ‘STAR --runThreadN 4 --outSAMattributes NH AS nM NM MD --outSAMtype BAM Unsorted --readFilesCommand gunzip -c’ plus options to specify the index directory and read 1 and 2 filenames. If any process issued recommendations for larger values, such as for parameters limitGenomeGenerateRAM or genomeSAindexNbases, then the program was re-run with the recommended values. For alignments to transcripts, the options ‘--alignIntronMin 100000 --alignIntronMax 0’ were added to preclude splicing. We extracted the primary alignment per read using BAM flags. Although STAR has an option to report multiple primary alignments, the option was not used because the documentation says no alignments would be reported if the observed number exceeded the given number.

Using Salmon (4) 1.9.0 in mapper mode (i.e. not in conjunction with an aligner), references were indexed with command line ‘salmon index’ and mapped with command line ‘salmon quant --index <index_dir> --libType A --threads 4 --output salmon_out --writeMappings=Aligned.sam’ plus options to specify the read filenames. Recent versions of Salmon recommend using ‘decoys’ to allow Salmon to identify RNA-seq from isoforms missing from the transcriptome, but decoys were not generated or used here. Salmon outputs were not subjected to machine learning because they lacked alignment fields required by our process.

The SAM/BAM format output files were rendered to text with samtools (46) and parsed by a custom Python script for feature extraction. The script relied on the following fields which are optional in BAM files: AS, XM, XO, XG, NM, MD. The script counted events such as mismatches and indels by parsing the CIGAR and MD strings. The script distinguished between such events by whether or not their base call quality score was maximal. The script accepted the maximum base call quality score encoding (e.g., ‘F’ or ‘K’) as a parameter, and the parameter value was selected by visual inspection of scores in the BAM files in each read set. The feature extractor ignored soft clipping; a cigar string like ‘1S99M’ with one soft-clipped base was treated as 100 aligned bases. Thus, the script could report more mismatches than given in the ‘NM’ field of the BAM file. Alignment spans and read lengths were not used as features since their distributions could be specific to an RNA-seq library or run.

Traditional machine learning was implemented with scikit-learn (76) version 1.2.2. The RandomForestClassifier class was used for random forest models (47). Feature ranking used the mean decrease in impurity (MDI) method. The GradientBoostingClassifier class was used as a gradient boosting model (77). The SVC class was used as a support vector machine (78). The multi-layer perceptron was built with Keras (79).

Reads were aligned in the order they appeared in the FASTQ files, which is essentially random. Models were trained using the first N read pairs that aligned to both parent references, with the first 80% used for training and the remaining 20% used for evaluation. Each experiment used one pair of RNA-seq FASTQ files, and N was adjusted according to data availability. N was two million for Equus where reads were most copious; N was one million for Arabidopsis, Brassica, and Mus; N had to be reduced to 400,000 for the experiment on STAR alignments to Brassica genomes, where STAR aligned fewer reads for unknown reasons. During training, models saw equal numbers of alignments from each parent, and alignments were interleaved such that even and odd samples came from different parents.

We employ several statistics to measure performance. Using TP, FP, TN, FN respectively to represent true positive, false positive, true negative, and false negative rates, we define: accuracy = 100*(TP+TN)/(TP+FP+TN+FN), sensitivity = 100*TP/(TP+FN), specificity = 100*TN/(TN+FP), F1 = 200*TP/(2*TP+FP+FN), precision = 100*TP/(TP+FP), and Matthews correlation coefficient or MCC = (TP*TN-FP*FN)/sqrt((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)). AUPRC is the area under the precision-recall curve, a plot of precision vs recall as the classifier’s score threshold varies from 0 to 1. AUROC is the area under the receiver-operator characteristic, a plot of sensitivity vs 1-specificity as the threshold varies. For all these statistics, a higher value is better. We show the map bias as preference for the positive class, “Pos Pref” = (TP+FP)/(TP+FP+TN+FN), for which 50% means no bias. Only the classifiers based on alignment score generated ties, and in these cases, ties were counted, but reads with tie scores were classified randomly for the purpose of generating comparable statistics.

Ethics approval and consent to participate: Not applicable.
Consent for publication: Not applicable.
Availability of data and materials: The text provides public database accessions for all the sequence data used in this study. An archived code repository (https://doi.org/10.5281/zenodo.8199652) contains the scripts and notebooks used for this study.
Competing interests: All the authors declare that they have no competing interest.
Funding: This work is supported in part by the US National Science Foundation (NSF), Award #1920920 and #2125872.
Authors' contributions: JRM designed and executed the experiments and wrote the manuscript. DAA managed the project.
Acknowledgements: The authors thank Granger G Sutton, Karina S Hornslien, and Paul E Grini for comments on the manuscript, and the SIGMA2 Norwegian Research Infrastructure Services for the use of the Saga computing cluster.
Authors' information: JRM is Associate Professor of Computer Science at Shepherd University. DAA is Professor of Computer Science at West Virginia University.

Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet. 2019;20(11):631–56.
Conesa A, Madrigal P, Tarazona S, Gomez-Cabrero D, Cervera A, McPherson A, et al. A survey of best practices for RNA-seq data analysis. Genome Biol. 2016;17(1):13.
Costa-Silva J, Domingues D, Lopes FM. RNA-Seq differential expression analysis: An extended review and a software tool. PLoS ONE. 2017;12(12):e0190152.
Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14(4):417–9.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Dobin A, Gingeras TR. Optimizing RNA-Seq Mapping with STAR. Methods Mol Biol. 2016;1415:245–62.
Srivastava A, Malik L, Sarkar H, Zakeri M, Almodaresi F, Soneson C, et al. Alignment and mapping methodology influence transcript abundance estimation. Genome Biol. 2020;21(1):239.
Lin H-N, Hsu W-L. DART: a fast and accurate RNA-seq mapper with a partitioning strategy. Bioinformatics. 2018;34(2):190–7.
Musich R, Cadle-Davidson L, Osier MV. Comparison of Short-Read Sequence Aligners Indicates Strengths and Weaknesses for Biologists to Consider. Front Plant Sci. 2021;12:657240.
Baruzzo G, Hayer KE, Kim EJ, Di Camillo B, FitzGerald GA, Grant GR. Simulation-based comprehensive benchmarking of RNA-seq aligners. Nat Methods. 2017;14(2):135–9.
Payá-Milans M, Olmstead JW, Nunez G, Rinehart TA, Staton M. Comprehensive evaluation of RNA-seq analysis pipelines in diploid and polyploid species. Gigascience. 2018;7(12).
Nodehi HM, Tabatabaiefar MA, Sehhati M. Selection of optimal bioinformatic tools and proper reference for reducing the alignment error in targeted sequencing data. J Med Signals Sens. 2021;11(1):37–44.
Zanti M, Michailidou K, Loizidou MA, Machattou C, Pirpa P, Christodoulou K, et al. Performance evaluation of pipelines for mapping, variant calling and interval padding, for the analysis of NGS germline panels. BMC Bioinformatics. 2021;22(1):218.
Donato L, Scimone C, Rinaldi C, D’Angelo R, Sidoti A. New evaluation methods of read mapping by 17 aligners on simulated and empirical NGS data: an updated comparison of DNA- and RNA-Seq data from Illumina and Ion Torrent technologies. Neural Comput Appl. 2021;33(22):15669–92.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.
Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15.
Huang H-C, Niu Y, Qin L-X. Differential Expression Analysis for RNA-Seq: An Overview of Statistical Methods and Computational Software. Cancer Inf. 2015;14(Suppl 1):57–67.
Wang T, Li B, Nelson CE, Nabavi S. Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data. BMC Bioinformatics. 2019;20(1):40.
Cleary S, Seoighe C. Perspectives on Allele-Specific Expression. Annu Rev Biomed Data Sci. 2021;4:101–22.
Feil R, Berger F. Convergent evolution of genomic imprinting in plants and mammals. Trends Genet. 2007;23(4):192–9.
Wolf JB, Oakey RJ, Feil R. Imprinted gene expression in hybrids: perturbed mechanisms and evolutionary implications. Heredity. 2014;113(2):167–75.
McCoy RC, Wakefield J, Akey JM. Impacts of Neanderthal-Introgressed Sequences on the Landscape of Human Gene Expression. Cell. 2017;168(5):916–927e12.
van Ekelenburg YS, Hornslien KS, Van Hautegem T, Fendrych M, Van Isterdael G, Bjerkan KN, et al. Spatial and temporal regulation of parent-of-origin allelic expression in the endosperm. Plant Physiol. 2023;191(2):986–1001.
Xu Q, Wu L, Luo Z, Zhang M, Lai J, Li L, et al. DNA demethylation affects imprinted gene expression in maize endosperm. Genome Biol. 2022;23(1):77.
Wei Y, Su J, Liu H, Lv J, Wang F, Yan H, et al. MetaImprint: an information repository of mammalian imprinted genes. Development. 2014;141(12):2516–23.
Picard CL, Gehring M. Identification and comparison of imprinted genes across plant species. Methods Mol Biol. 2020;2093:173–201.
Zhou Q, Guan P, Zhu Z, Cheng S, Zhou C, Wang H, et al. ASMdb: a comprehensive database for allele-specific DNA methylation in diverse organisms. Nucleic Acids Res. 2022;50(D1):D60–71.
He F, Steige KA, Kovacova V, Göbel U, Bouzid M, Keightley PD, et al. Cis-regulatory evolution spotlights species differences in the adaptive potential of gene expression plasticity. Nat Commun. 2021;12(1):3376.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Hornslien KS, Miller JR, Grini PE. Regulation of Parent-of-Origin Allelic Expression in the Endosperm. Plant Physiol. 2019;180(3):1498–519.
Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE. 2014;9(11):e112963.
Gehring M, Missirian V, Henikoff S. Genomic analysis of parent-of-origin allelic expression in Arabidopsis thaliana seeds. PLoS ONE. 2011;6(8):e23687.
Wolff P, Weinhofer I, Seguin J, Roszak P, Beisel C, Donoghue MTA, et al. High-resolution analysis of parent-of-origin allelic expression in the Arabidopsis Endosperm. PLoS Genet. 2011;7(6):e1002126.
Hasin-Brumshtein Y, Hormozdiari F, Martin L, van Nas A, Eskin E, Lusis AJ, et al. Allele-specific expression and eQTL analysis in mouse adipose tissue. BMC Genomics. 2014;15(1):471.
Krueger F, Andrews SR, SNPsplit. Allele-specific splitting of alignments between genomes with known SNP genotypes. [version 2; peer review: 3 approved]. F1000Res. 2016;5:1479.
Duchemin W, Dupont P-Y, Campbell MA, Ganley ARD, Cox MP. HyLiTE: accurate and flexible analysis of gene expression in hybrid and allopolyploid species. BMC Bioinformatics. 2015;16(1):8.
Reed E, Ferrari E, Soloviev M. Quality Control of Gene Expression Data Allows Accurate Quantification of Differentially Expressed Biological Pathways. Curr Bioinform. 2023;18(5):409–27.
Zhang X, Hassan MA, Prendergast JGD. Using machine learning to detect the differential usage of novel gene isoforms. BMC Bioinformatics. 2022;23(1):45.
He J, Sun M, Wang Z, Wang Q, Li Q, Xie H. Characterization and machine learning prediction of allele-specific DNA methylation. Genomics. 2015;106(6):331–9.
Zhang Z, van Dijk F, de Klein N, van Gijn ME, Franke LH, Sinke RJ, et al. Feasibility of predicting allele specific expression from DNA sequencing using machine learning. Sci Rep. 2021;11(1):10606.
Koren S, Rhie A, Walenz BP, Dilthey AT, Bickhart DM, Kingan SB et al. De novo assembly of haplotype-resolved genomes with trio binning. Nat Biotechnol. 2018.
Arita M, Karsch-Mizrachi I, Cochrane G. The international nucleotide sequence database collaboration. Nucleic Acids Res. 2021;49(D1):D121–4.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Breiman L. Random Forests. Springer Science and Business Media LLC; 2001.
Wang C, Li H, Guo Y, Huang J, Sun Y, Min J, et al. Donkey genomes provide new insights into domestication and selection for coat color. Nat Commun. 2020;11(1):6014.
Wyder S, Raissig MT, Grossniklaus U. Consistent Reanalysis of Genome-wide Imprinting Studies in Plants Using Generalized Linear Models Increases Concordance across Datasets. Sci Rep. 2019;9(1):1320.
Edwards CA, Watkinson WMD, Telerman SB, Hulsmann LC, Hamilton RS, Ferguson-Smith AC. Reassessment of weak parent-of-origin expression bias shows it rarely exists outside of known imprinted regions. eLife. 2023;12.
Degner JF, Marioni JC, Pai AA, Pickrell JK, Nkadori E, Gilad Y, et al. Effect of read-mapping biases on detecting allele-specific expression from RNA-sequencing data. Bioinformatics. 2009;25(24):3207–12.
Williams CR, Baccarella A, Parrish JZ, Kim CC. Trimming of sequence reads alters RNA-Seq gene expression estimates. BMC Bioinformatics. 2016;17:103.
King G, Zeng L. Logistic Regression in Rare Events Data. Political Anal. 2001;9(2):137–63.
Cunningham F, Allen JE, Allen J, Alvarez-Jarreta J, Amode MR, Armean IM, et al. Ensembl 2022. Nucleic Acids Res. 2022;50(D1):D988–95.
Sayers EW, Bolton EE, Brister JR, Canese K, Chan J, Comeau DC, et al. Database resources of the National Center for Biotechnology Information in 2023. Nucleic Acids Res. 2023;51(D1):D29–38.
Hu TT, Pattyn P, Bakker EG, Cao J, Cheng J-F, Clark RM, et al. The Arabidopsis lyrata genome sequence and the basis of rapid genome size change. Nat Genet. 2011;43(5):476–81.
Bramsiepe J, Krabberød AK, Bjerkan KN, Alling RM, Johannessen IM, Hornslien KS et al. Structural evidence for MADS-box type I family expansion seen in new assemblies of A. arenosa and A. lyrata [Internet]. bioRxiv. 2023 [cited 2023 Jul 28]. Available from: https://www.biorxiv.org/content/10.1101/2023.05.30.542816v1.
Kawabe A, Furihata HY, Tsujino Y, Kawanabe T, Fujii S, Yoshida T. Divergence of RNA editing among Arabidopsis species. Plant Sci. 2019;280:241–7.
Wang X, Wang H, Wang J, Sun R, Wu J, Liu S, et al. The genome of the mesopolyploid crop species Brassica rapa. Nat Genet. 2011;43(10):1035–9.
Zhang L, Cai X, Wu J, Liu M, Grob S, Cheng F, et al. Improved Brassica rapa reference genome by single-molecule sequencing and chromosome conformation capture technologies. Hortic Res. 2018;5:50.
Parkin IAP, Koh C, Tang H, Robinson SJ, Kagale S, Clarke WE, et al. Transcriptome and methylome profiling reveals relics of genome dominance in the mesopolyploid Brassica oleracea. Genome Biol. 2014;15(6):R77.
Li R, Nie S, Zhang N, Tian M, Zhang L. Transcriptome analysis reveals a major gene expression pattern and important metabolic pathways in the control of heterosis in chinese cabbage. Plants. 2023;12(5).
Zheng W, Shi J, Zhu Z-Y, Jin P, Chen J-H, Zhang L, et al. Transcriptomic analysis of succulent stem development of Chinese kale (Brassica oleracea var. alboglabra Bailey) and its synthetic allotetraploid via RNA sequencing. Front Plant Sci. 2022;13:1004590.
Church DM, Schneider VA, Graves T, Auger K, Cunningham F, Bouk N, et al. Modernizing reference genome assemblies. PLoS Biol. 2011;9(7):e1001091.
Lilue J, Doran AG, Fiddes IT, Abrudan M, Armstrong J, Bennett R, et al. Sixteen diverse laboratory mouse reference genomes define strain-specific haplotypes and novel functional loci. Nat Genet. 2018;50(11):1574–83.
Wang J, Geisert EE, Struebing FL. RNA sequencing profiling of the retina in C57BL/6J and DBA/2J mice: Enhancing the retinal microarray data sets from GeneNetwork. Mol Vis. 2019;25:345–58.
Wade CM, Giulotto E, Sigurdsson S, Zoli M, Gnerre S, Imsland F, et al. Genome sequence, comparative analysis, and population genetics of the domestic horse. Science. 2009;326(5954):865–7.
Krueger F. Trim Galore [Internet]. GitHub. 2023 [cited 2023 Jul 28]. Available from: https://github.com/FelixKrueger/TrimGalore.
Fu Y, Wu P-H, Beane T, Zamore PD, Weng Z. Elimination of PCR duplicates in RNA-seq and small RNA-seq using unique molecular identifiers. BMC Genomics. 2018;19(1):531.
Ferragina P, Manzini G. Opportunistic data structures with applications. Proceedings 41st Annual Symposium on Foundations of Computer Science. IEEE Comput. Soc; 2000. p. 390–8.
Ferragina P, Manzini G. Indexing compressed text. J ACM (JACM). 2005;52(4):552–81.
Burrows M, Wheeler DJ. A block-sorting lossless compression algorithm. Volume 10. Palo Alto, CA: Digital Equipment Corporation Systems Research Center; 1994.
Adjeroh D, Bell T, Mukherjee A. The Burrows-Wheeler Transform: Data Compression, Suffix Arrays, and Pattern Matching. Boston, MA: Springer US; 2008.
Manber U, Myers G. Suffix Arrays: A New Method for On-Line String Searches. SIAM J Comput. 1993;22(5):935–48.
Wu TD, Reeder J, Lawrence M, Becker G, Brauer MJ. GMAP and GSNAP for genomic sequence alignment: enhancements to speed, accuracy, and functionality. Methods Mol Biol. 2016;1418:283–334.
Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine Learning in Python. J Mach Learn Res. 2011;12:2825–30.
Friedman JH. Greedy function approximation: a gradient boosting machine. Ann Statist. 2001;29(5):1189–232.
Platt J. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Adv large margin classifiers. 1999;10(3):61–74.
Chollet F. Others. Keras [Internet]. 2015 [cited July 2023]. Available from: https://github.com/fchollet/keras.

No competing interests reported.

AlignmentSupplement.v16.docx

Download PDF

Journal Publication

published 12 Mar, 2024

Read the published version in BMC Bioinformatics →

Editorial decision: Revision requested
31 Oct, 2023
Editorial decision: Major revision
31 Oct, 2023
Reviews received at journal
14 Oct, 2023
Reviewers agreed at journal
11 Oct, 2023
Reviewers agreed at journal
04 Oct, 2023
Reviewers invited by journal
14 Aug, 2023
Editor invited by journal
07 Aug, 2023
Editor assigned by journal
07 Aug, 2023
Submission checks completed at journal
03 Aug, 2023
First submitted to journal
28 Jul, 2023

You are reading this latest preprint version

RNA-seq Parent-of-Origin Classification with Machine Learning applied to Alignment Features

Status:

Journal Publication

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Results

Discussion

Conclusion

Methods

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1