<jupyter_start><jupyter_text>A1196 Conservation Analysis*<NAME>**August 13th, 2019***Note to users of this code**: Let me know! I would definitely be interested in whatever it is being used for.Also, if you publish anything using it, please cite me!Question: **Is A1196/A1427 (*S. cerevisiae* / *E. coli* numbering) of the 18S/16S rRNA well conserved within RSCB deposited ribosome-structures? Is this also true of C1054?** BackgroundThe Weir Lab has been conducting MD simulations of the region of the ribosome within a 35 angstrom sphere surrounding G530 of the 16S/18S rRNA. The structures for these simulations were obtained from the five cryo-EM translocation intermediate structures obtained from [this paper](https://elifesciences.org/articles/14874), out of the Korostelev research group out of UMass Medical School.Within these simulations (presumably), it has been noticed that A1196 (which base stacks with C1054) transiently interacts via Hoogsteen/Watson-Crick base pairing with the second nucleotide of the incoming A-site codon, raising the question of a model in which the codon:anticodon base pairing is extended to 5 base pairs.The goals of this assessment is to find out whether or not the A1427 residue is well-conserved in the [RCSB](http://www.rcsb.org/)-deposited ribosome structures available. Additionally, as C1054 has been shown to [interact with the mRNA](https://academic.oup.com/nar/article/44/13/6036/2457638) during molecular dynamics simulations (as well as during our own Korostelev structure simulations), it would be interesting to look at conservation of this base within structures as well.One would expect that conservation of a residue, either on an intraspecies or interspecies level, would suggest mechanistic importance of the residue to the ribosome's function; so, getting a sense of a particular residue's level of conservation will allow us, by proxy, to get at a residue's functional importance.If we find intraspecies conservation but not interspecies conservation, we might expect that the residue is important for the ribosome's function at a species-specific level. If no conservation is found at all, it is unlikely that the residue has high functional relevance.With all of this background stated before I taint my thoughts with the analysis at hand, let us move forward to the task at hand.  Data Collection Data Source DescriptionThe data I am collecting to conduct my analysis is sourced from the [RCSB](http://www.rcsb.org/), an online database of Protein Data Bank (PDB) files representing the 3D shapes of proteins, nucleic acids, and complex assemblies. RCSB API implementationI'm reusing a RCSB API implementation I made last year to fetch data from the website, defined within `rcsb_api.py`.The benefit of this is that I can get any data I please, if I implement the API logic.The drawback is that any new queries must be implemented. See [PyPDB](http://www.wgilpin.com/pypdb_docs/html/) for a more mature project with more limited functionality.<jupyter_code># All dependencies, including the rcsb_api.py dependencies (of `requests`, `lxml`, `pandas`),
# should be installed for this to work properly.

# Additionally, you should install Clustal Omega to do the alignment later on
# Example:
# 1) Download the executable version from clustal.org/omega/
# 2) Run `sudo mv [CLUSTAL_OMEGA_PATH] /usr/bin/clustalo`
# 3) Run `sudo chmod +x /usr/bin/clustalo` to make it executable

# First, let's import the RCSB Python API I defined
import rcsb_api as rcsb
import pandas as pd

# And let's set the pandas display variables
pd.set_option('display.max_colwidth', 80)
pd.set_option('display.max_rows', 50)
# Then, let's query for all PDB IDs which have 'ribosome' in their title
pdb_ids = rcsb.query_for_pdb_ids(queryTitle="ribosome", verbose=False)
# And let's get all associated metadata with these IDs so we can see 
pdb_ids_w_info = rcsb.get_info_from_ids(pdb_ids)

# # The following lines print the DataFrame in full
# with pd.option_context('display.max_rows', None, 
#                        'display.max_columns', None,
#                        'display.max_colwidth', None):
#     display(pdb_ids_w_info)

# I don't want to keep scrolling past this, so let's actually only print some rows
display(pdb_ids_w_info)<jupyter_output><empty_output><jupyter_text>PDB ID FilteringIt is obvious if you look through the titles that only some of these structures are even of the ribosome, let alone structures containing the 16S/18S rRNA. So, let's try to get a list of what elements these structures contain, such that we can obtain only the structures we want, we need to apply a filter.We only want structures with the 16S/18S rRNA, so let's get the descriptions of the chains of each structure,and use it as a filtering condition.As an example, let's see what the chain descriptions of one structure look like: Let's collect this information for all of our PDB IDs,and try to see how many possible chain names there are.<jupyter_code># Sets containing all unique chain names and descriptions in all structures
# TODO: factor fasta join out of this function, instead do it within notebook logic
# (or add an option <--- this, as jeez it slows the queries down)
chain_name_set = set()
chain_description_set = set()
counter = 0

for (index, pdb_id) in enumerate(pdb_ids):
    # Prints out string to update on progress
    if index % 50 == 0:
        print("Getting chains for %s (%d/%d)" % (pdb_id, index+1, len(pdb_ids)))
    
    # Attempts to get chain descriptions from the API 
    # (HTTP request could fail due to flakiness; will print out if so)
    try:
        chain_df = rcsb.describe_chains(pdb_id)
        for chain_name in chain_df.Name:
            chain_name_set.add(chain_name)
        for chain_description in chain_df.Description:
            chain_description_set.add(chain_description)
    except:
        print("Request for '%s' failed", pdb_id)
# Let's see how many unique chain names and chain descriptions there are
print("There are %d unique chain names" % len(chain_name_set))
print("There are %d unique chain descriptions" % len(chain_description_set))
# What do chain names tend to look like?
display(list(chain_name_set)[1:20])
# What do chain descriptions tend to look like?
display(list(chain_description_set)[1:20])<jupyter_output><empty_output><jupyter_text>It seems at first glance that chain names and chain descriptions are used mostly interchangeably.<jupyter_code># Let's try printing all chain names with "16" or "18" in their names,
# to give us an idea of what kind of filters would work to get all of the 16S/18S rRNA Structures
# (Note: For display purposes, only printing first 10 names)
display(list(filter(lambda x: ("16" in x) or ("18" in x), chain_name_set ))[0:10] + ["etc..."])
# Do the same for chain descriptions
# (Note: For display purposes, only printing first 10 descriptions)
display(list(filter(lambda x: ("16" in x) or ("18" in x), chain_description_set ))[0:10]+ ["etc..."])<jupyter_output><empty_output><jupyter_text>Having looked through all of these chain names and decriptions, I'm officially convinced that the vast majority of 16S/18S rRNA candidates will have `16S` or `18S` in the their descriptionsSo, let's filter our set of structures to just contain those that have 16S or 18S in their descriptions.<jupyter_code># This function returns a bool representing whether or not the
# PDB structure has any of the given substrings in any of its chain names,
# chain descriptions, or both (depending on value of query_fields).
#
# Input:
#   pdb_id : PDB ID of structure we want to evaluate
#   substring_list : List of substrings to query on
#   query_fields : Which pdb chain column(s) we want to filter on
#                  (must be sublist of ["Name", "Description"])
#   num_attempts : The number of times to try the HTTP request to the RCSB API before
#                  giving up (and returning False)
# 
# Output:
#   Whether or not the PDB structure contains any chains satisfying these conditions.
# 
# EG  `chainFieldsContainSubstrings(["16S", "18S"], "5JUP", ["Description"])`
# should return `True`, as it has 18S chain description
# 
# And `chainFieldsContainSubstrings(["16S", "18S"], "5JUP", ["Name", "Description"])`
# should return `True`, as it has 18S chain description, but no 18S chain name
# 
# But `chainFieldsContainSubstrings(["16S", "18S"], "5JUP", ["Name"])`
# should return `False`, as it has no 18S chain description
def chainFieldsContainSubstrings(pdb_id, substring_list, 
                                 query_fields=["Name", "Description"],
                                 num_attempts=3):    
    # Retries the RCSB query multiple times (to reduce flakiness)
    num_attempted = 0
    success = False
    while (num_attempted < num_attempts):
        try:
            chain_df = rcsb.describe_chains(pdb_id)
            success = True
            break
        except:
            num_attempted += 1
    if not success:
        print("Failed to obtain chains from '%s'" % pdb_id)
        return False
    
    # Returns `True` if any of the substrings are in any of the chain names
    if "Name" in query_fields:
        for chain_name in chain_df.Name:
            for substring in substring_list:
                if substring in chain_name:
                    return True
    # Returns `True` if any of the substrings are in any of the chain descriptions
    if "Description" in query_fields:
        for chain_description in chain_df.Description:
            for substring in substring_list:
                if substring in chain_description:
                    return True
                
    # Otherwise, return `False`
    return False
# Let's test this function works properly
(chainFieldsContainSubstrings("5JUP", ["16S","18S"], ["Name", "Description"]) and
 not chainFieldsContainSubstrings("5JUP", ["16S","18S"], ["Name"]) and
 chainFieldsContainSubstrings("5JUP", ["16S","18S"], ["Description"]))<jupyter_output><empty_output><jupyter_text>Now I'm convinced the function actually works as it is supposed to, let's applythe filter to our PDB IDs, such that we only have IDs with 16S or 18S.<jupyter_code>desired_pdb_ids = []
for (index, pdb_id) in enumerate(pdb_ids):
    # Prints out string to update on progress
    if index % 50 == 0:
        print("Filtering structure ''%s' (%d/%d)" % (pdb_id, index+1, len(pdb_ids)))

    if chainFieldsContainSubstrings(pdb_id, substring_list=["16S", "18S"],
                                    query_fields=["Description"]):
        desired_pdb_ids.append(pdb_id)
print(len(desired_pdb_ids))
# Let's see how many 16S/18S rRNA candidates we have in all those structures
# (and which structures, if any, have repeats)
chain_list = []
id_seen_once = set()
id_seen_twice = set()

for (index, pdb_id) in enumerate(desired_pdb_ids):
    # Prints out string to update on progress
    if index % 50 == 0:
        print("Getting chains for %s (%d/%d)" % (pdb_id, index+1, len(desired_pdb_ids)))
    
    # Attempts to get chain descriptions from the API 
    # (HTTP request could fail due to flakiness; will print out if so)
    try:
        chain_df = rcsb.describe_chains(pdb_id)
        for chain_description in chain_df.Description:
            if "16S" in chain_description or "18S" in chain_description:
                chain_list.append((pdb_id, chain_description))
                if pdb_id not in id_seen_once:
                    id_seen_once.add(pdb_id)
                else:
                    id_seen_twice.add(pdb_id)
    except:
        print("Failed to get chains for structure %s" % pdb_id)
print(len(chain_list))<jupyter_output>507<jupyter_text>Notably, 507 != 458. So, let's see which structures have repeats.<jupyter_code>print(list(id_seen_twice))<jupyter_output>458
['4V7T', '2OM7', '4V7S', '3J0E', '6I7O', '1ZC8', '4V6F', '4CXH', '6BOH', '6OLE', '4V6C', '4WR6', '6OLI', '2NOQ', '6OM7', '1TRJ', '5W4K', '6HCF', '4V88', '6B4V', '6BOK', '3J0D', '4WU1', '4CXG', '5J4D', '486D', '5MYJ', '1MVR', '6OM0', '4V6D', '4V6E', '4WQ1', '6HCM', '4V7V', '6OLF']
35<jupyter_text>After a brief look at the FASTA files of those structures, it seems that they mostly have repeated 16S rRNA chains with identical sequences. We'll treat this more rigorously later.<jupyter_code># Now, let's get the associated DataFrame for all ribosome structures with 16S rRNA
# (Additionally, we cut off all structures that aren't at a resolution of at least 3.5 Angstroms)
# -- I got this cutoff from the 'Gene Machine' book by <NAME> on pg. 82
ribosome_df = rcsb.get_info_from_ids(desired_pdb_ids)
ribosome_df = ribosome_df[ribosome_df["resolution"] <= 3.5] # amino acid resolution cutoff
display(ribosome_df)<jupyter_output><empty_output><jupyter_text>Nice! So we now have our 257 candidate ribosome structures! Now, let's attempt to get the FASTA sequences for the 16S rRNA from this, and append them to one big data frame of the chains, the sequences, and their PDB accession IDs!<jupyter_code>acc_chain_data = pd.DataFrame()
for (index, pdb_id) in enumerate(ribosome_df["structureId"]):
    # Prints out string to update on progress
    if index % 50 == 0:
        print("Getting chains for %s (%d/%d)" % (pdb_id, index+1, len(ribosome_df["structureId"])))
    # Fetches chain data (with FASTA sequence) from RCSB
    df = rcsb.describe_chains(pdb_id, get_fasta=True)
    # Adds a PDB_ID column (so we can identify where each chain is from)
    df["PDB_ID"] = pdb_id
    # Filters out all chains that aren't 16S or 18S rRNA
    df = df[df["Description"].str.contains("16S") | df["Description"].str.contains("18S")]
    # Appends the desired rows to accumulated DataFrame
    acc_chain_data = acc_chain_data.append(df, ignore_index = True)
# Drop any duplicate rows
acc_chain_data.drop_duplicates(subset =["Taxonomy","Sequence","PDB_ID"], inplace = True) 
# # Let's take a look at this data
# with pd.option_context('display.max_rows', None, 
#                        'display.max_columns', None,
#                        'display.max_colwidth', 1000):
#     columnsTitles = ['PDB_ID', 'Description', 'Taxonomy', "Name", "Sequence"]
#     frame = acc_chain_data.reindex(columns=columnsTitles)
#     display(frame)<jupyter_output><empty_output><jupyter_text>Having looked through this data, it's obvious that some of it is mislabelled.A chain with a description of `16S rRNA` should not have any protein in it.So, let's apply another round of filtering.<jupyter_code># 16S/18S rRNA should not have amino acids in them
acc_chain_data = acc_chain_data[~acc_chain_data["Sequence"].str.contains("B")]
# Their names shouldn't include proteins
acc_chain_data = acc_chain_data[~acc_chain_data["Name"].str.contains("protein")]
# In case of duplicate chains with different sequences within same PDB_ID, drop all but the first one
acc_chain_data.drop_duplicates(subset="PDB_ID", inplace=True)
# Let's take a look at this data
# with pd.option_context('display.max_rows', None, 
#                        'display.max_columns', None,
#                        'display.max_colwidth', 1000):
#     acc_chain_data["Length"] = acc_chain_data["Sequence"].str.len()
#     columnsTitles = ['PDB_ID', 'Description', 'Taxonomy', "Name", "Length",  "Sequence"]
#     frame = acc_chain_data.reindex(columns=columnsTitles)
#     display(frame)
len(acc_chain_data) # Now we have 247 final sequences
acc_chain_data.head()<jupyter_output><empty_output><jupyter_text>Converting Chain Data to FASTA formatNow we have the data we want to feed into Clustal Omega! The problem is that Clustal Omega takes FASTA files as input. So, let's convert our Pandas DataFrame into that format.<jupyter_code># TODO: Annotate this function
def convert_df_to_fasta(df, file_name):
    fasta_string = ""
    for index, row in df.iterrows():
        sequence_header = ("> PDB: %s | Taxonomy: %s | Name: %s | Descr: %s \n" %
                        (row["PDB_ID"], row["Taxonomy"], row["Name"], row["Description"]))
        fasta_string += sequence_header
        fasta_string += row["Sequence"] + "\n"
        
    file = open(file_name,"w") 
    file.write(fasta_string)
    file.close()
    
convert_df_to_fasta(acc_chain_data, "test_file.txt")
%%bash
head test_file.txt<jupyter_output>> PDB: 1VVJ | Taxonomy: Thermus thermophilus | Name:  | Descr: 16S rRNA 
UUUGUUGGAGAGUUUGAUCCUGGCUCAGGGUGAACGCUGGCGGCGUGCCUAAGACAUGCAAGUCGUGCGGGCCGCGGGGUUUUACUCCGUGGUCAGCGGCGGACGGGUGAGUAACGCGUGGGUGACCUACCCGGAAGAGGGGGACAACCCGGGGAAACUCGGGCUAAUCCCCCAUGUGGACCCGCCCCUUGGGGUGUGUCCAAAGGGCUUUGCCCGCUUCCGGAUGGGCCCGCGUCCCAUCAGCUAGUUGGUGGGGUAAUGGCCCACCAAGGCGACGACGGGUAGCCGGUCUGAGAGGAUGGCCGGCCACAGGGGCACUGAGACACGGGCCCCACUCCUACGGGAGGCAGCAGUUAGGAAUCUUCCGCAAUGGGCGCAAGCCUGACGGAGCGACGCCGCUUGGAGGAAGAAGCCCUUCGGGGUGUAAACUCCUGAACCCGGGACGAAACCCCCGACGAGGGGACUGACGGUACCGGGGUAAUAGCGCCGGCCAACUCCGUGCCAGCAGCCGCGGUAAUACGGAGGGCGCGAGCGUUACCCGGAUUCACUGGGCGUAAAGGGCGUGUAGGCGGCCUGGGGCGUCCCAUGUGAAAGACCACGGCUCAACCGUGGGGGAGCGUGGGAUACGCUCAGGCUAGACGGUGGGAGAGGGUGGUGGAAUUCCCGGAGUAGCGGUGAAAUGCGCAGAUACCGGGAGGAACGCCGAUGGCGAAGGCAGCCACCUGGUCCACCCGUGACGCUGAGGCGCGAAAGCGUGGGGAGCAAACCGGAUUAGAUACCCGGGUAGUCCACGCCCUAAACGAUGCGCGCUAGGUCUCUGGGUCUCCUGGGGGCCGAAGCUAACGCGUUAAGCGCGCCGCCUGGGGAGUACGGCCGCAAGGCUGAAACUCAAAGGAAUUGACGGGGGCCCGCACAAGCGGUGGAGCA[...]<jupyter_text>Success!!! Piping this FASTA file into Clustal OmegaWe're going to use the `subprocess` Python module to allow us to run Clustal Omega in BASH and have its output as accessible within the Python environment.<jupyter_code># TODO: Implement this last part, and then data analysis of the results<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0865236930428352, 'fraction_numerical': 0.018317310427214288, 'mean_word_length': 4.340218712029162, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 46, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '480255', 'license': 'Apache-2.0', 'n_tokens_mistral': 5862, 'n_tokens_neox': 5391, 'n_words': 1930, 'path': 'A1196_Conservation_Analysis.ipynb', 'repo_name': 'lacoperon/A1196_Conservation_Analysis', 'size': '17600'}