drug-discovery-prompts/upstream/mims-harvard-ToolUniverse/skills/tooluniverse-rnaseq-deseq2/references/troubleshooting.md

title, task, lineage_type, upstream_source, upstream_sha, imported_at, prompt_class, upstream_changes, author, validated

title	task	lineage_type	upstream_source	upstream_sha	imported_at	prompt_class	upstream_changes	author	validated
Troubleshooting Guide		import	https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-rnaseq-deseq2/references/troubleshooting.md	e2520a96	2026-06-26	prompt	accepted	upstream	false

Troubleshooting Guide

Common issues and solutions for RNA-seq analysis with PyDESeq2.

Data Loading Issues

"No matching samples between counts and metadata"

Cause: Sample names don't match between count matrix and metadata.

Solutions:

# Check sample names
print("Count samples:", list(counts.index)[:5])
print("Metadata samples:", list(metadata.index)[:5])

# Try case-insensitive matching
counts.index = counts.index.str.lower()
metadata.index = metadata.index.str.lower()

# Remove whitespace
counts.index = counts.index.str.strip()
metadata.index = metadata.index.str.strip()

# Try transpose
if set(counts.columns) & set(metadata.index):
    counts = counts.T

"Non-integer counts detected"

Cause: Data is normalized (FPKM, TPM) or has floating point values.

Solutions:

# Option 1: Round to integers (acceptable if close to integers)
counts = counts.round().astype(int)

# Option 2: Use t-test instead of DESeq2 for normalized data
from scipy import stats
stat, pval = stats.ttest_ind(group1, group2)

# Check if data is truly raw counts
print(counts.head())
print("Min value:", counts.min().min())
print("Has decimals:", (counts % 1 != 0).any().any())

Matrix orientation wrong

Cause: Genes are rows instead of columns.

Solution:

# Check shape
print(f"Counts shape: {counts.shape}")  # Should be (n_samples, n_genes)
print(f"Samples >> genes? {counts.shape[0] > counts.shape[1] * 5}")

# Transpose if needed
if counts.shape[0] > counts.shape[1] * 5:
    counts = counts.T

DESeq2 Execution Issues

"Dispersion trend did not converge"

Cause: Small sample size or low variation.

Solution:

# Use mean fit instead of parametric
dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~condition",
    fit_type='mean',  # Instead of 'parametric'
    quiet=True
)
dds.deseq2()

"Contrast not found"

Cause: Wrong factor/level names in contrast.

Solution:

# Check available levels
print("Available factors:", metadata.columns)
print("Condition levels:", metadata['condition'].unique())

# Verify exact names (case-sensitive)
contrast = ['condition', 'treatment', 'control']  # Must match exactly

"All genes filtered out"

Cause: Too strict pre-filtering or all counts are zero.

Solution:

# Check data quality
print("Zero genes:", (counts.sum(axis=0) == 0).sum())
print("Low count genes:", (counts.sum(axis=0) < 10).sum())

# Remove only zero genes
nonzero = counts.sum(axis=0) > 0
counts = counts.loc[:, nonzero]

# Don't pre-filter too aggressively
# DESeq2 handles low counts internally

"Single replicate per condition"

Cause: Only 1 sample per condition - cannot estimate dispersion.

Solution:

# Cannot run DESeq2 with single replicates
# Use fold-change only
mean_A = counts[samples_A].mean(axis=1)
mean_B = counts[samples_B].mean(axis=1)
log2fc = np.log2((mean_B + 1) / (mean_A + 1))

# Or pool replicates from similar conditions

Reference Level Issues

Wrong reference level

Cause: Reference level not set correctly.

Solution:

# In PyDESeq2, FIRST category is reference
metadata['condition'] = pd.Categorical(
    metadata['condition'],
    categories=['control', 'treatment']  # Control first = reference
)

# Verify
print("Categories:", metadata['condition'].cat.categories)
print("First category (reference):", metadata['condition'].cat.categories[0])

LFC Shrinkage Issues

"Coefficient not found for shrinkage"

Cause: Wrong coefficient name format.

Solution:

# Check available coefficients
print("Available coefficients:", list(dds.varm['LFC'].columns))

# Standard format: factor[T.level]
coeff = 'condition[T.treatment]'

# Verify before shrinking
if coeff in dds.varm['LFC'].columns:
    stat_res.lfc_shrink(coeff=coeff)
else:
    print(f"ERROR: '{coeff}' not found")
    print("Skipping shrinkage")

Result Extraction Issues

"NaN in padj column"

Cause: Independent filtering removed genes with insufficient evidence.

Solution:

# This is EXPECTED behavior
# Remove NaN before counting DEGs
sig_genes = results.dropna(subset=['padj'])
sig_genes = sig_genes[sig_genes['padj'] < 0.05]

# Don't include NaN genes in counts
answer = len(sig_genes)  # Correct
# NOT: len(results[results['padj'] < 0.05])  # Wrong, includes NaN

"Gene not found in results"

Cause: Gene name case mismatch or gene filtered out.

Solution:

# Case-insensitive search
def find_gene(results_df, gene_name):
    # Exact match
    if gene_name in results_df.index:
        return gene_name

    # Case-insensitive
    idx_lower = {g.lower(): g for g in results_df.index}
    if gene_name.lower() in idx_lower:
        return idx_lower[gene_name.lower()]

    # Partial match
    matches = [g for g in results_df.index if gene_name.lower() in g.lower()]
    if len(matches) == 1:
        return matches[0]
    elif len(matches) > 1:
        print(f"Multiple matches: {matches}")

    return None

gene = find_gene(results, "TP53")
if gene:
    lfc = results.loc[gene, 'log2FoldChange']

Enrichment Analysis Issues

"No enrichment results"

Cause: Gene list too small or wrong organism/library.

Solution:

# Check gene list
print(f"Gene list size: {len(gene_list)}")
print(f"Sample genes: {gene_list[:5]}")

# Try different libraries
libraries = [
    'GO_Biological_Process_2023',
    'GO_Biological_Process_2021',
    'KEGG_2021_Human',
    'KEGG_2019_Mouse'
]

for lib in libraries:
    try:
        enr = gp.enrich(gene_list=gene_list, gene_sets=lib, outdir=None, no_plot=True)
        if len(enr.results) > 0:
            print(f"✓ {lib}: {len(enr.results)} results")
            break
    except Exception as e:
        print(f"✗ {lib}: {e}")

"gseapy organism parameter deprecated"

Cause: Using old gseapy syntax with organism parameter.

Solution:

# OLD (doesn't work in gseapy >= 1.1)
# enr = gp.enrich(gene_list=gene_list, organism='human', gene_sets='GO_Biological_Process')

# NEW (correct)
enr = gp.enrich(
    gene_list=gene_list,
    gene_sets='GO_Biological_Process_2023',  # Organism in library name
    outdir=None
)

Memory Issues

"Memory error with large gene set"

Cause: Too many genes (60K+ genes with many samples).

Solution:

# Pre-filter low-count genes
min_count = 10
keep = counts.sum(axis=0) >= min_count
counts_filtered = counts.loc[:, keep]
print(f"Kept {keep.sum()} / {len(keep)} genes")

# Use sparse matrices (if applicable)
from scipy.sparse import csr_matrix
counts_sparse = csr_matrix(counts.values)

Multi-Factor Design Issues

"Singular design matrix"

Cause: Confounded factors (e.g., batch perfectly correlated with condition).

Solution:

# Check confounding
print(pd.crosstab(metadata['batch'], metadata['condition']))

# If confounded, cannot separate effects
# Remove confounded factor from design
design = "~condition"  # Remove batch if confounded

Performance Issues

"DESeq2 takes too long"

Cause: Large dataset or many factors.

Solutions:

# Use parallel processing (if available)
import multiprocessing
n_cpus = multiprocessing.cpu_count()
# Note: PyDESeq2 doesn't directly support n_cpus, but numpy may use multiple cores

# Pre-filter more aggressively
min_count = 10
min_samples = 3
keep = (counts >= min_count).sum(axis=0) >= min_samples
counts_filtered = counts.loc[:, keep]

# Use mean fit instead of parametric (faster)
dds = DeseqDataSet(counts=counts, metadata=metadata, design="~condition",
                   fit_type='mean', quiet=True)

Validation Issues

Results don't match expected answer

Checks:

# 1. Check reference level
print("Reference:", metadata['condition'].cat.categories[0])

# 2. Check contrast direction
print("Contrast:", contrast)  # [factor, numerator, denominator]

# 3. Check filtering thresholds
print(f"padj < {padj_threshold}")
print(f"|log2FC| > {lfc_threshold}")

# 4. Check if shrinkage was applied
print("Shrinkage applied:", 'lfcSE' in results.columns)

# 5. Check for NaN handling
print("Genes with NaN padj:", results['padj'].isna().sum())

Debugging Checklist

Before reporting issues, verify:

• Count matrix has samples as rows, genes as columns
• Counts are non-negative integers
• Metadata index matches count matrix index exactly
• Design formula references valid column names in metadata
• Reference level is set correctly (first category in Categorical)
• Contrast factor and levels exist in metadata
• LFC shrinkage coefficient name matches pydeseq2 format
• Filtering thresholds match question exactly
• NaN values in padj are excluded from DEG counts

Getting Help

# Print diagnostic information
print("=" * 50)
print("DIAGNOSTIC INFORMATION")
print("=" * 50)
print(f"Counts shape: {counts.shape}")
print(f"Metadata shape: {metadata.shape}")
print(f"Design: {design}")
print(f"Reference level: {metadata['condition'].cat.categories[0]}")
print(f"Contrast: {contrast}")
print(f"PyDESeq2 version: {pydeseq2.__version__}")
print(f"\nSample counts:\n{counts.iloc[:3, :3]}")
print(f"\nMetadata:\n{metadata.head(3)}")
print(f"\nResults:\n{results.head(3)}")
print("=" * 50)