# RNA-seq Machine Learning Workflow in Python: Complete Guide

An RNA-seq machine learning Python workflow is a scripted pipeline that converts raw [RNA sequencing](/blog/guides/rna-sequencing) reads into a normalized [gene expression](/blog/guides/gene-expression) matrix and then trains, validates, and interprets a machine learning model on that matrix to predict labels such as cell type, disease state, or treatment response. The same workflow structure covers bulk tissue samples, single-cell data, and spatial [transcriptomics](/knowledge/bioinformatics/modern-transcriptomics-bulk-single-cell-spatial), with Python libraries such as scanpy, scikit-learn, and PyTorch handling the sequence processing, modeling, and interpretation stages.

This workflow matters because RNA-seq data is extreme in shape. A single experiment measures roughly 20,000 protein-coding genes per sample, so the number of features usually exceeds the number of samples by a wide margin. Patterns that predict disease or cell identity are distributed across hundreds of genes and are invisible to any pairwise test a human can run by hand. Machine learning finds those patterns, and Python matters because it keeps every stage, from reading a FASTQ file to explaining a trained model, in one reproducible environment. The difference between a correct pipeline and a leaky one is the difference between a model that generalizes to new patients and a model that only memorized its own training data.

## What Is an RNA-seq Machine Learning Python Workflow?

Three terms make up the name, and each one is worth defining precisely.

**RNA-seq** (ribonucleic acid sequencing) measures how much RNA from each gene is present in a sample, which serves as an estimate of gene activity. Bulk RNA-seq averages expression across a whole tissue. Single-cell RNA-seq (scRNA-seq) measures expression cell by cell, resolving the heterogeneity that bulk data averages away.

**Machine learning (ML)** refers to algorithms that learn a mapping from input features to an output label, or that learn structure in unlabeled data, rather than following hand-written rules. A model trained on expression profiles can predict a class (tumor versus normal, T cell versus B cell) or a continuous value (a cell-type proportion).

**Workflow** means an ordered set of scripted steps with defined inputs and outputs at every stage. In practice this is a set of Python scripts or a pipeline manager configuration that runs the same way every time, on the same inputs, with pinned software versions.

A typical workflow takes FASTQ files (raw reads plus base-quality scores) as input and produces three outputs: a fitted model, an honest estimate of its performance on unseen data, and an explanation of which genes drove the predictions. Common prediction tasks include cell-type annotation, which reviews of cardiovascular single-cell research describe as the step that determines whether downstream biological inference is trustworthy at all [1], disease classification from bulk expression, deconvolution of bulk samples into cell-type proportions [2], and inference of protein activity and master regulators from expression data [3]. The same design pattern appears across omics. MLMarker, a tissue-inference tool for proteomics, pairs a [random forest](/blog/data-analysis/what-is-random-forest) trained on 34 healthy tissues with SHAP-based explanations and a penalty factor for missing measurements, a structure that transfers directly to sparse expression data [4].

Students most often confuse this workflow with three neighboring analyses:

- **Differential expression analysis** asks which genes differ between two predefined groups and tests each gene separately. The machine learning workflow asks which combination of features predicts a label and how well that prediction generalizes. The two are complementary, and most solid projects run both.
- **Functional enrichment analysis** interprets a gene list after an analysis is done. It is a downstream interpretation step, not a modeling step.
- **[Variant calling pipelines](/knowledge/bioinformatics/variant-calling-pipelines-gatk-deepvariant)** work on DNA sequence differences between samples. They share tooling habits with RNA-seq pipelines but operate on different data.

Here is the workflow at a glance.

| Stage | Goal | Common tools | Output |
|--|--|--|--|
| Design and metadata | Define the label, replicates, and batch variables | pandas | sample sheet |
| Read quality control | Find failing samples, adapters, low-quality bases | Biopython, fastp, MultiQC | trimmed FASTQ, QC report |
| Quantification | Convert reads to per-gene counts | STAR, Salmon, kallisto, Cell Ranger | count matrix |
| Normalization | Remove library size and composition bias | scanpy, numpy | normalized, log-transformed matrix |
| Feature selection | Reduce 20,000+ genes to informative features | scanpy, scikit-learn | highly variable genes or principal components |
| Training | Fit the model | scikit-learn, XGBoost, PyTorch | fitted model |
| Validation | Estimate performance on unseen data | scikit-learn | metrics, calibration curves |
| Interpretation | Explain predictions biologically | SHAP, scanpy, pyVIPER | feature rankings, gene signatures |

## Python for Bioinformatics and Where Biopython Fits

Python for bioinformatics earned its position because one interpreter covers the entire arc of an analysis. The scverse ecosystem, built around scanpy and the AnnData object, is the de facto standard for single-cell analysis, and new tools are designed to plug into it. pyVIPER, for example, infers protein activity from single-cell expression data and integrates with scanpy, scverse, and widely adopted machine learning libraries [3]. mmContext stores learned omics embeddings in the AnnData.obsm slot so that any text encoder from Hugging Face can be paired with any numeric representation of an experiment [5]. When your data lives in AnnData, every tool in the ecosystem can read it.

**Biopython** is the sequence-handling layer. It parses FASTA and FASTQ files, queries NCBI databases through the Entrez module, runs BLAST searches, and manipulates alignments and phylogenetic trees. It does not train models. Its job is the sequence end of the workflow, and it does that job well.

This short Biopython tutorial fragment checks read lengths in a FASTQ file, which is a fast sanity check before committing to a long QC run:

```python
from Bio import SeqIO

lengths = [len(record.seq) for record in SeqIO.parse("sample_R1.fastq", "fastq")]
print("min length:", min(lengths))
print("max length:", max(lengths))
```

A second fragment searches NCBI's Sequence Read Archive for public RNA-seq data:

```python
from Bio import Entrez

Entrez.email = "you@institution.edu"
handle = Entrez.esearch(db="sra", term="Homo sapiens[Organism] AND RNA-seq")
print(Entrez.read(handle)["IdList"][:5])
```

For compute-heavy steps such as adapter trimming and quantification, the convention is to call dedicated command-line tools (fastp, Salmon, STAR) from scripts or a pipeline manager, and to use Python for everything around them.

## The RNA-seq Machine Learning Python Workflow Explained Step by Step

The nine steps below describe the standard arc. Bulk and single-cell variants differ mainly in steps 2 through 4, which the next sections cover.

### Step 1: Define the Question and Lock the Metadata

Write down the prediction target before touching data. "Classify tumor versus normal" and "predict the proportion of fibroblasts in a bulk sample" lead to different designs. Record the unit of replication. If the label belongs to a patient, the split between training and test data must be made at the patient level, never at the cell level, or the model will learn each patient's background rather than the biology. Build a sample sheet in CSV with sample identifier, condition, batch, library preparation, and FASTQ file paths. With small cohorts, expect wide uncertainty in performance estimates and say so in the analysis plan.

### Step 2: Quality Control on the Raw Reads

A FASTQ record contains the [nucleotide sequence](/knowledge/molecular-biology/nucleotide-sequence) and a Phred quality score per base, where Q30 means a 1 in 1,000 chance that the base call is wrong. QC examines per-base quality profiles, adapter contamination, and duplication. fastp handles trimming and reporting in one pass:

```bash
fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trimmed_R1.fastq.gz -O trimmed_R2.fastq.gz \
  -h fastp_report.html -w 8
```

Aggregate per-sample reports with MultiQC so one page reveals an outlier library. Single-cell data adds two QC steps that annotation reviews list as standard: filtering empty droplets that capture ambient RNA, and detecting doublets, which are droplets containing two cells [1].

### Step 3: Quantification, From Reads to a Count Matrix

Two routes convert reads into counts. Alignment-based quantification maps reads to a reference genome with STAR or HISAT2 and then assigns reads to genes with featureCounts or HTSeq. Pseudoalignment and lightweight quantification with Salmon or kallisto map reads to a transcriptome index directly and run much faster on ordinary hardware:

```bash
salmon quant -i grch38_salmon_index -l A \
  -1 trimmed_R1.fastq.gz -2 trimmed_R2.fastq.gz \
  -p 8 -o salmon_out/sample01
```

For single-cell data, Cell Ranger, STARsolo, or alevin-fry produce a genes-by-cells count matrix that loads directly into an AnnData object. Whatever route you take, record the reference genome version and annotation release, because gene identifiers change between releases.

### Step 4: Normalization and Transformation

Raw counts are not comparable across samples because [sequencing depth](/knowledge/diagnostics/molecular/how-to-calculate-sequencing-depth-and-coverage-for-your-ngs-run) differs by orders of magnitude, and highly expressed genes distort the counts of everything else (composition bias). Bulk workflows typically use size-factor normalization (the median-of-ratios method popularized by DESeq2) or trimmed mean of M values, or convert counts to TPM (transcripts per million, which sums to one million per sample) followed by a log transform such as log2(TPM + 1). Single-cell workflows in scanpy scale each cell to a fixed total, commonly 10,000 counts, then apply log1p:

```python
import scanpy as sc

adata = sc.read_h5ad("counts.h5ad")
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
```

Batch effects, which are technical differences between processing groups, are handled with methods such as Harmony, ComBat, or scVI. No algorithm rescues a design in which every treated sample was processed in one batch and every control in another, because batch and condition are then mathematically confounded. Check the design table before modeling, not after.

### Step 5: Feature Selection and Dimensionality Reduction

Genes with near-zero expression in almost every sample carry noise, not signal. The standard reduction is to keep a few thousand highly variable genes (HVGs, the genes with the greatest biologically meaningful variance), scale them, and summarize them with [principal component analysis](/blog/data-analysis/principal-component-analysis-explained) (PCA), which compresses correlated genes into 10 to 50 orthogonal components:

```python
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)
sc.pp.scale(adata, max_value=10)
sc.pp.pca(adata, n_comps=30)
```

UMAP and t-SNE produce two-dimensional plots that are excellent for seeing structure and useless as model features, because they distort global distances. Feed PCA coordinates or normalized expression to the model, and use UMAP only for figures. One rule governs this entire step: any selection that looks at the labels or at all samples at once must happen inside cross-validation folds, which is the leakage guard explained in step 7.

### Step 6: Model Training

Start with models that are hard to beat on tabular expression data. Regularized [logistic regression](/blog/data-analysis/logistic-regression-definition-formula) (including elastic net, which blends L1 and L2 penalties to handle correlated genes), random forests, and gradient boosting are the workhorses. One clinical-genomic study built random forest, logistic regression, and XGBoost models in Python 3.9.0 on 836 pneumonia cases with genomic and clinical features, split 7 to 3 into training and validation sets, and compared them on accuracy, precision, recall, F1 score, and ROC-AUC [6]. That comparison-first habit is the template.

Three practical problems dominate training:

- **Class imbalance.** Rare cell types are the classic case. sc-SynO addresses it with synthetic oversampling based on the LoRAS algorithm, generating synthetic minority-class profiles from real ones, and demonstrated the approach on cardiac glial cells that made up 17 nuclei out of 8,635, an imbalance of roughly 1 to 500 [7]. Simpler tools first: class_weight="balanced" in scikit-learn, then resampling if needed.
- **Very large datasets.** Training on millions of cells is wasteful when a well-chosen subsample carries the same signal. scValue ranks cells by data value using out-of-bag estimates from a random forest, gives more representation to cell types with high variability, and matched full-data performance on annotation tasks across four large datasets [8].
- **Borrowing strength from other cohorts.** Transfer learning formalizes this. LEARNER, a statistical transfer learning method for low-rank matrix estimation, penalizes differences between the latent row and column spaces of a source and target population and uses cross-validation to decide how much to borrow, often beating target-only estimation as source signal quality rises [9].

A minimal training loop with honest comparison:

```python
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_validate
from sklearn.pipeline import make_pipeline
from sklearn.[preprocessing](/blog/data-analysis/what-is-data-preprocessing) import StandardScaler

X = adata.obsm["X_pca"]
y = adata.obs["cell_type"].astype(str)

models = {
    "elastic_logreg": make_pipeline(
        StandardScaler(),
        LogisticRegression(penalty="elasticnet", l1_ratio=0.5,
                           solver="saga", max_iter=5000),
    ),
    "random_forest": RandomForestClassifier(
        n_estimators=500, class_weight="balanced",
        n_jobs=-1, random_state=0),
}

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
for name, model in models.items():
    scores = cross_validate(model, X, y, cv=cv,
                            scoring=["f1_macro", "balanced_accuracy"])
    print(name, round(scores["test_f1_macro"].mean(), 3))
```

Deep learning earns its complexity when data is abundant or when the model must learn representations, as in the autoencoder and multimodal models described later. For a matrix of a few hundred samples, tuned classical models are usually the right first answer.

### Step 7: Validation

Validation is where most student projects fail, so treat it as a design decision made in step 1, not a report written in step 9.

- **Split correctly.** Stratified k-fold cross-validation keeps class proportions equal in every fold, which is why the SNP-based schizophrenia study used stratified five-fold cross-validation across its 4,693 participants and reported results separately for each ethnicity and gender group [10]. For patient-level labels, use grouped folds so no patient appears in both training and test data.
- **Report the right metrics.** Accuracy collapses when classes are imbalanced. Report precision, recall, F1, and ROC-AUC (the probability that a random positive case ranks above a random negative one), and report sensitivity and specificity separately. The schizophrenia models illustrate why: several groups showed high sensitivity with low specificity, a pattern that overall accuracy would have hidden [10].
- **Check calibration.** For any model that outputs a probability meant to guide decisions, calibration curves show whether predicted probabilities match observed frequencies, and decision curve analysis evaluates net benefit across probability thresholds. The pneumonia study used both alongside its discrimination metrics.
- **Validate externally.** Cross-validation estimates performance on data from the same cohort. An independent cohort with its own batch effects is the real test.
- **Tune honestly.** Nested cross-validation, with an inner loop for hyperparameter tuning and an outer loop for performance, gives an unbiased estimate. Tuning on the same folds you report is a quiet form of leakage.

### Step 8: Interpretation

A model that cannot be inspected cannot be trusted or published usefully. Linear models expose coefficients directly. Tree ensembles support permutation importance, which measures how much a metric drops when one feature is shuffled. SHAP (SHapley Additive exPlanations) attributes each prediction to features using game-theoretic fair-division rules, and MLMarker demonstrates the design at production quality, pairing its random forest with SHAP-based explanations at the protein level [4].

Map top features back to biology rather than stopping at the ranking. Gene loadings from PCA connect components back to genes through adata.varm["PCs"]. For single-cell data, gene regulatory network (GRN) analysis moves the interpretation from correlated genes to regulatory proteins: pyVIPER infers protein activity and identifies master regulators from single-cell expression data inside the scanpy ecosystem, turning a cell-type classifier's confusion into a hypothesis about which transcription factors drive each state [3]. One caution applies to every method in this section: correlated features split importance between themselves unpredictably, so treat a ranking as a shortlist for experiments, never as a causal statement.

### Step 9: Scale, Automate, and Reproduce

Three forces shape the scaling end of the workflow. GPU acceleration turns intractable analyses routine, and pyVIPER's PyTorch-based implementation cut large single-cell analyses from hours to minutes, orders of magnitude faster than the original R implementation [3]. Subsampling keeps big-data training affordable without giving up accuracy, as the scValue results show [8]. Automation keeps the whole thing honest: AESTETIK, a spatial transcriptomics model, ships a Snakemake pipeline so its published results can be regenerated exactly [11], and Agentomics, an autonomous agentic system for biomedical machine learning, runs end-to-end experimentation with strict validation checkpoints between steps, evaluated across 20 datasets in [protein engineering](/knowledge/molecular-biology/protein-engineering), drug discovery, and regulatory genomics [12].

Reproducibility is a set of habits, not a tool: pin every package version in a conda environment file, set every random seed explicitly, log the reference genome and annotation release, and commit the sample sheet alongside the code.

## Worked Example: Cell-Type Classification With scanpy and scikit-learn

This script is one of the simplest RNA-seq machine learning Python workflow examples that runs end to end on a laptop. It assumes an AnnData file with raw counts in the .X slot and cell labels in adata.obs["cell_type"].

```python
import scanpy as sc
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix, f1_score
from sklearn.model_selection import train_test_split

# 1. Load, filter, and normalize
adata = sc.read_h5ad("pbmc_counts.h5ad")
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs["cell_type"].notna()].copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 2. Features
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)
sc.pp.scale(adata, max_value=10)
sc.pp.pca(adata, n_comps=30)
X = adata.obsm["X_pca"]
y = adata.obs["cell_type"].astype(str).values

# 3. Split. For patient-level labels, split on patient ID instead.
X_tr, X_te, y_tr, y_te = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=0)

# 4. Fit and evaluate
clf = RandomForestClassifier(n_estimators=500,
                             class_weight="balanced",
                             n_jobs=-1, random_state=0)
clf.fit(X_tr, y_tr)
pred = clf.predict(X_te)
print("macro F1:", round(f1_score(y_te, pred, average="macro"), 3))
print(confusion_matrix(y_te, pred))
```

Reading the output correctly matters more than the code. The [confusion matrix](/blog/data-analysis/confusion-matrix-definition-examples) shows which cell types the model merges, and those merges usually reflect real biological similarity or underrepresentation of a class. To interpret the model, compute permutation importance on the test fold, then trace important principal components back to genes through adata.varm["PCs"] and check the top genes against known markers. One honest caveat about this script: normalization and feature selection here run on the full object for readability. In a rigorous analysis, wrap scaling, HVG selection, and PCA in a scikit-learn Pipeline or refit them within training folds, so no information from the test cells shapes the features.

## Bulk, Single-Cell, and Spatial Versions of the Workflow

The skeleton is identical across data types. The differences concentrate in preprocessing and in what the label means.

**Bulk RNA-seq** works at the sample level, with labels such as diagnosis, subtype, or outcome. Because bulk expression is a mixture of cell types, a growing family of models deconvolves it, meaning they estimate the proportion of each cell type from the mixture. ACSCeND combines a supervised classifier that assigns single cells to stem-like states with an attention-guided autoencoder that deconvolves those states from bulk RNA-seq, and applied the framework to more than 25,000 tumor profiles from TCGA, PRECOG, relapse, and checkpoint inhibitor studies [2].

**Single-cell RNA-seq** works at the cell level for annotation tasks and at the sample level for patient-level tasks. The standard annotation pipeline described in recent reviews runs empty droplet filtering, doublet detection, normalization, and batch correction, then dimensionality reduction and graph-based clustering, then label assignment through marker-based curation, reference-correlation methods such as SingleR, or supervised classifiers such as CellTypist [1]. Rare populations make annotation an imbalanced classification problem, which is where synthetic oversampling methods such as sc-SynO enter [7]. Protein activity methods such as pyVIPER add a regulatory layer on top of the same AnnData objects [3].

**Spatial transcriptomics** preserves tissue location and adds a second data modality, histology. AESTETIK learns spot representations by jointly modeling spatial coordinates, expression, and tissue morphology, and reported a 21 percent increase in median adjusted Rand index (ARI, a clustering agreement score where 1 is a perfect match) on brain datasets, a 2-fold increase in breast cancer, and 79 percent in melanoma over previous methods [11].

**Multimodal extensions** pair expression with text or other data types. mmContext trains joint embeddings of omics and biological text on top of the Sentence Transformers library, accepts any numeric representation stored in AnnData.obsm, and demonstrated zero-shot classification (predicting labels the model never saw as labeled examples) of cell types and diseases across four independent datasets [5].

## Common Mistakes and Limitations

1. **Leakage through preprocessing.** Normalization parameters, HVG selection, PCA rotation, and feature rankings computed on the full dataset before splitting leak test information into training. Fix: put every learned transform inside the cross-validation loop, and keep a truly untouched external test set when the project warrants one.
2. **Cell-level splits for patient-level labels.** Cells from one patient are near-duplicates. Splitting by cell inflates accuracy that evaporates on new patients. Fix: GroupKFold or manual splits on patient or sample identifiers.
3. **Confounded batch and condition.** If all treated samples ran on one day and all controls on another, no algorithm can separate treatment from batch. Fix: randomize processing at the bench and inspect the design table before modeling.
4. **Using UMAP or t-SNE coordinates as features.** These methods preserve local neighborhoods and distort global geometry. Fix: model on PCA coordinates or normalized expression, and reserve embeddings for plots.
5. **Feeding raw counts to linear or distance-based models.** Sequencing depth differences dominate everything else. Fix: normalize and log-transform first, every time.
6. **Ignoring class imbalance.** With a 1-to-500 imbalance, a model that predicts the majority class every time scores high accuracy and zero usefulness. Fix: balanced metrics, class weights, and oversampling methods built for the setting [7].
7. **Ignoring missingness and platform shift.** Dropout in single-cell data and feature gaps across platforms behave like systematic missing data. Multi-output imputation models have been validated with up to 50 percent induced missingness and improved downstream age-prediction performance in methylation data, a result that doubles as a warning about how much missingness expression models can silently absorb [13].
8. **Reading importance as causation.** Feature importance ranks associations within one dataset and cohort. Experimental validation is a separate project.
9. **Stopping at cross-validation.** Good cross-validated numbers on one cohort say nothing about a hospital, a platform, or a population the model has never seen. Fix: external validation with calibration checks [6].
10. **Reproducibility gaps.** Unpinned versions, unset seeds, and untracked reference genomes make results unrepeatable, which reviewers and collaborators treat as a defect.

The ceiling on any of these models is the quality of its labels. Mis-annotated cell types propagate systematic bias into every downstream inference, which is why annotation reviews treat label accuracy as the [limiting factor](/blog/guides/limiting-factor-definition-types-and-examples) for the entire field [1]. Individual datasets also need expert biological review, and model predictions should be treated as hypotheses that require experimental follow-up, not conclusions.

### Misconceptions Students Often Hold

- **"More genes make a better model."** Past a few thousand features, noise genes dilute signal and slow training. Careful HVG selection usually improves both performance and interpretability.
- **"Deep learning always beats classical machine learning."** On expression matrices with modest sample counts, tuned random forests and gradient boosting remain competitive, faster, and easier to explain, which is why they anchor so many published omics tools [8], [4], [6].
- **"Machine learning replaces differential expression analysis."** The two answer different questions. A classifier says which features predict a label, a differential expression test says which genes differ between defined groups, and solid papers report both.
- **"Good cross-validation numbers mean the model is ready."** Cohort shift, batch effects, and poor calibration all appear only outside the training cohort [6].
- **"Normalization is optional because the model is flexible."** Depth and composition effects are technical, and a flexible model will happily learn them as if they were biology.

## Quick Review

- The pipeline order is fixed: metadata, read QC, quantification, normalization, feature selection, training, validation, interpretation.
- Split data by the unit of replication. Patient-level labels require patient-level splits.
- Every data-dependent transform, including scaling, HVG selection, and PCA, belongs inside the training folds.
- Single-cell normalization is per-cell scaling to a fixed total followed by log1p. Bulk normalization uses size factors or TPM plus a log transform.
- Report precision, recall, F1, ROC-AUC, sensitivity, and specificity, plus calibration for any probability output. Accuracy alone hides imbalance.
- UMAP and t-SNE are for figures. PCA coordinates are for models.
- Pin versions, set seeds, record the reference genome, and script every step from FASTQ to final figure.

## Frequently Asked Questions

### What is an RNA-seq machine learning Python workflow?

It is a scripted pipeline that converts raw RNA sequencing reads into a normalized gene expression matrix and then trains, validates, and interprets a machine learning model on that matrix to predict labels such as cell type, disease state, or treatment response.

### Which Python libraries does the workflow need?

scanpy with AnnData for single-cell data, scikit-learn for classical models, PyTorch for deep learning, Biopython for sequence files, and pandas for metadata. Quantification usually runs in Salmon, STAR, or Cell Ranger outside Python, called from scripts.

### Can I skip alignment and start from a count matrix?

Yes. Most machine learning work starts from a count matrix, provided you verify how it was generated, which reference genome and annotation were used, and that quality control was actually performed on the reads.

### How do I avoid data leakage in this workflow?

Split the data by patient or sample before anything else, fit every learned transform such as scaling, feature selection, and PCA inside the training folds only, tune with nested cross-validation, and keep an external test set when possible.

### How do I handle rare cell types or other rare classes?

![Calculation showing 17 cardiac glial nuclei out of 8,635 total nuclei gives a minority share of about 0.2 percent, an imbalance near 1 to 500.](/article-images/illustrations/rna-seq-machine-learning-python-workflow-imbalance-ratio.webp)

Use balanced metrics and per-class recall instead of accuracy, apply class weights, and consider resampling, including synthetic oversampling methods built for single-cell data and value-based subsampling for very large datasets.

### Does the same workflow apply to bulk, single-cell, and spatial data?

The skeleton is the same, while preprocessing, normalization, and label structure differ by data type. Spatial data adds image features and requires methods that model tissue location alongside expression.

<script type="application/ld+json">
{
  "@context": "https://schema.org",
  "@type": "FAQPage",
  "mainEntity": [
    {
      "@type": "Question",
      "name": "What is an RNA-seq machine learning Python workflow?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "It is a scripted pipeline that converts raw RNA sequencing reads into a normalized gene expression matrix and then trains, validates, and interprets a machine learning model on that matrix to predict labels such as cell type, disease state, or treatment response."
      }
    },
    {
      "@type": "Question",
      "name": "Which Python libraries does the workflow need?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "scanpy with AnnData for single-cell data, scikit-learn for classical models, PyTorch for deep learning, Biopython for sequence files, and pandas for metadata. Quantification usually runs in Salmon, STAR, or Cell Ranger outside Python, called from scripts."
      }
    },
    {
      "@type": "Question",
      "name": "Can I skip alignment and start from a count matrix?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Yes. Most machine learning work starts from a count matrix, provided you verify how it was generated, which reference genome and annotation were used, and that quality control was actually performed on the reads."
      }
    },
    {
      "@type": "Question",
      "name": "How do I avoid data leakage in this workflow?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Split the data by patient or sample before anything else, fit every learned transform such as scaling, feature selection, and PCA inside the training folds only, tune with nested cross-validation, and keep an external test set when possible."
      }
    },
    {
      "@type": "Question",
      "name": "How do I handle rare cell types or other rare classes?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Use balanced metrics and per-class recall instead of accuracy, apply class weights, and consider resampling, including synthetic oversampling methods built for single-cell data and value-based subsampling for very large datasets."
      }
    },
    {
      "@type": "Question",
      "name": "Does the same workflow apply to bulk, single-cell, and spatial data?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "The skeleton is the same, while preprocessing, normalization, and label structure differ by data type. Spatial data adds image features and requires methods that model tissue location alongside expression."
      }
    }
  ]
}
</script>

## Sources

1. [Annotation of cell types in single-cell sequencing for cardiovascular disease: Concepts, workflows, challenges, and best practices.](https://pubmed.ncbi.nlm.nih.gov/42716373/)
2. [Cross-modal mapping of cancer stem-like cell plasticity using deep learning.](https://pubmed.ncbi.nlm.nih.gov/42516214/)
3. [pyVIPER: a fast and scalable Python package for protein activity estimation and master regulator analysis of single-cell RNA sequencing data.](https://pubmed.ncbi.nlm.nih.gov/42332576/)
4. [MLMarker: a machine learning framework for tissue inference and biomarker discovery.](https://pubmed.ncbi.nlm.nih.gov/42343371/)
5. [mmContext: an open framework for multimodal contrastive learning of omics and text data.](https://pubmed.ncbi.nlm.nih.gov/42191659/)
6. [Machine learning prediction of antibiotic resistance in Gram-negative pneumonia: a clinical-genomic integrated model.](https://pubmed.ncbi.nlm.nih.gov/41957589/)
7. [Automated annotation of rare-cell types from single-cell RNA-sequencing data through synthetic oversampling.](https://pubmed.ncbi.nlm.nih.gov/34798805/)
8. [scValue: value-based subsampling of large-scale single-cell transcriptomic data for machine and deep learning tasks.](https://pubmed.ncbi.nlm.nih.gov/40515392/)
9. [LEARNER: A Transfer Learning Method for Low-Rank Matrix Estimation.](https://pubmed.ncbi.nlm.nih.gov/42438120/)
10. [SNP-based prediction of schizophrenia using machine learning.](https://pubmed.ncbi.nlm.nih.gov/41695265/)
11. [Representation learning for multi-modal spatially resolved transcriptomics data.](https://pubmed.ncbi.nlm.nih.gov/42179160/)
12. [Agentomics: an agentic system that autonomously develops novel state-of-the-art solutions for biomedical machine learning tasks.](https://pubmed.ncbi.nlm.nih.gov/42412783/)
13. [Multi-output learning for systematic missing value imputation in DNA methylation arrays.](https://pubmed.ncbi.nlm.nih.gov/41782682/)