Sequence-based mutation analysis TSD

There was only one catch and that was Catch-22, which specified that a concern for one's own safety in the face of dangers that were real and immediate was the process of a rational mind. Orr was crazy and could be grounded. All he had to do was ask; and as soon as he did, he would no longer be crazy and would have to fly more missions. Orr would be crazy to fly more missions and sane if he didn't, but if he was sane, he had to fly them. If he flew them, he was crazy and didn't have to; but if he didn't want to, he was sane and had to. Yossarian was moved very deeply by the absolute simplicity of this clause of Catch-22 and let out a respectful whistle.

"That's some catch, that Catch-22," he observed.

"It's the best there is," Doc Daneeka agreed.

-Catch 22

The journal for this task can be found here.

Mutations

Dataset

<figure id="fig:snpsOnstr">

SNPs highlighted on the HexA subunit structure

</figure>

The following SNPs, selected by an unbiased source, will be analysed: M1V, L39R, C58Y, L127R, R170W, R178H, S210F, D258H, L451V and E482K.

Chemical properties

The biochemical properties of the wildtype and mutant amino acids of the chosen SNPs ale listed in <xr id="tab:biochem"/>. Displayed are the hydrophobicity in form of the hydropathy index (Kyte Doolittle) and the according category, the volume with the matching characterisation, the charge and the grantham score. All are taken from <ref name="snpdbe">http://www.rostlab.org/services/snpdbe/</ref>.
The Grantham scores predicts the effect of substitutions between amino acids based on chemical properties, including polarity and molecular volume. It categorizes codon replacements into classes of increasing chemical dissimilarity, and it ranges from 5 to 215<ref name="grantham">Grantham R. Amino acid difference formula to help explain protein evolution. Science 1974; 185: 862-864 </ref>.
The highest Grantham score has the SNP C58Y where both amino acids have a high difference in volume. The change from glutamic acid to lysine on position 482 seems very severy as the charge changes from negative to positive. Altogether in 6 out of 10 cases the charge is different in wildtype and mutant. In 4 mutations the hydrohpobicity is changed from polar to non-polar or the other way round.

TODO

<figtable id="tab:biochem">

Table 1: Biochemical properties

Mutation	Wildtype			Mutant			Grantham score
	Hydrophpbicity	Volume	Charge	Hydrophpbicity	Volume	Charge
M1V	1.9 (nonpolar)	162.9 (bulky)	neutral	4.2 (nonpolar)	140.0 (small)	neutral	21
L39R	3.8 (nonpolar)	166.7 (bulky)	neutral	-4.5 (polar)	173.4 (bulky)	positive	102
C58Y	2.5 (polar)	108.5 (small)	neutral	-1.3 (polar)	193.6 (bulky)	neutral	194
L127R	3.8 (nonpolar)	166.7 (bulky)	neutral	-4.5 (polar)	173.4 (bulky)	positive	102
R170W	-4.5 (polar)	173.4 (bulky)	positive	-0.9 (nonpolar)	227.8 (bulky)	neutral	101
R178H	-4.5 (polar)	173.4 (bulky)	positive	-3.2 (polar)	153.2 (bulky)	neutral	29
S210F	-0.8 (polar)	89.0 (tiny)	neutral	2.8 (nonpolar)	189.9 (bulky)	neutral	155
D258H	-3.5 (polar)	111.1 (small)	negative	-3.2 (polar)	153.2 (bulky)	neutral	81
L451V	3.8 (nonpolar)	166.7 (bulky)	neutral	4.2 (nonpolar)	140.0 (small)	neutral	32
E482K	-3.5 (polar)	138.4 (bulky)	negative	-3.9 (polar)	168.6 (bulky)	positive	56

</figtable>

Structural observations

<figtable id="tbl:pymolstruc">

L39R	C58Y	L127R	R170W	R178H
S210F	D258H	L451V	E482K

</figtable>

 Using your secondary structure predictions from the previous tasks, investigate whether the mutations are inside secondary structure elements (Helix, Strand) or not.

Substitution matrices

BLOSUM62 and PAM

<xr id="tab:pamnlos"/> shows the substitution scores of all mutations for the BLOSUM62, PAM1 and PAM250 matrix. PAM1 and PAM250 are the extremes for this type of substitution matrix and likely not used often in reality <ref name="asd">TODO</ref>. BLOSUM62 however is one of the standard matrices used by default in large tools like BLAST. The BLOSUM62 matrix used here is the incorrectly calculated version, since it generally preforms better <ref name="wrongblosum">Styczynski, M. P., Jensen, K. L., Rigoutsos, I., & Stephanopoulos, G. (2008). BLOSUM62 miscalculations improve search performance. Nature biotechnology, 26(3), 274–275. Nature Publishing Group. Retrieved from http://www.nature.com/nbt/journal/v26/n3/full/nbt0308-274.html</ref>.

An entry in the PAM matrix describes the probability that the wt amino acid is replaced by the mt amino acid in a given timeframe, e.g. 1PAM or 250PAM. TODO Scores within a BLOSUM are log-odds scores that measure, in an alignment, the logarithm for the ratio of the likelihood of two amino acids appearing with a biological sense and the likelihood of the same amino acids appearing by chance.[2]

<figtable id="tab:pamnlos">

Table TODO: Substitution scores for the SNPs based on the matrices BLOSUM62, PAM1 and PAM250. Values for the PAM matrices have been multiplied by 10000 for better readability. Note that the values for BLOSUM62 are those from the incorrectly calculated matrix <ref name="wrongblosum"/>.

Mutation	BLOSUM62 <ref name="blosumvalsrc">http://www.ncbi.nlm.nih.gov/Class/BLAST/BLOSUM62.txt</ref>	PAM (1<ref name="pamvalsrc1">http://www.icp.ucl.ac.be/~opperd/private/pam1.html</ref>/250<ref name="pamvalsrc250">http://www.icp.ucl.ac.be/~opperd/private/pam250.html</ref>)
M1V	1	17/400
L39R	-2	1/200
C58Y	-2	3/300
L127R	-2	1/200
R170W	-3	2/200
R178H	0	8/500
S210F	-2	2/200
D258H	-1	3/400
L451V	1	11/1500
E482K	1	7/800

</figtable>

In a BLOSUM62 matrix, the values range between -4, a very unlikely substitution to 11, the most likely substitution. With this in mind it becomes apparent that most SNPs describe rare changes. The most likely ones only have values of 1 which does not seem much, compared to the maximum of 11, however the 75% quantile of all values in the matrix is only 0, so a value of 1 is not as low as it might seem (TODO I'm writing strange stuff here, can I say that and what do the values truly mean?). On the other hand a value of 1 seems a lot considering that E482K is a change from a negatively charged amino acid to a positively charged one.

In PAM1 and PAM250, the minimum and maximum values are 0,9976 and 9,72 respectively. The high maximum especially for PAM1, stems from the diagonal of the matrix; the 90% quantile is only 22. As to be expected, given more time for the changes to occur (PAM1 to PAM250), the probability for the substitutions to occur rises. However the general relation between the SNPs mostly stays comparable. The biggest change occurs for M1V which is seemed much more likely given only one unit of time.

From this small example the two matrices seem to correlate well. All values that are high in BLOSUM62 are comparably high in both PAM matrices as well, the only exception being M1V which becomes significantly less probable in PAM250. Concerning the phenotype it would be rather far fetched to make an assumption based only on these values, especially since they are mostly close together. However both matrices seem to agree that if there are disease causing SNPs among the set, they are M1V, L451V and E482K, possibly R178H as well.

   Look at the BLOSUM62 and PAM(1/250) matrix. What are the scores for the amino acid substitutions? Is it the worst possible substitution or not? Can we say anything about phenotype from this?

PSSM

PSI-Blast was used to create a PSSM that should contain more refined mutation scores since they are calculated separately for each position. The resulting scores are shown in <xr id="tab:pssm"/>. Comparing the values to the scores from position-independent matrices used above one can see that the general relations between the SNPs mostly stay the same. Exception are M1V, L451V and E482K becoming less likely and R178H becoming more likely.

The change in M1V is to be expected considering that the first position in the homologous sequences should mostly be occupied by methionine due to the start codon. The fact that R178H becomes more likely is surprising since this position is one of those originally identified to be potentially important for catalysis. This might be a small indicator that the original assumption was wrong or that it at least is not important in all homologs.

L451 is found in a single residue loop that could also be seen as an extension of the preceding alpha-helix and is directly followed by a beta-strand. In any case the side-chain at this position is pointing outside of the protein and with only beta-sheets and alpha-helices nearby it does not seem to fulfil a deeper purpose. Therefore exchanging the leucine against the even comparably similar valine should not be a problem and the decreased score is a surprise that cannot be explained for now. Apparently the residue does serve a purpose in other seqeunces. E482K is again an expected change, since there is not only the change in charge but structural analysis also suggests that glutamate forms hydrogen bonds to neighbouring residues.

<figtable id="tab:pssm">

Table TODO: Substitution scores for the SNPs based on a PSSM created with a PSI-Blast query for P06865 on BIG and five iterations. Parameters for the maximum number of returned hits were chosen as in Task 2, i.e. 3800 hits.

Mutation	PSI-Blast PSSM
M1V	-6
L39R	-4
C58Y	-2
L127R	-7
R170W	-5
R178H	3
S210F	1
D258H	1
L451V	-2
E482K	-2

</figtable>

   Getting a bit closer to evolution you will have to create a PSSM (position specific scoring matrix) for your protein sequence using PSI-BLAST (5 iterations). How conserved are the WT residues in your mutant positions? How is the frequency of occurrence (conservation) for the mutant residue type? Anything interesting?

Multiple sequence alignments

   And another step close to evolution: Identify all mammalian homologous sequences. Create a multiple sequence alignment for them with a method of your choice. Using this you can now calculate conservation for WT and mutant residues again. Compare this to the matrix- and PSSM-derived results.

Prediction

SIFT

The score does not seem to take factors like the number of sequences into account, however the webserver warns, that there is low confidence in the prediction for M1V, since there is only one sequence (the query) for this SNP.

The scaled probabilities from the full output suggest that for each position of interest, there are no mutations tolerated at all, and any substitution amino acid should exert an effect. Apparently the positions in question are well conserved in the homologous sequences retrieved by SIFT.

The number of sequences at this position are factored into this value <ref name="siftpub">Ng,P. (2001) Predicting deleterious amino acid substitutions. Genome research, 863-874.</ref>.

<figtable id="tab:sift">

Table TODO: Prediction is the effect prediction, which is always damaging when the score is below 0.05 <ref name="siftdoc">http://sift.jcvi.org/www/SIFT_help.html#SIFT_PROCEDURE</ref>. Score is said score, however the exact calculation remains unclear from both the online documentation as well as the original publication. Sequences are the number of sequences with amino acid information at the current position. Median sequence conservation is 3.00 +- 0.02 in all cases but M1V where it is 4.32.

Mutation	Prediction	Score	Sequences
M1V	AFFECT PROTEIN FUNCTION	0.00	1
L39R	AFFECT PROTEIN FUNCTION	0.00	54
C58Y	AFFECT PROTEIN FUNCTION	0.00	55
L127R	AFFECT PROTEIN FUNCTION	0.00	67
R170W	AFFECT PROTEIN FUNCTION	0.00	66
R178H	AFFECT PROTEIN FUNCTION	0.00	65
S210F	AFFECT PROTEIN FUNCTION	0.00	68
D258H	AFFECT PROTEIN FUNCTION	0.00	68
L451V	AFFECT PROTEIN FUNCTION	0.00	66
E482K	AFFECT PROTEIN FUNCTION	0.00	65

</figtable>

PolyPhen2

The results of the effect prediction by PolyPhen2 are shown in <xr id="tab:pph2"/>. As can be seen all SNPs but M1V are predicted to have a damaging effect with very high confidence. In fact the false positive rates for these probabilities are below 0.1. Obserbing the full output it can be seen that PolyPhen2 used structural information from PDB entry 2gk1 for all mutations but M1V, where no structure is available. In addition UniProtKB sequence annotations were used, marking C58 as being part of a disulfide bridge as well as noting that M1V is part of the signal peptide and L39 and C58 are part of a propeptide. Whether these predictions are therefore more reliable remains to be seen at the end when checking the real effect in HGMD, however only from the scores given by PolyPhen2 this would not be apparent, since many other mutations without these additional features have an equally high probability as well as FDR, FPR, TPR.

<figtable id="tab:pph2">

Table 1: Prediction results for all SNPs by PolyPhen2. Prediction is the predicted effect for the specific SNP. With 'probably damaging' there are less than 10% false positives <ref name="pph2">Peshkin,L. et al. (2010) A method and server for predicting damaging missense mutations. Nature methods, 7, 248-249.</ref>. pph2_class is the output of a binary effect classification and pph2_prob the raw probability of the mutation having a damaging effect. The column based_on present in the output file as well is not shown since it is a remnant of PolyPhen1 and obsolote <ref name="pph2_output>http://genetics.bwh.harvard.edu/pph2/dokuwiki/appendix_a</ref>.

Mutation	Prediction	pph2_class	pph2_prob
M1V	benign	neutral	0.063
L39R	probably damaging	deleterious	1
C58Y	probably damaging	deleterious	1
L127R	probably damaging	deleterious	1
R170W	probably damaging	deleterious	1
R178H	probably damaging	deleterious	1
S210F	probably damaging	deleterious	1
D258H	probably damaging	deleterious	1
L451V	probably damaging	deleterious	0.994
E482K	probably damaging	deleterious	1

</figtable>

SNAP

   Finally, we use three different approaches to score our mutants.
       SIFT
       Polyphen2
       SNAP is installed on the student cluster and should be used command-line only. You will need to create your own ~/.snapfunrc (unless Tim will change the default one) to point to the correct paths. -- As blast is the bottleneck of SNAP, and you are doing that anyway, we might as well look at all possible substitutions in the position of our mutations. This way we can learn much more about the nature of the given mutation: Is our mutation problematic because we introduce an unwanted effect, or because the WT residue is essential and by mutating we remove that?

SNAP was developed by Bromberg et al. in 2007 <ref name="yana">Yana Bromberg and Burkhard Rost. Snap: predict effect of non-synonymous polymorphisms on function. Nucleic Acids Res, 35(11):3823–3835, 2007.</ref>. The method is based on a neural network with many input features such as SWISS-PROT annotations, biochemical properties, sequence information, PSI-BLAST profiles and predicted 1D-structure. As input SNAP needs only sequence information but will also consider and benefit of available functional and structural annotations. SNAP is stated to have a 77% accuracy of non neutral predictions and 76% of the neutral substitutions at 80% accuracy. The reliability index RI indicates the confidence of a certain prediction and should be considered together with the expected accuracy.

<figtable id="tab:pph2">

Table : Prediction results for all SNPs by SNAP. Prediction is the predicted effect for the specific SNP. RI is the reliability index wth a range from 0-9. The expected accuracy ranges from 0 to 100%.

Mutation	Prediction	RI	Accuracy	Non-neutral at this position
M1V	neutral	9	91%	0
L39R	non-neutral	2	63%	10
C58Y	non-neutral	4	71%	16
L127R	non-neutral	2	63%	7
R170W	non-neutral	7	80%	15
R178H	non-neutral	1	60%	16
S210F	neutral	0	51%	1
D258H	non neutral	3	67%	16
L451V	neutral	6	79%	9
E482K	neutral	2	59%	2

</figtable>

Consensus

   Compare ALL results and create an overview table.
   Try to come up with a consensus between all the findings requested above.

Evaluation

   Check whether you are right in the HGMD – were you able to predict a change? 

For this task it is very important to us that you properly interpret and discuss your results. The production of the data should not take that long – so you have more time to do real science!

References

<references/>