301 lines
8.4 KiB
Markdown
301 lines
8.4 KiB
Markdown
---
|
||
title: "Data Loading and Validation"
|
||
task: ""
|
||
lineage_type: import
|
||
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-rnaseq-deseq2/references/data_loading.md
|
||
upstream_sha: e2520a96
|
||
imported_at: 2026-06-26
|
||
prompt_class: prompt
|
||
upstream_changes: accepted
|
||
author: upstream
|
||
validated: false
|
||
---
|
||
|
||
# Data Loading and Validation
|
||
|
||
Detailed guide for loading count matrices and metadata.
|
||
|
||
## Handling R Data Files (RDS)
|
||
|
||
If your data is in R format (`.rds`, `.RData`), convert to CSV first:
|
||
|
||
### Option 1: Convert in R (if available)
|
||
|
||
```r
|
||
# In R
|
||
result <- readRDS("deseq2_results.rds")
|
||
result_df <- as.data.frame(result)
|
||
write.csv(result_df, "deseq2_results.csv", row.names=TRUE)
|
||
```
|
||
|
||
Then load the CSV in Python as usual.
|
||
|
||
### Option 2: Use rpy2 in Python (requires R installed)
|
||
|
||
```python
|
||
import rpy2.robjects as ro
|
||
from rpy2.robjects import pandas2ri
|
||
pandas2ri.activate()
|
||
|
||
# Read RDS file
|
||
ro.r(f'df <- as.data.frame(readRDS("{rds_file}"))')
|
||
|
||
# Convert to pandas
|
||
df = pandas2ri.rpy2py(ro.r['df'])
|
||
```
|
||
|
||
**Note**: This requires R and the relevant R packages (e.g., DESeq2) installed. If unavailable, use Option 1 on a machine with R.
|
||
|
||
### Working with Pre-computed DESeq2 Results
|
||
|
||
If you have an RDS file containing DESeq2 results (not count data):
|
||
|
||
1. **Convert to CSV** (Option 1 or 2 above)
|
||
2. **Filter genes** directly in pandas:
|
||
|
||
```python
|
||
# Read pre-computed results
|
||
results = pd.read_csv("deseq2_results.csv", index_col=0)
|
||
|
||
# Filter upregulated genes
|
||
upregulated = results[
|
||
(results['log2FoldChange'] > 0) & # Positive = upregulated
|
||
(results['padj'] < 0.05) # Significant
|
||
]
|
||
|
||
# Get gene list
|
||
genes = upregulated.index.tolist()
|
||
|
||
# Use with /tooluniverse-gene-enrichment skill
|
||
```
|
||
|
||
3. **Skip to enrichment** - Use gene list with enrichment analysis (Step 5 in workflow)
|
||
|
||
## Load Count Matrix
|
||
|
||
```python
|
||
import pandas as pd
|
||
import numpy as np
|
||
import os
|
||
|
||
def load_count_matrix(file_path, **kwargs):
|
||
"""Load count matrix from various formats.
|
||
|
||
Expects: genes as rows/columns, samples as rows/columns.
|
||
PyDESeq2 requires: samples as rows, genes as columns.
|
||
"""
|
||
ext = os.path.splitext(file_path)[1].lower()
|
||
|
||
if ext in ['.csv']:
|
||
df = pd.read_csv(file_path, index_col=0, **kwargs)
|
||
elif ext in ['.tsv', '.txt']:
|
||
df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
|
||
elif ext in ['.h5ad']:
|
||
import anndata
|
||
adata = anndata.read_h5ad(file_path)
|
||
df = pd.DataFrame(
|
||
adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
|
||
index=adata.obs_names,
|
||
columns=adata.var_names
|
||
)
|
||
return df, adata.obs # Return metadata too if available
|
||
else:
|
||
# Try tab-separated as default
|
||
df = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
|
||
|
||
return df
|
||
```
|
||
|
||
## Orient Matrix
|
||
|
||
**CRITICAL**: PyDESeq2 expects samples as rows, genes as columns.
|
||
|
||
```python
|
||
def orient_count_matrix(df, metadata_samples=None):
|
||
"""Ensure samples are rows and genes are columns.
|
||
|
||
Heuristic: if column count >> row count, genes are likely columns (correct).
|
||
If row count >> column count, genes are likely rows (need transpose).
|
||
If metadata_samples provided, match against index and columns.
|
||
"""
|
||
if metadata_samples is not None:
|
||
# Check if samples match rows or columns
|
||
row_match = len(set(df.index) & set(metadata_samples))
|
||
col_match = len(set(df.columns) & set(metadata_samples))
|
||
if col_match > row_match:
|
||
df = df.T
|
||
return df
|
||
|
||
# Heuristic: typical RNA-seq has 10-1000 samples and 10000-60000 genes
|
||
if df.shape[0] > df.shape[1] * 5: # Many more rows than columns
|
||
df = df.T # Transpose: genes were rows
|
||
|
||
return df
|
||
```
|
||
|
||
## Load Metadata
|
||
|
||
```python
|
||
def load_metadata(file_path, **kwargs):
|
||
"""Load sample metadata (colData in R)."""
|
||
ext = os.path.splitext(file_path)[1].lower()
|
||
|
||
if ext in ['.csv']:
|
||
meta = pd.read_csv(file_path, index_col=0, **kwargs)
|
||
elif ext in ['.tsv', '.txt']:
|
||
meta = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
|
||
else:
|
||
meta = pd.read_csv(file_path, sep='\t', index_col=0, **kwargs)
|
||
|
||
return meta
|
||
```
|
||
|
||
## Validate and Align
|
||
|
||
```python
|
||
def validate_inputs(counts, metadata):
|
||
"""Validate count matrix and metadata alignment."""
|
||
issues = []
|
||
|
||
# Check sample alignment
|
||
count_samples = set(counts.index)
|
||
meta_samples = set(metadata.index)
|
||
|
||
if count_samples != meta_samples:
|
||
common = count_samples & meta_samples
|
||
if len(common) == 0:
|
||
# Try matching columns
|
||
if set(counts.columns) & meta_samples:
|
||
counts = counts.T
|
||
count_samples = set(counts.index)
|
||
common = count_samples & meta_samples
|
||
|
||
if len(common) > 0:
|
||
counts = counts.loc[sorted(common)]
|
||
metadata = metadata.loc[sorted(common)]
|
||
issues.append(f"Aligned to {len(common)} common samples")
|
||
else:
|
||
issues.append("ERROR: No matching samples between counts and metadata")
|
||
return None, None, issues
|
||
|
||
# Ensure integer counts
|
||
if counts.dtypes.apply(lambda x: x == float).any():
|
||
if (counts % 1 == 0).all().all():
|
||
counts = counts.astype(int)
|
||
else:
|
||
# Might be normalized data - round to integers for DESeq2
|
||
issues.append("WARNING: Non-integer counts detected. Rounding to integers.")
|
||
counts = counts.round().astype(int)
|
||
|
||
# Remove genes with zero counts across all samples
|
||
nonzero_mask = counts.sum(axis=0) > 0
|
||
n_removed = (~nonzero_mask).sum()
|
||
if n_removed > 0:
|
||
counts = counts.loc[:, nonzero_mask]
|
||
issues.append(f"Removed {n_removed} genes with zero counts across all samples")
|
||
|
||
# Remove negative values
|
||
if (counts < 0).any().any():
|
||
issues.append("WARNING: Negative counts detected. Setting to 0.")
|
||
counts = counts.clip(lower=0)
|
||
|
||
return counts, metadata, issues
|
||
```
|
||
|
||
## Subset Samples
|
||
|
||
```python
|
||
def subset_samples(counts, metadata, condition_col, values=None, exclude=None):
|
||
"""Subset samples based on metadata conditions."""
|
||
if values is not None:
|
||
mask = metadata[condition_col].isin(values)
|
||
elif exclude is not None:
|
||
mask = ~metadata[condition_col].isin(exclude)
|
||
else:
|
||
return counts, metadata
|
||
|
||
metadata = metadata[mask]
|
||
counts = counts.loc[metadata.index]
|
||
return counts, metadata
|
||
```
|
||
|
||
## Example: Load and Validate
|
||
|
||
```python
|
||
# Load data
|
||
counts = load_count_matrix("counts.csv")
|
||
metadata = load_metadata("metadata.csv")
|
||
|
||
# Orient if needed
|
||
counts = orient_count_matrix(counts, metadata.index)
|
||
|
||
# Validate and align
|
||
counts, metadata, issues = validate_inputs(counts, metadata)
|
||
|
||
# Print validation issues
|
||
for issue in issues:
|
||
print(f" {issue}")
|
||
|
||
# Subset if needed
|
||
counts, metadata = subset_samples(
|
||
counts, metadata,
|
||
condition_col='treatment',
|
||
values=['control', 'treated']
|
||
)
|
||
|
||
print(f"\nFinal dimensions:")
|
||
print(f" Counts: {counts.shape[0]} samples × {counts.shape[1]} genes")
|
||
print(f" Metadata: {metadata.shape[0]} samples")
|
||
```
|
||
|
||
## Handling Different Input Formats
|
||
|
||
### CSV with genes as columns (correct)
|
||
```
|
||
,Gene1,Gene2,Gene3
|
||
Sample1,100,50,200
|
||
Sample2,120,60,180
|
||
```
|
||
No action needed.
|
||
|
||
### CSV with genes as rows (needs transpose)
|
||
```
|
||
,Sample1,Sample2,Sample3
|
||
Gene1,100,120,90
|
||
Gene2,50,60,45
|
||
```
|
||
Call `orient_count_matrix()`.
|
||
|
||
### H5AD (AnnData)
|
||
```python
|
||
import anndata
|
||
adata = anndata.read_h5ad("data.h5ad")
|
||
counts = pd.DataFrame(
|
||
adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
|
||
index=adata.obs_names,
|
||
columns=adata.var_names
|
||
)
|
||
metadata = adata.obs
|
||
```
|
||
|
||
### Pre-normalized data (FPKM, TPM)
|
||
If data is NOT raw counts:
|
||
```python
|
||
# Option 1: Round to integers (acceptable for DESeq2)
|
||
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)
|
||
```
|
||
|
||
## Common Issues
|
||
|
||
| Issue | Symptom | Solution |
|
||
|-------|---------|----------|
|
||
| Samples don't match | "No matching samples" error | Check if transpose needed, strip whitespace |
|
||
| Float counts | "Non-integer counts" warning | Round to integers or use t-test |
|
||
| Sample name mismatch | Different # of samples | Use `set(counts.index) & set(metadata.index)` |
|
||
| Zero-count genes | All-zero columns | Pre-filter before DESeq2 |
|
||
| Negative counts | Impossible biology | Set to 0, investigate source |
|