# Comparing early- and later-onset colorectal cancer: finding EOCRC-specific expression candidates

> Run paired DESeq2 analyses on 21 early-onset colorectal cancer pairs from GSE196006 and 22 later-onset pairs from GSE251845, then compare Tumor−Normal changes and the paper's 8-gene signature.

## 1. What you will practice

This exercise asks one question.

> **How similar are expression changes between tumor and adjacent normal tissue in early- and later-onset colorectal cancer, and which genes show more pronounced changes in younger patients?**

**Early-onset colorectal cancer** (EOCRC) generally means colorectal cancer diagnosed before age 50. The comparison group, **later-onset colorectal cancer** (LOCRC), consists in this study of patients diagnosed after age 50. These are not different cancer types. They are study groups divided by age at diagnosis.

As in the [later-onset colorectal cancer exercise](/en/practice/colorectal-deg-gse251845/), you first compare Tumor and Normal samples within each patient. This time, run the analysis separately in the two age groups, then place the `log2FoldChange` values for the same genes side by side.

The data is split across two GEO accessions.

| Age group | GEO dataset | Patients | RNA-seq samples | Samples from one patient |
| --- | --- | ---: | ---: | --- |
| EOCRC | [GSE196006](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE196006) | 21 | 42 | 1 Tumor and 1 Normal |
| LOCRC | [GSE251845](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE251845) | 22 | 44 | 1 Tumor and 1 Normal |

Both datasets use a pair of surgically resected tumor tissue and adjacent non-tumor tissue from each patient. Because this is bulk RNA-seq from a mixture of cells rather than isolated cancer cells, the differences reflect both cancer-cell expression and changes in the cellular composition of the tissue. [Tumor tissue and adjacent normal tissue](/en/reference/tumor-adjacent-normal/) explains the meaning of this control group.

### What do the two count tables look like?

The first two patient pairs in the EOCRC file `GSE196006_raw_counts.csv.gz` look like this. These are not values invented for illustration. They come from the first three genes in the actual file.

| gene ID | `X15.018_L0_G821_htseq.out` | `X15.018_L7_G821_htseq.out` | `X15.130_B0_G821_htseq.out` | `X15.130_B7_G821_htseq.out` |
| --- | ---: | ---: | ---: | ---: |
| `ENSG00000000003` | 2,935 | 3,326 | 2,620 | 331 |
| `ENSG00000000005` | 33 | 189 | 39 | 1 |
| `ENSG00000000419` | 1,124 | 4,156 | 817 | 1,126 |

`15.018` and `15.130` are patient IDs. In the final part of each filename, `0` means Normal and `7` means Tumor. A letter such as `L` or `B` in the middle encodes tissue-location information and must not be interpreted as the condition.

The LOCRC file uses a different naming convention, such as `24C_htseq.out` and `24N_htseq.out`. Here, `C` means cancer, or Tumor, and `N` means Normal. In both files, rows are Ensembl gene IDs, columns are samples, and values are integer [HTSeq raw counts](/en/reference/htseq-counts/).

### Why not combine the two files directly?

EOCRC and LOCRC come from different GEO experiments. If you join `GSE196006` and `GSE251845` into one count matrix, **age group and experimental batch overlap completely**. Even if the two age groups separate in a plot, you cannot tell whether age caused the separation or whether it came from batch differences such as the time of library preparation or sequencing conditions.

This exercise does not combine the two count tables directly. It first calculates the `Tumor − Normal` change within each cohort, then connects those changes by gene ID.

```text
21 EOCRC pairs → paired DESeq2 → EOCRC log2FC for each gene
22 LOCRC pairs → paired DESeq2 → LOCRC log2FC for each gene
join both result tables by gene ID → compare changes for the same gene
```

Within-patient comparisons help avoid some large batch differences, but they do not eliminate experimental differences between the two cohorts. The result is therefore a set of **EOCRC-specific candidates**, not definitive evidence that age caused the difference.

### Order of checks in this exercise

1. Can all 21 EOCRC pairs and 22 LOCRC pairs be reconstructed correctly from the filenames?
2. Within each cohort, do the overall expression patterns of Tumor and Normal samples separate?
3. How similar are the directions of Tumor−Normal change for the same genes in the two cohorts?
4. Under the current reproduction conditions, do the eight genes from the paper still show larger changes in EOCRC?
5. Does a difference visible in the cohort average repeat across multiple patients?

### Scope of interpretation

This page does not reproduce the full analysis in the [original paper](https://pmc.ncbi.nlm.nih.gov/articles/PMC11047122/). It compares **gene-level expression changes** using the two public HTSeq count tables.

:::caution[A difference between significant and nonsignificant is not automatically a significant difference]
If `padj < 0.05` in EOCRC and `padj > 0.2` in LOCRC, that does not automatically make the effect difference between the two age groups statistically significant. When sample sizes and variances differ, one group can cross a threshold even if the actual changes are similar. The set filter in this exercise follows the paper's candidate-discovery procedure. It is not a direct test of the difference between the two effects.
:::

## 2. Ask an agent to prepare the exercise environment

Run Codex CLI or Claude Code, then copy and submit the prompt below. The agent will manage both cohorts together while keeping their analyses separate.

### Exercise environment setup prompt

```text
Prepare an environment for comparing EOCRC and LOCRC RNA-seq data.

Requirements:
1. Check whether the current folder is already the eocrc-locrc project. If not, create eocrc-locrc under the current folder and work inside it.
2. Check whether uv is installed. If it is not, do not install it yourself. Show only the official installation method and stop.
3. Initialize the project with uv init --bare --python 3.11.
4. Run uv add "pydeseq2==0.5.4" pandas matplotlib seaborn scikit-learn adjustText gprofiler-official.
5. Create analysis.py and an outputs directory. Do not add analysis code to analysis.py yet.
6. Add .venv/, data/, and __pycache__/ to .gitignore. Do not exclude outputs, so the results remain available for inspection.
7. Use uv run python to print the versions of Python, PyDESeq2, pandas, and matplotlib and verify the installation.
8. Report the files created or modified and the main package versions, then stop.

Do not modify the global Python environment or other files in the repository.
```

When preparation is complete, `data/` does not yet exist and `analysis.py` is empty. Starting with the next prompt, preserve the code from previous stages and continue adding to it.

## 3. Compare the two cohorts one step at a time with prompts

:::tip[How to compare your output with the example]
The values below are examples from a run using Python 3.13, PyDESeq2 0.5.4, and GEO files downloaded in August 2026. The exercise prompt pins Python 3.11, so exact values may differ slightly. Do not force the numbers to match. Check whether the sample structure, directions of change, and interpretation agree.
:::

### 3-1. Check patient pairs in both files

The first step converts filenames into metadata and confirms that all 43 patients have both Tumor and Normal samples. EOCRC and LOCRC use different naming conventions, so do not force both through one regular expression.

### Step 1 prompt: download data and inspect its structure

```text
Perform only the step that inspects the EOCRC and LOCRC data structures.

1. Create a data directory in the project root.
2. Download the two files below and validate their gzip archives. If a file already exists and passes validation, do not download it again.
   - data/GSE196006_raw_counts.csv.gz
     https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE196006&format=file&file=GSE196006_raw_counts.csv.gz
   - data/GSE251845_htseq_raw_counts.csv.gz
     https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE251845&format=file&file=GSE251845_htseq_raw_counts.csv.gz
3. Read each file with the first column as the index, separate HTSeq special counters that begin with __, then transpose so samples are rows.
4. EOCRC sample names look like X15.018_L0_G821_htseq.out. Interpret 15.018 as the patient ID, 0 as Normal, and 7 as Tumor. Do not use the letter in the middle as the condition.
5. LOCRC sample names look like 24C_htseq.out. Interpret C as Tumor and N as Normal, and normalize lowercase names such as 35c to uppercase.
6. Create cohort-specific metadata with cohort, patient, and condition columns.
7. Verify the following with assertions:
   - EOCRC has 42 samples from 21 patients, with 21 Tumor and 21 Normal
   - LOCRC has 44 samples from 22 patients, with 22 Tumor and 22 Normal
   - every patient has exactly one Tumor and one Normal
   - the row order of each count matrix matches its metadata
8. Print the original dimensions, special counters, the first rows of metadata, and sample counts by patient for both files. Explain the results in English, then stop.

Do not combine the count tables or run low-expression filtering, DESeq2, or PCA yet.
```

When read correctly, both files have 60,605 rows including special rows. Their sample axes contain 42 EOCRC and 44 LOCRC samples. The important check at this stage is that `0/7` and `C/N` were not handled in the same way.

### 3-2. Paired DESeq2 and PCA by cohort

Include the patient-specific baseline in each cohort's model. Both models use `~ patient + condition` as the [design formula](/en/reference/design-formula/), but patient IDs, size factors, and dispersions are calculated separately within each cohort.

### Step 2 prompt: paired DESeq2 and PCA by cohort

```text
Keep the previous data-validation code and add only paired DESeq2 and PCA for each cohort.

1. In each cohort, retain genes with a count of at least 10 in at least 3 samples from that cohort. Calculate the filters separately for the two cohorts.
2. Create separate DeseqDataSet objects for EOCRC and LOCRC with design="~ patient + condition" and n_cpus=1, then run deseq2 on each.
3. In both models, calculate a Tumor-versus-Normal contrast with Normal as the reference condition.
4. Save the results to outputs/eocrc-tumor-vs-normal.csv and outputs/locrc-tumor-vs-normal.csv. Preserve Ensembl gene IDs as the index.
5. Calculate VST with use_design=False from each DeseqDataSet, then calculate separate PCAs for the two cohorts. Draw them in left and right panels, but do not place samples from one cohort into the coordinates of the other cohort's PCA.
6. Show Normal samples as teal circles and Tumor samples as orange triangles, and save the plot to outputs/cohort-pca.png at 180 dpi.
7. Print the number of genes passing the filter in each cohort, the number satisfying both padj < 0.05 and abs(log2FoldChange) >= 1, and the variance explained by PC1 and PC2.
8. Run uv run python analysis.py and inspect the figure. Explain in English the separation and overlap of Tumor and Normal within each cohort, along with any distant samples, then stop.

Do not combine the two count matrices into one PCA or DESeq2 model.
```

[![Two PCA panels made by transforming the 42 EOCRC and 44 LOCRC samples separately with VST](/images/practice/eocrc-locrc/cohort-pca.png)](/images/practice/eocrc-locrc/cohort-pca.png)

*The current experiment's `outputs/cohort-pca.png`. The left and right PCA panels were calculated separately, so their coordinate values should not be compared directly.*

In this run, 21,650 genes passed the low-expression filter in EOCRC and 24,116 passed it in LOCRC. The numbers satisfying both `padj < 0.05` and `|log2FC| >= 1` were 4,570 and 8,256, respectively.

In EOCRC, PC1 explained 24.73% of the variance and PC2 explained 14.58%. In LOCRC, PC1 explained 34.03% and PC2 explained 9.46%. Tumor and Normal separate mainly along PC1 in both cohorts, but EOCRC also contains samples near the opposite condition and samples far from the rest.

PCA is not given the condition as an answer. It finds the major directions of variation among samples. The separation therefore means that tumor status corresponds to a large part of the overall expression difference. Imperfect separation does not mean the analysis failed, and distant samples should not be removed before checking their quality and clinical information.

### 3-3. Put gene-level changes on the same coordinates

An `inner join` of the two result tables by Ensembl gene ID places changes for the same gene in one row. The x-axis is the LOCRC Tumor-versus-Normal `log2FC`, and the y-axis is the same value from EOCRC.

### Step 3 prompt: compare EOCRC and LOCRC log2FC

```text
Keep the previous analysis and add only a gene-level comparison of log2FoldChange between the two cohorts.

1. Take baseMean, log2FoldChange, lfcSE, pvalue, and padj from the eocrc and locrc result tables and inner join them by Ensembl gene ID.
2. Prefix each column with eo_ or lo_, then calculate lfc_difference = eo_log2FoldChange - lo_log2FoldChange.
3. Exclude rows with missing or infinite log2FoldChange from the plot, but preserve the original values in the result file.
4. Draw a scatter plot with LOCRC log2FoldChange on the x-axis and EOCRC log2FoldChange on the y-axis. Each point represents one shared gene.
5. Draw a dashed y=x line and reference lines at x=0 and y=0. A point above the dashed line has a more positive change in EOCRC, while a point below it has a more negative change.
6. Calculate the Pearson correlation between the two log2FoldChange columns.
7. Save the comparison table to outputs/eo-vs-lo-lfc-comparison.csv and the figure to outputs/lfc-comparison.png at 180 dpi.
8. Print the number of shared genes, the correlation coefficient, and the meanings of the four quadrants and the y=x line.
9. Run uv run python analysis.py and explain in English how closely the overall point cloud follows the diagonal, then stop.

Do not expand a high correlation into a claim that the two cohorts are identical or that EOCRC-specific genes do not exist.
```

[![Gene-level scatter plot with LOCRC log2FC on the x-axis and EOCRC log2FC on the y-axis. Most points follow the diagonal.](/images/practice/eocrc-locrc/lfc-comparison.png)](/images/practice/eocrc-locrc/lfc-comparison.png)

*The current experiment's `outputs/lfc-comparison.png`. Gray points are shared genes, and orange points are the eight genes reported in the original paper.*

There were 21,420 shared genes with finite `log2FC`, and the Pearson correlation between the two changes was **0.930**. The point cloud running from the bottom left to the top right means that the large increases and decreases in expression in cancer tissue are generally shared by the two age groups.

Read the graph as follows:

| Position | Meaning |
| --- | --- |
| Upper right | Increased in Tumor in both cohorts |
| Lower left | Decreased in Tumor in both cohorts |
| Upper left | Decreased in LOCRC, increased in EOCRC |
| Lower right | Increased in LOCRC, decreased in EOCRC |
| Near `y=x` | Similar changes in the two cohorts |
| Above `y=x` | EOCRC change is more positive than LOCRC change |
| Below `y=x` | EOCRC change is more negative than LOCRC change |

Points far from the diagonal are a starting point for finding EOCRC candidates. Do not select candidates by distance alone. Check each cohort's `baseMean`, `padj`, and patient-level values as well.

### 3-4. Recheck the paper's 8-gene signature

The original paper narrowed genes that change in EOCRC but barely change in LOCRC through set-based conditions and presented an eight-gene signature: `ALDOB`, `FBXL16`, `IL1RN`, `MSLN`, `RAC3`, `SLC38A11`, `WBSCR27`, and `WNT11`.

The official symbol for `WBSCR27` has since changed to `METTL27`, but the figure preserves the name used at the time so it can be compared with the paper. Its Ensembl gene ID is `ENSG00000165171`.

### Step 4 prompt: compare the paper's criteria and eight genes

```text
Keep the previous analysis and add only the paper's EOCRC candidate rule and a check of the 8-gene signature.

1. Define the following symbol-to-Ensembl-gene-ID mapping explicitly in the code. Preserve the name WBSCR27 for comparison with the original paper, but add a comment that its current symbol is METTL27.
   ALDOB=ENSG00000136872, FBXL16=ENSG00000127585,
   IL1RN=ENSG00000136689, MSLN=ENSG00000102854,
   RAC3=ENSG00000169750, SLC38A11=ENSG00000169507,
   WBSCR27=ENSG00000165171, WNT11=ENSG00000085741
2. For this exercise, define EOCRC significant regardless of direction as eo_padj < 0.05, abs(eo_log2FoldChange) > 1, and eo_baseMean > 50.
3. Define LOCRC unchanged as lo_padj > 0.2, abs(lo_log2FoldChange) < 0.7, and lo_baseMean > 50.
4. Mark rows satisfying both conditions above and abs(lfc_difference) > 1.5 as paper_like_candidate.
5. On the Step 3 scatter plot, overlay paper_like_candidate in teal and the eight published genes in orange, and label them.
6. For each of the eight genes, show the LOCRC and EOCRC log2FoldChange as two points connected by a line, and save the figure to outputs/signature-eight-genes.png.
7. Print a table containing the two log2FoldChange values, two padj values, lfc_difference, and paper_like_candidate status for each of the eight genes.
8. Print the number and names of candidates passing the current package and filter settings.
9. Run uv run python analysis.py and check whether the paper's eight genes are selected again exactly. If they differ, do not hide the discrepancy. Explain that differences in filters, DESeq2 implementation, versions, and preprocessing can change the candidate set, then stop.

Do not secretly alter the current criteria to force the eight genes from the paper to pass.
```

[![Comparison of LOCRC and EOCRC log2FoldChange for the 8-gene EOCRC signature reported in the paper](/images/practice/eocrc-locrc/signature-eight-genes.png)](/images/practice/eocrc-locrc/signature-eight-genes.png)

*The current experiment's `outputs/signature-eight-genes.png`. Gray represents the LOCRC Tumor-versus-Normal change, and orange represents the EOCRC change.*

For all eight genes, the EOCRC point lies farther to the right than the LOCRC point, or farther to the left in the case of `SLC38A11`. In other words, the changes show the same direction reported by the original paper, with **a tendency to be more pronounced in EOCRC**.

| gene | EOCRC log2FC | EOCRC padj | LOCRC log2FC | LOCRC padj | EO−LO difference |
| --- | ---: | ---: | ---: | ---: | ---: |
| `ALDOB` | 1.87 | 0.00024 | 0.72 | 0.133 | 1.16 |
| `FBXL16` | 0.81 | 0.0084 | -0.23 | 0.442 | 1.05 |
| `IL1RN` | 1.18 | 0.0025 | 0.45 | 0.259 | 0.73 |
| `MSLN` | 1.85 | 0.000025 | 0.09 | 0.851 | 1.76 |
| `RAC3` | 0.58 | 0.031 | -0.09 | 0.738 | 0.68 |
| `SLC38A11` | -1.55 | 0.000042 | -0.61 | 0.182 | -0.94 |
| `WBSCR27` | 0.96 | 0.0018 | 0.26 | 0.351 | 0.70 |
| `WNT11` | 1.13 | 0.0045 | 0.36 | 0.353 | 0.77 |

Under the unified set conditions in this exercise, however, **only `MSLN`** passed every filter. The other seven showed similar directions but failed to cross one boundary, such as `|EO−LO| > 1.5`, EOCRC `|log2FC| > 1`, or LOCRC `|log2FC| < 0.7`.

This does not mean the paper was wrong, nor that the code should be changed to produce the desired answer. The original paper used R DESeq2 results, a separate paired Wilcoxon test, and the selection procedure recorded in the paper. This page uses PyDESeq2 0.5.4 and the stated low-expression filter. Threshold-based candidate lists are sensitive to implementation, version, filtering, and missing-value handling, so **the full values and reproduction environment must be recorded together**.

### 3-5. Expand the eight genes by patient

One cohort-level `log2FC` summarizes changes across 21 or 22 patients. To check whether the same direction repeats across patients behind the average difference, expand normalized counts back into patient pairs.

### Step 5 prompt: patient-level change heatmap for the eight genes

```text
Keep the previous analysis and add only a patient-level paired-change heatmap for the 8-gene signature.

1. Retrieve normed_counts from each cohort's DeseqDataSet.
2. For each of the eight genes, pivot values from the same patient by condition and calculate paired_log2_delta = log2(Tumor + 1) - log2(Normal + 1).
3. Save cohort, patient, symbol, and paired_log2_delta to outputs/signature-patient-deltas.csv.
4. Draw EOCRC and LOCRC heatmaps in upper and lower panels, with the eight genes as rows and patients as columns.
5. Use the same color range centered at zero. Show positive values in red, negative values in blue, and values near zero in white.
6. Save the figure to outputs/signature-patient-heatmap.png at 180 dpi.
7. For each gene, print the number of patients with positive and negative values and the median paired_log2_delta in each cohort.
8. Run uv run python analysis.py and inspect the figure. Explain in English an example that shows the cohort average and patient-level heterogeneity together, then stop.

Do not describe paired_log2_delta as the same value as the DESeq2 model's log2FoldChange.
```

[![Heatmap of Tumor−Normal paired log2 changes for the paper's eight genes across 21 EOCRC and 22 LOCRC patients](/images/practice/eocrc-locrc/signature-patient-heatmap.png)](/images/practice/eocrc-locrc/signature-patient-heatmap.png)

*The current experiment's `outputs/signature-patient-heatmap.png`. Each column is one patient. Red means higher expression in that patient's Tumor, and blue means higher expression in Normal.*

The `MSLN` row in EOCRC contains many red cells, showing that the positive cohort-level `log2FC` repeats across multiple patients. The same row still contains white and blue cells, however, and the magnitude and direction of change for the other genes also differ among patients. This is why the eight names should not be read as one fixed EOCRC type.

The heatmap's `paired_log2_delta` directly transforms and subtracts two normalized counts from one patient. DESeq2's `log2FoldChange` is a model estimate that uses all patients, the count distribution, and dispersion together, so the two values are not exactly the same.

## 4. What else is needed for the rest of the paper's analysis?

The original paper did not stop at comparing gene counts. It further narrowed the differences between EOCRC and LOCRC with these analyses.

| Analysis | Question asked by the paper | Possible with HTSeq gene counts alone? |
| --- | --- | --- |
| GO enrichment of 48 EOCRC-specific genes | Which functions contain the candidates? | Possible with the candidate list |
| TCGA survival analysis | Is the 8-gene score associated with survival? | No, additional clinical survival data is required |
| xCell deconvolution | Do cell-composition scores differ in bulk tissue? | Possible, but requires separate signatures and validation |
| Alternative splicing | Are there EOCRC-specific transcript isoforms? | No, alignment files and exon/junction information are required |
| Confirming HOTAIRM1 splicing | Does a particular splice event differ by age group? | No, one gene-count row contains no isoform information |
| Neoantigen prediction | Can an abnormal splice peptide bind MHC-I? | No, transcript sequence, translation, HLA type, and binding prediction are required |

[HTSeq gene counts](/en/reference/htseq-counts/) combine multiple transcripts of one gene into a single value. Even if you have the total count for `HOTAIRM1`, you cannot tell which exon was included or whether a new peptide was produced. Do not extend the result of a gene-level DEG exercise into conclusions about splicing or neoantigens.

In addition to the 8-gene signature, the original paper reported 48 genes specifically altered in EOCRC, immune-related functions, differences in RNA splicing factors, and seven potential neoantigen events. These findings are starting points for a later exercise, not conclusions reconfirmed by the current count comparison alone.

### Questions to answer after viewing all the graphs

1. If the two count matrices are combined directly, what overlaps with age group?
2. Why use `~ patient + condition` in each cohort?
3. What do positions above and below `y=x` mean in the log2FC scatter plot?
4. Why can EOCRC candidates exist even when the correlation of log2FC across all genes is high?
5. Why can you not conclude that the two effects differ solely because a gene is significant in EOCRC and nonsignificant in LOCRC?
6. How should you report that only one of the original paper's eight genes passes the current set filter?
7. Why can alternative splicing and neoantigens not be validated directly with gene-level HTSeq counts?

## 5. Sources and reproduction information

- [GSE196006: 42 samples from 21 EOCRC patients](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE196006)
- [GSE251845: 44 samples from 22 LOCRC patients](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE251845)
- [Identification of differentially expressed genes and splicing events in early-onset colorectal cancer](https://pmc.ncbi.nlm.nih.gov/articles/PMC11047122/)
- [PyDESeq2 reference](/en/reference/pydeseq2/)
- [Reading DESeq2 and PyDESeq2 result tables](/en/reference/deseq2-results/)

For reproduction, record the download dates of both GEO files, Python and PyDESeq2 versions, low-expression filter, `design`, contrast, candidate thresholds, and date of gene-symbol mapping. For genes with renamed symbols such as `WBSCR27`, save the Ensembl gene ID together with the symbol.