1.4 Running RNA-Seq-Pop on the Busia data¶

Theme: Tools & technology

In Module 1.3 you built a three-rule Snakemake workflow. Now you configure and launch a real one. RNA-Seq-Pop [1] is a Snakemake workflow that takes RNA-seq reads all the way to differential expression and population genetics: allele frequencies of resistance mutations, Fst and PBS selection scans, genetic diversity, chromosomal inversion karyotypes and species ancestry.

We run it on the study at the centre of this course: a pyrethroid-resistant Anopheles gambiae colony from Busia, Uganda, sequenced together with the susceptible Kisumu strain (ENA PRJNA748581). By the end of the module your own copy of the workflow will be running on a subsampled “mini” version of the data in GitHub Codespaces. On Day 2 we explore the results of the same workflow run on the full dataset.

Learning objectives¶

At the end of this module you will be able to:

• Describe what RNA-Seq-Pop does, and why it is worth doing population genomics on RNA-seq data.

• Find your way around the RNA-Seq-Pop repository: config/, resources/, workflow/ and results/.

• Write a sample sheet (samples.tsv) and configuration file (config.yaml) for an experiment, including contrasts, ploidy and PBS comparisons.

• Dry-run the workflow, read its job list and rule graph, and launch it in a GitHub Codespace.

Setup¶

The Python cells need only pandas and PyYAML, which Colab already has.

What RNA-Seq-Pop does¶

Most RNA-seq pipelines stop at differential expression. RNA-Seq-Pop also uses the sequence in the reads. It has five groups of modules, each switched on or off in config.yaml:

At the end it builds a results book (a small Jupyter Book, like this course) at results/rna-seq-pop-results/_build/html/index.html, with one page per analysis.

Why do population genomics on RNA-seq?¶

Transcriptomic studies of insecticide resistance usually sequence pools of mosquitoes from a resistant and a susceptible group. The reads cover every expressed gene, often at high depth, so the data already hold a lot of genetic information. Using it has several benefits:

• It comes free. No extra sequencing is needed. In the Busia paper, 73% of the ~734,000 SNPs called were in exons: RNA-seq behaves much like exome sequencing.

• It finds target-site resistance. Mutations such as Vgsc-995S (kdr) change the protein, not its expression. In Busia the Vgsc gene was not differentially expressed between the parental and selected lines, but Fst in Vgsc was in the top 5% of genes and 995S rose from ~25% to 100%. An expression-only analysis would have missed it.

• It explains expression. A gene may be overexpressed because it sits on a haplotype that is being selected, for example inside a chromosomal inversion. Fst, PBS and karyotypes put the expression results in context.

• It checks your samples. PCA, ancestry markers and karyotypes can reveal mislabelled samples, contamination or hybridisation. RNA-Seq-Pop found that the Kisumu lab strain carries a lot of putative An. coluzzii ancestry on its autosomes.

There are caveats (allele-specific expression, uneven coverage, pooled samples), which we cover in Module 2.2.

The Busia experiment¶

A colony was founded from An. gambiae collected in Busia, Uganda, in November 2018 [1]. After 24 generations in the insectary it had lost most of its pyrethroid resistance. Part of the colony was then selected with deltamethrin for four generations, and the G28 survivors of a final deltamethrin exposure were sequenced. RNA was extracted from pools of five 4-day-old females and sequenced as 150 bp paired-end reads on an Illumina NovaSeq 6000.

The ENA names say “Sus” and “Res”, but the parental line is not a fully susceptible strain: it is the unselected parent of the selected line. We use clearer treatment names in the workflow.

The full data is ~113 GB of FASTQ. For the course we subsampled each sample to 1 million read pairs (the “Busia mini” dataset, ~1.5 GB), which is enough to exercise every step of the workflow in a few hours. The script that made it is data/prepare_busia_mini.sh in the course repository.

How RNA-Seq-Pop is organised¶

RNA-Seq-Pop follows the standard Snakemake layout. You do not install it; you clone the repository, add your configuration and data, and run snakemake from the top folder.

rna-seq-pop/
├── config/
│   ├── config.yaml           ← you write this: which modules, contrasts, reference files
│   └── samples.tsv           ← you write this: one row per sample
├── resources/
│   ├── reads/                ← FASTQs: {sampleID}_1.fastq.gz, {sampleID}_2.fastq.gz
│   ├── reference/            ← genome, transcriptome, GFF, gene map, GO annotations
│   ├── exampleMutations.tsv  ← variants of interest (Vgsc kdr, Cyp4j5, Gste2, Ace1, Rdl ...)
│   ├── gamb_vs_colu.zarr     ← ancestry-informative markers
│   └── karyotype_tag_snps.csv
├── workflow/
│   ├── Snakefile             ← the entry point: includes the rule files, defines rule all
│   ├── rules/*.smk           ← qc, diffexp, hisat2-freebayes, variantAnalysis, ...
│   ├── envs/*.yaml           ← one conda environment per group of rules
│   ├── scripts/              ← R and Python scripts called by rules
│   └── notebooks/            ← notebooks run by papermill to make the results book
├── .test/                    ← a tiny test dataset (X chromosome only) used by CI
└── results/  logs/           ← created when you run the workflow

Let’s look at the Snakefile of the pinned release.

Note three things:

1. configfile: "config/config.yaml": Snakemake reads the configuration from config/, relative to the folder you run it from. That is why we always run from the rna-seq-pop/ folder.

2. include: pulls in the rule files. The variants-of-interest rules are only included if that module is switched on.

3. rule all asks for GetDesiredOutputs, a Python function in rules/common.smk that builds the list of target files from the config. Switch a module off, and its outputs disappear from the list, so its rules never run.

The sample sheet¶

samples.tsv has one row per sample. RNA-Seq-Pop needs four columns:

• sampleID must match the FASTQ names (Kis1 → resources/reads/Kis1_1.fastq.gz and Kis1_2.fastq.gz).

• treatment defines the groups that are compared. No underscores, because contrasts are written as groupA_groupB.

• species is used for sorting and in some plots.

• strain groups samples into populations for freebayes. We set it to the treatment, so each group is treated as its own population during variant calling.

Here is the sample sheet for Busia, from workflows/rna-seq-pop/ in the course repository:

The configuration file¶

config.yaml controls everything else. The most important settings for Busia are:

dataset: 'Ag_Busia'
fastq:
      auto: True          # reads are resources/reads/{sampleID}_1.fastq.gz / _2.fastq.gz
      paired: True
contigs: ['2L', '2R', '3L', '3R', 'X']
contrasts:                # control_case: log2FC > 0 means higher in the case group
      - 'Kisumu_BusiaParental'
      - 'BusiaParental_BusiaSelected'
VariantAnalysis:
      ploidy: 10          # pools of 5 diploid females = 10 genome copies
      selection:
            population-branch-statistic:
                  contrasts:
                        - 'BusiaSelected_BusiaParental_Kisumu'   # pop1_pop2_outgroup

• Contrasts follow the design of the paper. Kisumu vs parental asks what differs between a lab-susceptible strain and a field-derived colony; parental vs selected asks what changed during deltamethrin selection.

• Ploidy matters for pooled data. Freebayes calls each pool as if it were one organism with 10 copies of the genome, so allele frequencies can take values 0, 0.1, 0.2 ... 1. For individual diploid mosquitoes you would use ploidy: 2.

• PBS (population branch statistic) needs three groups: two closely related ones and an outgroup. BusiaSelected_BusiaParental_Kisumu measures how much the selected line has changed since it split from the parental line, using Kisumu to polarise the comparison.

• References. RNA-Seq-Pop does not download a genome. We provide AgamP4 from Ensembl Metazoa (release 63), prepared by data/prepare_busia_mini.sh --reference and shipped with the mini dataset.

Let’s load the full file and run the same consistency checks that RNA-Seq-Pop’s first rule (CheckInputs) runs.

The DAG¶

With 16 samples and every module on, the dry run lists 535 jobs from 39 rules. The full job DAG is too big to read, so we look at the rule graph instead: one node per rule. This is a simplified version (the notebook and results-book rules are grouped):

Two branches run side by side from the trimmed reads: pseudoalignment to the transcriptome (kallisto → DESeq2) for expression, and alignment to the genome (HISAT2 → freebayes) for variants. Freebayes is split into 8 chunks for each of the 5 chromosome arms (40 jobs), and allele counts at the 22 variants of interest are taken separately for each of the 16 samples (352 small jobs). These are independent, so Snakemake runs them in parallel.

Practical: run RNA-Seq-Pop in a Codespace¶

Step 0: before you start¶

GitHub stops an idle codespace after 30 minutes by default, even if a job is running. Go to GitHub → Settings → Codespaces → Default idle timeout and set it to 240 minutes. If your codespace does stop, restart it and run the same Snakemake command again: finished jobs are not repeated.

Remember that codespaces run on your own free allowance (see the course home page). Delete your codespace at the end of the course.

Step 1: fork the course repository¶

1. Go to https://github.com/sanjaynagi/anopheles-omics-training.

2. Click Fork → Create fork. You now have https://github.com/<your-username>/anopheles-omics-training.

Step 2: open a codespace¶

1. On your fork, click Code → Codespaces → ⋯ → New with options.

2. Choose 4-core (16 GB RAM) as the machine type, then Create codespace.

3. Wait for the setup to finish (the first start installs the Snakemake environment and downloads the workflows; it can take several minutes). VS Code opens in the browser, with the course repository at /workspaces/anopheles-omics-training.

4. Open a terminal: ☰ menu → Terminal → New Terminal.

Check that the tools are there:

which snakemake          # a path inside the course's pixi environment
snakemake --version      # should print 7.32.4
conda --version          # still needed, see below
echo $SNAKEMAKE_PROFILE  # the course's Snakemake profile

Two kinds of environment. The codespace uses two package managers, for two different jobs:

• pixi manages the environment that runs Snakemake (Snakemake 7.32.4, Python, pandas, graphviz). It is defined in workflows/pixi.toml in the course repository, with a lock file that pins every package, and it is activated automatically in every terminal.

• conda is still needed because RNA-Seq-Pop gives each group of rules its own environment (workflow/envs/*.yaml), and Snakemake’s --use-conda builds those with conda. A Snakemake profile ($SNAKEMAKE_PROFILE) sets conda-frontend: conda and rerun-incomplete: true, so you do not need to type those options.

For your own projects, we recommend pixi: one pixi.toml per project, a lock file for exact reproducibility, and pixi run / pixi shell instead of activating environments by hand (Module 1.2).

Step 3: find RNA-Seq-Pop, pinned to a release¶

The codespace setup has already cloned the v2.3.0 release of RNA-Seq-Pop next to the course repository, at /workspaces/rna-seq-pop. You do not need to clone it again. Check which version you have:

cd /workspaces/rna-seq-pop
git describe --tags      # v2.3.0
git log -1 --oneline     # 38869ef

Always use a tagged release rather than master: the main branch keeps changing (rules have been renamed and output files changed since v2.3.0), and a pinned version means your results can be reproduced. (On another computer you would get the same version with git clone --branch v2.3.0 --depth 1 https://github.com/sanjaynagi/rna-seq-pop.git.)

To browse the files in VS Code, run code -a /workspaces/rna-seq-pop (adds the folder to the explorer).

Step 4: explore the repository¶

ls
ls config resources workflow workflow/rules workflow/envs
less workflow/rules/diffexp.smk          # press q to quit
grep -n "^rule" workflow/rules/*.smk | wc -l

Find the rule KallistoQuant in workflow/rules/diffexp.smk. What are its input, output and params? (Kallisto is run with -b 100: 100 bootstrap samples, which sleuth needs for isoform-level testing.)

Step 5: add the Busia configuration¶

The config and sample sheet are already written, in the course repository. Copy them into place and read them:

cp /workspaces/anopheles-omics-training/workflows/rna-seq-pop/config.yaml config/config.yaml
cp /workspaces/anopheles-omics-training/workflows/rna-seq-pop/samples.tsv config/samples.tsv
cat config/samples.tsv
less config/config.yaml

Step 6: get the Busia mini dataset¶

The mini dataset contains reads/ (32 FASTQ files, 1 M read pairs per sample) and reference/ (the AgamP4 files named in config.yaml). Unpack it into resources/:

BUSIA_MINI_URL="..."   # TODO(data): Zenodo record of the Busia mini dataset (busia-mini.tar.gz)
curl -L -o busia-mini.tar.gz "$BUSIA_MINI_URL"
tar -xzf busia-mini.tar.gz -C resources
ls -lh resources/reads | head
ls -lh resources/reference
zcat resources/reads/Kis1_1.fastq.gz | head -4

Step 7: dry run¶

Always dry-run first. It checks the config, finds every input and prints the jobs that would run:

snakemake --cores 4 --use-conda -n

The first thing printed is RNA-Seq-Pop’s banner (version, dataset). Then comes the list of jobs, and at the end a Job stats table like the one above (535 jobs). If something is wrong, for example a FASTQ file is missing or a contrast names a treatment that does not exist, the dry run fails here with an error message, before anything expensive happens.

To see why each job would run, add --reason.

Step 8: draw the rule graph¶

snakemake --rulegraph | sed -n '/^digraph/,$p' | dot -Tsvg > rulegraph.svg

(RNA-Seq-Pop prints its banner to the same output, so sed keeps only the graph.) Open rulegraph.svg from the VS Code explorer. Compare it with the diagram above. --dag draws the full job graph instead; try it with | dot -Tpng > dag.png and see why we prefer the rule graph.

Step 9: launch¶

Creating the ~10 conda environments takes the longest part of the first run. You can build them on their own first:

snakemake --cores 4 --use-conda --conda-create-envs-only

Then launch the run in the background, so it keeps going if you close the terminal, and follow the log:

nohup snakemake --cores 4 --use-conda --keep-going > snakemake-run.log 2>&1 &
tail -f snakemake-run.log      # Ctrl+C stops following the log, not the run

• --use-conda creates each rule’s environment from workflow/envs/*.yaml (or a wrapper).

• --keep-going carries on with independent jobs if one fails.

• The profile adds --rerun-incomplete, which redoes any job that was interrupted (for example if the codespace stopped).

While it runs, watch the outputs appear:

ls results/ results/counts/ | head
ls logs/
grep -c "Finished job" snakemake-run.log

Expect the mini run to take a few hours on 4 cores. You do not need to wait: the Day 2 notebooks use results from the full dataset, prepared in advance.

Step 10: save your work to git¶

Keep a record of what you ran in your fork of the course repository:

mkdir -p /workspaces/anopheles-omics-training/course-work/rna-seq-pop
cp rulegraph.svg config/config.yaml config/samples.tsv /workspaces/anopheles-omics-training/course-work/rna-seq-pop/
cd /workspaces/anopheles-omics-training
git add course-work/rna-seq-pop
git commit -m "Busia RNA-Seq-Pop run: config and rule graph"
git push

Troubleshooting¶

Exercise 1¶

Without running anything, use the sample sheet and config to work out how many KallistoQuant, VariantCallingFreebayes and mpileupVariantsOfInterest jobs there will be. Then check your answer against the dry run’s Job stats table.

Exercise 2¶

You want a third comparison: Kisumu vs the selected line. Add the contrast to config/config.yaml and dry-run again. Which new jobs appear, and which existing jobs must be rerun? Why?

Exercise 3¶

A colleague sequenced 12 individual (not pooled) mosquitoes from Siaya, Kenya: 6 that survived a deltamethrin net and 6 that died. What would you change in samples.tsv and config.yaml? Would you switch on PBS?

Exercise 4¶

Which rule writes results/genediff/BusiaParental_BusiaSelected.csv, and which script does the work? Use grep in the terminal.

Quiz¶

Check your understanding. The quiz runs in Colab and in the course book.

Summary¶

• RNA-Seq-Pop analyses both expression (kallisto, DESeq2, sleuth, GO enrichment) and sequence (freebayes, allele frequencies, Fst, PBS, diversity, ancestry, karyotypes) from the same RNA-seq reads.

• The sequence data can find things expression cannot, such as the rise of Vgsc-995S in Busia.

• You configure it with two files: samples.tsv (samples and treatments) and config.yaml (modules, contrasts, references, ploidy). Treatment names must not contain underscores.

• Use a tagged release (v2.3.0), dry-run first, look at the rule graph, then launch with --use-conda in the background.

Well done!¶

Your Busia run is under way. In Module 2.1 we explore the full-dataset results, starting with differential expression.

References¶

1. Nagi SC, Oruni A, Weetman D, Donnelly MJ (2023). RNA-Seq-Pop: Exploiting the sequence in RNA sequencing—A Snakemake workflow reveals patterns of insecticide resistance in the malaria vector Anopheles gambiae. Molecular Ecology Resources 23:946–961. Nagi et al. (2023)

2. Mölder F et al. (2021). Sustainable data analysis with Snakemake. F1000Research 10:33. https://doi.org/10.12688/f1000research.29032.2

3. Bray NL, Pimentel H, Melsted P, Pachter L (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34:525–527. https://doi.org/10.1038/nbt.3519

4. Garrison E, Marth G (2012). Haplotype-based variant detection from short-read sequencing. arXiv 1207.3907. Garrison & Marth (2012)

5. Kim D et al. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 37:907–915. Kim et al. (2019)

RNA-Seq-Pop documentation: https://sanjaynagi.github.io/rna-seq-pop/. GitHub Codespaces documentation: https://docs.github.com/en/codespaces.