search

neherlab / ncov-simple

2 stars 1 forks source link

feat: reconstruct emergence of insertions and place on tree as branch labels #26

Closed corneliusroemer closed 2 years ago

corneliusroemer commented 2 years ago

Resolves #5

corneliusroemer commented 2 years ago

Encoding nextclade insertions output into binary vector works:

[Running]  /usr/local/Caskroom/mambaforge/base/envs/nextstrain/bin/python "/Users/cr/code/ncov-simple/scripts/reconstruct_insertions.py"
{'28257:CTGT': 0, '28251:CTG': 1, '19298:GGTTGTG': 2, '21992:ACT': 3, '29786:TAT': 4, '28263:AACA': 5, '15:C': 6, '3677:T': 7, '13344:ACATGGAA': 8, '19298:GGTTGTGATGGTGGCAGTTTGTGGTTGTG': 9, '21751:TGCAGATGAAAC': 10, '22206:CGGCAGGCT': 11, '29685:T': 12, '26734:T': 13, '5286:ATCGCCA': 14, '3939:G': 15, '29754:G': 16, '10388:T': 17, '29774:GAGAGCTGC': 18, '24359:N': 19, '25555:T': 20}
strain
Chile/LI-115265/2021                 [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
USA/OR-TRACE-BENT-032821-444/2021    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
USA/OR-OHSU-212100289/2021           [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
USA/MA-CDC-ASC210033067/2021         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
DominicanRepublic/Yale-6231/2021     [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
                                                           ...                        
Mexico/ROO-INMEGEN-06-01-333/2021    [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
Peru/ARE-UPCH-0729/2021              [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
USA/TX-DSHS-6055/2021                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
USA/WI-WSLH-218057/2021              [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
Canada/BC-BCCDC-143600/2021          [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
Name: insertions_vector, Length: 65, dtype: object
corneliusroemer commented 2 years ago

@rneher I've set up the sequences correctly I think but treetime can't match my alignment strain names with the tree tip strain names. What am I doing wrong?

Alignment with 476 rows and 21 columns
AAAAAAAAAAAAAAAAAAAAA USA/NJ-CDC-2-4122345/2020
AAAAAAAAAAAAAAAAAAAAA USA/WA-UW-21042246599/2021
AAAAAAAAAAAAAAAAAAAAA Japan/RIMD00399/2021
AAAAAAAAAAAAAAGAAAAAA Chile/LI-115265/2021
AAAAAAAAAAAAAAAAAAAAA India/WB-1931300247197/2021
AAAAAAAAAAAAAAAAAAAAA Philippines/PH-PGC-47340/2021
AAAAAAAAAAAAAAAAAAAAA USA/MN-CDC-IBX178570093152/2021
AAAAAAAAAAAAAAAAAAAAA USA/MO-CDC-ASC210332079/2021
AAAAAAAAAAGAAAAAAAAAA USA/OR-TRACE-BENT-032821-444/2021
AAAAAAAAAAAAAAAAAAAAA Belgium/rega-9437/2021
AAAAAAAAAAAAAAAAAAAAA USA/GA-CDC-STM-000055342/2021
AAAAAAAAAAAAAAAAAAAAA USA/OR-CDC-ASC210047521/2021
AAAAAAAAAAAAAAAAAAAAA India/WB-1931300256989/2021
AAAAAAAAAAAAAAAAAAAAA Denmark/DCGC-129279/2021
AAAAAAAAAAAAAAAAAAAAA Netherlands/ZH-RIVM-12316/2021
AAAAAAAAAAAAAAAAAAAAA CzechRepublic/NRL_10270/2021
AAAAAAAAAAAAAAAAAAAAA USA/NH-CDC-QDX24155544/2021
AAAAAAAAAAAAAAGAAAAAA USA/OR-OHSU-212100289/2021
...
AAAAAAAAAAAAAAGAAAAAA Canada/BC-BCCDC-143600/2021
augur ancestral is using TreeTime version 0.8.4

0.07    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        Wuhan/Hu-1/2019

0.08    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        Philippines/PGC01/2020

0.08    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        Wuhan/WH01/2019

0.08    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        Canada/AB-ABPHL-08857/2021

0.08    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        CostaRica/INC-0049/2020

0.09    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        USA/WA-QDX-2176/2020

...

0.35    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        USA/MO-CDC-ASC210332079/2021

0.35    ***WARNING: TreeAnc._attach_sequences_to_nodes: NO SEQUENCE FOR LEAF:
        USA/CA-CDC-FG-054156/2021
Traceback (most recent call last):
  File "/Users/cr/code/ncov-simple/scripts/reconstruct_insertions.py", line 69, in <module>
    main()
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/click/core.py", line 1128, in __call__
    return self.main(*args, **kwargs)
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/click/core.py", line 1053, in main
    rv = self.invoke(ctx)
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/click/core.py", line 1395, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/click/core.py", line 754, in invoke
    return __callback(*args, **kwargs)
  File "/Users/cr/code/ncov-simple/scripts/reconstruct_insertions.py", line 65, in main
    print(ancestral.ancestral_sequence_inference(tree,alignment))
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/augur/ancestral.py", line 48, in ancestral_sequence_inference
    tt = TreeAnc(tree=tree, aln=aln, ref=ref, gtr='JC69',
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/treetime/treeanc.py", line 162, in __init__
    self._check_alignment_tree_gtr_consistency()
  File "/usr/local/Caskroom/mambaforge/base/envs/nextstrain/lib/python3.8/site-packages/treetime/treeanc.py", line 380, in _check_alignment_tree_gtr_consistency
    raise MissingDataError("TreeAnc._check_alignment_tree_gtr_consistency: At least 30\\% terminal nodes cannot be assigned a sequence!\n"
treetime.MissingDataError: TreeAnc._check_alignment_tree_gtr_consistency: At least 30\% terminal nodes cannot be assigned a sequence!
Are you sure the alignment belongs to the tree?
corneliusroemer commented 2 years ago

It works, can be tested here: https://nextstrain.org/groups/neherlab/ncov/europe/2021-11-11?branchLabel=insertions%20%28internal%29&m=div image

corneliusroemer commented 2 years ago

Closed in favour of #27 as this branch had a messed up source