File size: 2,919 Bytes
18410f0 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 | """Inspect the E-MTAB-5009 count table and sample metadata."""
from pathlib import Path
import pandas as pd
PROJECT_ROOT = Path(__file__).resolve().parents[1]
DATA_DIR = PROJECT_ROOT / "data" / "raw" / "E-MTAB-5009"
COUNTS_PATH = DATA_DIR / "WP2_100acc.htseq.txt"
METADATA_PATH = DATA_DIR / "E-MTAB-5009.sdrf.txt"
def inspect_counts() -> None:
# The first field in each data row is a gene ID. Its header cell is absent,
# so index_col=0 prevents the gene IDs from being mistaken for counts.
counts = pd.read_csv(COUNTS_PATH, sep="\t", index_col=0)
counts.index.name = "gene_id"
library_sizes = counts.sum(axis=0)
print("COUNT TABLE")
print(f"Path: {COUNTS_PATH}")
print(f"Shape: {counts.shape[0]:,} genes x {counts.shape[1]:,} count columns")
print(f"Duplicate gene IDs: {counts.index.duplicated().sum():,}")
print(f"Duplicate column names: {counts.columns.duplicated().sum():,}")
print(f"Missing values: {counts.isna().sum().sum():,}")
print(f"Smallest count: {counts.min().min():,}")
print(f"Largest count: {counts.max().max():,}")
print(f"Genes with zero counts in every column: {(counts.sum(axis=1) == 0).sum():,}")
print(
"Library sizes (min / median / max): "
f"{library_sizes.min():,.0f} / "
f"{library_sizes.median():,.0f} / "
f"{library_sizes.max():,.0f}"
)
print("\nFirst 5 genes and first 5 count columns:")
print(counts.iloc[:5, :5].to_string())
def inspect_metadata() -> None:
metadata = pd.read_csv(METADATA_PATH, sep="\t")
# One assay can occupy multiple SDRF rows because paired-end reads and
# sequencing lanes are described separately. Deduplicate before counting
# biological assays and conditions.
assays = metadata.drop_duplicates(subset="Assay Name")
print("\nMETADATA")
print(f"Path: {METADATA_PATH}")
print(f"Raw metadata shape: {metadata.shape[0]:,} rows x {metadata.shape[1]:,} columns")
print(f"Unique assays: {assays['Assay Name'].nunique():,}")
print(f"Unique ecotypes: {assays['Factor Value[ecotype]'].nunique():,}")
print("Unique assays by condition:")
print(assays["Factor Value[environmental stress]"].value_counts().to_string())
print(f"Missing cells in metadata: {metadata.isna().sum().sum():,}")
print("\nSelected metadata columns for the first 5 unique assays:")
selected = [
"Assay Name",
"Factor Value[ecotype]",
"Factor Value[environmental stress]",
"Factor Value[block]",
]
print(assays[selected].head().to_string(index=False))
def main() -> None:
for required_path in (COUNTS_PATH, METADATA_PATH):
if not required_path.exists():
raise FileNotFoundError(
f"Missing {required_path}. Run scripts/download_data.ps1 first."
)
inspect_counts()
inspect_metadata()
if __name__ == "__main__":
main()
|