# How to Perform Differential Expression Analysis in Python By the end of this guide you will be able to take a raw RNA-seq count matrix, run a full differential expression analysis in Python, produce a statistically defensible results table with adjusted p-values and shrunk log2 fold changes, and generate the standard QC and interpretation plots (sample clustering, MA plot, volcano plot, heatmap). If your data are single-cell, you will also learn the pseudobulk strategy that brings cell-level data into the same statistical framework. Everything runs on a laptop or a free Google Colab session, and the whole pipeline stays in Python from start to finish. Before you start, have the following on hand: - Python 3.10 or newer, ideally through conda or mamba - A gene-by-sample count matrix (raw integer counts, not TPM or FPKM) - A sample metadata table with at least one biological condition and, ideally, batch information - Roughly 8 GB of RAM for datasets up to a few hundred samples ## Why Python Works for Differential Expression Python for bioinformatics has matured to the point where the heavy lifting of [transcriptomics](/knowledge/bioinformatics/modern-transcriptomics-bulk-single-cell-spatial) now happens natively in Python. The scverse ecosystem (AnnData, Scanpy, PyDESeq2) gives you data containers, single-cell tooling, and statistical testing that interoperate cleanly. PyDESeq2 is a Python reimplementation of the widely used DESeq2 workflow for bulk [RNA-seq differential expression](/knowledge/bioinformatics/rna-seq-differential-expression-deseq2-edger), released under the MIT license and documented at pydeseq2.readthedocs.io. It reaches model likelihoods similar to or higher than the original R implementation, scales better on large datasets such as TCGA cohorts, and interfaces directly with pandas DataFrames and modern Python data science tools [1]. That said, differential expression analysis is a statistics problem wearing a sequencing costume. The computational steps of a typical RNA-seq workflow are read trimming, alignment, quantification, normalization, and differential expression testing [2]. Python covers all five stages, but the statistical core only behaves correctly if you feed it the right input. The single most common error in the entire field is handing a statistical model normalized values or TPM instead of raw counts. This guide is structured so that mistake cannot happen. ## The Differential Expression Analysis Python Workflow at a Glance The end-to-end differential expression analysis python workflow looks like this: 1. Set up a reproducible environment (conda environment, pinned versions) 2. Move from FASTQ files to a gene-by-sample count matrix 3. Load counts and metadata into pandas or AnnData 4. Filter genes with near-zero counts 5. Normalize (size factors) and fit a negative binomial model with PyDESeq2 6. Extract results, correct for multiple testing, shrink fold changes 7. Visualize and sanity-check 8. Run functional enrichment on the gene list Steps 2 through 4 are data engineering. Steps 5 and 6 are the statistics. Steps 7 and 8 are where biology comes back in. Most of your debugging time will go into steps 2 and 3, and most of your scientific judgment goes into steps 6 and 8. | Step | Tool | Input | Output | Key rule | |--|--|--|--|--| | QC and trimming | FastQC, Trimmomatic | FASTQ | clean FASTQ | never trim after alignment | | Alignment or pseudoalignment | HISAT2, Subread, Salmon | FASTQ + reference | BAM or transcript abundances | keep the genome annotation version fixed | | Counting | htseq-count, featureCounts | BAM + GTF | raw count matrix | counts must be integers | | Transcript-to-gene | pytximport | transcript-level estimates | gene-level counts [3] | use the same transcript-to-gene map for every sample | | Filtering | PyDESeq2, scanpy | count matrix | filtered matrix | filter before normalization | | Normalization and testing | PyDESeq2 | raw counts | size factors, DE results [1] | never feed TPM, FPKM, or log-transformed values | | Enrichment | g:Profiler | gene list | enriched terms [4] | use the ranked or thresholded list consistently | ## Step 1: Set Up a Reproducible Environment Create a dedicated environment. Version drift between PyDESeq2, scanpy, and pandas is a real source of unreproducible results. ```bash mamba create -n dea -c conda-forge -c bioconda \ python=3.11 pydeseq2 scanpy pytximport htseq \ hisat2 samtools subread fastqc jupyterlab mamba activate dea ``` If you have no local computing resources, the full workflow, including the alignment and counting steps, can be run on free cloud platforms such as Google Colab using SRA Toolkit, FastQC, Trimmomatic, HISAT2, samtools, and Subread, with Python handling the downstream analysis [4]. Verify the install: ```bash python -c "import pydeseq2; print(pydeseq2.__version__)" ``` Expected output: ``` 0.4.4 ``` Pin whatever version prints in your environment file so collaborators get identical numbers. ## Step 2: From FASTQ Files to a Count Matrix ![Flowchart splitting FASTQ input into an alignment and counting route or a Salmon and pytximport route, both ending in an integer count matrix.](/article-images/illustrations/differential-expression-analysis-python-workflow-counts-routes.webp) PyDESeq2, like DESeq2 before it, models raw counts with a [negative binomial distribution](/blog/data-analysis/negative-binomial-distribution-formula-examples). Two routes produce those counts. ### Route A: Alignment and gene-level counting The classical route is splice-aware alignment followed by counting reads that overlap exons. A per-sample loop looks like this: ```bash hisat2 -x GRCh38_genome -1 sample1_R1.fastq.gz -2 sample1_R2.fastq.gz \ | samtools sort -o sample1.bam samtools index sample1.bam ``` Then count reads per gene with htseq-count, which preprocesses RNA-seq data for differential expression analysis by counting the overlap of reads with genes, built on the HTSeq Python framework [5]: ```bash htseq-count -f bam -r pos -t exon -i gene_id \ sample1.bam GRCh38.gtf > sample1_counts.tsv ``` The two flags that matter most: - `-t exon` tells the counter which annotation feature type defines a gene's extent - `-i gene_id` names the GTF attribute used as the gene identifier Run the identical command for every sample, then concatenate the per-sample files into one matrix in Python: ```python import pandas as pd samples = ["sample1", "sample2", "sample3", "sample4"] counts = pd.concat( [pd.read_csv(f"{s}_counts.tsv", sep="\t", header=None, names=["gene", s]) for s in samples] ).groupby("gene").sum() counts.to_csv("counts_matrix.csv") ``` Expected shape check: ```python counts.shape ``` ``` (60746, 4) ``` The exact row count depends on your annotation. Around 60,000 rows for a full GENCODE human release is normal, and most will be filtered in step 4. ### Route B: Transcript quantification followed by pytximport Modern pipelines often quantify transcripts directly with tools such as Salmon rather than aligning. The output is transcript-level estimated counts, which statistical DE tools cannot consume directly. pytximport solves this by providing a Python implementation of the tximport approach, supporting multiple input formats, bias correction modes, inferential replicates, gene-level summarization, and optional filtering by transcript biotype [3]. In Python: ```python from pytximport import tximport txi = tximport( file_paths=[f"salmon/{s}/quant.sf" for s in samples], sample_names=samples, transcript_to_gene_mapping="gencode_v44_tx2gene.tsv", counts_from_abundance="length_scaled_tpm", output_format="wide", ) counts = txi.to_df() counts.round().astype(int).to_csv("counts_matrix_tximport.csv") ``` Note the cast to integers at the end. PyDESeq2 expects integer-like counts, and estimated counts arrive as floats. Rounding estimated counts to integers for use with the DESeq2-style workflow is the standard practice that tximport itself recommends. ## Step 3: Load Counts and Metadata ![Decision path checking whether counts are raw integers, sample names match metadata, and batch is in the design before testing.](/article-images/illustrations/differential-expression-analysis-python-workflow-raw-counts-check.webp) Keep two artifacts strictly separate: the count matrix (genes as rows, samples as columns) and the metadata table (samples as rows, covariates as columns). Sample names must match exactly between the two. ```python import pandas as pd counts = pd.read_csv("counts_matrix.csv", index_col=0) metadata = pd.read_csv("metadata.csv", index_col=0) ``` A minimal metadata file: ``` condition batch sample1 control A sample2 control A sample3 treated B sample4 treated B ``` Run the first integrity check now, not later: ```python assert set(counts.columns) == set(metadata.index), "Sample name mismatch" assert counts.dtypes.map(lambda x: pd.api.types.is_integer_dtype(x)).all(), \ "Counts are not integers" ``` If the assertion fails, you have either a whitespace problem in one file or float counts from a quantification tool. Fix it here, because a silent sample-order mismatch is the error most likely to survive into a published figure. A concrete real-world case shows the stakes of a well-organized multi-tissue design: the Burmese python digestion study built a five-organ, three-time-point transcriptome and used it to find condition-specific gene sets in heart, liver, stomach, pancreas, and intestine [6], which is only possible when sample metadata is unambiguous. ## Step 4: Filter Low-Count Genes Genes with almost no counts carry no information, cost degrees of freedom, and hurt multiple-testing correction. The DESeq2 convention is to keep genes with a minimum number of counts across samples. In PyDESeq2 you can filter before constructing the dataset: ```python genes_to_keep = counts.columns[counts.sum(axis=0) >= 10] # Note: PyDESeq2 expects genes as columns and samples as rows, # or you can transpose here for consistency with the API ``` A common, defensible rule is to require roughly 10 counts total across all samples, or counts in at least a minimum number of samples. Being more aggressive here means fewer tests and therefore more power in step 6, at the cost of losing genes expressed in a tiny fraction of cells or samples. Document the exact rule you use, because reviewers will ask. ## Step 5: Normalize and Fit the Model with PyDESeq2 This is the heart of the differential expression analysis python workflow. PyDESeq2 implements the DESeq2 pipeline in three phases: median-of-ratios size-factor normalization, dispersion estimation, and a negative binomial generalized linear model per gene [1]. ```python from pydeseq2.dds import DeseqDataSet from pydeseq2.ds import DeseqStats counts_t = counts.T # samples as rows, genes as columns counts_t = counts_t.loc[metadata.index] # enforce matching order dds = DeseqDataSet( counts=counts_t, metadata=metadata, design="~ batch + condition", refit_cooks=True, n_cpus=8, ) dds.deseq2() ``` Parameters that actually matter: - `design` is an R-style formula. `~ batch + condition` adjusts the model for batch while testing condition. Put nuisance variables (batch, sex, sequencing lane) before the variable of interest. If you have paired samples (for example, the same patient before and after treatment), include the patient identifier: `~ patient + condition`. - `refit_cooks=True` triggers Cooks distance-based outlier refitting. With fewer than 7 replicates per group, outlier detection is disabled by the underlying method, so do not count on it for small designs. - `n_cpus` controls parallelism. For large cohorts such as TCGA, PyDESeq2's speed advantage over the original implementation becomes substantial [1]. Expected console output during `dds.deseq2()`: ``` Fitting size factors... ... done in 1.20 seconds. Fitting dispersions... ... done in 24.31 seconds. Fitting dispersions MAP... ... done in 19.87 seconds. Fitting LFCs... ... done in 22.14 seconds. Replacing outliers... ... done in 4.10 seconds. ``` A critical detail: `DeseqDataSet` expects samples as rows and genes as columns, which is transposed relative to how most count matrices are stored. The `.T` transpose above prevents the classic column-name error. If you work outside this framework, normalization is still the non-negotiable step. Alternative methods exist for other data types. CONSTANd, for example, is a fast normalization method for relative quantification in omics experiments, with implementations in both R and Python, designed to fit automated pipelines of any scale [7]. Whatever you choose, normalize once, before testing, and never mix normalized and raw values in the same design matrix. ## Step 6: Extract Results, Correct, and Shrink ```python stat_res = DeseqStats( dds, contrast=("condition", "treated", "control"), alpha=0.05, n_cpus=8, ) stat_res.summary() stat_res.lfc_shrink(coeffs=["condition[T.treated]"]) res = stat_res.results_df ``` What each piece does: - `contrast` sets which two factor levels are compared and in which direction. Here positive log2 fold change means higher in treated. - `alpha=0.05` sets the target for independent filtering and the Benjamini-Hochberg false discovery rate (FDR) correction, which adjusts p-values so that, among genes you call significant, the expected fraction of false positives is at most 5 percent. - `lfc_shrink` applies fold-change shrinkage to reduce the inflation of log2 fold changes for low-count genes. Shrunken estimates are what you should use for ranking, visualization, and effect-size claims. Expected results table: ``` baseMean log2FoldChange lfcSE stat pvalue padj gene_A 842.1 2.310 0.191 12.09 1.1e-33 3.4e-30 gene_B 1290.7 -1.104 0.098 -11.27 1.7e-29 2.6e-26 gene_C 3.2 5.119 1.302 3.93 8.4e-05 0.081 ``` Read it as follows. `baseMean` is average normalized expression. `log2FoldChange` is the effect size. `lfcSE` is its standard error. `stat` is the Wald test statistic. `pvalue` is uncorrected. `padj` is the FDR-adjusted value you filter on. Gene C illustrates why: it has an impressive raw p-value but only 3 counts on average, so the adjusted value is not significant. Never call a gene DE based on `pvalue` alone. ## Step 7: Visualize and Sanity-Check Before believing any result table, check three plots. **Sample-level QC with PCA.** Compute the normalized and regularized log transform, then project samples: ```python dds.layers["normalized_counts"] dds.obs["condition"] = metadata["condition"] import scanpy as sc dds.layers["log1p"] = sc.pp.normalize_total # see note below ``` In practice, the cleanest path is to use the PyDESeq2 built-in: ```python dds.deseq2() # if not already run from pydeseq2.dds import DeseqDataSet dds_vst = DeseqDataSet( counts=counts_t, metadata=metadata, design="~ batch + condition" ) dds_vst.deseq2() dds_vst.obsm["normalized_counts"] import numpy as np log_counts = np.log2(counts_t[keep] + 1) ``` Then run PCA on the log-transformed normalized matrix and color points by condition and by batch. Samples should separate by condition if you expect a strong treatment effect, and must not separate by batch. If batch dominates, the batch term in your design becomes essential, and if batch is confounded with condition (all controls in lane A, all treated in lane B), no statistical method can fully rescue the design. **MA plot.** Mean expression versus log2 fold change. Shrunken fold changes should show a funnel shape, with variance decreasing as expression increases. A tilted or banana-shaped cloud signals a normalization problem. **Volcano plot.** ```python import matplotlib.pyplot as plt sig = res[(res.padj < 0.05) & (res.log2FoldChange.abs() > 1)] plt.scatter(res.log2FoldChange, -np.log10(res.padj), s=3, alpha=0.4) plt.scatter(sig.log2FoldChange, -np.log10(sig.padj), s=3, color="crimson") plt.axhline(-np.log10(0.05), ls="dashed", lw=0.7) plt.axvline(1, ls="dashed", lw=0.7) plt.axvline(-1, ls="dashed", lw=0.7) plt.xlabel("log2 fold change") plt.ylabel("-log10 adjusted p-value") plt.savefig("volcano.png", dpi=300) ``` Final sanity checks on the table itself: - Known positive-control genes for your perturbation should appear among the top hits. For a heat-shock experiment, heat-shock protein genes should rank highly. - The number of significant genes should be plausible. Zero significant genes suggests low power or a design problem. Ten thousand suggests a threshold or contrast error. - Check the direction. A common sign error is a swapped contrast producing the exact mirror image of the biology you expect. ## Step 8: Functional Enrichment A ranked list of 300 genes means little until it is mapped to pathways. The g:Profiler service is a standard choice and is what published minimal cloud-based RNA-seq protocols pair with pyDESeq2 for functional enrichment [4]. A typical Python call: ```python from gprofiler import GProfiler gp = GProfiler(return_dataframe=True) enrich = gp.profile( organism="hsapiens", query=list(sig.index), sources=["GO:BP", "KEGG", "REAC"], ) print(enrich[["name", "p_value", "intersection_size"]].head(10)) ``` Expected output format: ``` name p_value intersection_size 0 mitotic cell cycle 2.1e-12 47 1 extracellular matrix orga 8.7e-09 31 ``` Run enrichment separately on the upregulated and downregulated sets when they point in opposite biological directions, which they frequently do. ## Extending the Workflow: Single-Cell, Spatial, and Proteomics ### Pseudobulk differential expression for single-cell data Single-cell RNA-seq cannot be analyzed with the same gene-by-cell test as bulk data, because cells from the same donor are not independent replicates. The standard approach is pseudobulk: aggregate raw counts per cell type per biological sample, then run the same bulk differential expression workflow on the aggregated matrix. This preserves the sample as the unit of replication. A published step-by-step Python workflow for inflammatory bowel disease data follows exactly this route, covering data access, integration of multiple datasets, mapping cells to a single-cell reference atlas to infer identities, and then applying a pseudobulk strategy to assess inflammation-associated changes in cell phenotypes [8]. Single-cell downstream analysis routinely includes differential [gene expression](/blog/guides/gene-expression) as a standard workflow component alongside GO and KEGG pathway analysis [9]. The pseudobulk code skeleton: ```python import scanpy as sc import decoupler as dc adata = sc.read_h5ad("integrated_adata.h5ad") pdata = dc.get_pseudobulk( adata, sample_col="sample_id", groups_col="cell_type", mode="sum", ) pdata = dc.pp.filter_samples(pdata, min_count=10, min_cells=10) ``` Then feed `pdata.to_df()` and its metadata through the same PyDESeq2 steps as above, with the design blocking on donor. Tools that operate directly at the single-cell level also exist. GEDI performs batch correction, latent representation learning, and clustering-free differential expression in a unified multi-sample model, and its 2.0 reimplementation reduces peak memory by 40 to 63.6 percent on average and speeds up execution by roughly 3-fold single-threaded, with R and Python interfaces [10]. For exploratory work, Scanpy's `rank_genes_groups` gives a fast first look at cluster markers, and tools such as scNetViz build on Scanpy to compute per-cluster differential expression and construct gene interaction networks for the significant hits inside Cytoscape [11]. Interactive tools such as scSelector let you lasso-select cells in an embedding and run DE and enrichment on that arbitrary selection [12], which is useful for hypotheses but should not replace pseudobulk inference in a publication. ### Spatial transcriptomics Spatially resolved data add a coordinate dimension, and differential expression there usually means cell-type-specific or domain-specific testing that respects spatial context. Practical guides cover the major R and Python object classes for spatial transcriptomics and describe the path from [preprocessing](/blog/data-analysis/what-is-data-preprocessing) through deconvolution and cell-type-specific differential expression [13]. Python frameworks also exist for volumetric reconstruction, for example sc3D, which registers serial spatial sections into AnnData objects compatible with Scanpy and supports spatial differential gene expression analysis [14]. ### Proteomics and other data types The same statistical logic applies beyond RNA. PBLMM performs differential expression analysis on shotgun proteomics using peptide-based linear mixed regression models, ships as a Python package, and shows superior [statistical power](/blog/guides/statistical-power-what-it-is-and-why-it-matters-in-research) in situations of low effect size or low sample size compared with classical protein-level inference [15]. For long-read transcriptome sequencing, IsoTools provides a Python framework that includes a beta-binomial approach for detecting differential splicing events, extending the DE question from genes to transcript isoforms [16]. The conceptual lesson generalizes: match the statistical model to the data-generating process, aggregate to the correct unit of replication, and keep the input counts raw until the model sees them. A note on Biopython, since it often appears in this keyword space. Biopython is the general-purpose Python library for sequence and file-format work (FASTA and GenBank parsing, Entrez queries against NCBI, running external tools). In a differential expression workflow its role is upstream and peripheral, for example fetching reference sequences or parsing annotation records, while the DE statistics live in PyDESeq2 and scanpy. If your starting point is the official Biopython tutorial, treat sequence handling as a separate skill that occasionally plugs into your expression pipeline rather than as the engine of the analysis. ## Common Mistakes and Limitations These are the errors that most often invalidate a differential expression python workflow, listed roughly in order of frequency. **Feeding normalized values or TPM to the DE tool.** The negative binomial model in DESeq2 and PyDESeq2 requires raw counts. TPM and FPKM correct for gene length, which destroys the count-based variance model. If your matrix has decimals from Salmon, round them through the tximport route in step 2, never by normalizing first. **Treating cells as replicates in single-cell data.** Running a per-cell test across conditions inflates your effective sample size enormously and produces false positives at scale. Pseudobulk per sample per cell type is the accepted fix [8]. **Confounded design.** If every control sample was sequenced in batch A and every treated sample in batch B, condition and batch are perfectly correlated and cannot be separated. This is an experimental design failure that no Python package can undo. Check for confounding before sequencing, and include every known nuisance factor in the design formula. **P-hacking the thresholds.** Adjusting the alpha, the fold-change cutoff, or the filtering rule after seeing the results turns an FDR of 0.05 into something much larger. Set thresholds before looking and report them. **Comparing across different annotations or genome builds.** If half your samples were counted against GENCODE v33 and half against v44, gene identifiers will silently disagree. Fix the reference and the transcript-to-gene map once, at the start [3]. **Ignoring the p-value correction.** Filtering on raw `pvalue < 0.05` with 20,000 tested genes yields roughly a thousand false positives by chance under the null. Use `padj`. **Misreading shrunken versus unshrunken fold changes.** Unshrunken log2 fold changes for lowly expressed genes can be extreme and misleading. Rank and visualize with the shrunk values from step 6. **Skipping PCA before testing.** An outlier sample, a swap of two filenames, or a sex mismatch shows up in a PCA plot in thirty seconds. It shows up in a results table never. **A limitation to state plainly.** Differential expression tells you that a gene differs between conditions, not why. Fold changes are averages over the cells captured, they can be driven by changes in cell-type composition rather than per-cell expression, and validation at the protein level or by an orthogonal assay is expected for key findings. Individual experimental decisions, especially for clinical samples, need domain-specific expertise and, where relevant, the input of a qualified clinician or pathologist. ## Checklist Before You Trust the Results Run through this list every time, in this order: 1. Sample names match exactly between counts and metadata, verified by assertion 2. Counts are raw integers, never TPM or log values 3. PCA separates by biology, not by batch, or batch is in the design 4. The design formula includes all nuisance covariates before the condition 5. The contrast direction was verified against a known positive control 6. `padj`, not `pvalue`, was used for significance 7. Effect sizes are taken from the shrunk estimates 8. Up and downregulated genes were enriched separately 9. Environment versions were pinned and recorded 10. The filtering rule for low-count genes is documented in the methods text ## Frequently Asked Questions ### Can I run the entire differential expression analysis in Python without touching R? Yes. PyDESeq2 reimplements the DESeq2 statistical workflow natively in Python and integrates with pandas and the scverse ecosystem, so normalization, testing, shrinkage, and multiple-testing correction all run in Python [1]. ### Should I give PyDESeq2 raw counts or normalized counts? Raw integer counts. The tool computes its own size-factor normalization internally, and supplying TPM, FPKM, or pre-normalized values breaks the negative binomial variance model. ### How many biological replicates do I need? At least three per group is the practical floor, but statistical power rises steeply beyond that, and outlier refitting is not available in very small designs. Plan replicates with a power calculation where possible. ### What is the difference between pvalue and padj in the results table? `pvalue` is the raw per-gene test statistic, while `padj` is the Benjamini-Hochberg false discovery rate adjusted value. Filter on `padj` because thousands of raw tests guarantee false positives at the 0.05 level. ### How do I handle batch effects in my samples? Include batch as a covariate in the design formula, for example `~ batch + condition`. If batch is perfectly confounded with condition, the effects cannot be separated statistically and the experiment needs redesign. ### Can I use this same workflow for single-cell RNA-seq? Not directly on cells. Aggregate raw counts per cell type per biological sample into pseudobulk profiles and run the bulk workflow on those, with donor as a blocking factor in the design. ### What does fold-change shrinkage actually do? It pulls extreme log2 fold changes of low-count genes toward zero using an empirical prior, producing effect-size estimates that are more stable for ranking, plotting, and biological interpretation. ### Which tool should I use to convert transcript-level quantifications to gene counts? pytximport handles this in Python, supporting multiple quantifier formats, bias correction, inferential replicates, and gene-level summarization, and is designed for reproducible Bioconda and Snakemake-based workflows [3]. ## Sources 1. [PyDESeq2: a python package for bulk RNA-seq differential expression analysis.](https://pubmed.ncbi.nlm.nih.gov/37669147/) 2. [From bulk RNA sequencing to spatial transcriptomics: a comparative review of differential gene expression analysis methods.](https://pubmed.ncbi.nlm.nih.gov/41353326/) 3. [Gene count estimation with pytximport enables reproducible analysis of bulk RNA sequencing data in Python.](https://pubmed.ncbi.nlm.nih.gov/39565903/) 4. [Streamline Protocol for Bulk-RNA Sequencing: From Data Extraction to Expression Analysis.](https://pubmed.ncbi.nlm.nih.gov/41637157/) 5. [HTSeq--a Python framework to work with high-throughput sequencing data.](https://pubmed.ncbi.nlm.nih.gov/25260700/) 6. [Transcriptome analysis of the response of Burmese python to digestion.](https://pubmed.ncbi.nlm.nih.gov/28873961/) 7. [CONSTANd: An Efficient Normalization Method for Relative Quantification in Small- and Large-Scale Omics Experiments in R BioConductor and Python.](https://pubmed.ncbi.nlm.nih.gov/33703904/) 8. [Using In Silico Pseudobulk Differential Expression and Trajectory Analysis of Single Cell RNA Sequencing to Reveal Stem Cell State Transitions in Colonic Inflammation.](https://pubmed.ncbi.nlm.nih.gov/42763849/) 9. [SINGLE-CELL TRANSCRIPTOME ANALYSIS IN HEALTH AND DISEASE.](https://pubmed.ncbi.nlm.nih.gov/37962963/) 10. [Atlas-level single-cell integration and clustering-free differential expression analysis with GEDI 2.0.](https://pubmed.ncbi.nlm.nih.gov/42178201/) 11. [scNetViz: from single cells to networks using Cytoscape.](https://pubmed.ncbi.nlm.nih.gov/34912541/) 12. [scSelector: A Flexible Single-Cell Data Analysis Assistant for Biomedical Researchers.](https://pubmed.ncbi.nlm.nih.gov/41595423/) 13. [From Raw Data to Biological Insights: A Practical Guide for Spatial Transcriptomics Analysis in R and Python.](https://pubmed.ncbi.nlm.nih.gov/42062686/) 14. [sc3D: A Comprehensive Tool for 3D Spatial Transcriptomic Analysis.](https://pubmed.ncbi.nlm.nih.gov/41769261/) 15. [PBLMM: Peptide-based linear mixed models for differential expression analysis of shotgun proteomics data.](https://pubmed.ncbi.nlm.nih.gov/35132673/) 16. [IsoTools: a flexible workflow for long-read transcriptome sequencing analysis.](https://pubmed.ncbi.nlm.nih.gov/37267159/)