Measuring evolutionary rates of proteins in a structural context

1. Collect amino acid sequences

One of the most popular methods to collect sequences is BLAST^32,33. To search for orthologous sequences in the NCBI database, determine a query sequence, on which BLAST will base its homology search. The BLAST output will specify the number and percent of sites with matches, near matches, and no matches. BLAST refers to these as identities, positives, and gaps, respectfully. We recommend the specific algorithm PSI-BLAST³³ if one is interested in collected amino-acid sequences (as opposed to nucleotide). For either data type, we recommend that users specify that BLAST query NCBI “RefSeq” (reference sequence)³⁴, which has been heavily curated to contain only nonredundant and reliably annotated sequences. Aside from BLAST, many other approaches, including collecting orthologs from databases such as ENSEMBL³⁵ or UniProt³⁶, are also suitable for this step. Regardless of the approach taken, we recommend that a final dataset contain at least 20 sequences, with more being preferable, to achieve reliable evolutionary rate estimates.

2. Align sequences with MAFFT (see Note 4)

Store all of the sequences you wish to align into one FASTA-formatted file. The FASTA format contains two pieces of information for each sequence: the sequence ID preceded by a ">" sign and followed by a new line, and then the sequence itself. We will use the FASTA file HRH1_unaligned.fasta that contains homologous sequences that are not aligned. We align them with the command:

mafft --auto --inputorder \
       HRH1_unaligned.fasta > \
       HRH1_aligned.fasta

Arguments above correspond to the following:

The output file HRH1_aligned.fasta will contain the aligned sequences.

3. Infer tree with RAxML (see Notes 2, 5)

Using the file with the alignment HRH1_aligned.fasta, run RAxML with the following command:

raxmlHPC-SSE3 -s HRH1_aligned.fasta \
                -n HRH1_tree \
                -m PROTCATLG \
                -p 12345

Arguments above correspond to the following:

The desired tree file is RAxML_bestTree.HRH1_tree.

4. Infer site-wise rates with HyPhy (see Note 6)

We calculate rates with the LEISR method in HyPhy³⁹. To run this method, the file runLEISR.bf must be edited to specify the directories and file names that will be used in the analysis. Edit these two lines of runLEISR.bf

"0": "/path/to/HRH1_aligned.fasta",
"1": "/path/to/RAxML_bestTree.HRH1_tree",

Here, "0" should specify the full path to the alignment file HRH1_aligned.fasta, and "1" should specify the full path to the tree file RAxML_bestTree.HRH1_tree.

Run HyPhy with the command

HYPHYMP runLEISR.bf

An output file HRH1_aligned.fasta.LEISR.json is written to the folder that contains the alignment.

5. Parse HyPhy output (see Note 3)

For further downstream processing, the HyPhy output file in JSON format needs to be converted to CSV format. The custom python script parse_LEISR.py will extract each site’s position, rate, and other site-specific inference information from the JSON output file. Parse the JSON file with the command

python parse_LEISR.py \
-j HRH1_aligned.fasta.LEISR.json \
-r extracted_HRH1_rates.csv

Arguments above correspond to the following:

6. Calculate relative site-wise rates (see Note 7)

As discussed by Jack et al.¹⁰, we recommend calculating relative evolutionary rates by normalizing inferred site-specific rates by their average. In other words, to compute the relative amino-acid rates, calculate the mean rate of the entire sequence and divide each site’s rate by this mean rate. Once normalized, a rate below 1 will indicate a site that evolves more slowly than average. For example, a rate of 0.5 implies that the corresponding site evolves half as quickly as does the average. Similarly, a rate above 1 will indicate a site that evolve more quickly than average. For example, a rate of 2 implies that the corresponding site evolves twice as quickly as does the average.

1. Collect nucleotide sequences

Collect nucleotide sequence using step 1 in Protocol 1, using a nucleotide sequence as the query. Alternatively, the Ensembl database’s Biomart tool³⁵ may represent a more reliable approach for collecting strictly protein-coding sequences. By contrast, the UniProt database³⁶ should not be used, as it contains only protein sequences and lacks clear cross-references to the corresponding nucleotide sequences.

2. Translate codon sequences (see Note 3)

In this section, both codon and amino-acid sequences are required to perform site-wise rate calculations. Store all of the desired nucleotide sequences into one FASTA file. Use our custom script to convert a codon FASTA file to an amino acid FASTA files. We use the FASTA file HRH1_unaligned_codon.fasta that contains homologous nucleotide sequences we wish to translate. Translate with the command:

python translate_aln_codon_to_aa.py \
        -n HRH1_unaligned_codon.fasta \
        -o HRH1_unaligned_aa.fasta

Arguments above correspond to the following:

3. Align amino acid sequences with MAFFT

Align amino acid sequences using step 2 in Protocol 1.

4. Back-translate the amino acid alignment into a codon alignment (see Note 3)

This step requires the original codon sequences and the amino acid alignment. Note that the amino acid alignment is retained, and the script simply inserts corresponding codons in place of amino acids at each column of the alignment. Use this command to back-translate the sequences:

python translate_aln_aa_to_codon.py \
         -a HRH1_aligned_aa.fasta \
         -n HRH1_unaligned_codon.fasta \
         -o HRH1_aligned_codon.fasta
                  

Arguments above correspond to the following:

5. Infer tree with RAxML

The following step is the same as step 3 in Protocol 1. Use the amino-acid alignment file HRH1_aligned_aa.fasta to infer the tree.

6. Infer site-wise rates with HyPhy (see Note 8)

To calculate site-wise dN/dS, we use the Fixed Effects Likelihood (FEL) method in HyPhy²¹. To run FEL in HyPhy, the file runFEL.bf must be edited to specify the directories and file names that will be used in the analysis. Edit the following two lines of the runFEL.bf script:

"1": "/path/to/HRH1_aligned_codon.fasta",
"2":  "/path/to/RAxML_bestTree.HRH1_tree",

Here, "1" should specify the full path to the align- ment file HRH1_aligned_codon.fasta, and "2" should specify the full path to the tree file RAxML_bestTree.HRH1_tree.

Run HyPhy with the following command:

HYPHYMP  runFEL.bf

An output file HRH1_aligned_codon.fasta.FEL.json is written to the folder that contains the alignment file.

7. Parse HyPhy output (see Note 3)

For further downstream processing, the HyPhy output file in JSON format needs to be converted to CSV format. The custom python script parse_FEL.py will extract the site’s position, dN (referred to as ‘beta’ in HyPhy output), dS (referred to as ‘alpha’ in HyPhy output), and other site information outputted by HyPhy:

python parse_FEL.py \
-j HRH1_aligned_codon.fasta.FEL.json \
-r extracted_HRH1_dNdS.csv

Arguments above correspond to the following:

8. Calculate site-specific dN/dS (see Notes 3)

FEL will calculate dS and dN values for all informative sites. However, the FEL method will assign dS = 0 and dN = 0 to sites without any synonymous or non-synonymous substitutions, respectively. When calculating dN/dS at these entirely conserved sites, we recommend to use the value dN/dS = 0. For sites with only one non-gap residue, FEL will similarly assign both dN and dS a value of 0. For those sites, we also recommend to use the value dN/dS = 0.

We provide a custom script that will calculate site-wise dN/dS and will assign dN/dS = 0 to such sites. For simplicity, this script checks for conserved sites in amino acid alignments. That the script considers amino acid rather than codon conservation does not influence rate assignments, as dN/dS = 0 in both cases of fully conserved amino acids and codons. For sites with substitutions, HyPhy’s FEL method, used specifically as recommended here, assigns dS = 1. At those sites, our script calculates site-wise dN/dS by simply dividing site’s dN by site’s dS. The original format of extracted_HRH1_dNdS.csv will not be changed.

python calc_dNdS.py \
                  -a HRH1_aligned_aa.fasta \
                  -r extracted_HRH1_dNdS.csv \
                  -o   processed_HRH1_dNdS.csv

Arguments above correspond to the following:

1. Calculate relative solvent accessibility (RSA) from the PDB file (see Note 3)

We provide a custom script calc_rsa.py that will run the software mkDSSP^27,28, extract absolute solvent accessibilities from the mkDSSP output , and calculate relative solvent accessibilities²³. The first argument is the PDB file, and the second op- tional argument (-o 3rze) is the prefix used for the output files.

python calc_rsa.py 3rze.pdb -o 3rze

This command will generate two output files: 3rze.asa.txt containing the raw mkDSSP output, and 3rze.rsa.csv containing RSA values and secondary structure classifications.

2. Calculate weighted contact numbers (WCN) from the PDB file (see Note 3)

WCN measures amino acid packing density and may be calculated with respect to either the α-carbon or the geometric center of the side-chain^19,40. We provide a custom script that will calculate both types of WCN values, although we strongly recommend using the side-chain WCN values¹⁹. The command line arguments follow the same format as the calc_rsa.py script.

python calc_wcn.py 3rze.pdb -o 3rze

The above command will produce an output file 3rze.wcn.csv that contains both side-chain and α-carbon WCN values for each position in the input PDB file.

1. Generate sequence alignment map (see Note 3)

To map site specific evolutionary rates to residues in a PDB structure, we first align the sequence of amino acids extracted from the PDB structure to the multiple sequence alignment used for rate inference. We provide a script that calls MAFFT to align a PDB sequence to a multiple sequence alignment and reformat the output.

python make_map.py \
                   HRH1_aligned.fasta  3rze.pdb

Running the above command produces a CSV file 3rze.map.csv with four columns. The first and second columns contain the numbered position of a given residue in the alignment used for rate inference and the numbered position of a given residue in the PDB structure, respectively. The numbered positions in the second column are obtained directly from the PDB input file and may therefore include PDB insertion codes (see Note 10). The third and fourth columns contain the single-letter amino acid present in the PDB structure and the PDB chain, respectively. If an amino acid is in the alignment but not in the PDB structure, the PDB position is assigned a value of NA. Likewise, if an amino acid is in the PDB structure but not the alignment, the alignment position is assigned NA.

2. Map rates to structural features (see Note 3)

After mapping the alignment used for rate inference to the sequence of the PDB structure, we merge rates with structural features for each residue. We provide a script that uses the map generated above to combine rates and structural features into a single CSV.

python map_features.py 3rze.map.csv \
-r processed_HRH1_dNdS.csv \ extracted_HRH1_rates.csv \
-f 3rze.rsa.csv 3rze.wcn.csv

Arguments above correspond to the following: