Model selection - Revision history

Gorm at 18:27, 15 April 2026

2026-04-15T18:27:05Z

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T15:09:47Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 17:09, 15 April 2026
Line 240:		Line 240:
	: tells the program to keep positions with gaps		: tells the program to keep positions with gaps

	: '''Important:''' codeml requires a phylogenetic tree as input. The tree supplies the topology for the codon-model analysis, while codeml estimates the branch lengths and the remaining model parameters itself. In this exercise, we use the tree produced by ~~'''~~IQ-TREE 3~~'''~~ (~~'''~~gp120.nexus.treefile~~'''~~) as input.		: '''Important:''' codeml requires a phylogenetic tree as input. The tree supplies the topology for the codon-model analysis, while codeml estimates the branch lengths and the remaining model parameters itself. In this exercise, we use the tree produced by IQ-TREE 3 (gp120.nexus.treefile) as input.

	: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.		: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T15:08:58Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 17:08, 15 April 2026
Line 240:		Line 240:
	: tells the program to keep positions with gaps		: tells the program to keep positions with gaps

	: '''Important:''' codeml requires a phylogenetic tree as input. The tree supplies the topology for the codon-model analysis, while codeml estimates the branch lengths and the remaining model parameters itself. In this exercise, use the tree produced by '''IQ-TREE 3''' ('''gp120.nexus.treefile''') as input.		: '''Important:''' codeml requires a phylogenetic tree as input. The tree supplies the topology for the codon-model analysis, while codeml estimates the branch lengths and the remaining model parameters itself. In this exercise, we use the tree produced by '''IQ-TREE 3''' ('''gp120.nexus.treefile''') as input.

	: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.		: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T10:08:22Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 12:08, 15 April 2026
Line 310:		Line 310:
	----		----

	'''Question 12: '''Is M2 better than M1?		'''Question 12: '''Is M2 better than M1? Based on this comparison, is there evidence that some codon sites in gp120 have dN/dS > 1?

	----		----
Line 317:		Line 317:

	'''Examine list of positively selected sites'''		'''Examine list of positively selected sites'''
	: If your M2 is clearly better than M1, then you have evidence for the existence of positively selected sites in the gp120 gene. Now scroll down to the end of the result file and locate a list similar to the one below. Note: this is the "Bayes Empirical Bayes" table, not the "Naive Empirical Bayes" table.

	~~<pre>~~		: If your M2 is clearly better than M1, then you have evidence for the existence of positively selected sites in the gp120 gene. Now scroll down to the end of the result file and locate the "Bayes Empirical Bayes (BEB) analysis" table. Note: use the BEB table, not the Naive Empirical Bayes (NEB) table. (It is not important what the distinction is in this context, but if you are interested: NEB ignores the fact that there is uncertainty about maximum likelihood estimates, especially for smaller data sets - for instance w for some codon is perhaps not exactly 3.046, but could lie in a region around that value - whereas [https://pubmed.ncbi.nlm.nih.gov/15689528/ BEB accounts for that uncertainty].)
	Bayes Empirical Bayes (BEB) analysis
	~~Positively selected sites~~

	~~Prob~~(w>1) mean w		: This table lists codon positions, the posterior probability that each site belongs to the positively selected class (that is, the class with dN/dS > 1), and the estimated mean value of w for the site.

	~~25 A 0~~.~~959* 3.133 +- 0.769~~		'''Question: '''List all sites having more than 95% probability of belonging to the positively selected class.
	~~27 P 0.906 2.990 +- 0.877~~
	~~56 K 0.987* 3.197 +- 0.687~~		----
	~~59 V 0.915 3.032 +~~- ~~0.873~~
	~~78 R 0.637 2.351 +~~- ~~1.129~~
	~~88 K 0.573 2.148 +~~- ~~1.077~~
	~~95 V 0.925 3.046 +~~- ~~0.843~~
	~~...~~
	~~</pre>~~

	: It is not important what the distinction is in this context, but very briefly: NEB ignores the fact that there is uncertainty about maximum likelihood estimates, especially for smaller data sets (for instance w for some codon is perhaps not exactly 3.046, but could lie in a region around that value), whereas [https://pubmed.ncbi.nlm.nih.gov/15689528/ BEB accounts for that uncertainty].		'''Question 14'''

	: ~~This gives you~~ a ~~list~~ of ~~which residues (if any) were found to~~ belong to the positively selected ~~dN/dS~~ class~~. Also listed is the probability that the site really is in the codon class where~~ dN/dS > 1~~, and a weighted average of the w at the site~~.		: Consider your results from Questions 10–13. In model M2, only a fraction of sites belong to the positively selected class, while many sites belong to classes with dN/dS ≤ 1.

	'''Question: '''~~List all sites having more than 95% probability~~ of ~~belonging to the~~ positively selected ~~class~~.		'''Question: '''What does this suggest about the overall pattern of selection acting on gp120? Explain briefly why one might expect a mixture of negatively selected, neutral, and positively selected sites in this gene.

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T09:48:45Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 11:48, 15 April 2026
Line 249:		Line 249:
	codeml		codeml
	: This will start the codeml program using the settings in the file codeml.ctl. Depending on your computer, this will take some minutes to finish. (You may be able to see how the optimization procedure results in progressively better fits: the likelihood increases, meaning that negative log-likelihood decreases, as the fit improves.)		: This will start the codeml program using the settings in the file codeml.ctl. Depending on your computer, this will take some minutes to finish. (You may be able to see how the optimization procedure results in progressively better fits: the likelihood increases, meaning that negative log-likelihood decreases, as the fit improves.)

			: '''Note:''' On some systems the current version of codeml (4.10.10) may print an error message near the end of the run, even though the relevant results have already been written to '''selection.results'''. If this happens, inspect '''selection.results''' before assuming that the run failed completely.

	'''Inspect result file:'''		'''Inspect result file:'''
Line 254:		Line 256:
	nedit selection.results &		nedit selection.results &
	: This file contains a wealth of information concerning your analysis. You may want to turn off line-wrapping to more clearly see the structure of the file (there will be some long lines that will otherwise wrap around). The top part of the file gives an overview of your sequences, codon usage and nucleotide frequencies. You can ignore this information for now, and move on to the interesting part, namely the model likelihoods and parameter values.		: This file contains a wealth of information concerning your analysis. You may want to turn off line-wrapping to more clearly see the structure of the file (there will be some long lines that will otherwise wrap around). The top part of the file gives an overview of your sequences, codon usage and nucleotide frequencies. You can ignore this information for now, and move on to the interesting part, namely the model likelihoods and parameter values.

	: '''Note:''' On some systems codeml may print an error message near the end of the run, for instance during the BEB step, even though the relevant results have already been written to '''selection.results'''. If this happens, inspect '''selection.results''' before assuming that the run failed completely.

	'''Find likelihood, and number of free parameters for model M1'''		'''Find likelihood, and number of free parameters for model M1'''

Gorm at 09:42, 15 April 2026

2026-04-15T09:42:36Z

← Older revision		Revision as of 11:42, 15 April 2026
Line 240:		Line 240:
	: tells the program to keep positions with gaps		: tells the program to keep positions with gaps

	: '''Important ~~note~~:''' in this exercise, use the tree ~~from~~ '''IQ-TREE 3''' ~~for the codeml analysis~~. ~~On some systems, codeml may have problems with the PAUP tree file or may stop with a confusing error message even though the downstream results are still written to the output file~~. ~~Using the IQ-TREE tree avoids this issue~~.		: '''Important:''' codeml requires a phylogenetic tree as input. The tree supplies the topology for the codon-model analysis, while codeml estimates the branch lengths and the remaining model parameters itself. In this exercise, use the tree produced by '''IQ-TREE 3''' ('''gp120.nexus.treefile''') as input.

	: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.		: The settings entered by us will cause codeml to analyze two hypotheses about dN/dS ratios. M1 says there are two classes of codons with different dN/dS ratios in the sequence: one class with dN/dS < 1 (codons under purifying or negative selection), and one class with dN/dS = 1 (no selection - neutrally evolving sites). M2 says there are 3 distinct dN/dS ratios for different sites in the sequence: one class with dN/dS < 1, one class with dN/dS = 1 (these are the same type of classes as for M1), and one class with dN/dS > 1 (corresponding to sites under positive selection). The value of the dN/dS ratios (for those classes that have dN/dS < 1 or dN/dS > 1), the fraction of sites belonging to each class, and the position of sites belonging to each class, are unknown at first and will be determined during the analysis.

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T09:28:30Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 11:28, 15 April 2026
Line 223:		Line 223:
	'''Edit the parameter file'''		'''Edit the parameter file'''
	nedit codeml.ctl &		nedit codeml.ctl &
	: The file "codeml.ctl" contains several settings that control how '''codeml''' is run. Before continuing, edit the file and save it so that the following ~~arguments~~ have ~~exactly~~ these values:		: The file "codeml.ctl" contains several settings that control how '''codeml''' is run. Before continuing, edit the file and save it so that the following settings have these values. Any text appearing after an asterisk (*) on a line is only a comment, and you do not need to change it.
	'''seqfile = gp120align.fasta''': name of alignment file
	'''treefile = gp120.nexus.treefile''': name of tree file produced by IQ-TREE 3		'''seqfile = gp120align.fasta'''
	'''seqtype = 1''': tells the program that ~~our~~ data ~~consists~~ of coding DNA		: name of alignment file
	'''NSsites = 1 2''': tells the program to analyze models M1 and M2
	'''cleandata = 0''': tells the program to keep positions with gaps		'''treefile = gp120.nexus.treefile'''
			: name of tree file produced by IQ-TREE 3

			'''seqtype = 1'''
			: tells the program that the data consist of coding DNA

			'''NSsites = 1 2'''
			: tells the program to analyze models M1 and M2

			'''cleandata = 0'''
			: tells the program to keep positions with gaps

	: '''Important note:''' in this exercise, use the tree from '''IQ-TREE 3''' for the codeml analysis. On some systems, codeml may have problems with the PAUP tree file or may stop with a confusing error message even though the downstream results are still written to the output file. Using the IQ-TREE tree avoids this issue.		: '''Important note:''' in this exercise, use the tree from '''IQ-TREE 3''' for the codeml analysis. On some systems, codeml may have problems with the PAUP tree file or may stop with a confusing error message even though the downstream results are still written to the output file. Using the IQ-TREE tree avoids this issue.

Gorm: /* Detection of positively selected sites in gp120 */

2026-04-15T09:25:14Z

Detection of positively selected sites in gp120

← Older revision		Revision as of 11:25, 15 April 2026
Line 221:		Line 221:
	:* If M2 is better: identify the positively selected codons		:* If M2 is better: identify the positively selected codons

	'''~~Inspect~~ the parameter file'''		'''Edit the parameter file'''
	nedit codeml.ctl &		nedit codeml.ctl &
	: The file "codeml.ctl" contains several settings that ~~are relevant for running the program~~ '''codeml'''. ~~Find~~ the ~~following lines~~ and ~~ensure~~ that the ~~file contains~~ these values:		: The file "codeml.ctl" contains several settings that control how '''codeml''' is run. Before continuing, edit the file and save it so that the following arguments have exactly these values:
	'''seqfile = gp120align.fasta''': name of alignment file		'''seqfile = gp120align.fasta''': name of alignment file
	'''treefile = gp120.nexus.treefile''': name of tree file produced by IQ-TREE 3		'''treefile = gp120.nexus.treefile''': name of tree file produced by IQ-TREE 3

Gorm: /* Tree inference using the model selected by jModelTest2 (PAUP) */

2026-04-15T09:18:50Z

Tree inference using the model selected by jModelTest2 (PAUP)

← Older revision		Revision as of 11:18, 15 April 2026
Line 148:		Line 148:
	'''Question 5'''		'''Question 5'''

	: Above, you used jModelTest2 to select a substitution model based on AIC. You will now use that model to construct a maximum likelihood tree in '''PAUP'''. The purpose of doing this in PAUP is mainly pedagogical: it makes the individual steps of model-based phylogenetic inference more explicit. In practice, one would often use a ~~more integrated~~ program such as IQ-TREE 3 (see next section).		: Above, you used jModelTest2 to select a substitution model based on AIC. You will now use that model to construct a maximum likelihood tree in '''PAUP'''. The purpose of doing this in PAUP is mainly pedagogical: it makes the individual steps of model-based phylogenetic inference more explicit. In practice, one would often use a program such as IQ-TREE 3 (see next section).

	: Close the results table. In the main window you should now scroll up to the lines giving PAUP commands that implement the selected model. The command is enclosed between "BEGIN PAUP" and "END;" and should look something like this:		: Close the results table. In the main window you should now scroll up to the lines giving PAUP commands that implement the selected model. The command is enclosed between "BEGIN PAUP" and "END;" and should look something like this:

Gorm at 09:18, 15 April 2026

2026-04-15T09:18:00Z

Show changes

← Older revision		Revision as of 20:27, 15 April 2026
Line 10:		Line 10:
	: The role gp120 plays in infection and the fact that it is situated on the surface of the HIV particle, means it is an obvious target for the immune response. That means that there may be a considerable selective pressure on gp120 for creating immune-escape mutants, where amino acids in the gp120 epitopes have been substituted. In this exercise you will construct a maximum likelihood tree that we will subsequently use to investigate whether you can detect such a selective pressure on parts of gp120, again using maximum likelihood methods.		: The role gp120 plays in infection and the fact that it is situated on the surface of the HIV particle, means it is an obvious target for the immune response. That means that there may be a considerable selective pressure on gp120 for creating immune-escape mutants, where amino acids in the gp120 epitopes have been substituted. In this exercise you will construct a maximum likelihood tree that we will subsequently use to investigate whether you can detect such a selective pressure on parts of gp120, again using maximum likelihood methods.

	: One major goal with the exercise is to introduce you to statistically based methods for assessing the strength of evidence for a set of alternative hypotheses about some biological system of interest. The model selection method we will use is AIC (Akaike Information Criterion), based on which you will compute model probabilities. A second goal is to make you aware that phylogenetic analysis is not only about constructing trees, but that it is also a useful framework for analyzing biological questions more generally.		: One major goal with the exercise is to introduce you to statistically based methods for assessing the strength of evidence for a set of alternative hypotheses about some biological system of interest. The model selection method we will use is AIC (Akaike Information Criterion), based on which you will compute Akaike weights, which can be interpreted as relative support for the models, or approximately as relative model probabilities within the candidate set. A second goal is to make you aware that phylogenetic analysis is not only about constructing trees, but that it is also a useful framework for analyzing biological questions more generally.

	: Specifically, you will		: Specifically, you will
Line 21:		Line 21:
	----		----

	== Recipe for computing AIC values and model probabilities ==		== Recipe for computing AIC values and Akaike weights (interpretable as relative model probabilities) ==

	: Later in ~~today's~~ exercise you will be asked to compute AIC values and ~~model probabilities~~. Return to this section and follow the instructions when you need to do so.		: Later in today’s exercise you will be asked to compute AIC values and Akaike weights. Return to this section and follow the instructions when you need to do so.

	:# Fit a set of models to your data, note the maximized log likelihoods (lnL) and the number of free parameters (K) for each model in the investigated set. The models you fit should represent a plausible and comprehensive set of hypotheses about your data.		:# Fit a set of models to your data, note the maximized log likelihoods (lnL) and the number of free parameters (K) for each model in the investigated set. The models you fit should represent a plausible and comprehensive set of hypotheses about your data.
Line 30:		Line 30:
	:# Compute the "ΔAIC" values for each model: '''ΔAIC = AIC - AICmin''' <br>For each model subtract the minimum AIC value. The best model will have a ΔAIC of zero. The rest of the models will have positive ΔAICs.		:# Compute the "ΔAIC" values for each model: '''ΔAIC = AIC - AICmin''' <br>For each model subtract the minimum AIC value. The best model will have a ΔAIC of zero. The rest of the models will have positive ΔAICs.
	:# For each model compute the following quantity: '''numerator = exp(-0.5 x ΔAIC)''' <br>For example, a model with ΔAIC=4.2 will have numerator = exp(-0.5 x 4.2) = exp(-2.1) = 0.1225. Also compute the '''sum''' of the numerator values for all models.		:# For each model compute the following quantity: '''numerator = exp(-0.5 x ΔAIC)''' <br>For example, a model with ΔAIC=4.2 will have numerator = exp(-0.5 x 4.2) = exp(-2.1) = 0.1225. Also compute the '''sum''' of the numerator values for all models.
	:# Finally, the ~~model probabilities (also called model weights)~~ for each model ~~are~~ found as:<br> '''P(model) = numerator / sum''' <br>For example, if sum = 3.75 and a model has numerator = 1.3, then it has P(model) = 1.3 / 3.75 = 0.35		:# Finally, the Akaike weight for each model is found as:<br> '''weight(model) = numerator / sum''' <br>For example, if sum = 3.75 and a model has numerator = 1.3, then it has P(model) = 1.3 / 3.75 = 0.35

	: You may want to keep track of the computations by constructing a table along the following lines:		: You may want to keep track of the computations by constructing a table along the following lines:
Line 304:		Line 304:

	'''Assess strength of evidence for models M1 and M2'''		'''Assess strength of evidence for models M1 and M2'''
	: M2 will always have a better (higher) log-likelihood than model M1 because M2 has more free parameters, and M1 is nested within M2. You should now use the recipe given above to compute AIC values and ~~model probabilities~~ for M1 and M2.		: M2 will always have a better (higher) log-likelihood than model M1 because M2 has more free parameters, and M1 is nested within M2. You should now use the recipe given above to compute AIC values and Akaike weights for M1 and M2.

	'''Question: '''Report: AIC, ΔAIC, weight (model probability) for M1 and M2		'''Question: '''Report: AIC, ΔAIC, weight (model probability) for M1 and M2