File size: 4,997 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 | """Combine technical counts, filter genes, and normalize E-MTAB-5009."""
from math import ceil
from pathlib import Path
import numpy as np
import pandas as pd
PROJECT_ROOT = Path(__file__).resolve().parents[1]
RAW_COUNTS_PATH = (
PROJECT_ROOT / "data" / "raw" / "E-MTAB-5009" / "WP2_100acc.htseq.txt"
)
SAMPLE_SHEET_PATH = PROJECT_ROOT / "data" / "processed" / "sample_sheet.csv"
OUTPUT_DIR = PROJECT_ROOT / "data" / "processed"
MIN_CPM = 1.0
MIN_SAMPLE_FRACTION = 0.10
def load_inputs() -> tuple[pd.DataFrame, pd.DataFrame]:
counts = pd.read_csv(RAW_COUNTS_PATH, sep="\t", index_col=0)
counts.index.name = "gene_id"
sample_sheet = pd.read_csv(SAMPLE_SHEET_PATH, keep_default_na=False)
if counts.columns.duplicated().any():
raise ValueError("Raw count columns are not unique.")
if sample_sheet["count_column"].duplicated().any():
raise ValueError("Sample-sheet count columns are not unique.")
if set(counts.columns) != set(sample_sheet["count_column"]):
raise ValueError("Raw counts and sample sheet contain different columns.")
if (counts < 0).any().any():
raise ValueError("Raw counts contain negative values.")
# Put sample-sheet rows in exactly the same order as the count columns.
sample_sheet = sample_sheet.set_index("count_column").loc[counts.columns]
sample_sheet.index.name = "count_column"
sample_sheet = sample_sheet.reset_index()
return counts, sample_sheet
def combine_technical_counts(
counts: pd.DataFrame, sample_sheet: pd.DataFrame
) -> tuple[pd.DataFrame, pd.DataFrame]:
column_to_group = sample_sheet.set_index("count_column")["sample_group_id"]
# Raw sequencing counts from technical measurements are additive.
aggregated = (
counts.T.join(column_to_group)
.groupby("sample_group_id", sort=False)
.sum()
.T
)
aggregated.index.name = "gene_id"
group_sheet = (
sample_sheet.groupby("sample_group_id", sort=False)
.agg(
ecotype=("ecotype", "first"),
condition=("condition", "first"),
block=("block", "first"),
technical_count_columns=("count_column", "count"),
source_count_columns=("count_column", ";".join),
)
.reset_index()
)
if list(aggregated.columns) != group_sheet["sample_group_id"].tolist():
raise ValueError("Aggregated counts and group sheet have different orders.")
return aggregated, group_sheet
def normalize_and_filter(
aggregated: pd.DataFrame,
) -> tuple[pd.DataFrame, pd.DataFrame, int]:
library_sizes = aggregated.sum(axis=0)
if (library_sizes == 0).any():
raise ValueError("At least one sample group has a library size of zero.")
cpm = aggregated.div(library_sizes, axis="columns") * 1_000_000
minimum_samples = ceil(MIN_SAMPLE_FRACTION * aggregated.shape[1])
keep_gene = (cpm >= MIN_CPM).sum(axis=1) >= minimum_samples
filtered_counts = aggregated.loc[keep_gene]
log2_cpm = np.log2(cpm.loc[keep_gene] + 1).astype("float32")
return filtered_counts, log2_cpm, minimum_samples
def main() -> None:
counts, sample_sheet = load_inputs()
aggregated, group_sheet = combine_technical_counts(counts, sample_sheet)
filtered_counts, log2_cpm, minimum_samples = normalize_and_filter(aggregated)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
aggregated.to_csv(
OUTPUT_DIR / "counts_aggregated.tsv.gz", sep="\t", compression="gzip"
)
filtered_counts.to_csv(
OUTPUT_DIR / "counts_filtered.tsv.gz", sep="\t", compression="gzip"
)
log2_cpm.to_csv(
OUTPUT_DIR / "expression_log2_cpm.tsv.gz", sep="\t", compression="gzip"
)
group_sheet.to_csv(OUTPUT_DIR / "sample_groups.csv", index=False)
original_libraries = counts.sum(axis=0)
combined_libraries = aggregated.sum(axis=0)
print(f"Original count columns: {counts.shape[1]:,}")
print(f"Sample groups after technical summation: {aggregated.shape[1]:,}")
print(f"Genes before filtering: {aggregated.shape[0]:,}")
print(
f"Filter: CPM >= {MIN_CPM:g} in at least "
f"{minimum_samples} of {aggregated.shape[1]} sample groups"
)
print(f"Genes retained: {filtered_counts.shape[0]:,}")
print(f"Genes removed: {aggregated.shape[0] - filtered_counts.shape[0]:,}")
print(
"Original library sizes (min / median / max): "
f"{original_libraries.min():,.0f} / "
f"{original_libraries.median():,.0f} / "
f"{original_libraries.max():,.0f}"
)
print(
"Combined library sizes (min / median / max): "
f"{combined_libraries.min():,.0f} / "
f"{combined_libraries.median():,.0f} / "
f"{combined_libraries.max():,.0f}"
)
print("Wrote counts_aggregated.tsv.gz")
print("Wrote counts_filtered.tsv.gz")
print("Wrote expression_log2_cpm.tsv.gz")
print("Wrote sample_groups.csv")
if __name__ == "__main__":
main()
|