anopheles-genomic-surveillance / selection-atlas

Here be dragons.

selection-atlas

Here be dragons.

Development docs:

• Implementation plan
• UI design
• Workflow design

Conda environment management

Assuming you have a recent version of mamba installed.

The file requirements.yml has the dependencies required to build the site. To ensure reproducibility we currently also maintain the file environment.yml which contains an export of the solved environment.

To create and activate an environment on your own computer:

mamba env create --force --file environment.yml
mamba activate selection-atlas

To create and activate an environment on datalab-bespin:

mamba env create --force --prefix=${HOME}/envs/selection-atlas --file environment.yml
conda activate ${HOME}/envs/selection-atlas

If you need to add or upgrade a package, edit requirements.yml. Do not edit environment.yml.

To upgrade environment.yml:

mamba env create --force --file requirements.yml
mamba env export -f environment.yml -n selection-atlas-requirements --override-channels --channel conda-forge --channel bioconda
sed -i "s/selection-atlas-requirements/selection-atlas/" environment.yml

Running the workflow

If running on your local system with GCS caching enabled, you'll need to run the build without any parallelisation:

snakemake -c1

If running on Google Cloud and GCS caching is disabled, you can try running with parallelisation, e.g.:

snakemake -c4

selection-atlas's Issues

Country page, code cell showing

Run a first build for all of Ag3.0

Run a first build of the site for all samples in Ag3.0 and all chromosomes.

Cohort pages add number of collection sites (locations)

In the preamble text on the cohort pages, insert the number of distinct collection sites (locations). This can be found from distinct values of either the "location" column or the "latitude" and "longitude" columns in the sample metadata.

Lock/pin conda environment

Suggest we use conda-lock so we make sure everyone is using exactly the same environment.

Country pages and home page, use representative point for markers in cohort maps

In the map on the home page and the maps in the country pages, use the representative points for plotting markers, which should then bring together all cohorts which share a common admin2 unit.

I.e., there should be exactly one marker for each admin2 unit.

Plot H12 and G123 in cohort pages with Hampel filter applied

We are using the Hampel filter to clear up noise in the H12 scans before we fit signals. But in the cohort pages we show the unfiltered H12 and G123 data. It might be a good idea to show the filtered data, as it will have less isolated outliers and help to focus attention on the stronger signals.

Would probably require some work upstream in malariagen_data -- malariagen/malariagen-data-python#351

Cohort pages add H12 window size calibration plots

In the diagnostics section, output window size calibration plot.

Final tidy up for alpha1

Do an editorial pass over all pages, adding or refining some text where needed, removing darts videos, and generally making the site look presentable to demo.

Slow build caused by pathological xarray concatenation

I'm finding the build is slow, especially when trying to work with joined arm contigs. I've pinned it down to some pathological behaviour within xarray's concat() function I think. I don't fully understand what's happening there, but have hacked a workaround in malariagen/malariagen-data-python#395.

script to copy build to GCS bucket and download back

So that users can re-build the selection-atlas without having to re-run all the analyses, we want to create a script that copies the build to a GCS bucket, and another script to download the data.

Loading cohort pages is slow and causes page unresponsive

Probably the GWSS plots are the main factor here, plotting lots of points.

Cohort page, add hover text on map markers

Small suggestion, on the cohort pages, we could should the number of samples for each collection site as hover text over map markers.

Fix overview plots in chromosome pages

Currently the overview plots showing all the signals on the chromosome pages don't seem to be working.

When we have got them working, do any visual styling and tweaking required so they work OK with the numbers of signals we're getting from the full Ag3.0 build.

Chromosome page, signals plot, some tools getting duplicated

Like the X wheel zoom tool, there seem to be two of them? And wheel zoom not actually working?

Cohort page, fix styling of headers for chromosomes

Currently the sub-headings for the chromosomes within the "Selection scans" section on the cohort page is supposed to be coming out as h3 but is coming out as p for some reason...

ipyleaflet icons should not be draggable

Cohort page, suggestion for how to present collection site and date information

Currently on the cohort pages we have a map of collection locations and a bar chart of number of samples by month.

The bar chart is visually a bit odd because it's big but doesn't convey much information. Also we are not showing a breakdown of numbers of samples by location and month.

Suggest to replace the bar chart with a pivot table with collection locations as rows and months as columns, showing the numbers of samples.

Dodgy signals

I'm seeing a few dodgy looking signals, inferred where there is no sign of a peak. E.g.:

Decide GWSS methods

Issue to discuss which GWSS methods to run.

Compact stacking of selection signal bars in chromosome pages

Stack up the bars only when they overlap, to reduce the amount of vertical space needed to show signals from all cohorts.

Initial zoom on home page map

Set the initial zoom on the home page map to zoom out a little, so that all Africa is visible. Currently some of west africa gets cut off...

Cohort pages add study and contributor information

In the cohort pages preamble text, add information on the study (and investigator?) who contributed the specimens.

Reorder selection alerts page

1. Signals
2. Cohorts affected
3. Candidate genes
   3a. information on candidates

"CONTIG" KeyError in H12 calibration

When trying to build, I'm getting an error from within malariagen_data - not quite sure whats going on...

I am using the new environment. Have tried removing the previous cache and the old build/ folder. The error also occurs with both joined contigs (2RL) and not joined (2R).

CacheMiss                                 Traceback (most recent call last)
File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:5846, in AnophelesDataResource.h12_calibration(self, contig, analysis, sample_query, sample_sets, cohort_size, min_cohort_size, max_cohort_size, window_sizes, random_seed)
   5845 try:
-> 5846     calibration_runs = self.results_cache_get(name=name, params=params)
   5848 except CacheMiss:

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:853, in AnophelesDataResource.results_cache_get(self, name, params)
    852 if not results_path.exists():
--> 853     raise CacheMiss
    854 results = np.load(results_path)

CacheMiss: 

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
Cell In[13], line 1
----> 1 ag3.plot_h12_calibration(
      2     contig=h12_calibration_contig,
      3     analysis=phasing_analysis,
      4     sample_sets=sample_sets,
      5     sample_query=sample_query,
      6     min_cohort_size=min_cohort_size,
      7     max_cohort_size=max_cohort_size,
      8     window_sizes=window_sizes,
      9 );

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:5915, in AnophelesDataResource.plot_h12_calibration(self, contig, analysis, sample_query, sample_sets, cohort_size, min_cohort_size, max_cohort_size, window_sizes, random_seed, title, show)
   5889 @doc(
   5890     summary="Plot h12 GWSS calibration data for different window sizes.",
   5891     parameters=dict(
   (...)
   5913 ) -> gplt_params.figure:
   5914     # get H12 values
-> 5915     calibration_runs = self.h12_calibration(
   5916         contig=contig,
   5917         analysis=analysis,
   5918         sample_query=sample_query,
   5919         sample_sets=sample_sets,
   5920         window_sizes=window_sizes,
   5921         cohort_size=cohort_size,
   5922         min_cohort_size=min_cohort_size,
   5923         max_cohort_size=max_cohort_size,
   5924         random_seed=random_seed,
   5925     )
   5927     # compute summaries
   5928     q50 = [np.median(calibration_runs[str(window)]) for window in window_sizes]

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:5849, in AnophelesDataResource.h12_calibration(self, contig, analysis, sample_query, sample_sets, cohort_size, min_cohort_size, max_cohort_size, window_sizes, random_seed)
   5846     calibration_runs = self.results_cache_get(name=name, params=params)
   5848 except CacheMiss:
-> 5849     calibration_runs = self._h12_calibration(**params)
   5850     self.results_cache_set(name=name, params=params, results=calibration_runs)
   5852 return calibration_runs

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:5867, in AnophelesDataResource._h12_calibration(self, contig, analysis, sample_query, sample_sets, cohort_size, min_cohort_size, max_cohort_size, window_sizes, random_seed)
   5854 def _h12_calibration(
   5855     self,
   5856     contig,
   (...)
   5865 ):
   5866     # access haplotypes
-> 5867     ds_haps = self.haplotypes(
   5868         region=contig,
   5869         sample_sets=sample_sets,
   5870         sample_query=sample_query,
   5871         analysis=analysis,
   5872         cohort_size=cohort_size,
   5873         min_cohort_size=min_cohort_size,
   5874         max_cohort_size=max_cohort_size,
   5875         random_seed=random_seed,
   5876     )
   5878     gt = allel.GenotypeDaskArray(ds_haps["call_genotype"].data)
   5879     with self._dask_progress(desc="Load haplotypes"):

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:5717, in AnophelesDataResource.haplotypes(self, region, analysis, sample_sets, sample_query, inline_array, chunks, cohort_size, min_cohort_size, max_cohort_size, random_seed)
   5715 debug("normalise parameters")
   5716 sample_sets = self._prep_sample_sets_param(sample_sets=sample_sets)
-> 5717 resolved_region = self.resolve_region(region)
   5718 del region
   5720 if isinstance(resolved_region, Region):

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anopheles.py:1736, in AnophelesDataResource.resolve_region(self, region)
   1731 @doc(
   1732     summary="Convert a genome region into a standard data structure.",
   1733     returns="An instance of the `Region` class.",
   1734 )
   1735 def resolve_region(self, region: base_params.region) -> Region:
-> 1736     return resolve_region(self, region)

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/util.py:436, in resolve_region(resource, region)
    433     raise TypeError("The region parameter must be a string or Region object.")
    435 # check if region is a whole contig
--> 436 if region in _valid_contigs(resource):
    437     return Region(region, None, None)
    439 # check if region is a region string providing coordinates

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/util.py:406, in _valid_contigs(resource)
    404 def _valid_contigs(resource):
    405     """Determine which contig identifiers are valid for the given data resource."""
--> 406     valid_contigs = resource.contigs
    407     # allow for optional virtual contigs
    408     valid_contigs += getattr(resource, "virtual_contigs", ())

File ~/projects/selection-atlas/.snakemake/conda/becd6af5994f9e79c1329307d32adf4a_/lib/python3.10/site-packages/malariagen_data/anoph/genome_sequence.py:23, in AnophelesGenomeSequenceData.contigs(self)
     21 @property
     22 def contigs(self) -> Tuple[str, ...]:
---> 23     return tuple(self.config["CONTIGS"])

KeyError: 'CONTIGS'

cc @alimanfoo

Selection alert for Vgsc locus

Create a selection alert page for the Vgsc locus.

As part of this, create any necessary utility functions that will be used on all selection alert pages, e.g., a function to show all selection signals overlapping the locus of interest.

Cohort page, some code cells are showing

Need to remove input on a few code cells, e.g.:

Cohort page, add overlay of inferred signals to H12 scan plots

In the cohort pages, on the H12 scan plots, overlay the inferred selection signals somehow. E.g., similar to the original selection atlas prototype, could use color shading to indicate focus and span:

...perhaps with some stats like delta I in hover text.

Chromosome page, styling tweaks to the plot of signals by chromosome

For the plot of selection signals, suggest a few styling tweaks:

• Reduce row height, to allow for situations where lots of signals get stacked up at the same location.
• Add some kind of color legend (we are coloring by taxon)
• Consider using diamonds instead of rectangles maybe?
• Currently possible to zoom out beyond end of chromosome...

Use joined arm contigs

Update config to use contigs 2RL, 3RL, X.

Error in cohort page for CI-LG_Agneby-Tiassa_colu_2012

The CI-LG_Agneby-Tiassa_colu_2012 cohort seems to have a -1 for the month column in the metadata (based on the error at least, I haven't checked it yet). This gives an error when doing pandas datetime stuff on the cohort page.

The error:

cols = pd.MultiIndex.from_tuples(
     13     [("Location", "Name"), ("Location", "Longitude"), ("Location", "Latitude")] + 
---> 14     [("Month", pd.to_datetime(x, format="%m").month_name())

ValueError: time data "-1" doesn't match format "%m", at position 0. You might want to try:
    - passing `format` if your strings have a consistent format;
    - passing `format='ISO8601'` if your strings are all ISO8601 but not necessarily in exactly the same format;
    - passing `format='mixed'`, and the format will be inferred for each element individually. You might want to use `dayfirst` alongside this.

Upgrade malariagen_data to 7.4.0

When it's released, upgrade malariagen_data to 7.4.0 in order to bring in support for joined arm virtual contigs, and IHS and G123.

Cohort labels are long for some countries

E.g.:

Maybe we should use the ISO alpha-2 country codes in the cohort labels instead of the full country names.

Deploy alpha1

After having resolved all other issues in the alpha1 milestone, do a full build of the site and deploy to github pages.

Deploy to a sub-url:

https://anopheles-genomic-surveillance.github.io/selection-atlas/alpha1/

Selection alert for Cyp6aa/p locus

Create a selection alert page for the Cyp6aa/p locus.

Set up to support contribution of selection alerts

Make it possible for someone to add a new selection alert page without having to run a full build of all the underlying selection scans and signal detection workflow.

Error in cohort pages when no signals were found by h12-signal-detection

It seems that the h12-signal-detection notebook saves an empty file with no column names when no signals are found for a given cohort's contig.

This seems to cause an error when running the cohort page notebook for that cohort. For example, in the MZ-I_Morrumbene_gamb_2004_Q1 cohort, we get the following error:

----> 3 df_signals = [
      4     pd.read_csv(here() / "build/h12-signal-detection/" / f"{cohort_id}_{contig}.csv")
      5     for contig in contigs
      6 ]

EmptyDataError: No columns to parse from file

Limit zoom out on selection signals plots

Configure limits and zooming on selection signals plots so you can't zoom out beyond ends of the contig.

Add caching of GCS data in github actions

Could speed up CI potentially by caching data from GCS. See here for an example.

Chromosome pages, improve table of selection signals

Consider which columns to show in the table of selection signals.

Also potentially add hyperlinks to cohort pages.

Add favicon

Add a pretty favicon - either a mosquito or DNA probably

Add workflow rules for running G123 and IHS GWSS

Currently for H12 we have a dedicated notebook and a task in the workflow to run the GWSS for all contigs and cohorts.

However, we don't have equivalent tasks for G123 or IHS.

Later, when building the cohort pages, IHS and G123 scans get run because they are requested when building the plots. This is OK, but it might be nice to have dedicated tasks so we can see when those scans are actually being run during a snakemake build.

Intermittent "Address already in use" errors when running the workflow

I'm getting occasional errors when running the snakemake workflow with multiple cores, I believe it is the same as this: nteract/papermill#511

Seems fine to workaround by relaunching the pipeline, but thought I'd record here just in case we're wondering about this in future.

Add G123 selection scans

Run G123 window size calibration and selection scans, and add plots to the cohort pages.

Use full width for GWSS plots in cohorts page

Add G123 window size calibration

Currently the size runs G123 GWSS using the window size chosen from H12 calibration. However, the G123 calibration curves generally look a bit different, and indicate a different window size should be used for G123. Suggest to add a separate G123 window size calibration step to the workflow.

Btw, it's probably safe to assume that if a cohort passes H12 calibration then it will also pass G123 calibration, and therefore filtering of final cohorts may still only need to inspect the results of H12 calibration. Something to consider though, as there could be an edge case where H12 calibration passes but G123 calibration fails, in which case the build would break.

Error in `fit_exponential_peak()` in CD-NU_Gbadolite_gamb_2015_Q3-3RL H12-signal-detection

The offending line:

loc_max = np.argmax(y_peak1)

• Potentially occurs when there are no signals found.

Add IHS selection scans

Run IHS scans and add them to the cohort pages.

Cohort pages add G123 window size calibration plot

Maybe add G123 window size calibration plot in the diagnostics section.

(Or can we use the H12 window size for G123 as well, do we need to calibrate separately?)

Improve table of cohorts in country pages

In the table of cohorts, put a hyperlink under the cohort label column, and remove the cohort page (ID) column.

Also add quarter column.

Preparatory work in malariagen_data

In order to run GWSS for this project, we'll need some work upstream in the malariagen_data package. Creating this issue to collect ideas for things we'll need in malariagen_data.

• Ability to run GWSS across whole chromosomes, i.e., 2RL and 3RL. This is particularly useful to see the signals associated with Vgsc, which usually span the join between 2R and 2L. Can be achieved by concatenating data from the two chromosome arms then offsetting coordinates for the left arm. Will need this support at least in all GWSS functions we'll want to run. Raised upstream: malariagen/malariagen-data-python#332
• Implementation of G123. Raised upstream: malariagen/malariagen-data-python#305
• Implementation of IHS. Raised upstream: malariagen/malariagen-data-python#310

We probably don't want to use XPEHH or PBS for this project, because of the manual effort required to select suitable comparison populations. However, PBS is coming upstream anyway, via malariagen/malariagen-data-python#295.