Ingest

[j23414]

Description of proposed changes

See commit messages.

Related issue(s)

Testing

• Checks pass

If reviewer wants to run a local check as documented in ingest/README.md and README.md:

# Checkout branch
git clone https://github.com/nextstrain/ebola.git
cd ebola
git checkout ingest

# Test Ingest
cd ingest
nextstrain build .
ls -ltr data

# Test build
cd ..
nextstrain build .
nextstrain view auspice

[huddlej]

Thank you for putting this ingest PR together so quickly, @j23414. I'm sorry it took me so long to review this more closely. Let me know if you want to chat synchronously about any of the comments. Even though there are some bigger topics that fall outside the scope of this PR (e.g., How to best reuse ingest logic across pathogen repos?), there seem to be simple local solutions we can make in the short term.

[huddlej]

I tested this approach using sequences we have in fauna.

# From `master` of ebola repo:
snakemake -j 1 -p data/ebola.fasta

# From ingest directory of ingest branch:
seqkit sample -n 50 ../data/ebola.fasta > local_ebola.fasta

wget https://raw.githubusercontent.com/nextstrain/monkeypox/master/ingest/bin/fasta-to-ndjson -O bin/fasta-to-ndjson
chmod 755 bin/fasta-to-ndjson

./bin/fasta-to-ndjson \
        --fasta local_ebola.fasta \
        --fields strain virus accession collection_date region country division location source locus authors url title journal puburl \
        --separator "|" > data/local_ebola.ndjson

These steps all ran as expected. When I tried to dry-run the ingest after modifying the config file to include this new source, I got the following error:

$ nextstrain build .  -np
Building DAG of jobs...
Relative file path './source-data/geolocation-rules.tsv' starts with './'. This is redundant and strongly discouraged. It can also lead to inconsistent results of the file-matching approach used by Snakemake. You can simply omit the './' for relative file paths.
MissingInputException in line 49 of /Users/jlhudd/projects/nextstrain/ebola/ingest/workflow/snakemake_rules/fetch_sequences.smk:
Missing input files for rule fetch_all_sequences:
data/local_ebola_all.ndjson

I tried renaming the NDJSON file I created above to data/local_ebola_all.ndjson and ran the workflow again. It ran as expected. In cases where the strain names were different between fauna and GenBank for the same accession (e.g., Ebolavirus/H_sapiens_wt/SLE/2014/Makona_J0005 and J0005 for accession KP759628), both records appeared in the output metadata and sequences. While we can disambiguate the metadata records by strain name, the sequences use only the accession, so we end up with duplicate records in sequences.

These duplicate sequences cause augur filter to exit with an error (with a list of duplicate strains) in the top-level workflow. This outcome is probably close to what we would expect, since we want users to know duplicates exist and handle these accordingly. It might be nicer for users for the ingest to warn them about duplicates, but I don't think we have any tools to check for dups currently, do we?

Regarding the issue with the NDJSON filename mismatch above, we could probably address this problem by removing the serotype wildcard from the sources, since we plan to collect all Ebola serotypes into a single pair of metadata and sequence files and split by serotype in the Nextstrain workflow anyway.

[j23414]

Regarding the issue with the NDJSON filename mismatch above,

Done with commit: fc12467

"I don't think we have any tools to check for dups currently, do we?

Would be wonderful if we did...I tend to use smof uniq , uniq_merge.py, or write up a one-off perl scripts to remove duplicate sequences. There are many possible ways to remove sequence duplicates. Was there a de-dups rule you had in mind? Maybe from a nextstrain repo?

While we can disambiguate the metadata records by strain name

I've seen some gnarly strain names in my time (and don't get me started on resequenced Isolates, or sequenced segments from the same strain)... so highly recommend keeping with the Monkeypox method of IDs by Accession. I see you're checking old fauna datasets.... My method of checking fauna was writing a quick script to replace

> strain_name

with its associated:

>Accession_number

I can spec out a quick perl/python script to do the swap where it (1) reads in a metadata file to generate strain_name-to-accession_number hash and (2) read in sequence file to replace strain_name with Accession_number.

In cases where the strain names were different between fauna and GenBank for the same accession (e.g., Ebolavirus/H_sapiens_wt/SLE/2014/Makona_J0005 and J0005 for accession KP759628), both records appeared in the output metadata and sequences.

Ah this is exactly where I use uniq_merge.py to flag accessions with multiple strain names and create a fix-stain-name.pl script. I'll think about it more. Thank you for the feedback.