How to use this course (self-study with a .md file)
You can read the course on this page. You can also give the .md file below (a plain-text teaching file that opens on a computer or a smartphone) to a conversational AI, which will then act as your tutor and take you through the course one session at a time.
- Download the file below.
- Give the file to a conversational AI (attach it, or paste in its contents) and ask: “Please start Session 1.”
- Do one session a day, about 15 minutes each (up to 10 sessions).
- Ask the AI about anything unclear as you go. When a session involves a program, ask it to “explain it line by line” and write the program together.
- After the final session, send us your feedback (subject line below).
AI answers can be wrong. Check numbers and facts against the course text on this page. Do not enter personal information or other people's data into the AI.
If you use this course, please write to the address on our Contact page with the subject “Course feedback: Read which genes are at work in single cells”, and tell us how many sessions you completed, what was unclear, what worked well, and how long each session took.
Read which genes are at work in single cells (10 sessions, 15 minutes each, v1.2)
Goal of this course: Use public data to see, cell by cell, which genes are switched on.
This course is for undergraduates who have completed high-school biology and mathematics. We use public data in which cells were measured one at a time while hormone-producing cells (such as the insulin-producing β cells) were developing in the pancreas of mouse embryos. With your own hands you will open the data, decide which cells to discard, give the cells names, and go as far as checking "is that difference real?". At the end, you will write one page on "what I can now judge".
If you would like to learn in English, first tell the AI: "Please teach this course in English."
Version: v1.2 (November 2026). This version was revised based on feedback from people who actually used it (see "What was changed in this version" below).
Introduction (how to use this course)
What this is
This is a teaching material published for students' self-study by the Division of Neurophysiology, Kobe University Graduate School of Medicine. It is a conversational course of 10 sessions, 15 minutes each. In our laboratory, medical students and researchers followed this path "15 minutes a day"; we have rebuilt it so that it can be completed using public data only.
Who it is for
- High-school biology (DNA, genes, cell division) and mathematics (proportions, logarithms) are enough.
- It does not matter if you have never written a program. In each session you copy and run 10–25 lines of Python (a language for writing out the steps of a calculation), and change only 1–2 places to your own values.
- The point where beginners get stuck most is losing track of "what in the table is connected to what, and where". Session 1 deals with exactly that, using a small table.
What you need
- A conversational AI (one that answers your questions in a chat). You load this file into it and use it.
- Somewhere to run Python. A free notebook environment that runs in the browser (no installation needed) is fine. In the first cell of Session 1 you install the tools (about 1 minute). If you run it on your own computer, type pip install scanpy igraph leidenalg statsmodels once.
- You download one public dataset (Session 2, about 52MB). From Session 3 onward you reuse it.
- Notes (on paper or on screen). At the end of each session you write a 1–3 line "record". These records become the material for the final page.
How to proceed
- Load this entire file into the conversational AI.
- Type "Please start Session 1 of this course."
- Answer the AI's questions in your own words (it does not matter if you are wrong). Copy and run the code yourself, and paste the results into the AI.
- Stop after 15 minutes. The next day, load the same file again and type "Please start Session 2 of this course." - If you do not use an AI, read each session from the top, answer the check questions yourself, and then open "Answers and reasoning". - The course is written so that you can follow the storyline without reading the code. Reading only the "Explanation" and the "Expected results" of each session, you can still follow what was checked and what was learned. The code is for those who want to reproduce the same results with their own hands. - In some notebook environments, files disappear when you close them. In that case, running the code of Sessions 2, 3 and 4 in a row will rebuild them (about 1 minute of computation in total).
Overview of the 10 sessions:
| Session | Title | Hands-on task |
|---|---|---|
| 1 | Building a small container — the three rooms and their order | Assemble a 3-cell × 4-gene table yourself (fill in the blanks) |
| 2 | Opening real public data | Load a table of 3,696 cells and check the counts by a different route |
| 3 | Which cells to discard | Write down your quality thresholds first, then look at the breakdown of discarded cells |
| 4 | A map of cells and their names | Name the clusters from the marker table alone, before looking at the answers |
| 5 | Picking cells with a set of genes | Look at the table for all cell types and split the errors into "boundary" and "specificity" |
| 6 | Can shallow data be read? | Write down your pass/fail threshold first, then judge six practice samples |
| 7 | Does "not detected" mean "absent"? | Compare with the number expected from the sequencing depth |
| 8 | What counts as "one"? | Tell apart females and males within a single sample, and build a pseudobulk |
| 9 | Is that "significant" result real? | Check by shuffling the labels (overall / within depth strata) |
| 10 | Placing cells on a map | Using the same table as for the markers, see where cell types absent from the map end up |
Glossary (please read through it once at the start)
In sessions where many new terms appear at once (Sessions 3–5 and 7–9), the numbers from this table are listed at the start of the session as "terms that appear for the first time in this session". If you lose track, come back here.
| No. | Term | Meaning (short) | First appears in session |
|---|---|---|---|
| G1 | Gene expression | How much of a gene's copy (mRNA: messenger RNA) is being made | 1 |
| G2 | NGS | Next-generation sequencing. A device that reads the sequences of tens of millions of small fragments at once | 1 |
| G3 | Single-cell RNA sequencing (scRNA-seq) | A method that separates cells one by one and counts the mRNA in each | 1 |
| G4 | Molecule counts (UMI counts) | The number of mRNA molecules, counted using a tag that makes copies of the same molecule count as one | 2 |
| G5 | Container (AnnData) / the three rooms | One object bundling the count table X (cells × genes), the cell description table obs, and the gene description table var | 1 |
| G6 | Linkage by order | The convention that the cell in row 1 of X and row 1 of obs are the same cell. They are tied together by labels (row names) | 1 |
| G7 | Sequencing depth (depth) | The total number of molecules counted in one cell. The deeper, the more low-abundance genes are detected | 3 |
| G8 | Quality control (QC) | Finding broken cells and empty droplets by their numbers and removing them | 3 |
| G9 | Mitochondrial fraction | The fraction of all molecules taken up by mitochondrial genes. Tends to be high in broken cells | 3 |
| G10 | Normalization and log transformation | Scaling each cell's depth to the same total (10,000) and taming large values with log(1 + x) | 4 |
| G11 | Highly variable genes | Genes whose values differ greatly between cells. Only these are used to build the map | 4 |
| G12 | Principal component analysis (PCA) | Compressing the information of many genes into a small number of "directions" | 4 |
| G13 | Clustering and resolution | Grouping similar cells together. Raising the resolution produces many small clusters | 4 |
| G14 | UMAP | A method to draw the closeness of cells as a 2D picture. Distances and sizes of the islands are not reliable | 4 |
| G15 | Marker gene | A gene used especially in one cell type that serves as a clue to identify it | 4 |
| G16 | Lineage and stage | The path along which the same cells develop, and each stage along the way (ductal cells → Ngn3-low progenitors → Ngn3-high progenitors → cells just before producing hormones → split into four types: α, β, δ and ε) | 4 |
| G17 | Adjacent stage | On the same lineage, the stage immediately before or immediately after the stage of interest. Where the path branches, all of the immediately following branches are adjacent (the branches are not adjacent to each other). Cells that merely share a marker gene but are far apart on the lineage are not "adjacent" | 5 |
| G18 | Boundary problem / specificity problem | Adjacent stages mixing in = boundary problem (hard to avoid because development is continuous). Non-adjacent cells mixing in = specificity problem (the chosen genes are not unique to that cell type) | 5 |
| G19 | Score and threshold | The score given to each cell by the chosen gene set, and the value decided as "pick everything above this" | 5 |
| G20 | Positive predictive value (PPV) and sensitivity | The proportion of picked cells that were correct / the proportion of true cells that were picked | 5 |
| G21 | Base rate (baseline proportion) | What fraction of all cells the type you want to pick makes up. A clue for setting the threshold | 5 |
| G22 | Downsampling | Randomly discarding molecules to make the data deliberately shallower. Used when comparing at matched depth | 6 |
| G23 | Dropout | Failing to count even one molecule although mRNA is actually present. More frequent the shallower the data | 7 |
| G24 | Number expected from sequencing depth | If the gene were made at the same rate as in the reference group, in how many cells it should be detected at this depth | 7 |
| G25 | Independent experimental unit | One unit taken without influencing the others (usually one animal). The "one case" when comparing | 8 |
| G26 | Pseudobulk | A table in which the molecules of the cells from the same unit are summed, giving one column per unit | 8 |
| G27 | Log fold change | The ratio of the values of two groups expressed as log2. +1 means twofold, −1 means half | 8 |
| G28 | p-value and correction | "The probability of a difference this large or larger arising by chance when there is actually no difference", and the calculation that corrects the inflation when many genes are examined at once | 9 |
| G29 | Pseudoreplication | Counting cells from the same animal as if they had been sampled separately. Makes differences look exaggerated | 9 |
| G30 | Confounding | When two factors (e.g. group and measurement batch) always change together, so that their effects cannot be separated | 9 |
| G31 | Permutation check (permutation test) | Randomly shuffling the group labels to directly generate how large a difference arises by chance | 9 |
| G32 | Permutation within strata | Permuting only among cells of similar depth. Prevents differences in depth from appearing as differences between groups | 9 |
| G33 | Placing on a map (label transfer) | Placing your own cells on a map built from other cells and borrowing the name of the nearest cluster | 10 |
| G34 | GEO | The public repository of gene expression data run by the US National Center for Biotechnology Information (NCBI) | 2 |
| G35 | Circularity (double dipping: using the same data twice) | Comparing clusters using the very genes that were used to split them. A difference always appears, so it is not evidence | 9 |
What was changed in this version (v1.1 → v1.2)
- The dataset was replaced with a single public dataset measuring endocrine cells (hormone-producing cells) during their development in the mouse pancreas. Downloading one file of about 52MB once lets you complete all 10 sessions.
- The four fixes made in response to comments on the previous version have been kept as they are: ① next to the answer-checking table, the marker values of each cluster are shown (Sessions 4 and 10); ② a table of how often the chosen genes are detected in every cell type is shown first (Session 5); ③ "adjacent stage" was defined in one sentence (glossary G17 and Session 5); ④ the glossary was placed at the beginning so that sessions with many terms can refer to it by number.
- Session 10 was changed into a form where you yourself try out "where do cell types absent from the map end up".
- The text was written so that the storyline connects through the "Explanation" and "Expected results" of each session alone, so that people who do not read the code can follow it.
- The format in which every line of the final page must include a "because" is unchanged from the previous version.
Cautions
- The AI can make mistakes. Check important points against textbooks or original sources (papers, official documentation). The "Expected results" in this course may also vary slightly depending on tool versions and random numbers.
- Do not enter your own or other people's personal information, clinical information, or other people's data into the AI. The only things you may paste into the AI in this course are results you yourself produced from public data.
- This course is not intended as medical advice. The pancreas and insulin appear, but there is no discussion of the diagnosis or treatment of disease.
Request for feedback (required)
This course is revised every time based on the voices of those who used it. If you used it, please be sure to send us feedback.
Send to: the address on the "Contact" page of the laboratory homepage Subject: "Course feedback: NGS gene expression analysis" Four things we would like you to write: 1. How many sessions you completed 2. What was hard to understand (which table or which term, in which session) 3. What was good 4. How long one session took
(If any) please also attach figures or examples you would like added, and your final page. Feedback from people who stopped partway is just as useful. The four fixes in this course came from feedback that was actually sent to us.
Terms of use
- Use for personal study and in classes is free.
- If you distribute a modified version, please contact the laboratory.
- For the public data used, cite the paper and the tool paper, and follow GEO's terms of use (Session 2 and "Data and references used" at the end).
- The content is current as of November 2026.
Instructions to the AI (tutor role)
This section contains the rules that the AI that has loaded this file must follow.
- Keep each session to 15 minutes. Respect the "Time guide" at the start of each session. If it looks likely to run over, stop at the next natural break, have the learner write the record, and move the rest to the next session.
- Ask only one question at a time and wait for the learner's answer. Keep each turn within 200 Japanese characters (within 100 words in English). Go longer only when asked "more detail".
- Do not give the answer first. For thresholds (quality thresholds, score thresholds) and gene sets, first have the learner write down their own values, and only then present the AI's proposal. Use the collapsed "Answers and reasoning" only after the learner has answered or when they say "show me the answer".
- Have numbers checked by a different route. When the AI states a number (threshold, count, fold change), prompt the learner to compare it with their own output or to check it by another calculation (.shape and len(), table and figure, etc.). Biological claims too should be checked against genes in the data, not against the AI's words.
- Do not fabricate. Do not invent dataset IDs, gene names, p-values or papers that are not in this file. When recommending a paper, give its DOI or PubMed ID and add "please check that it exists and what it says".
- Do not say "I can see" about things you have not seen. If the learner's output has not been pasted, ask them to paste it.
- The learner copies and runs the code. Do not extensively rewrite each session's code or add steps. When an error is pasted, explain one by one which line did what. Have the learner compare their results with the "Expected results" of each session (small differences due to versions are fine; if there is a large difference, review the thresholds and settings of the previous session together). For learners who do not run the code, read the "Explanation" and "Expected results" tables together and have them answer the same questions.
- At the start of each session, ask one review question based on the key points of the previous session. If there was something the learner got wrong in the previous session, give it priority. In sessions where three or more new terms appear, first show the relevant glossary numbers and briefly restate each one on the spot.
- Show each session's "Figure" as is, in its box, at the point in the explanation where it belongs.
- From Session 5 onward, when the discussion turns to wrongly picked cells, have the learner judge using the definition in glossary G17 (only the stages immediately before or after on the same lineage are "adjacent"; after a branch point, all branches are adjacent). Do not call cells that merely share a marker but are far apart on the lineage "adjacent".
- In each session's "What this session does, in one sentence", check whether the learner's sentence includes "which room (X, obs, var) and which column"; if not, prompt them to add just one (practice in writing instructions that ask an AI to do work).
- At the end of each session, give three key points and two mini check questions. Then help tidy up the 1–3 line record, and finally ask how long it took.
- Do not move on to the next session on your own. Start only when the learner types "Please start Session N of this course."
- If the learner writes a medical consultation, their own or others' personal information, or other people's data, decline to handle it and return to general explanation.
- At the very start of Session 1, ask one question at a time: "What do you want to be able to do through this course?", "Have you ever run Python?", "Have you ever handled rows and columns of a table in a spreadsheet?". Adjust your analogies and pace to the answers.
- After Session 10 (or when the learner says they will stop partway), help finish the completion page (with a "because" in every line), ask the four feedback items one at a time, and hand them over in a form that can be pasted directly into an e-mail. The destination is the "Contact" page of the laboratory homepage.
- Talk about tables, thresholds and reasons, not about the learner's personality. Give one good point and at most two points to fix.
- If the learner asks for English, proceed in English. Give the English for technical terms alongside.
Session 1 Building a small container — the three rooms and their "order"
Question: Where in the count table is the information about cell type (for example, "β cell") written? ↓ It is nowhere in the count table. It is written in a separate table, and the two tables are connected only by their "order".
Time guide: reading 5 min, typing 8 min, record 2 min
Goals of this session
- Be able to say that the container for single-cell data consists of three rooms (X, obs, var).
- Be able to say that the rooms are connected by their order (the row labels).
- See with your own eyes that when the order is broken, no error appears and a different value is returned.
Explanation
The cells of the body carry almost the same DNA. Yet a β cell of the pancreas differs from a skin cell because it differs in which genes it copies and how much — that is, in gene expression (how much of a gene's copy, mRNA, is being made). With NGS (next-generation sequencing: a device that reads the sequences of tens of millions of small fragments at once), this mRNA can be counted. In single-cell RNA sequencing (a method that separates cells one by one and counts the mRNA in each), you get, for each individual cell, "how many molecules of which gene were present".
What comes out is, in the end, one big table. The rows are cells, the columns are genes, and the entries are the numbers of molecules counted. But this table alone does not tell you "what kind of cell the cell in this row is". So the Python container (AnnData) has three rooms.
The count table itself is X, the cell description table is obs (one row per cell), and the gene description table is var (one row per gene). The only thing connecting the three rooms is their order (the row labels). There is nothing but the convention that the cell in row 1 of X and row 1 of obs are the same cell.
If this convention is broken, nobody tells you. Today we check this with a small table of 3 cells × 4 genes. What people handling this kind of data for the first time struggled with most was "where and how the attribute data are connected". Settling this point first makes the remaining nine sessions much easier.
Figure
Hands-on task (Python, fill in the blanks)
In the first cell of the notebook, install the tools (about 1 minute; it is fine as long as ERROR does not appear).
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 1 with me and explain it line by line.”
Fill in the three ____ in the next cell. Ins1 is the insulin gene and a marker of β cells, Gcg is the glucagon gene and a marker of α cells, and Sox9 is a marker of ductal cells (the ducts that carry pancreatic juice). Decide each cell's type from the gene with the largest count in that cell.
(There is another program here. As above, write it together with a conversational AI and ask it to explain each line as you run it.)
Expected results: AnnData object with n_obs × n_vars = 3 × 4, the retrieved value [[8.]], and finally OK: the order matches.
Next, deliberately get the order wrong. Change the cell_type line to ["alpha cell", "beta cell", "duct"] and run it again. Before running, predict what will be returned, then run it.
What this session does, in one sentence: Write today's task in one sentence, as an instruction asking an AI to do it. Example form: "toy is an AnnData with rows = cells. Put the cell types into the cell_type column of obs, in the order of the rows of X." Naming which room and which column makes it more likely that you get the intended answer in one go.
Record (example): Session 1: three rooms = X, obs, var. The only link is the order. When the order was broken, [[0.]] was returned and no error appeared.
Check questions
- What does one row of obs represent?
- When you got the order wrong, why did no error appear?
- If you cut the genes in half (if the rows of var decreased), would the number of rows of obs change?
Answers and reasoning (for the AI; learners open after answering)
- The blanks in the task: ["beta cell", "alpha cell", "duct"] (the maximum is Ins1 = 5 for cellA, Gcg = 8 for cellB, and Sox9 = 7 for cellC).
- The version with the wrong order: [[0.]] is returned, and finally NG: the order is shifted... appears (confirmed in the author's environment). This is because Gcg in cellA, which is labeled "alpha cell", is 0.
- Q1: One cell. Descriptions of that cell, such as its type, are lined up.
- Q2: Because the program simply runs trusting the convention that "row 1 and row 1 are the same cell", and does not know whether the contents are correct. That is why we put in a step to check with our own eyes (today's OK/NG cell).
- Q3: No. Reducing the genes (columns) does not break the order of the cells (rows). So the cell_type in obs can be used as it is.
Summary
- The container has three rooms: the count table X, the cell descriptions obs, and the gene descriptions var.
- The only link between the rooms is their order. Even if it breaks, no error appears.
- Instructions to the AI get across when you name "which room and which column".
Session 2 Opening real public data
Question: If the authors of a paper write "we measured 3,696 cells", how can you check that the table you loaded really matches that? ↓ Count the same number twice by different routes. If they agree, that is evidence that you loaded it correctly.
Time guide: reading 4 min, typing 9 min, record 2 min (you download about 52MB; around 1 minute)
Goals of this session
- Be able to obtain the public data file and describe the contents of the three rooms.
- Be able to check the number of cells by counting again by a different route.
- Be able to remove unneeded shelves and move the raw counts to a backup shelf.
Explanation
Today we use public data that measured hormone-producing cells during their development in the pancreas of mouse embryos (embryonic day 15.5 after fertilization) (Bastidas-Ponce et al. 2019, Development 146(12):dev173849, doi:10.1242/dev.173849). The original data are deposited in GEO (the public repository of gene expression data run by the US National Center for Biotechnology Information) as GSE132188. Today we use a single file (about 52MB) in which the developers of an analysis tool (scVelo) packaged the day-15.5 portion for practice.
In the pancreas, progenitor cells that have started to use a gene called Ngn3 (Neurogenin 3) arise from the family of ductal cells (the ducts that carry pancreatic juice), pass through a stage just before producing hormones, and split into four types of hormone cells: insulin β cells, glucagon α cells, somatostatin δ cells and ghrelin ε cells. The authors of this file have given each individual cell the name of its stage.
The contents are a table of molecule counts (UMI counts: numbers counted using a tag that makes copies of the same molecule count as one). This file is in the container (AnnData) format from the start, so just loading it gives you all three rooms. However, it also comes with the authors' map and many shelves for other analyses. Today we keep only the authors' names, in a column called cell_type, and remove the other extras. We keep the names as "the answers" until Session 4 and do not look at them beforehand.
The idea of the check is simple. The number of rows of the cell description table obs, the number of rows of the count table X, and the sum of the cell counts per author name are counts of the same cells by different routes. If they do not match, something was missed or mixed up somewhere during loading.
Figure
Hands-on task (Python, one blank)
Fill in ____ and then run. Hint: of the three rooms from Session 1, which room contains the raw counts?
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 2 with me and explain it line by line.”
Expected results (confirmed in the author's environment. When loading, a warning about moving where things are stored may appear, but you need not worry about it):
Run it with the AI and you will see numbers like these (example)
There are eight author names. Ductal is ductal cells; Ngn3 low EP and Ngn3 high EP are endocrine progenitors using a little / a lot of Ngn3 (EP stands for endocrine progenitor); Pre-endocrine is cells just before producing hormones; and Alpha, Beta, Delta and Epsilon are the four types of hormone cells.
What this session does, in one sentence: Try writing it naming the room and the column, for example: "Count the number of cells for each name in the cell_type column of obs in adata, and check whether the total equals the number of rows of obs."
Record (example): Session 2: pancreas day-15.5 data, 3696 cells × 27998 genes. 8 names (ductal 916 to delta 70). Counting by 3 routes agreed.
Check questions
- In 3696 × 27998, which are the cells and which are the genes? Answer without looking at the code.
- How many cells does the smallest type have? What trouble might small types cause in later sessions?
- Look at the total number of molecules in the first three cells with (a line of program). Why do you think the counts differ between cells?
Answers and reasoning (for the AI)
- The blank: X ((a line of program)). In this file, X contains the integers as counted (confirmed in the author's environment that it contains integers only).
- Q1: 3,696 are cells (rows) and 27,998 are genes (columns).
- Q2: Delta, with 70 cells. In clusters with few cells, "not detected" easily happens by chance (Session 7), and the cluster proportions also tend to fluctuate.
- Q3: In the author's environment, 4954, 7071 and 4070. The total number of molecules counted (depth) differs depending on cell size, how damaged the cell is, and how well it was read. This is used in Session 3.
Summary
- Handle public data as a set: the paper, the accession number (GSE132188), and the terms of use.
- Count numbers again by a different route and check that they agree.
- Keep the authors' names as "the answers" and do not look at them beforehand.
Data source and conditions: The original data are GSE132188 (NCBI GEO). When using them, cite the paper by Bastidas-Ponce et al. 2019 and the paper of the tool that distributes the practice file (Bergen et al. 2020, Nature Biotechnology 38(12):1408–1414, doi:10.1038/s41587-020-0591-3), and follow GEO's terms of use (NCBI places no restrictions, but submitters may claim rights).
Session 3: Which cells to discard: draw the quality threshold first
Question: Can't we just decide the threshold for discarding poor-quality cells after looking at the data? ↓ Write the threshold first, then look at the "breakdown" of the discarded cells. If one cell type alone disappears in large numbers, what you are discarding is not noise but biology.
Time guide: reading 3 min · typing 10 min · record 2 min
Terms introduced in this session: G7 sequencing depth · G8 quality control · G9 mitochondrial fraction (check the glossary before going on)
Goals of this session
- You can name the three per-cell numbers (sequencing depth, number of genes detected, mitochondrial fraction).
- You can write your own threshold and count how many cells it removes.
- You can look at the fraction of removed cells per cell type and decide whether to revise the threshold.
Explanation
When cells are separated one by one, some cells break during the process, and some droplets contain no cell at all. Finding them is quality control (QC: identifying broken cells and empty droplets by numbers and removing them). Three numbers are used: sequencing depth (the total number of molecules counted in one cell), the number of genes detected, and the mitochondrial fraction (the proportion of all molecules that come from mitochondrial genes). When the cell membrane ruptures, cytoplasmic mRNA leaks out, but mRNA inside the mitochondria tends to remain, so this fraction rises.
What matters here is the order. First, write your own threshold. Next, count how many cells that threshold removes. Finally, look at which cell types the removed cells are concentrated in. If some cell type naturally has little mRNA, a uniform threshold will cut away that type alone.
This dataset holds one more clue. Look at the minimum sequencing depth and the maximum mitochondrial fraction. There is not a single very shallow cell, nor a single cell with high mitochondrial content. This means the authors had already discarded them before depositing the data.
Figure
Hands-on task (Python)
Write your own values into MIN_GENES and MAX_MT before looking at the results, then run it (the example below uses 1000 and 5).
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 3 with me and explain it line by line.”
Expected result (author's environment, thresholds 1000 and 5): the minimum sequencing depth (total_counts) is 3020, the minimum number of genes detected is 1473, and the maximum mitochondrial fraction is 4.3%. Cells discarded: 0 of 3696.
Next, change only MIN_GENES to 2000 and run it again (do not run the final write_h5ad line; if you continue to the next session with the 1000-and-5 version, it will be easier to compare with the later "Expected result" sections). In the author's environment, 374 cells were discarded, and the breakdown was Pre-endocrine (cells just before they begin producing hormones) 23.8%, Ngn3 high EP 15.6%, and Alpha 10.8%. The earlier a cell is in its development, the fewer genes are detected, and the more easily it is cut away by a uniform threshold.
What this session does, in one sentence: try writing something like "Select cells using n_genes_by_counts and pct_counts_mt in obs, and show the removed cells as fractions per cell_type."
Record (example): Session 3: Threshold = genes ≥1000, mitochondria ≤5% → 0 cells. Minimum depth 3020, mitochondria max 4.3% = authors already discarded cells. With genes 2000, 23.8% of Pre-endocrine disappears.
Check questions
- What can you infer from the fact that the minimum sequencing depth is 3020 and the maximum mitochondrial fraction is 4.3%?
- When the gene threshold was raised to 2000, nearly a quarter of the cells just before hormone production disappeared. What would you do with the threshold? Give one reason.
- Answer in one sentence the opinion that "we can just let the AI decide the threshold."
Answers and reasoning (for the AI)
- Q1: Before depositing the data, the authors had already discarded shallow cells and broken cells. That is why almost nothing is removed even when you draw your own threshold.
- Q2: There is no single correct answer. Examples: do not use 2000 (because what disappears is not broken cells but cells at a stage that naturally has few genes) / use a separate threshold for each cell type. What matters is having "decided after looking at the breakdown of the removed cells" and writing the reason in the log.
- Q3: Example: ask for the AI's proposal after writing your own value, and decide with data — the number and breakdown of removed cells — which is better.
Summary
- There are three quality numbers: sequencing depth, number of genes detected, and mitochondrial fraction.
- Write the threshold first. Always look at the fraction of removed cells per cell type.
- In public data, someone may already have drawn a threshold once.
Session 4: Build a cell map and name the clusters before seeing the answer
Question: Why shouldn't we overlay the authors' labels (the answer) on the computationally derived clusters from the start? ↓ If you look first, you get no practice reading with your own eyes. Name the clusters using only the marker table, then check your answers against the right-hand end of the same table.
Time guide: reading 3 min · typing 10 min (about 30 s of computation) · record 2 min
Terms introduced in this session: G10 normalization and log transformation · G11 highly variable genes · G12 principal component analysis · G13 clustering · G14 UMAP · G15 marker genes · G16 lineage and stage (this session has many terms; open the glossary first)
Goals of this session
- You can state the standard workflow (normalization → highly variable genes → principal component analysis → nearest-neighbour graph → map → clustering), together with the purpose of each step.
- You can name clusters from a table of mean marker gene values.
- You can state, in order, the stages through which hormone-producing cells develop in the pancreas.
Explanation
Because sequencing depth differs from cell to cell, we first perform normalization (scaling each cell to the same total of 10,000 and taming large values with log(1 + x)). Next, we use only the 2,000 "highly variable genes" that differ most between cells, and compress them into 30 directions by principal component analysis (compressing the information of many genes into a small number of directions). From there, each cell is connected to its nearest cells, and the cells are grouped (clustering: grouping similar cells together). UMAP (a method that draws the proximity of cells as a two-dimensional picture) is a picture for viewing that result by eye. The distances between islands and the sizes of islands are unreliable, so do not draw conclusions from the picture alone.
Names are assigned with marker genes (genes particularly used by a given cell type). The cells in this dataset develop by passing through the following path in order. This path is called the lineage.
Ductal cells (high Sox9) → progenitors that have started using Ngn3 (high Neurog3, the gene encoding Ngn3; dividing cells also have high Mki67) → cells just before hormone production (high Fev) → split into 4 types (β = Ins1, α = Gcg, δ = Sst, ε = Ghrl).
The answer-checking table is a single table with the authors' labels placed at the right-hand end of the marker value table. Going back and forth between two tables makes it easy to mix up cluster numbers (a previous user told us, "I got the names wrong while switching between tables").
Figure
Hands-on task (Python)
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 4 with me and explain it line by line.”
When table (1) appears, before looking at (2), write down the name you assign next to each cluster number.
Expected result (author's environment, scanpy 1.11.5; the number of clusters and their numbering vary by version): 10 clusters.
Clusters in which one marker stands out are named correctly as is (cluster 7 with Ins1 at 2.7 = β, fraction 0.99; cluster 9 with Ghrl at 5.5 = ε). On the other hand, one of the authors' labels, Ngn3 low EP, was not the most frequent label in any cluster. These cells are mixed into ductal clusters 1 and 2 (this is why the fraction for cluster 1 is as low as 0.71). Cluster 6, with high Sst, also had a fraction of only 0.50, and most of the rest were cells just before hormone production.
What this session does, in one sentence: try writing something like "For each leiden cluster in obs, tabulate the mean of the markers genes, and add at the right-hand end the most frequent cell_type label and its fraction."
Record (example): Session 4: 10 clusters. Named using markers only … of … clusters matched. Ngn3 low EP has no cluster of its own and is mixed into ductal clusters. δ cluster has fraction 0.50 (mixed with pre-endocrine cells).
Check questions
- Why do you think Ngn3 low EP did not have a cluster of its own and was mixed into the ductal clusters?
- What name would you give the cluster with a fraction of 0.50 (cluster 6)?
- What do you think happens to the number of clusters if resolution=0.5 is changed to 1.0?
Answers and reasoning (for the AI)
- Q1: The transition from ductal cells to cells that have started using Ngn3 is continuous. Cells that have just started still share most of their genes with ductal cells, so when grouped by proximity they gather on the ductal side. Placing labeled boxes on a boundary makes the boxes mix with each other.
- Q2: There is no single correct answer. Example: "δ cells and cells just before becoming δ." If you add the marker values (Sst 2.7, Fev 1.3) and the fraction, readers will understand that the cluster is mixed. Do not force it into a single name.
- Q3: The number of clusters increases and they split into smaller ones (the author has not tried this; the learner may be allowed to try it).
Summary
- Each step — normalization → highly variable genes → principal component analysis → neighbours → clusters — has a purpose.
- Name clusters first from the marker table, then check the answers at the right-hand end of the same table.
- The cells in this dataset develop along a branching path (lineage) from ductal cells to 4 types of hormone-producing cells.
Session 5: Calling cells with a gene set: separate "boundary" from "specificity"
Question: To call the progenitors that strongly use Ngn3 (Ngn3 high EP), is it enough to choose 5 familiar textbook "endocrine cell markers"? ↓ If you line up the 5 chosen genes across all cell types, you will find genes that are detected far more often in the stage right after the target than in the target stage itself. Looking at them one at a time, you will not notice.
Time guide: reading 3 min · typing 9 min · record 3 min
Terms introduced in this session: G17 adjacent stage · G18 boundary problem / specificity problem · G19 score and threshold · G20 positive predictive value and sensitivity · G21 base rate
Goals of this session
- You can see, in a single table, in what percentage of cells of every type each chosen gene is detected.
- You can describe the calling result with two numbers: positive predictive value and sensitivity.
- You can classify wrongly called cells into "adjacent stage" (boundary problem) and "non-adjacent cells" (specificity problem).
Explanation
The target is Ngn3 high EP. Choose 5 genes, give each cell a score for that set (a score for how high the chosen genes are), and treat cells above the threshold as "called."
The calling result is described with two numbers: positive predictive value (PPV; the proportion of called cells that were correct) and sensitivity (the proportion of true cells that were called). Improving only one of them is easy. Setting the threshold very high raises the positive predictive value, and setting it low raises sensitivity. A clue for the threshold is the base rate (the proportion of all cells that belong to the type you want to call). Ngn3 high EP makes up about 17% of all cells in this dataset, so starting with a threshold that "calls the top 17%" makes it easier to balance the two numbers.
Wrongly called cells are considered in two kinds. Here, "adjacent" has a clear rule.
An adjacent stage is the stage immediately before or immediately after the target stage on the same lineage. Where the path branches, all of the branches immediately after are adjacent. The 4 types after branching (β and α, etc.) are neither before nor after one another, so they are not adjacent to each other.
For Ngn3 high EP, only Ngn3 low EP immediately before and Pre-endocrine immediately after are adjacent. Mixing in of adjacent stages is a boundary problem. Because development is continuous, it is unavoidable to some extent. On the other hand, mixing in of stages two or more steps away (Ductal, or post-branching Beta, Alpha, etc.) is a specificity problem. It means the chosen genes are not unique to that stage, so you choose the genes again.
This was a point on which a previous user hesitated: "Can cells that share markers be called adjacent?" The answer is "no." Whether a stage is adjacent is decided by its position on the lineage, not by marker similarity. β cells share several markers with the progenitors that strongly use Ngn3, but there is one stage (Pre-endocrine) between them, so they are not adjacent. If they are mixed in, it is a specificity problem.
Figure
Hands-on task (Python)
Replace my_genes with your own 5 genes before running (the example is a well-known set of "endocrine cell markers").
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 5 with me and explain it line by line.”
Expected result (author's environment, the example 5 genes, threshold 83):
Table (1) (excerpt, %):
Chga, Fev, and Neurod1 are detected far more often in Pre-endocrine, immediately after the target Ngn3 high EP, and in β and α further on, than in the target itself.
The "true" cells in (2) are the cells the authors labeled Ngn3 high EP (642 cells). Because the true cells make up about 17% of all cells (642 ÷ 3,696), the threshold was set at the top 17% of scores. Called 629 cells, true 642 cells: positive predictive value 0.17 (about 17% of the 629 called cells were true) and sensitivity 0.16 (about 16% of the 642 true cells were called). In (3), the largest group of errors was Pre-endocrine (416 cells), which is adjacent (boundary), but non-adjacent cells such as Beta (73 cells), Alpha (29 cells), and Epsilon (7 cells) were also mixed in.
Next, look at table (1) and choose genes again that are frequent in the target stage and rare elsewhere. With the set the author tried, ["Neurog3", "Ppp1r14a", "Cbfa2t3", "Cdc14b", "Btbd17"] (chosen by detection rate in the target stage and difference from other cell types), positive predictive value was 0.85 and sensitivity 0.83; the errors were almost all from adjacent stages (Pre-endocrine 79 cells, Ngn3 low EP 6 cells), and non-adjacent cells totalled 8. The remaining errors are a boundary problem, so the proper course is to write "the boundary is continuous" in the record rather than to change the genes.
What this session does, in one sentence: try writing something like "Using layers["counts"], tabulate, for every cell type by cell_type in obs, the percentage of cells with at least 1 molecule of each gene in my_genes."
Record (example): Session 5: Example 5 genes = PPV 0.17, sensitivity 0.16. Chga, Fev, Neurod1 are more frequent in Pre-endocrine and β (noticed in the all-types table). Reselected = 0.85/0.83; remaining errors are adjacent stages = boundary.
Check questions
- State the rule for "adjacent stage" in one sentence, in your own words.
- If β cells are called, is that a boundary problem or a specificity problem?
- Why is it not enough to look only at the "most confusable cell type" for each gene individually?
Answers and reasoning (for the AI)
- Q1: On the same lineage, only the stage immediately before or immediately after the target stage is adjacent (all branches right after a branch point are adjacent). Cells that merely share markers but are far apart on the lineage are not adjacent.
- Q2: A specificity problem. β cells are two steps after Ngn3 high EP (with Pre-endocrine in between). Some of the chosen genes (Chga, Fev, etc.) are also used in β cells.
- Q3: If 3 of the 5 genes each leak a little into the same cell type (here, Pre-endocrine), this is not conspicuous while you look at genes one at a time, but in the set score the leaks add up into a large error. A table of all cell types × all genes reveals it at a glance.
Summary
- After choosing genes, look in a single table at what percentage of every cell type each gene is detected in.
- Describe the calling result with two numbers: positive predictive value and sensitivity. A clue for the threshold is the base rate.
- Adjacent stage = immediately before or after on the same lineage. If adjacent stages are mixed in, it is a boundary problem; if anything else is mixed in, it is a specificity problem.
Session 6: Can shallow data be read? Write the pass/fail threshold first
Question: When very little mRNA can be captured from a single cell, at how few molecules does it become impossible to tell "what cell this is"? ↓ You can measure this by deliberately discarding molecules from healthy cells to make them shallow. And if you write the pass/fail threshold before measuring, you cut off the temptation to move the threshold after seeing the results.
Time guide: reading 3 min · typing 9 min · record 3 min
Terms introduced in this session: G22 downsampling (also review G7 sequencing depth)
Goals of this session
- You can describe, using a table, how the number of genes detected decreases as sequencing depth is reduced.
- You can explain why data of different sequencing depths should be compared after matching depth.
- You can write the pass/fail threshold before the data and judge the 6 practice libraries.
Explanation
In experiments that collect and read only a very small number of cells, each piece of data (library: one sample submitted for sequencing) contains very few molecules. By downsampling deep cells (randomly discarding molecules to make them deliberately shallow), you can create "how the same cell would look if it were shallower."
There are two things to look at. One is the number of genes detected. The other is the proportion of cells in which that cell type's marker (here, Pdx1 for β cells) is detected. The shallower the data, the more genes that are truly present fail to be counted even once. Therefore "the number of genes detected" is a number dragged along by sequencing depth, and data of different depths cannot be compared directly. When comparing, downsample everything to the same depth before counting.
In research, before sequencing deeply (which costs money and time), one sometimes does a pilot with a little sequencing to decide whether to proceed. The most important thing then is to write the threshold before the results come out. If you decide the threshold after seeing the results, it almost always becomes "a threshold that passes what you want to pass." Today, you will try this once with "6 practice libraries": 6 β cells from the public data deliberately reduced to different depths.
Figure
Hands-on task (Python)
Write your own values into MIN_DEPTH (minimum number of molecules) and MIN_GENES (minimum number of genes after matching depth) before running (the example uses 1000 and 650).
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 6 with me and explain it line by line.”
Expected result (author's environment, thresholds 1000 and 650):
Run it with the AI and you will see numbers like these (example)
Even for the same β cells, downsampling to 1000 molecules reduced the proportion of cells in which the marker Pdx1 was detected from 95.3% to 51.9%. Some libraries, like L3, "have sufficient depth but fall slightly short of genes after matching." If you now feel like lowering the threshold to 645, that is the temptation to "move the threshold after seeing the results."
What this session does, in one sentence: "Using layers["counts"], downsample only the cells whose cell_type in obs is Beta to 2000, 1000, 500, and 200 molecules per cell, and tabulate the median number of genes detected."
Record (example): Session 6: Genes 2528→1183→694→392→174 (all→200 molecules). Pdx1 95.3%→13.9%. Wrote threshold (1000 molecules, 650 genes) first; 2 of 6 passed. L3 failed on gene count.
Check questions
- If Pdx1 is not detected in a β cell with 200 molecules, can you say "this cell does not use Pdx1"?
- Why must you not directly compare "the number of genes detected" between data with depth 3000 and data with depth 900?
- With your threshold, how many of the 6 practice libraries passed? State the reason for your threshold in one sentence, without using the results.
Answers and reasoning (for the AI)
- Q1: No. At full depth it is detected in 95.3% of cells, but at 200 molecules in 13.9%. A 0 may mean "could not be counted" rather than "absent" (this will be checked with numbers in the next session).
- Q2: The number of genes detected is largely determined by depth. Downsample to the same depth before counting.
- Q3: The answer differs from person to person. The reason should ideally come from the upper half of today's table, e.g. "in table (1), about 700 genes are detected at 1000 molecules = cell types can still be distinguished," or from how many errors you can tolerate. "Because I want L3 to pass" is not a reason.
Summary
- Making data shallower hides genes that are truly present. A 0 does not necessarily mean "absent."
- Compare data of different depths after matching them to the same depth.
- Write the pass/fail threshold before the data.
Session 7: Does "not detected" mean "absent"? — the number expected from sequencing depth
Question: In the cells just before they start making hormones (Pre-endocrine), very few cells have the mature β-cell gene G6pc2 detected. But these precursor cells are sequenced slightly less deeply, so isn't this just dropout? ↓ If you calculate "how many cells would have it detected at this depth if it were made at the same proportion as in the reference group," you can tell in numbers whether the shallower depth explains it.
Time guide: reading 4 min · typing 6 min · hand calculation 3 min · record 2 min
Terms introduced in this session: G23 dropout · G24 number expected from sequencing depth
Aims of this session
- You can state that dropout (failing to count even one molecule of a gene that is actually present) arises from sequencing depth.
- You can calculate the number expected from sequencing depth and compare it with the number actually detected.
- When there are few cells, you can calculate by hand how often "0 cells" occurs purely by chance.
Explanation
Suppose a gene accounts for a proportion f of all molecules in the reference group (here, β cells). In a cell with sequencing depth N, the probability that at least one molecule of that gene is detected is approximately 1 − e^(−N×f) (N×f is the number of molecules expected in that cell). Summing this over all cells of the group being examined (here, Pre-endocrine) gives "the number of cells in which the gene is expected to be detected, given their sequencing depth."
If the actual number is about the same as expected, the difference can be explained by depth. If the actual number is far below expectation, shallowness cannot explain it, which gives grounds to conclude that the gene is genuinely low at that stage. The trick is to place this comparison side by side with a gene that every cell uses (here, Actb, a gene for the cell's cytoskeleton). If expected and actual agree for Actb, that confirms the comparison itself is properly correcting for depth.
One more caution, for when there are few cells. Even a gene that is truly detected in 5% of cells may easily go undetected in every cell if there are only 40 cells. Multiplying 0.95 by itself 40 times gives about 0.13; that is, "0 cells" happens about once in 8 times. Treat a "0" in a small group with suspicion first.
Figure
Hands-on task (Python and hand calculation)
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 7 with me and explain it line by line.”
Expected result (in the author's environment):
Run it with the AI and you will see numbers like these (example)
(Depending on the environment, the expected number may be printed with a long decimal, such as 296.600006.)
Hand calculation (a calculator is fine): Suppose a group has 40 cells and the true detection proportion of a gene is 5%. The probability that not a single molecule is detected in any of the 40 cells is 0.95 to the 40th power. What is it? And for 20 cells? For 200 cells?
What to do in this session, in one sentence: "Using layers["counts"], take the cells whose cell_type in obs is Beta as the reference, calculate in how many Pre-endocrine cells G6pc2 should be detected given their sequencing depth, and put that next to the actual number."
Log (example): Session 7: G6pc2 in Pre-endocrine: 4 detected vs 297 expected = shallowness cannot explain it. Actb is 591 vs 592, the same. Nkx6-1 is 495 vs 570 = mostly explained by depth, only slightly lower. With 40 cells at 5%, 0 cells occurs about 13% of the time.
Check questions
- What can you conclude from "4 vs 297" for G6pc2? Answer by comparing it with "591 vs 592" for Actb.
- How do you read "495 vs 570" for Nkx6-1?
- From your hand calculation, what do you think you should do before writing "this gene is 0 at this stage" for a group of 40 cells?
Answers and reasoning (for the AI)
- Hand calculation: 0.95^40 ≈ 0.13, 0.95^20 ≈ 0.36, 0.95^200 ≈ 0.00004 (3.5×10⁻⁵). The fewer the cells, the more easily "0 cells" happens by chance.
- Q1: G6pc2 is about 1/74 of the number that depth could explain = grounds to conclude it is genuinely low at the preceding stage (not yet in use). Actb matches expectation, so the comparison itself is properly correcting for depth.
- Q2: About 87% of expected. Most of the difference is explained by depth, and the remainder is small. Saying neither "used about as much as in β cells" nor "not used," but writing something like "somewhat lower," fits the numbers.
- Q3: State the number of cells and calculate the expected number (or the probability above). If possible, confirm with another sample or dataset that has more cells at the same stage.
- Supplementary (if asked): The reference group is where you place the assumption "if it were made in the same way as in the group being examined." If you use as reference a cell type that normally uses the gene differently, the expectation itself loses its meaning.
Summary
- 0 does not necessarily mean "absent." Compare with the number expected from sequencing depth.
- Place a gene that every cell uses alongside, to check the comparison itself.
- For a "0" in a small group, first calculate how likely it is by chance with (1 − proportion)^number of cells.
Session 8: What counts as "one"? — females and males within a single sample
Question: There are 591 β cells and 481 α cells. Can we say we have as many as 1,072 data points to compare? ↓ One "data point" in a comparison is one animal. These data come from a single sample in which the pancreases of several embryos were pooled and measured together, and there is no record of which cell came from which embryo.
Time guide: reading 4 min · typing 8 min · record 3 min
Terms introduced in this session: G25 independent experimental unit · G26 pseudobulk · G27 log fold change
Aims of this session
- You can state that the independent experimental unit (usually one animal) is the "one data point" in a comparison.
- You can confirm from sex genes that females and males are mixed in a single sample.
- You can build a pseudobulk table and check whether the same answer is obtained within females and within males.
Explanation
Cells taken from the same animal share everything: genotype, the mother's body in which they developed, the time of dissection, and sample processing. So no matter how many hundreds of cells you collect, you have only "one animal's worth" of information about differences between animals. One data point in a comparison is an independent experimental unit (one item taken without influencing the others).
Because these data have no per-cell record of "which embryo," the units cannot be counted correctly. But there is one clue: sex genes. Xist is the gene that silences one of the X chromosomes in female cells, and Ddx3y, Eif2s3y, Kdm5d and Uty are Y-chromosome genes. Cells in which only Xist is detected can be classed as female-like, and cells in which only Y genes are detected as male-like (cells in which neither or both are detected are classed as "undetermined," because of shallow depth or mixing). If there are both female-like and male-like cells, we can say that at least two embryos are mixed.
So we sum the molecules of female-like cells together and of male-like cells together, one column each. This is pseudobulk (a table with one column per unit). We calculate the difference between β and α separately within females and within males, and see whether the two agree. A difference that agrees is weak confirmation that it is at least not "a quirk of a single animal."
Always add the number of cells in each group and the number of units to tables and figures. Even if there are hundreds of cells, if they come from one or two animals, the reader needs to know that.
Figure
Hands-on task (Python)
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 8 with me and explain it line by line.”
Expected result (in the author's environment): In table (1), every cell type had both female-like and male-like cells (for example, β had 340 female-like, 134 male-like and 117 undetermined cells). At least two embryos are mixed. The pseudobulk in (2) has 4 columns (β·female / β·male / α·female / α·male). Comparing β and α by log fold change, the agreement (correlation) between within-females and within-males is 0.68. The correlation is how well "the β-vs-α fold change calculated within females" and "the same fold change calculated within males" line up when placed side by side gene by gene; 1 means they line up perfectly for every gene, and 0 means they are unrelated. 0.68 means roughly "they line up well for many genes, but some genes with small differences diverge." The genes with the largest differences went in the same direction within females and within males.
Run it with the AI and you will see numbers like these (example)
(Excerpt. + means higher in β, − means higher in α. +1 is 2-fold; +5.8 is about 56-fold (2 to the 5.8th power). Ins1 and Ins2 are insulin genes, and Gcg is the glucagon gene.)
What to do in this session, in one sentence: "Using Xist and the 4 Y-chromosome genes in layers["counts"], create a sex column in obs, and make a table of cell counts for each cell_type."
Log (example): Session 8: Both female-like and male-like cells in the single sample = at least 2 animals. No per-embryo record = no statistical test possible. β vs α agrees within females and within males (correlation 0.68; Ins and Gcg in the same direction).
Check questions
- When there is no per-embryo record, why is it difficult to attach a p-value to the difference between β and α?
- If you compared only "female-like β cells" with "male-like α cells," what kinds of genes would likely appear at the top of the fold-change ranking (assuming you had not removed the genes used for assignment)?
- If you wanted to compare by embryo, how should the data be collected again (how many animals, what to record)?
Answers and reasoning (for the AI)
- Q1: Because the units cannot be counted, the variability within a group (differences between embryos) cannot be measured. Without variability, you cannot judge whether a difference is within the range of chance. A list of fold changes can be produced, but a statistical test cannot be performed. The two columns, female and male, are not "a substitute for units"; they are only a clue that "there are at least 2 animals" (we do not even know how many animals the female-like cells come from).
- Q2: Xist and the Y-chromosome genes (Ddx3y, Eif2s3y, etc.). Because the two groups being compared differ not only in cell type but also in sex. When sex overlaps (is confounded) with what you want to compare, sex leaks into the top genes. Reading the identity of the top genes lets you notice the mix-up.
- Q3: Example = at least 3 animals per group, with each embryo as a separate sample, or measure with a tag on each cell that identifies which embryo it came from. Record sex and balance it across both groups. There is no single correct answer, but it is acceptable as long as it includes "the number of units" and "knowing which cell belongs to which unit."
Summary
- The "one data point" in a comparison is one animal. The number of cells is not the number of data points.
- Even in data without a record of units, sex genes tell you "at least how many" animals there are.
- Checking whether the same answer appears within subsets, such as within females and within males, lets you suspect a quirk of a single animal.
Session 9: Is that "significant" result real? — checking by permutation
Question: When female-like β cells and male-like β cells were compared with a statistical test treating each cell as one data point, "significant" genes appeared. When the labels are randomly shuffled, there are 0. So are the genes that appeared "differences between female and male β cells"? ↓ Look at the names of the genes that appeared and you will know the answer. And the permutation check protects you only "when cells are exchangeable." Differences between animals and differences in sequencing depth have to be checked with a different approach.
Time guide: reading 4 min · typing 8 min · record 3 min
Terms introduced in this session: G28 p-value · correction · G29 pseudoreplication · G30 confounding · G31 permutation check · G32 permutation within strata · G35 circularity (This session has many terms. Please open the glossary first.)
Aims of this session
- You can state what "significant genes vs 0 after permutation" shows and what it does not show.
- You can explain with examples "circularity" (comparing groups using the genes used to define the groups) and pseudoreplication (counting cells from the same animal as if they had been taken separately).
- When a difference in sequencing depth looks like a difference between groups, you can check it by permuting within depth strata.
Explanation
A p-value is "the probability that a difference this large or larger would arise by chance when there is actually no difference." When thousands of genes are tested at once, chance "hits" get mixed in, so you count them after correction (a calculation that removes the inflation). Today's statistical test compares by ranks (the Wilcoxon test).
In Cell 1, we compare the female-like β cells and male-like β cells from Session 8. The "significant" genes that come out are the very genes used to define the groups. This is circularity (using the same data twice). If you compare groups using what you used to make the groups, a difference will always appear. It is not a discovery; it is merely a confirmation of how the groups were split.
The permutation check (permutation test) directly generates how large a difference can arise by chance, by randomly shuffling the group labels. If shuffling gives 0, the statistical test tool itself is not broken. However, permutation mixes cells one at a time. If the female-like cells came from one female and the male-like cells from one male, "the difference between sexes" and "the difference between those two individuals" cannot be separated. Permutation cannot see this. This mistake, counting hundreds of cells as hundreds of independent pieces of evidence, is pseudoreplication.
Another pitfall is sequencing depth. In Cell 2, we compare the "detection proportion" of a gene that every cell uses (Gapdh, the gene for a sugar-degrading enzyme) between β cells and the cells just before they start making hormones. The precursor cells are shallower, and that alone lowers the detection proportion. Permuting over all cells generates chance with even the depth difference erased, so the apparent difference looks large. So we divide the cells into 10 strata by depth and shuffle the labels only within the same stratum. In other words, we compare the actual difference with the chance obtained by comparing cells of similar depth.
Finally, a word about sequencing runs. These data are one sample measured in one run, so the issue does not arise here, but if a sample rich in β cells had been measured in the first run and a sample rich in α cells in the second run, the difference in cell type and the difference in run would overlap completely: confounding (two factors always changing together so they cannot be separated). No statistics applied afterwards can separate them. Only planning before data collection can.
Figure
Hands-on task (Python, 2 cells)
Cell 1 (a cell-level statistical test, and shuffling the labels):
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 9 with me and explain it line by line.”
Cell 2 (permuting within depth strata):
(There is another program here. As above, write it together with a conversational AI and ask it to explain each line as you run it.)
Expected result (in the author's environment): In Cell 1, the cell-level statistical test gives 3 genes (Ddx3y, Eif2s3y, Xist), and 0 after shuffling. All 3 were the very genes used to split female from male. In Cell 2, median depth was 4876.5 for the precursor cells and 6172.0 for β. The difference in the proportion of cells with Gapdh detected was 0.145. Permuting over all cells, the mean permuted difference was 0.0 and p = 0.001. Next, when the cells were sorted by depth and divided into 10 groups (strata), and the labels were shuffled only within the same stratum, the mean permuted difference rose to 0.08. This means that even with random labels, a difference of 0.08 arises from the depth difference alone (p remains 0.001. This is the smallest value obtainable with 1000 permutations, and means that none of the 1000 permutations reached the actual difference). Of the 0.145 difference, about half (0.08) is an apparent difference produced by depth alone, and the remaining half persists even after matching depth.
One step further (optional, 3 min): You can check the same thing with an equation that explains "detected or not (0 or 1)" by "group" and "log depth" (logistic regression). Write this as a continuation of Cell 2.
(There is another program here. As above, write it together with a conversational AI and ask it to explain each line as you run it.)
In the author's environment, the group coefficient was 2.19 in the equation without depth and 1.49 in the equation with depth (p below 0.001 in both). Adding depth shrinks the group effect to about two thirds, but does not make it disappear. This is the same reading as the permutation result.
What to do in this session, in one sentence: "For cells of two cell_type values in obs, divide the depth from layers["counts"] into 10 strata, shuffle the cell_type labels 1000 times only within strata, and build the chance distribution of the difference in Gapdh detection proportion."
Log (example, 3 lines): Session 9: Female vs male β: 3 significant = the very genes used for assignment (circularity); 0 after shuffling. Permutation does not see differences between individuals = pseudoreplication. Of the 0.145 Gapdh difference, about 0.08 is an apparent effect of depth (permutation within strata); the rest persists even after matching depth.
Check questions
- What does "0 after shuffling" guarantee, and what does it not guarantee?
- Why can't the "3 significant genes" in Cell 1 be called a discovery?
- When planning an experiment that compares two groups, write one rule about sequencing runs (sample-processing batches) that should be decided before data collection.
Answers and reasoning (for the AI)
- Q1: Guarantees = if cells are exchangeable, the statistical test does not produce too many chance "hits." Does not guarantee = apparent differences due to differences between animals (individual, mother, time of dissection, etc.) or differences in depth. Permutation that mixes cells cannot see these.
- Q2: Because the groups were compared using the very genes used to define them (Xist and the Y genes) (circularity). A difference is bound to appear; it is merely confirmation that the split worked as intended. The fact that no other significant genes appeared is also not evidence that "there is no difference between sexes" (the number of units is unknown, and the number of cells is small).
- Q3: Example = include both groups in roughly equal numbers in every run (and, if possible, balance sex as well). Do not put all of group A in the first run and all of group B in the second.
Summary
- The permutation check is a powerful tool, but what it protects against changes depending on what you permute.
- Do not compare groups using what you used to define them (circularity). Always read the names of the genes that come out.
- Check depth differences by permuting within depth strata. Prevent confounding between runs and groups by planning before data collection.
Session 10: Placing cells on a map — where do cell types that are not on the map end up?
Question: If you place your own cells on a cell map that someone else built, names are assigned automatically. Can you say that those names are correct? ↓ A placement method always places each cell in some cluster. It never says "I don't know." Today, you will deliberately remove one stage from the map and see with your own eyes where the cells of that stage end up.
Time guide: reading 3 min, typing 10 min (about 30 s of computation), record 2 min
Terms introduced in this session: G33 placing on a map (label transfer)
Goals of this session
- You can say what placing on a map (label transfer) actually does.
- You can read the table of counts of cells that arrived alongside the marker values of the clusters they arrived at.
- You can name one assumption that the placement method relies on, and describe a situation in which it breaks down.
Explanation
In research, when naming your own cells, it is common to place them on a cell map (a reference) built in another study and borrow the name of the nearest cluster. Today, we try this within a single dataset. We build a map using only the female-like cells, and place the male-like cells on it. The cells that build the map and the cells that are placed come from different embryos.
However, we deliberately remove Pre-endocrine (cells just before they start making hormones) from the map. In a research setting, this corresponds to the situation "the reference map did not contain cells of the very stage I want to study." Where will cells of a type that is not on the map end up?
The map here is not the UMAP picture, but a computed space built by the same procedure as in Session 4 (highly variable genes → principal component analysis (PCA)). Each cell used to build the map has its position determined by 30 numbers (principal components), and these cells have been divided into clusters. The cells to be placed are also converted into 30 numbers using the same genes and the same reduction; then, in that space, the nearest map cells are found, and their cluster number is borrowed. Because this method "always places each cell somewhere," cell types that are not on the map are also forced into some cluster.
That is why the table for checking the answer is a single table with the counts of cells that arrived, and to its right the mean marker values of each cluster. As in Session 4, moving back and forth between two tables leads to mixing up cluster numbers (people who actually used this course asked for this format). Finally, we also line up the mean marker values on the side of the placed cells, and check whether they match the clusters where those cells arrived.
Figure
Hands-on task (Python)
The program for this session is not shown on this page. Give the file for the AI (above) to a conversational AI and ask: “Please write the program for Session 10 with me and explain it line by line.”
Expected result (author's environment, scanpy 1.11.5; cluster numbers and proportions may change with the version or random seed): 1,766 map cells and 886 cells placed. The map has 9 clusters. Table of counts arrived and markers (excerpt):
Run it with the AI and you will see numbers like these (example)
("Clu" is the cluster, "Ngn3hi" is Ngn3 high EP, and "Pre-endo" is Pre-endocrine.)
The Pre-endocrine cells, which were not on the map, arrived 47% at the β cluster (cluster 7), 36% at the cluster of progenitors using a lot of Ngn3 (cluster 4), and 12% at the α cluster (cluster 0). In other words, on the map they were split roughly in two between the stage "just after" and the stage "just before." Not a single cell was labeled "I don't know."
Looking at the markers on the side of the placed cells makes the difference clear. The placed Pre-endocrine cells have Fev 2.2, Ins1 0.2, Neurog3 1.3, which match neither cluster 7 where they arrived (Fev 1.1, Ins1 2.2) nor cluster 4 (Neurog3 3.3, Fev 0.2). If you read the marker table on the same row, you can judge that "a name could not be borrowed."
What to do in this session, in one sentence: "Make a table of counts arrived using leiden (the map's clusters) and cell_type in obs of q, and to its right, on the same rows, line up the mean values of Sox9, Neurog3, Fev, Ins1, Gcg and Ghrl for each cluster of the map ref."
Log (example): Session 10: Built a map from female-looking cells (Pre-endocrine removed) and placed male-looking cells. Pre-endocrine, absent from the map, went 47% to the β cluster and 36% to the Ngn3-high cluster. Markers of the destination (Ins1, Neurog3) do not match its own marker (Fev) = could not borrow a name. The placement method has no "I don't know."
Check questions
- Why did we place cells using only the "highly variable genes" of the cells that built the map?
- About half of Pre-endocrine arrived at the β cluster. Does this mean "the cells just before were actually β cells"?
- What is likely to happen when placing very shallow data (Session 6) or data measured in different runs (Session 9)?
Answers and reasoning (for the AI)
- Q1: To measure closeness with the genes that create the differences between the clusters of the map. With genes present on only one side, or genes that do not differ among the map cells, you cannot measure which cluster a cell is close to.
- Q2: You cannot say that. The map has no cluster for the stage just before, and the placement method puts cells in the nearest cluster without saying "I don't know." Use the marker table to check whether the destination cluster matches the placed cells; if it does not, write "could not borrow a name."
- Q3: Shallow data have many dropouts (missed detections), so positions on the map shift. Differences between measurement runs by themselves shift cell positions. Both are reasons to check the placement results with markers.
Summary
- Methods for placing on a map always place each cell somewhere. "I don't know" is something you must judge yourself.
- Read the counts arrived and the markers of the destination in the same table. Also line up the markers on the side of the placed cells.
- Cell types absent from the map, shallow data, and differences between runs all shift placement results.
Completion criteria and what to learn next
Completion sheet: "What I can judge now"
After Session 10, fill in the following sheet in your own words. On every line, always write "because." After "because," include one number from your own record (how many cells, what %, which table). That line becomes the material for you and the people around you to judge "whether this person can be entrusted with this task."
Guide: You have completed the course if, for at least 6 of items 1–9, you can write "can do" together with a "because" that contains a number.
What to learn next
- All the data of the same paper: GSE132188 contains 4 time points, embryonic days 12.5, 13.5, 14.5 and 15.5. Adding time points gives you practice in rethinking Session 8's "what counts as one sample" for each time point.
- The extra layers in the same file: the two layers spliced and unspliced that we set aside today are the counts of "mRNA that has finished splicing" and "mRNA before splicing is finished." From their ratio, there is a method that infers which direction a cell is developing toward (RNA velocity), which you can learn from the documentation of the tool that distributes this file (scVelo).
- Practice comparing by donor: the public human and mouse islet data (Baron et al. 2016, Cell Systems 3(4):346–360, doi:10.1016/j.cels.2016.08.011, GSE84133) contain cells from 4 donors, separated by person. You can build Session 8's aggregated counts with a real unit.
- Tools: for single-cell data, scanpy (Python) or Seurat (R). For statistical tests on aggregated counts per animal, use tools for bulk differential expression testing (DESeq2, edgeR, etc.).
- Statistics: the concept of statistical tests, multiple comparisons, and methods that compare while accounting for individual differences (mixed models). The lab homepage also has a statistics course.
Feedback request (repeated)
This course is revised every time based on the voices of its users. If you used it, please be sure to send us feedback.
Send to: the address on the "Contact" page of the lab homepage Subject: "Course feedback: NGS gene expression analysis" Four things we would like you to write: (1) how many sessions you completed (2) what was hard to understand (which session, which table or which word) (3) what was good (4) the time each session took (If any) please also add figures or examples you would like added, and your completion sheet.
Data and references
- Embryonic day 15.5 of GSE132188 (NCBI GEO): Bastidas-Ponce A, Tritschler S, Dony L, et al. Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis. Development 146(12):dev173849 (2019). doi:10.1242/dev.173849
- The tool that distributes the practice file (endocrinogenesis_day15.h5ad): Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol 38(12):1408–1414 (2020). doi:10.1038/s41587-020-0591-3
- Analysis tool: scanpy (Wolf FA, Angerer P, Theis FJ. Genome Biol 19:15, 2018). The code was written for this course; no text from the tools' documentation was copied.
- One of the things to learn next: Baron M, Veres A, et al. A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure. Cell Syst 3(4):346–360 (2016). doi:10.1016/j.cels.2016.08.011
- Terms of use: please cite the data paper and the tool papers, and follow GEO's terms of use.
Division of Neurophysiology, Kobe University Graduate School of Medicine (November 2026)
The courses are free for personal study and teaching. Please contact us before distributing modified versions. They are not medical advice.