Measuring evolutionary rates of proteins in a structural context

1. Align sequences with MAFFT

Store all of the sequences you wish to align into one file. The file must be in the FASTA format. 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.

2. Infer tree with RAxML (see Notes 2, 4)

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.

3. Infer site-wise rates with HyPhy (see Note 5)

To run HyPhy, the file runRelativeProtRates.bf must be edited to specify the directories and file names that will be used in the analysis. Edit these two lines of runRelativeProtRates.bf

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

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

Run HyPhy with the command

HYPHYMP runRelativeProtRates.bf

An output file HRH1_aligned.fasta.site-rates.json is written to the folder that contains the alignment.

4. 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_prot_rates.py will extract the site’s position, rate, and other site information outputted by HyPhy. Parse the JSON file with the command

python parse_prot_rates.py \
-j HRH1_aligned.fasta.site-rates.json \
-r HRH1_rates.csv

Arguments above correspond to the following:

5. Calculate relative site-wise rates

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 fast as 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 fast as the average.

1. 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:

2. Align amino acid sequences with MAFFT

Align amino acid sequences using step 1 in Protocol 1.

3. 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. The command to back-translate the sequences is:

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:

4. Infer tree with RAxML

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

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

To run 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" specifies the full path to the alignment file HRH1_aligned_codon.fasta, and "2" specifies 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.

6. 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/dS, 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:

7. Change dN/dS values for conserved sites (see Note 3)

The FEL method as implemented in HyPhy assigns dN/dS = 1 to sites without any synonymous and any non-synonymous substitutions. We recommend to express rate at entirely conserved sites with dN/dS = 0. We provide a custom script that will assign dN/dS = 0 to completely conserved sites. This script will not change extracted_HRH1_dNdS.csv’s original format.

python fix_dNdS_conserved_sites.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 mkdssp^27,28, extract absolute solvent accessibilities, and calculate relative solvent accessibilities²³. The first argument is the PDB file, and the second optional 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^33,19. We provide a custom script that will calculate both types of 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 WCN 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 with four columns. The first column contains, for each residue, the numbered position of that residue in the alignment used for rate inference. The second column similarly contains the numbered position of each residue in the PDB structure. Numbered positions are extracted directly from the PDB input file and may include PDB insertion codes (see Note 8). 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 align rates with structural features. 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 \
        HRH1_rates.csv \
        -f 3rze.rsa.csv 3rze.wcn.csv

Arguments above correspond to the following:

The output from this command provides all the data needed to compute correlations between rates and structural features and corresponding visualizations. See Figure 1 as an example.