<jupyter_start><jupyter_text>Finding diverging positions in a single spike protein across samples To recap:     1. I took a few GISAID samples, translated them to proteins        accoridng to open reading frames.     2. I then globally align all these sequences with a single         canonical protein.     3. Parsing the alignment file gives me a one-to-many sequence        and sequences that are >90% identical are selected and written         in a separate file, which is then passed on to this code.   And now, I do the following:         This code will turn each sequence into a numpy array and then calculate the conservation score via Shannon Entropy and plot the standardized normalized scoreon a seaborn line plot The code itselfIf you're looking for the code, my friend Viktor contributed by pulling thecode from the scripts folder. The code can be found here:    https://github.com/nadzhou/CoVPeek/blob/master/scripts/parse_variation.py<jupyter_code># Load the aligned FASTA file
aligned_path = "gisaid_results/mafft_aligned.fasta"
parser1 = DivergenceParser.retrieve_sequence(aligned_path)
fig1 = parser1.plot_variation()
fig1.show()<jupyter_output><empty_output><jupyter_text>Cross-check with another sampleThe plot here shows that for this given GISAID dataset (about 16k sequences in total), the divergence at the 613th position shows the highest amount of divergence. In fact, we can see that the position can have all the possible amino acids in that position. I've done this for another GISAID sample (a comparatively smaller one) that shows that this position does indeed show maximum divergence.<jupyter_code>aligned_path2 = "/home/nadzhou/SEQs/id_MAFFT0.fasta"
# we use the previous dataset, since GISAID data can not be shared
parser2 = DivergenceParser.retrieve_sequence(aligned_path)
fig2 = parser2.plot_variation()
fig2.show()<jupyter_output><empty_output><jupyter_text>ValidationSeems like this 614th position is not new. This has already been reported in the following Biorxiv paper:        https://www.biorxiv.org/content/10.1101/2020.04.29.069054v1.full.pdf         From the article's abstract:    "Here we describe a three-stage data pipeline (analysis of daily GISAID data, structural modeling   of sites of interest, and experimental evaluation) and the identification of several sites of   positive selection, including one (D614G) thatmay have originated either in China or Europe, but begin   to spread rapidly first in Europe, and then in otherparts of the world,and which is now the dominant   pandemic form in many countries."  Contribution The 614th psoition doesn't just show substition between 2 amino acids, but multiple. Here however, I present the result that this 614th position doesn't only contain aspartic acid / glutamine substitution, but this position can be substituted by ALL amino acids.  There's a strong bias against Histidine, Methionine and Tryptophan at the 614th position.As can be seen in the subsequent graph, there's a strong bias against the three amino acids. This seems to be an important result coming from the given samples from across the globe.<jupyter_code>aa1 = parser1.aminoacids_in_variable_positions()[614]
aa2 = parser2.aminoacids_in_variable_positions()[614]
countplots = make_countplots(aa1,aa2)
countplots.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.051475015051173995, 'fraction_numerical': 0.020770620108368453, 'mean_word_length': 4.860670194003528, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 14, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '544488', 'license': 'MIT', 'n_tokens_mistral': 943, 'n_tokens_neox': 851, 'n_words': 412, 'path': 'notebooks/variation_parser.ipynb', 'repo_name': 'nadzhou/covpeek', 'size': '3322'}