AI agents for bioinformatics
≈ 10 min
AI tools can now do bioinformatics for you: you say what you want, and an agent maps the reads and calls the variants. In this practical you use a real AI model in two ways – as an assistant that answers and as an agent that acts – on real sequencing reads, with real programs. Then you turn what the agent did into an analysis that anyone can repeat.
The question
The gene CYP2C19 makes an enzyme that turns the anti-clotting drug clopidogrel into its active form. One common variant, CYP2C19*2 (rs4244285, G>A), stops the enzyme from being made, and people who carry it respond less well to the drug. Your question for today:
Which genotype does the sample NA12878 have at CYP2C19*2?
You have the sequencing reads of that sample and a reference sequence. In this data the variant is at position 11,616 of the sequence called human_CYP2C19. You will first look at the data yourself, then answer the question twice – with an AI assistant, and with an AI agent – and in the end make the answer one that someone else can check.
| Terminal | A Linux-like terminal. The programs are the real ones – fastp, minimap2, Bowtie 2, samtools, bcftools and the GNU text tools – compiled to run in your browser, on this computer. The data is small, so each program takes seconds. |
|---|
| Files | Your folders, an editor for text files, and a viewer for reports. |
|---|
| AI | A real AI model, in two roles: Chat, an assistant that answers questions, and Agent, which carries out a task itself. |
|---|
Assistant or agent?
An assistant writes text. You ask, it answers; you run the commands and read what they print.
An agent is the same kind of model in a loop: it writes a command, the command is run, the model reads what it printed and writes the next command – until it decides the task is done. You say what you want, not how.
Your taskin words
→
AI modelreads everything so far, writes the next command
command ↓↑ what it printed
Your terminalthe page runs the command
→
Reportwhen the model decides it is done
Connect the AI
- Open the AI tab and read What is sent, and where.
- Paste an API key and press Connect. Your lecturer will say which key to use: one that they hand out, or one of your own – for Google’s Gemini API you make one at aistudio.google.com/apikey (sign in, Create API key, copy it).
- In the Terminal, see which programs are there:
help A show me button only types the command – press Enter to run it.
Three rules for today – and for every day after
- No personal or patient data in an AI tool. What you type, and what the agent’s commands print, goes to the company that runs the model. Today’s data is a public reference sample, so that is fine.
- The result is yours. An AI can be wrong and sound certain. Whatever you report, you must be able to show where it comes from.
- A key is a password. Do not share it or put it in a file. This page keeps it in this browser tab only.
Good to know
- A free key allows each model only a few requests per minute, and a limited number per day. The agent uses one request for each step.
- If a Gemini model is busy, over a limit or gives no answer, the page asks another one, and leaves the first alone for a while: an amber dot beside the model’s name at the top of the AI tab says so. Under each answer you can see which model replied.
- When no model can answer for the moment, the agent waits – a countdown is shown – and goes on by itself. If there is still no answer after three waits, the run ends with a message that says what each model said.
- Your work is kept in this browser. My answers (top right) downloads your answers; chapter 4 ends with a download of your analysis.
Suppose an AI agent tells you: “NA12878 is heterozygous for CYP2C19*2.” What would you want to see before you put that sentence in a report?
The sentence is a claim. The evidence would be: the commands that were run, and the versions of the programs; the line of the variant file for that position – the genotype, its quality, how many reads cover the position and how many support each allele; how well the reads mapped; and that the analysis can be run again and gives the same answer. By the end of this practical you will have all of these.
Look before you delegate
≈ 20 min
You cannot judge an agent’s work on data you have never looked at. So, before any AI: what is in the files, are they intact, how good are the reads, which programs are here? A few commands of your own.
The files
- List the data:
cd ~ && ls -l data Then read what it is: cat data/README.md
- Look at one read – four lines: a name, the bases, a plus sign, a quality for each base.
head -4 data/NA12878_R1.fastq
- Count the reads of the first file: its lines, divided by four.
echo $(( $(wc -l < data/NA12878_R1.fastq) / 4 ))
- The reference: which sequences does it have?
grep '>' data/reference.fa
How many read pairs are there? (The second file has the same number of reads: read 1 and read 2 of each pair.)
3,519 pairs – 14,076 lines in each file, four lines per read. Each read is 76 bases long.
Are the files what they should be?
- A checksum is a fingerprint of a file’s contents: change one byte and the checksum changes. The
data folder has a list of the expected checksums, MD5SUMS; check the files against it: (cd ~/data && md5sum -c MD5SUMS) The brackets make the cd last for this one line only.
Three times OK: these are exactly the files the practical was written for. You will use checksums again in chapter 4, to compare results.
Quality control with fastp
- fastp measures the quality of the reads, removes poor reads and trims adapters. It writes the cleaned reads and two reports.
fastp -i data/NA12878_R1.fastq -I data/NA12878_R2.fastq -o trimmed_R1.fastq -O trimmed_R2.fastq -h fastp.html -j fastp.json
- Open the report:
open fastp.html Read the Summary, then scroll to the quality curves: quality along the read, before and after filtering.
- The same numbers, for programs to read, are in the JSON report.
jq .filtering_result fastp.json jq .summary.before_filtering fastp.json
How many of the 7,038 reads passed fastp’s filter?
6,170 (87.7 %). 864 reads failed for low quality and 4 for too many unknown bases (N); none was too short. Adapters were trimmed from 42 reads.
fastp removed 868 reads – about one in eight. Are these bad data? What in the report helps you decide?
Not bad, but not new either. Before filtering, 82 % of the bases have a quality of at least 30 (one error in 1,000) and 90 % of at least 20; the curves show quality falling towards the end of the reads, more so in read 2 – usual for Illumina reads of this age. Only 42 reads carried adapter sequence, and the duplication rate is about 7 %. Nothing here says the data cannot answer the question; it does say that a careful analysis might use the filtered reads, and that a result should not hang on a few low-quality bases.
Which programs, which versions?
fastp --version; minimap2 --version; bowtie2 --version | head -1; samtools --version | head -1; bcftools --version | head -1
What you now know – and an agent does not, until it looks
3,519 pairs of 76-base reads; a reference with two sequences, human_CYP2C19 (30,000 bases) and human_CYP2C9 (26,000 bases); about one read in eight of low quality; and the programs at hand: fastp 0.20.1, minimap2 2.22, Bowtie 2 2.4.2, samtools 1.17, bcftools 1.10. With that you can tell whether what an AI says about this data makes sense.
The assistant answers – you act
≈ 20 min
An AI assistant is a language model that writes text; here, commands. It runs nothing, and it sees nothing of your files unless you show it. You are the one who acts – and the one who checks.
Ask for the mapping
- Make sure the terminal is in your home folder , then give the assistant this prompt – read it first: what does it tell the assistant, and what does it leave open?
A prompt for the assistant
I have paired-end Illumina reads in data/NA12878_R1.fastq and data/NA12878_R2.fastq, and a reference in data/reference.fa. Give me the commands to map the reads with minimap2 and to make a sorted, indexed BAM file called mapped.bam. Say in a few words what each option does.
- Read the answer. Under each block of commands is a button, Put command … in the terminal: it types one command for you. Read the command, press Enter, read what it prints – then the next one, until
mapped.bam and its index are there.
- Check the result yourself:
samtools flagstat mapped.bam
What percentage of the reads were mapped? (If the assistant gave you other options than -ax sr, your number may differ a little.)
With minimap2 -ax sr: 96.86 % – 6,817 of 7,038 reads; 94.77 % as proper pairs.
Ask for the variants
A prompt for the assistant
Now give me the commands to call the variants from mapped.bam with bcftools into a compressed, indexed file variants.vcf.gz, and a command that shows me the line for position 11616 of human_CYP2C19.
variants.vcf.gz is there.- If a command fails, do not guess: under the error message in the terminal press ✦ Ask the AI assistant about this error. The command and its message go to the assistant. (No error? Good – tick this step.)
- The line for the position should now be on your screen. If not:
bcftools view -H -r human_CYP2C19:11616 variants.vcf.gz
Which genotype does NA12878 have at position 11,616? Write it as it stands in the last column of the line (before the colon).
0/1 – heterozygous: one copy of the reference allele G and one of the alternative allele A. NA12878 carries one CYP2C19*2 allele.
Read the whole line. What is the evidence for this genotype?
With minimap2 and bcftools’ default settings: the reference base is G, the alternative A; the quality of the call (QUAL) is 222, about as high as bcftools goes; 106 reads cover the position (DP). DP4=40,4,46,2 counts the good-quality reads for each allele, forward and reverse: 44 for G and 48 for A – about half and half, as expected when one of the two chromosomes carries the variant. GT:PL 0/1:255,0,255: the heterozygous genotype is far more likely than either homozygous one.
Find the assistant’s limits
A prompt for the assistant
How many variants are there in my file variants.vcf.gz?
bcftools view -H variants.vcf.gz | wc -l A prompt for the assistant
Give me the bwa mem command to map these reads.
Why this assistant knows that bwa is not here
Before your first question, this page gave the model a description of the terminal: the programs and their versions, the data files, what is not installed. Such a text is called a system prompt; with every question the page also sends the name of your folder and of your files. A chatbot in another browser tab knows none of this. It would most likely give you a bwa mem command – correct in general, and useless here.
Did the assistant tell you how many variants there are? What can it know about your file, and what not?
It knows the file’s name (sent with your question) and how such files are made. It cannot open the file, and it did not see what your commands printed – unless you pasted that in. So the honest answer is “I cannot see it; run this command”, and the right number is the one you counted: 60 with minimap2 -ax sr and bcftools’ defaults. If an assistant gives a number anyway, it has made one up that looks right. That is what people call a hallucination, and it is why a number from an AI is worth nothing until you have seen the output it comes from.
Agents: four ways to let an AI act
≈ 45 min
An agent runs the commands itself. That saves you the typing – and takes away the moments in which you would have looked. How much you still see and decide depends on how you let it work.
An agent does what you did in the last chapter – run a command, read what it prints, decide what comes next – by itself. The page runs each of its commands in your terminal, in a folder of its own (~/runs/1-approve, ~/runs/2-plan …) that starts with a copy of the data. While it works the terminal is locked; switch to the Terminal tab to watch.
How much an agent may do without asking is a setting. There are four settings here. You will give the same task with each of them – with one of them twice, so five runs in all:
The task – the same for all five runs
Find the variants of sample NA12878. Map the paired reads data/NA12878_R1.fastq and data/NA12878_R2.fastq to data/reference.fa, call the variants and save them as variants.vcf.gz (compressed and indexed). Then tell me the genotype at position 11616 of human_CYP2C19, and the evidence for it.
1 · Approve each step
The agent proposes, you decide
- In the AI tab choose Agent and Approve each step, put the task in the box and press Send.
- For every step, read what the agent wants to run before you allow it. Would you have done the same? Press Run it if so. (For a step that only writes a file the button reads Write it.)
- At least once, do not just agree. Change the command in its box before you press Run it – or press Refuse, with your reason in the line for the agent. The agent is told what you did; watch how it reacts.
- Go on until the agent writes its report. (Tired of approving? Run, and stop asking lets it finish alone – many people press that button rather soon. If you do, remember at which step.)
- Now check the report against the file. On the run’s card press Go to its folder, then Enter – the prompt of the terminal now shows the run’s folder,
~/runs/1-approve. Then: bcftools view -H -r human_CYP2C19:11616 variants.vcf.gz
Compare the agent’s report with the line you just printed. Is every statement in the report backed by something a command printed? Is there anything you cannot find in the outputs?
The genotype (0/1), the quality and the read counts should match the line exactly, if the agent printed that line and read it. Look for the rest in the record (The record (RUN.md) on the card): how many reads mapped, how many variants. Statements that no command printed – “high confidence”, “a poor metaboliser”, a clinical recommendation – come from what the model knows in general, not from this analysis. They may be right (they often are), but nothing in this run shows it. And a carrier of one *2 allele is usually called an intermediate metaboliser, not a poor one: exactly the kind of detail to check.
2 · Plan first
Agree on the route, then let it drive
- Choose Plan first and send the same task. The agent answers with a plan and runs nothing yet.
- Read the plan as you would read a colleague’s. Then ask for one change – for example: “Check the quality of the reads with fastp first, and map the filtered reads.” or “Use Bowtie 2 instead of minimap2.” Press Ask for changes.
- When the plan is what you want, press Approve the plan. From here the agent works alone.
- When it has finished, open The record (RUN.md). It has the plan, and under What the agent did every command. Did the agent do what it said it would?
Did the agent keep to the plan? Where did it depart from it – and did it tell you?
Look for departures of these kinds: a command fails and is repaired; an extra check is added; an option or a file name is not what the plan said; a step is quietly dropped. A good agent says so in its report – but the report is the agent’s own account. The record shows what happened whether the agent mentions it or not. A plan is a promise; the record is what was done.
3 · Autonomous – twice
Let it run, and then let it run again
- Choose Autonomous and send the same task. Watch; you are not asked anything.
- When it has finished, send exactly the same task again, again autonomous.
- You now have four runs of one task. Compare them. This line prints, for each run, the model that answered and the main commands that worked:
cd ~/runs && for d in */; do echo "== $d"; jq -r '"model: " + (.models | join(", ")), (.steps[].commands[] | select(.code == 0) | .cmd)' $d/run.json | grep -E '^(model|fastp|minimap2|bowtie2|bcftools (mpileup|call|filter|norm|view -[a-zA-Z]*[ie]))'; done (jq reads run.json – the record of a run in a form for programs – and grep keeps the lines that start with one of the main programs. For everything a run did, open its RUN.md.)
- And this one prints the number of variants each run found, and what each found at the position:
cd ~/runs && for d in */; do echo "$d $(bcftools view -H $d/variants.vcf.gz | wc -l) variants, at 11616: $(bcftools query -f '%REF>%ALT QUAL=%QUAL [%GT]' -r human_CYP2C19:11616 $d/variants.vcf.gz)"; done
Error messages and 0 variants: the run has no variants.vcf.gz. Error messages and a number, but nothing after at 11616: – the file has no index (bcftools index). Both are results worth noticing.
Same task, same data – and, if the model: lines agree, the same model. Were your runs the same? List the differences that could change the result.
Most likely not – compare with your neighbours, too. A model need not give the same answer twice, and the task leaves choices open: which mapper (minimap2 or Bowtie 2); the raw reads or the reads filtered by fastp; default settings or others; a filter on the variants or none; which checks. (When a model is busy or over its limit the page asks another one, so even the model can differ between your runs – the model: lines show it.) Those choices change the list of variants: with default settings minimap2 and bcftools give 60 variants (27 of them with QUAL of 30 or more), Bowtie 2 and bcftools 50; after fastp it is 60 and 49. The answer to your question should be the same in all of them – G>A, 0/1, with QUAL 222 from minimap2 alignments and 178 from Bowtie 2’s – because that position is well covered. A result that survives different reasonable analyses is a robust one; but you only know that because you compared. And if your runs happened to be identical: nothing guaranteed it.
4 · Write a script
Ask for the analysis, not for the answer
- Choose Write a script and send the same task once more. The agent is now asked to write one script,
analysis.sh, to run it, and to repair it until it works.
- Open
analysis.sh: press Open next to its name on the run’s card (or find it in the Files tab under runs). Read it from top to bottom. Can you follow every step?
Name one thing in analysis.sh that you would change or check before you used it on other data.
Whatever you found – for example: file and sample names written into the commands, so that the script works for this sample only; a threshold (a minimum quality, a depth) that the agent chose without giving a reason; no test that the input files exist or that a step produced something; nothing that records the versions of the programs; an option you do not know and would look up with --help. The point of a script is that you can read and question all of this – a chat that ended with “the genotype is 0/1” gives you nothing to question.
The four ways, side by side
| Who decides each step | What you must do | Use it when |
|---|
| Approve each step | You, one step at a time | Read every command before it runs | the tools or the data are new to you; a command could destroy something |
|---|
| Plan first | You approve the route; the agent drives | Judge the plan; afterwards check that it was kept | the approach matters more than each single command |
|---|
| Autonomous | The agent | Check afterwards: the record and the results | routine work in a place where nothing can be damaged, and checking is cheap |
|---|
| Write a script | The agent writes; the script does | Read the script; run it again | you will repeat, share or publish the analysis |
|---|
Outside this page
Agent tools for the command line and for code editors have the same settings under other names: ask before each command, plan mode, run without asking. The difference is where they run. Here the agent works in a web page, in a folder of its own; what it changes anywhere else is put back. On your laptop or on a cluster an autonomous agent can delete or overwrite real data, and it will not always notice. So: keep raw data read-only, work in a copy or a container, keep your work under version control – and do not give an agent more access than the task needs.
Make it reproducible
≈ 25 min
The model that did your analysis will be replaced before long, and it might not write the same commands twice anyway. What can be kept is what was done: the data, the commands, the versions. This chapter turns one of your runs into an analysis that anyone can run again – and tests whether it really gives the same result.
Three things make an analysis repeatable: the data (and a way to know it is the same data), the commands (as a script, not as a memory), and the programs (with their versions). The agent’s runs have left you all three – in pieces. Now you put them together and test them.
The record: what was done
- Choose the run you trust most. In the AI tab, on its card, press The record (RUN.md). Find: the task as you gave it; the model; each command with its exit status; the programs and versions; the checksums of the input files.
What does RUN.md tell you that the agent’s report did not? And who wrote it?
Every command, also the ones that failed, with exit status and the end of what it printed; what you approved, changed or refused; the versions of the programs; the checksums of inputs and results; the model and the time. The report is the agent’s summary of its own work, and a summary leaves things out. RUN.md was written by the page, from what really ran – the agent cannot improve on it. Agent tools on your own computer keep such logs too (session transcripts). Keep them with the results: they are your lab notebook for work you did not do by hand.
From the record to a script
- On the same card press Make a script from this run. The page writes
rerun.sh: the commands that worked, in the order they ran. A command that failed and left nothing behind is there as a comment – and so is any command that reached outside the run’s folder. A command that ended with an error status but did write a file is kept (diff ends with status 1 when two files differ: no mistake, and the next step may read its output). (A Write a script run has its script already, analysis.sh, and no such button: open that one, and tick this step yourself. If the agent did not write it, the card says so – and has the button.)
- Read the script. Take out what only looked at things (
ls, head, cat …) and is not needed to make the result. Leave the comments that say what a step is for – or write better ones. Where the agent named its folder in full (/home/student/runs/…), the page has written "$RUN_FOLDER": the folder the script runs in, set near the top. A line that begins with # left out says why a command of the run is not in the script; a line that begins with # CHECK: stands above a command that the page kept but could not check – read those first. The script stops at the first command that fails and at a variable that has no value (set -eu); the agent’s own commands ran without that. Where its run went on all the same – a grep that finds nothing inside a loop, a variable that was never set – the page has put set +e or set +u before the command, with the reason, and switched it on again after it: decide whether the analysis should go on there. A command that ended the way a script ends – an exit inside it, the agent’s own set -e – stands in a subshell, ( … ), so that it ends alone, as it did in the run; the lines above it say so. Check, too, that no command reaches outside its folder (a path that starts with ~, / or .., or a cd that leaves the folder): in the agent’s run the page put such changes back; when you run the script, nothing does. Save with Ctrl+S.
Does it repeat?
- A script that works in the folder where everything already exists proves little. Test it in a new folder that has only the data – a fresh copy of the course data, not the files the agent has worked with – and the script. On the run’s card press Go to its folder and Enter; then:
mkdir -p ~/check1 && cp -rf ~/data rerun.sh RUN.md ~/check1/ && cd ~/check1 && bash rerun.sh For a Write a script run
Its script is called analysis.sh. This line copies it under the name rerun.sh, so that everything below is the same for you: mkdir -p ~/check1 && cp -rf ~/data RUN.md ~/check1/ && cp -f analysis.sh ~/check1/rerun.sh && cd ~/check1 && bash rerun.sh
- If the script stopped with an error: read the message, repair the script (
nano rerun.sh) and run it again with bash rerun.sh – or ask the Chat. To start afresh instead, remove the test folder (cd ~ && rm -rf ~/check1), go to the run’s folder again and repeat the line above.
- Once more, in a second new folder:
mkdir -p ~/check2 && cp -rf ~/data ~/check1/rerun.sh ~/check2/ && cd ~/check2 && bash rerun.sh
- Are the two results the same file?
md5sum ~/check1/variants.vcf.gz ~/check2/variants.vcf.gz
- Not the same! Find the difference: unpack both and compare them line by line.
zcat ~/check1/variants.vcf.gz > ~/v1.vcf; zcat ~/check2/variants.vcf.gz > ~/v2.vcf; diff ~/v1.vcf ~/v2.vcf
- Now compare what matters – the variants, without the header – for both test folders and for all the agent’s runs:
for d in ~/check1 ~/check2 ~/runs/*; do echo "$(bcftools view -H $d/variants.vcf.gz | md5sum) $d"; done
Which line differs between the two files, and why?
A line of the header: ##bcftools_callCommand=call …; Date=…. bcftools writes into the file when it was run. (If your script filters or converts the file with bcftools, there are more such lines.) One changed character gives another checksum, so two .vcf.gz files made at different times are not identical – although every variant in them is.
So: did the script give “the same result” twice? What do you compare to decide – and what did the last command show about the agent’s runs?
Yes: the variant lines (bcftools view -H) have the same checksum in check1 and check2 – and the same as the run the script was made from. “Reproducible” needs a definition: here, the same variants with the same genotypes and qualities, not the same bytes in a file that holds a date. Say which you mean, and give the checksum of the thing you compare. The other runs of the agent show other checksums wherever the agent chose other commands: identical checksums mean the same variant list; different ones, a different list.
Say what was done – the AI included
- A folder that someone else can use needs a note. Create this file in
~/check1 and fill in every ... – from RUN.md and from what you did yourself.~/check1/README.md
# CYP2C19*2 in NA12878
## Question
Which genotype does the sample NA12878 have at CYP2C19*2
(rs4244285; position 11,616 of the sequence human_CYP2C19)?
## Answer
Genotype: ...
Evidence: QUAL ..., ... reads at the position, ... support G and ... support A
(bcftools view -H -r human_CYP2C19:11616 variants.vcf.gz)
## How to repeat the analysis
In a folder with data/ and rerun.sh: bash rerun.sh
The same result is: bcftools view -H variants.vcf.gz | md5sum -> ...
(The .vcf.gz file itself differs from run to run: its header holds the date.)
## Data (md5)
d2b8ddb488975fa5ef18e6e59edb5681 data/NA12878_R1.fastq
aee82e1b82bedb037208f5ff5b9ff34b data/NA12878_R2.fastq
d9155e00b786f74b0837f6a675483563 data/reference.fa
## Programs
... (copy the list from RUN.md)
WebAssembly builds, run in a web browser.
## Use of AI
Tool: the AI agent of the practical "AI agents for bioinformatics"
Model and date: ...
What the AI did: chose and ran the commands (see RUN.md, which the page wrote, not the AI)
What I did: ... (for example: read every command; took out ...; ran the
script twice in new folders; compared the variants; read the
line for position 11,616 myself)
- Fill it in and save. Under Use of AI, be exact about who did what: it protects you, and it tells the reader what was checked by a person.
- Download the folder: data, script, record, results, note. That file is your analysis.
A colleague opens your .zip two years from now. The model you used no longer exists. What can they still do – and what not?
They can read what was asked and what was done (RUN.md), check that the data is the same (the checksums), run rerun.sh with the same versions of the programs and get the same variants, and read which parts a person checked. (With minimap2 the checksum in your README comes out the same on another computer. Bowtie 2 places a few reads that fit two places equally well by a pseudo-random number, which can differ between computers: there, compare positions and genotypes.) They cannot have the same conversation with the agent again – and might not get the same commands if they could. That no longer matters: the analysis does not depend on the model any more. What it does depend on is the programs, and old versions are not always easy to get. That is why they are written down; in real projects an environment file or a container fixes them.
Reference
A checklist, the commands, the words, and what to do when something does not work.
Before you believe an agent
- Do I know the data well enough to notice a wrong answer? (chapter 1)
- Is each statement in the report backed by something a command printed?
- Have I looked at the result myself – the line, the count, the file – and not only at the report?
- Do I know what the agent decided for me: the program, the options, the filters?
- Is there a record of what was run, written by the tool and not by the model?
- Is there a script, and have I run it somewhere fresh?
- Have I said what the AI did and what I checked?
Commands of this practical
fastp -i R1.fastq -I R2.fastq -o out1.fastq -O out2.fastq -h fastp.html -j fastp.json | quality control and filtering of paired reads |
minimap2 -ax sr ref.fa R1.fastq R2.fastq > mapped.sam | map short paired reads |
bowtie2-build ref.fa ref · bowtie2 -x ref -1 R1.fastq -2 R2.fastq -S mapped.sam | the same with Bowtie 2: index first, then map |
samtools sort -o mapped.bam mapped.sam · samtools index mapped.bam | sort and index the alignments |
samtools flagstat mapped.bam | how many reads mapped |
bcftools mpileup -f ref.fa mapped.bam | bcftools call -mv -Oz -o variants.vcf.gz | call variants |
bcftools index variants.vcf.gz | index the variants (needed for -r) |
bcftools view -H -r human_CYP2C19:11616 variants.vcf.gz | the line for one position |
bcftools view -H variants.vcf.gz | wc -l | count the variants |
bcftools view -H variants.vcf.gz | md5sum | a checksum of the variants, without the header |
md5sum -c MD5SUMS | check files against a list of checksums |
bash script.sh | run a script; set -euo pipefail at its top stops it at the first error |
Every program explains itself: samtools --help, samtools view --help, bcftools call --help, fastp --help. (Some commands of samtools and bcftools print their help and then end with an error status, as they do on Linux – the help is there all the same. One, bcftools mpileup, shows its help only when it is typed with nothing after it.) help lists what the terminal has.
Words
| Language model | A program that continues text. Asked a question, it writes a likely answer – from what it was trained on and from what is in front of it. |
|---|
| Assistant | A language model you talk to. It writes; you act. |
|---|
| Agent | A language model in a loop with tools – here it can run commands in the terminal and write text files. It acts, reads the result and acts again, until it decides it is done. |
|---|
| Prompt | The text the model is given. The system prompt is the part you do not see: what the tool told the model about its job and its surroundings. |
|---|
| Context | Everything the model can “see” at one moment: system prompt, your messages, its own earlier replies, outputs of commands. What is not in the context does not exist for the model. |
|---|
| Hallucination | Text that looks right and is not backed by anything: an invented number, option, file or reference. |
|---|
| Non-deterministic | The same prompt need not give the same answer twice. |
|---|
| Provenance | Where a result comes from: which data, which commands, which programs, who or what ran them. |
|---|
| Reproducible | Someone else, with the same data and the same instructions, gets the same result. |
|---|
A statement on the use of AI
Rules differ between courses, journals and employers – read yours. A useful statement answers four questions:
| Which tool and model, and when? | for example: the agent of this practical; the model’s name as it stands in RUN.md; the date |
|---|
| What did the AI do? | for example: chose and ran the commands for mapping and variant calling |
|---|
| What did I do? | for example: set the task; read each command; ran the script again in a new folder; checked the genotype in the VCF file myself |
|---|
| Where is the evidence? | for example: RUN.md (the record of the run), rerun.sh (the script) |
|---|
The note you wrote in chapter 4 (README.md, Use of AI) is such a statement.
If something does not work
| “Too many requests” (429) | Every model that the page can ask is over its limit per minute. The agent waits (up to a minute and a half, three times at most) and goes on by itself; in the chat, wait a minute. If an agent’s run ends with this message, many people may be using the same key: give the task again a little later. |
|---|
| “No requests are left for today” | Every model is over its limit per day for your key; the message says when the service takes requests again. Waiting does not help today, and a second key from the same Google project does not either (the limits are counted per project): work with a neighbour, or ask for another key. |
|---|
| “Busy”, “high demand” (503) | The service is overloaded. The page asks other models by itself and names the one that answered. If a run ends with this message, every model was busy: wait a few minutes and give the task again. |
|---|
| “… sent an empty reply” | A model sometimes ends a reply before it has begun. In the chat, ask again; the agent asks again by itself. |
|---|
| The key is not accepted | Copy it again, without spaces, and check that the service chosen in the form is the one the key belongs to. An old Gemini key may no longer be accepted by Google: make a new one. |
|---|
| A command never ends | Press Ctrl+C in the terminal; press it again to stop the program by force. The programs then start again. A file that a program wrote in the last seconds before may be lost, or come back as it was before – the terminal says which: make it again. |
|---|
| The agent goes in circles | Press Stop. Its folder and the record stay. Start a new run, and say in the task what it should avoid. |
|---|
| The page was reloaded | Your ticks, answers and files are still there. (A file that a program wrote in the last seconds before is not: the terminal names it.) An agent that was running cannot go on, and that run has no record – its folder holds a note instead: give the task again. |
|---|
| “Your browser storage is full” | The browser gives a site only a few megabytes for text files, and other practicals of the same site share them. Keep what you have with My answers → Save all my work to a file; then make room: remove large text files you do not need (a SAM file, say: rm mapped.sam), or use Start again in a practical that you have finished. |
|---|
| “This practical is open in another tab or window” | Work in one tab or window only. Two of them save to the same place, and each writes over what the other has saved. Close one. |
|---|
| Start from scratch | My answers → Start again clears everything in this browser. |
|---|
The data
3,519 pairs of 76-base Illumina reads of the public reference sample NA12878 (exome sequencing, run SRR098401 of the 1000 Genomes Project), from the regions of CYP2C19 and CYP2C9; and two pieces of the human reference genome hg19. Details: cat ~/data/README.md