Capstone - Expression Atlas

The finished program reads a raw count matrix and a sample table and writes a differential expression summary: it validates every count row, measures each library's size and detection rate, decides which libraries are fit to use, converts counts to counts per million, compares the two condition groups with a log2 fold change that survives a group mean of zero, sets aside the genes with too little evidence to rank, orders what remains by effect size, and ends with a verdict on the run as a whole.

The point of the exercise is the order of those steps. Normalising before you have decided which libraries to keep, or ranking before you have removed the low-count genes, produces a table that looks exactly as convincing and is wrong.

The data

The count matrix. The sandbox mounts counts.tsv, which is the shape you are parsing:

gene	control_1	control_2	treated_1	treated_2
TP53	1204	1150	2310	2402
BRCA1	455	470	441	462
EGFR	3021	2988	1502	1477
KRAS	880	912	1760	1690
MYC	150	143	900	945
GAPDH	15200	15410	15100	15330

The first line is a header of gene followed by the sample names, in the order the columns appear. Every later line is a gene name followed by one integer count per sample. Real matrices are not this clean: a row can carry the wrong number of fields, a value can be blank or negative or a float where an integer was promised, and a gene name can appear twice because two annotation sources were concatenated. The second appearance is not an error; it represents a second set of counts for the same gene, requiring summation across duplicate rows.

The sample table. A sample name and its condition group, one pair per line. Only control and treated take part in the comparison. Any other label, and any column of the matrix with no row in this table at all, is a library that is present in the file and absent from the analysis: a spike-in, a failed run, a sample from a different experiment left in the matrix by accident.

Note that the two tables agree on names and nothing else. The matrix decides the column order; the sample table decides what each column means.

Requirements

1. Read the two sections and the header. Treat a count row as malformed, count it, and drop it when its field count differs from the header's, when the gene name is empty, or when any value is not a non-negative integer.
2. Sum the rows of a gene that appears more than once, position by position, before anything else is computed. The library sizes depend on it.
3. Compute each library's size, the total of its column over the retained genes, and its detection rate, the fraction of retained genes with a non-zero count in that column. With no genes retained the detection rate is 0.0 rather than an error.
4. Give a library the status unassigned when its name is missing from the sample table or carries a group other than the two named above. An unassigned library is reported and takes no further part.
5. Give any other library the status flagged when its size is below half the median size of all libraries, or its detection rate is below 0.5, and ok otherwise. Only ok libraries are retained for the comparison. The median is taken over every library including the unassigned ones, because it is a description of the run, not of the analysis.
6. Convert each count to counts per million against its own library's size. A library of size zero yields 0.0 throughout and must not divide.
7. Set aside a gene whose total raw count across the retained libraries is below 10. Count these separately: they are not a failure, they are genes the experiment has no power to speak about.
8. For every other gene take the mean CPM within each group and report log2(mean treated + 1) - log2(mean control + 1). The added one is what keeps a gene absent from one group finite, and it is why the fold change of a very lowly expressed gene is pulled toward zero rather than toward infinity.
9. Rank by absolute log fold change descending, then by gene name ascending.
10. Emit the library count and the malformed count, one line per library in header order, the retained gene count and the low-count count, the ranked rows, and a verdict: qc: pass when nothing is flagged, qc: fail (<n> flagged) otherwise. Say no comparison possible when either group has no retained library, and no genes ranked when both groups do but nothing survives the low-count filter.

Milestones

Run these in order. They give you the QC measurements and the shape of the fold change, and leave the assembly to you.

Measure the libraries.

genes: 6 libraries: 4
control_1 20910 1.0
control_2 21073 1.0
treated_1 22013 1.0
treated_2 22306 1.0

The sizes differ by a few per cent here, which is normal. A library at a tenth of the others is not a library with less signal; it is a library whose CPM values are dominated by a handful of counts, and requirement 5 exists to keep it out of the means.

Normalise and compare, with the guard in place.

TP53 56076.2 106311.0 0.923
BRCA1 22031.7 20372.8 -0.113
EGFR 143134.6 67223.9 -1.09
KRAS 42681.6 77858.6 0.867
MYC 6979.8 41625.1 2.576
GAPDH 729096.2 686608.7 -0.087

counts / sizes divides each column by its own total because pandas aligns the Series on the column index. Doing it by hand, per gene, is where the wrong axis creeps in.

Apply the flagging rule to a library that ought to fail it.

median size: 21073
control_1 20910 ok
control_2 21073 ok
treated_1 22013 ok
treated_2 22306 ok
treated_3 40 flagged

Adding one small library moves the median very little, which is the property you want: a mean would be dragged down by the outlier until the outlier no longer looked like one.

What to watch for

• Summing the duplicate genes after computing the library sizes. The size is the total of the column, and the total changes when two rows become one.
• Dropping the unassigned libraries before taking the median. Requirement 5 takes the median over every library; excluding some of them changes which of the rest are flagged.
• Computing CPM against the retained libraries' mean size instead of each library's own size. That is not CPM, and it reintroduces exactly the depth difference the normalisation was there to remove.
• Dividing by a library size of zero. It happens whenever every count in a column is zero, which is what a failed run looks like.
• Leaving out the +1. A gene with a group mean of zero gives an infinite or undefined fold change, and it will be at the top of your ranking.
• Ranking on the signed fold change. The strongest results are at both ends; sorting without the absolute value silently discards half of them.
• Filtering low-count genes on the raw column total across all libraries rather than across the retained ones, so a gene supported only by a flagged library survives.

Exercise

Write the whole summary in one function, using the standard library, pandas or NumPy as you prefer. Every stage of the brief is present: library QC, CPM normalisation, the guarded log fold change, the low-count filter and the ranking.

Write atlas(text), returning a list of output lines. The input arrives in two sections separated by a line containing exactly --: the sample table (sample name and group, tab separated, where the only groups that count are control and treated) and the count matrix, whose first line is a header of gene followed by the sample names and whose later lines are a gene name followed by one integer per sample. Ignore blank lines. A count row is malformed, and is counted and dropped, when its field count differs from the header's, when its gene name is empty, or when any value is not a non-negative integer. A gene name that appears more than once is not an error: sum its rows position by position.

Compute each library's size, the total of its column across the retained genes, and its detection rate, the fraction of retained genes with a non-zero count in that column, or 0.0 when no genes remain. A library whose sample name is missing from the sample table, or carries any group other than the two above, has status unassigned and takes no part in the comparison. Any other library has status flagged when its size is below half the median size of all libraries or its detection rate is below 0.5, and ok otherwise. Only ok libraries are retained for the comparison.

Convert each count to counts per million using its own library's size, so a library of size zero gives 0.0 throughout. A gene whose total raw count across the retained libraries is below 10 is low-count: count it and drop it from the ranking. For every other gene take the mean CPM within each group and report log2(mean treated + 1) - log2(mean control + 1), the added one being what keeps a group mean of zero finite.

Print libraries=<n>, the number of columns in the header, then malformed=<n>. Follow those with one line per library in header order, holding the literal lib, the sample name, the group or unassigned, the size, the detection rate to exactly three decimal places, and the status, separated by tabs. Then print genes=<n>, the number of distinct genes retained, and low_count=<n>. Then one line per ranked gene, ordered by the absolute log fold change descending and then by gene name ascending, holding the literal gene, the name, the log fold change, the mean control CPM and the mean treated CPM, each to exactly three decimal places, separated by tabs. Print no comparison possible in place of those lines when either group has no retained library, and no genes ranked when both groups do but no gene survives the low-count filter. Close with qc: pass when no library is flagged, or qc: fail (<n> flagged) otherwise.