File size: 6,470 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 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 | """Map every processed count column to the E-MTAB-5009 metadata."""
from pathlib import Path
import re
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"
OUTPUT_PATH = PROJECT_ROOT / "data" / "processed" / "sample_sheet.csv"
COUNT_COLUMN_PATTERN = re.compile(
r"^(?P<ecotype>.+)_(?P<condition>[CS])_B(?P<block>.+?)"
r"(?:_T(?P<technical_number>\d+))?\.TAIR10_allchr\.htseq$"
)
def normalize_name(value: str) -> str:
"""Make names such as Ag_0 and Ag-0 comparable."""
return re.sub(r"[^a-z0-9]", "", value.lower())
def technical_sort_key(value: str) -> int:
match = re.search(r"(\d+)$", value.strip())
return int(match.group(1)) if match else 1
def parse_count_columns(columns: list[str]) -> pd.DataFrame:
records = []
for column in columns:
match = COUNT_COLUMN_PATTERN.fullmatch(column)
if match is None:
raise ValueError(f"Could not parse count column: {column}")
condition_code = match.group("condition")
technical_number = match.group("technical_number")
records.append(
{
"count_column": column,
"normalized_ecotype": normalize_name(match.group("ecotype")),
"condition": "control" if condition_code == "C" else "drought",
"block": match.group("block"),
"column_technical_number": int(technical_number or 1),
}
)
return pd.DataFrame(records)
def prepare_run_metadata() -> pd.DataFrame:
metadata = pd.read_csv(METADATA_PATH, sep="\t", keep_default_na=False)
# Paired-end read 1 and read 2 occupy separate SDRF rows but share one run.
runs = metadata.drop_duplicates(subset="Comment[ENA_RUN]").copy()
runs["normalized_ecotype"] = runs["Factor Value[ecotype]"].map(normalize_name)
runs["condition"] = runs["Factor Value[environmental stress]"].map(
{"control": "control", "Mild drought stress": "drought"}
)
runs["block"] = runs["Factor Value[block]"].str.strip()
runs["metadata_technical_number"] = runs[
"Characteristics[technical replicate]"
].map(technical_sort_key)
return runs
def build_sample_sheet() -> tuple[pd.DataFrame, pd.DataFrame]:
count_columns = pd.read_csv(COUNTS_PATH, sep="\t", nrows=0).columns.tolist()
parsed = parse_count_columns(count_columns)
runs = prepare_run_metadata()
keys = ["normalized_ecotype", "condition", "block"]
mapped_groups = []
for key, count_group in parsed.groupby(keys, sort=False):
run_group = runs
for key_name, key_value in zip(keys, key):
run_group = run_group.loc[run_group[key_name].eq(key_value)]
count_group = count_group.sort_values("column_technical_number")
run_group = run_group.sort_values("metadata_technical_number")
if len(count_group) != len(run_group):
raise ValueError(
f"Mapping mismatch for {key}: "
f"{len(count_group)} count columns and {len(run_group)} metadata runs"
)
count_group = count_group.reset_index(drop=True)
run_group = run_group.reset_index(drop=True)
mapped_groups.append(
pd.DataFrame(
{
"count_column": count_group["count_column"],
"sample_group_id": ["|".join(key)] * len(count_group),
"ecotype": run_group["Factor Value[ecotype]"],
"condition": count_group["condition"],
"block": count_group["block"],
"technical_replicate": run_group[
"Characteristics[technical replicate]"
].str.strip(),
"assay_name": run_group["Assay Name"],
"ena_sample": run_group["Comment[ENA_SAMPLE]"],
"ena_experiment": run_group["Comment[ENA_EXPERIMENT]"],
"ena_run": run_group["Comment[ENA_RUN]"],
}
)
)
sample_sheet = pd.concat(mapped_groups, ignore_index=True)
processed_keys = set(map(tuple, parsed[keys].drop_duplicates().to_numpy()))
metadata_groups = runs[keys + ["Factor Value[ecotype]"]].drop_duplicates()
missing_mask = [
tuple(row) not in processed_keys
for row in metadata_groups[keys].to_numpy()
]
missing_processed_groups = metadata_groups.loc[missing_mask].copy()
return sample_sheet, missing_processed_groups
def validate(sample_sheet: pd.DataFrame) -> None:
count_columns = pd.read_csv(COUNTS_PATH, sep="\t", nrows=0).columns
if len(sample_sheet) != len(count_columns):
raise ValueError("The sample sheet does not contain every count column.")
if sample_sheet["count_column"].duplicated().any():
raise ValueError("The sample sheet contains duplicate count columns.")
if set(sample_sheet["count_column"]) != set(count_columns):
raise ValueError("Sample-sheet and count-table column names differ.")
if sample_sheet[
["ecotype", "condition", "block", "ena_run"]
].eq("").any().any():
raise ValueError("Required sample-sheet fields contain blank values.")
def main() -> None:
sample_sheet, missing_groups = build_sample_sheet()
validate(sample_sheet)
OUTPUT_PATH.parent.mkdir(parents=True, exist_ok=True)
sample_sheet.to_csv(OUTPUT_PATH, index=False)
group_sizes = sample_sheet.groupby("sample_group_id").size()
condition_pairs = (
sample_sheet[["ecotype", "condition"]]
.drop_duplicates()
.groupby("ecotype")["condition"]
.nunique()
)
print(f"Wrote: {OUTPUT_PATH}")
print(f"Mapped count columns: {len(sample_sheet):,}")
print(f"Sample groups after combining technical replicates: {group_sizes.size:,}")
print(f"Groups with more than one technical count column: {(group_sizes > 1).sum():,}")
print(f"Ecotypes in processed counts: {sample_sheet['ecotype'].nunique():,}")
print(f"Ecotypes having both conditions: {(condition_pairs == 2).sum():,}")
print("Metadata groups absent from the processed count table:")
print(
missing_groups[["Factor Value[ecotype]", "condition", "block"]]
.to_string(index=False)
)
if __name__ == "__main__":
main()
|