SARS-Cov-2 Recombinant Finder for fasta sequences

Right now the output is good for interactive human analysis, but there's a lack of csv/tsv machine readable output for sharing/further analysis.

From my experience with Nextclade, main difficulty here is the design of the specs of the file, which columns to include etc, which separators to take if you need an intra-column separator etc.

Maybe best to discuss on this issue before implementing something as one will kind of get locked in to the format.

See this comment by @AngieHinrichs which even contains an alternative.

Thanks a lot for your detailed explanation! I'm trying to move this over here so it's easier to find for me.

(Also, if the comment thread over at pange-designation gets locked down after too many "off topic" comments, I won't be able to comment there at all. Already happened in other issues.)

There's a bit of a problem with using covSpectrum's current mutation API implementation: Ns in any sample is treated as reference.

This can cause confusion. For example, I thought that this intermission here within Spike was a bad sign: https://github.com/cov-lineages/pango-designation/issues/498

But it isn't! Both 22813 and 22882 are defining for both BA.1 and BA.2. However, both are apparently N in 40% of sequences in BA.1. Causing sc2rf to think that it's in fact not a defining mutation in BA.1 making spurious intermissions appear.

I'm not sure how to work around this best. Really, this should be fixed in covSpectrum: Ns should be left out of mutation proportion calculations - and not be treated as reference (implicitly).

@chaoran-chen can you think of a workaround? How can one get the share of Ns for a query? Could that maybe be supplied by a new API endpoint?

Usually, Ns don't make up 40% of a site, but sometimes they do and that can cause problems like here, where one falsely thinks there's a non-clean breakpoint.

This is a fantastic tool, and I've already put it to good use in Arkansas to research some strange lineages. Great work!

I do have to share the visuals, and I an wondering if there is a way to pipe the results to an outside file, such as png or txt. I am more of an applied researcher, so if I missed something, I would appreciate any directions.

Again, great tool already!

Thanks,

I'm analyzing one sequence and am wondering why you output all potential donors/parents, not just the two that seem most relevant here: BA.1/21J?

Are my arguments wrong? When I reduce parents to 0-5, I get not output which is weird. Don't quite understand what's going on here.

Originally posted by @Vjimenez-vasquez in https://github.com/lenaschimmel/sc2rf/issues/25#issuecomment-1089053922:

Hi there,

I ran the following command :

python3 sc2rf.py test2.fasta --unique 1

And got the following message :

Traceback (most recent call last):
  File "sc2rf.py", line 987, in <module>
    main()
  File "sc2rf.py", line 132, in main
    reference = read_fasta('reference.fasta', None)['MN908947 (Wuhan-Hu-1/2019)']
  File "sc2rf.py", line 476, in read_fasta
    with my_tqdm(total=os.stat(path).st_size, desc="Read " + path, unit_scale=True) as pbar:
  File "sc2rf.py", line 199, in my_tqdm
    return tqdm(*margs, delay=0.1, colour="green", disable=bool(args.hide_progress), **kwargs)
  File "/home/hp/anaconda3/lib/python3.7/site-packages/tqdm/std.py", line 922, in __init__
    TqdmKeyError("Unknown argument(s): " + str(kwargs)))
tqdm.std.TqdmKeyError: "Unknown argument(s): {'delay': 0.1, 'colour': 'green'}"

Do you have any suggestion, please ?

Shouldn't be difficult, you need a setup.py and account of Pypi

You can have a look at this repo of mine that can be installed via Pypi as a command line tool (if you install it, the command becomes automatically available in Path!) https://github.com/corneliusroemer/fasta_zstd_sqlite/blob/master/setup.py

Hey!

First of all, great tool to find the potential recombinants. Made my life easy. I needed to parse the output of sc2rf only to get the potential recombinant sequences and the breakpoints of it. I see the --csvfile option in the README. But, it must not have been included in the sc2rf python executable. I get this error.

sc2rf.py: error: unrecognized arguments: --csvfile output.csv

Any idea if I could get the ouput in the way I need?

Would be nice to see how things are going

tqdm makes this very easy with python

a bit more logging while the analysis is going would be cool too, just so that I know what's going on, instead of seeing nothing for a minute

Thanks for the tool.

Just had a quick note that I think Python 3.9 is required due to the | operator in dict.

I was getting an error before trying it with 3.9.

Getting this error while trying to run the program:

Reading reference genome, lineage definitions... Done. Reading actual input. Traceback (most recent call last): File "search_recombinants.py", line 539, in <module> main() File "search_recombinants.py", line 96, in main all_samples = all_samples | read_samples TypeError: unsupported operand type(s) for |: 'dict' and 'dict'

Hi there,

I just was wondering why I have no output and tried the second example from here: https://github.com/lenaschimmel/sc2rf#no-output--some-sequences-not-shown

So I added --clades all --force-all-parent to my call, but it seems that they can't be used both:

The number of allowed parents, the number of selected clades, and the --force-all-parents conflict so that the results must be empty.

Also, --clades all can't be used as the last argument (before the input) because the input won't be recognized

Input sequences must be provided, except when rebuilding the examples. Use --help for more info. Program exits.

I'm not sure if this is only my setup/input problem.

Would you suggest to use -c all or -f? My full command is

  python3 sc2rf.py --csvfile ../${name}_sc2rf.csv --parents 1-35 --breakpoints 1-2 \
                      --max-intermission-count 3 --max-intermission-length 1 \
                      --unique 1 --max-ambiguous 10000 --max-name-length 55 \
                      ### --clades all  --force-all-parents  \ ###
                      ../${fasta}

Best Marie

First, thanks to the authors for bringing the useful tool for us.

We have been using sc2rf to scan for recombinant sequences and determine breakpoint, but i found from the result to the Pangolin X* lineage calls there is a gap. I was wondering whether it is possible to bridge the gap by: 1. take in the lineage designation from Pangolin X* lineages, scan and store the profiles for each of the recombinant lineages; 2. for a new query sequence, if the breakpoint profile matches existing Pangolin X* lineages, in the result not just suggest the parent lineages and breakpoint, provide a possible X* lineage call as well. More or less in the way of how the Scorpio Constellation works.

I expect this would be a more accurate way of assigning recombinant lineages than the current UShER calls, where the breakpoint positions may not match.

Thanks for considering the suggestion.

Hi, I've noticed that sc2rf.py (version sc2rf-7427d2f94b69c965362034c2597b643c5dfaa1cf) could not find any recombination for XT samples available on GISAID python sc2rf.py nextclade.aligned_XT_Gisaid.fasta. Here are the available aligned sequences. nextclade.aligned_XT_Gisaid.txt

Nextclade: sc2rf:

Thanks for looking into this and other lineages that might be in the same situation.

• I've been experimenting with a flag --ignore-shared that ignores positions that are shared (have the exact same nucleotide) across all parents/examples.
• I like this option because it makes the breakpoints visually clearer, as there's a direct color change (red -> green) rather than having the intermediate shared positions (red -> white -> green)
• For testing, a nextclade fasta alignment of XM-like recombinants (public on genbank): XM.txt

1. Do you think this is scientifically sound for reporting? And if so,
2. Would you be interested in a PR if I tidy up the code?

Default Output:

python3 sc2rf.py XM.fasta --ansi --unique 1

Proposed Option:

python3 sc2rf.py XM.fasta --ansi --unique 1 --ignore-shared

While investigating https://github.com/cov-lineages/pango-designation/issues/590, I noticed that samples with the BA.2 S2M deletion (29734:29759) were being incorrectly visualized as having reference bases in sc2rf:

Consensus View:

sc2rf View:

I think this could be for a couple of reasons:

1. When --enable-deletions is used, perhaps deletions should not be considered missing data?

   missings_matches = ["N"]
   if not args.enable_deletions:
       missings_matches.append("-")
   

2. I think there is missing logic when detecting a run of Ns, to catch if that runs proceeds to the end of the genome?

   if s in missings_matches:
       # we've been tracking a run of N's, this base marks the end              
       if start_n == -1:
           start_n = i  # mark the start of possible run of N's
   elif start_n >= 0:
       missings.append((start_n, i-1))  # Python-style (closed, open) interval
       start_n = -1
   
   # Missing logic to catch missing data at the end of the genome?
   if i == len(reference) and s in missings_matches:
       missings.append((start_n, i-1))
   

With these changes, the sc2rf output more closely matches the consensus sequence/my expectation:

I think this is a bug, but if it's the intended behaviour for deletions, please let me know. Thanks!