Amelia Edmondson-Stait

3 minute read

First post in series.

/assets/images/tests-bug.png

I’ve been developing an R package that formats and plots output from a genetic software (called SuSiEx).

I wrote a few functions that I wanted to put into a package so that my colleagues and other researchers could use these functions if they use the SuSiEx software, and can format and plot their results quickly rather than duplicating the work.

I’ve also been learning about testing, so I have been using the R package testthat to write some unit tests for my functions.

The package is full of bugs and I want to fix these. I also want to refactor the code and get in a better habit of writing smaller functions that make up a larger main function. Then test those smaller functions in the unit tests.

Firstly, I’m just going to focus on one of these functions, plotAncestryCausal(). This function takes the summary results from SuSiEx (ie. the output from formatResults()) and the genetic ancestries used to create the summary results. It outputs a plot showing whether the genetic variant (x-axis) was causal for that ancestry (shape and colour of point) using a post-hoc probability statistic (y-axis). When I last used this function I found a bug, the points are labeled with the incorrect ancestries.

When looking up best practices for writing unit tests I read in several places that when you find a bug, you write a test for it, and then the function (ie. test-driven development). So this is what I aim to do here.

Some other useful advice from Hadley Wickham & Jennifer Bryan’s book on R packages:

“try to avoid putting too many expectations in one test - it’s better to have more smaller tests than fewer larger tests.”

So my aim is to 1) recreate the bug, 2) write & run a test that fails 3) update the function, 4) rerun test and check it passes.

This is the first blog post in a series of posts going through each of these steps.

1. Reproduce the bug 🪲

I’m going to run the script that produced the bug and see if I can figure out what bit of the code might be failing. I’m going to look out for any error messages or warnings. I may also run the function with the test data (that I think shouldn’t produce this bug), to see what the differences are. See the script susiexR/inst/scripts/fix_bug_with_tests.R for details on how I did this. But briefly:

  • Used the debugger, putting break points inside plotAncestryCausal().
  • Discovered how difficult it is to debug when you have very long lines of code. I need to refactor this function into smaller functions.
  • But I also saw that the ancestry columns were not what i was expecting in the data frame summary_results. And that there were fewer rows in the dataframe than I expected.
  • I then went back and checked the data output from SuSiEx and discovered it was not what I was expecting.
  • I recently re-ran SuSiEx and I think somehow the old results were still in the GitHub repo. (I have a suspision this happened during a merge conflict as I was careful to remove old files from the repo before re-running).
  • I have cleared the nextflow cache and am re-running SuSiEx again.
  • Once I am confident I have the correct results from SuSiEx I will re-run the script to see if the same plot occurs.
  • A different plot occured.
  • Inspected the folder with the SuSiEx output in it.

Found the bug! 🥳

  • The files in this folder had ancestry names in different orders, leading to inconsistencies when extracting ancestry information in plotAncestryCausal(). eg.:
fineMapping/output/SuSiEx.AFR-SAS-EUR.output.cs95_7:146868144:147068144.cs
fineMapping/output/SuSiEx.EUR-AFR-SAS.output.cs95_2:22470074:22670074.cs
fineMapping/output/SuSiEx.SAS-AFR-EUR.output.cs95_4:60234365:60434365.cs
  • Since ancestry orders varied, this caused errors when processing the files.
  • Nextflow processes data using channels, which do not preserve a fixed order.
  • This means that every time SuSiEx runs, the ancestries may be in a different order.
  • I thought I had accounted for this in Nextflow by ensuring that SuSiEx ran within a single combined process, but I need to double-check if that logic is correctly enforced.

Reflections so far

  • I discovered this bug when I compared a previously correct plot from an older SuSiEx run with a new output that looked different. Without that comparison, I might never have questioned the plot’s accuracy. Nextflow channels, by design, execute tasks in parallel and don’t guarantee order, so variations in the ancestry order were always a risk I anticipated during pipeline development. Had I implemented tests early on to verify consistent ancestry order, I likely would have caught this issue sooner. This whole process has highlighted the importance of writing tests, and test-driven development.

Resources:

Updated: