Protein S, human

The human variation (HumVar) and human divergence (HumDiv) data sets used to assess SIFT’s performance were obtained from UniProtKB (20 (link)). Adzhubei et al. (20 (link)) compiled the HumDiv deleterious list using mutations annotated to cause Mendelian diseases in humans. They created the HumDiv neutral data set by comparing human proteins to their homologs in closely related mammals, and identifying amino acids that are different. For the HumVar deleterious data set, the authors included any mutation annotated to cause human disease, regardless of whether they are Mendelian in origin or not. The HumVar neutral data set is made up of nonsynonymous polymorphisms not annotated as disease causing. We mapped the HumVar and HumDiv data to Ensembl, RefSeq and UCSC Known ids using the UniProtKB id mapping tool (http://www.uniprot.org/help/uniprotkb). Not all mutations from the data sets could be mapped. Hence, the final number of mutations used is less than that of the original dataset (Table 1). True positives (TP) are defined as disease-causing mutations correctly predicted to affect protein function, and false negatives (FN) are those incorrectly predicted to be tolerated. True negatives (TN) are neutral variations correctly predicted as tolerated and false positives (FP) are neutral variations incorrectly predicted to affect protein function.
Table 1.

Number of HumDiv and HumVar data points used to assess SIFT’s performance

Data set	Number of data points			Coverage** (%)
	From original dataset (20 (link))	Used in evaluating SIFT*	With SIFT predictions
HumDiv neutral	6027	5816	5582	96.0
HumDiv deleterious	3055	2893	2791	96.5
HumVar neutral	8638	7475	7178	96.0
HumVar deleterious	12 598	11 982	11 561	96.5

*Lookups to the SIFT database required Ensembl, RefSeq and UCSC Known protein identifiers and the chromosome associated with the given identifier. Not all data points could be mapped to these types of protein identifiers using UniProtKB’s ID mapping tool. Furthermore, we were not able to map some proteins to their chromosomes.

**Coverage = (Number with predictions/Number of data points tested)

The various statistics are computed as follows:

Sensitivity = TP/(TP + FN)

Specificity = TN/(TN + FP)

Accuracy = (TP + TN) / (TP + TN + FP + FN)

Precision = TP / (TP + FP)

Negative predictive value (NPV) = TN / (TN + FN)

Matthews correlation coefficient (MCC) = X / Y

where X = [(TP × TN) – (FP × FN)] and Y = SQRT[(TP + FP) (TP + FN) (TN + FP) (TN + FN)].
We generated receiver operating characteristic (ROC) curves for each protein database by computing the SIFT score for each substitution and categorizing them as tolerated or deleterious using different threshold values. For each threshold, the true positive rate (sensitivity) and false positive rate (1 – specificity) are then computed and plotted in Figure 1.
Figure 1.

Performance statistics of SIFT predictions on PolyPhen-2’s (a) HumVar and (b) HumDiv data sets when using various protein databases. ROC curves on the (c) HumVar and (d) HumDiv data sets. Although UniRef-100 shows slightly better performance than UniRef-90, it has lower coverage.

Coding-potential prediction is essentially a binary decision problem, which makes logistic regression a suitable approach. As an alignment-free method, all selected features (predictor variables) were calculated directly from the sequence. The first feature was the maximum length of the open reading frame (ORF). ORF length is one of the most fundamental features used to distinguish ncRNA from messenger RNA because a long putative ORF is unlikely to be observed by random chance in noncoding sequences. Despite the simplicity, ORF length has high concordance with more sophisticated discrimination methods and remains the primary criterion in almost all coding-potential prediction methods (21 (link)). The second feature was ORF coverage defined as the ratio of ORF to transcript lengths. This feature also has good classification power, and it is highly complementary to, and independent of, the ORF length (Supplementary Figures S1 and S3). Some large bona fide ncRNAs may contain putative long ORFs by random chance (5 (link)), and thus cannot be classified correctly by ORF length alone. Fortunately, those large ncRNAs usually have much lower ORF coverage than protein-coding RNAs (Figure 1B).
Figure 1.

Score distribution between coding (red) and noncoding (blue) transcripts for the four linguistic features selected to build the logistic regression model; training data set containing 10 000 coding and 10 000 noncoding transcripts were used. (A) ORF size. (B) ORF coverage. (C) Fickett score (TESTCODE statistic). (D) Hexamer usage bias measured by log-likelihood ratio.

The third feature we used was the Fickett TESTCODE score (termed ‘Fickett score’ hereafter), which is a simple linguistic feature that distinguishes protein-coding RNA and ncRNA according to the combinational effect of nucleotide composition and codon usage bias (22 (link)). Briefly, the Fickett score is obtained by computing four position values and four composition values (nucleotide content) from the DNA sequence. The position value reflects the degree to which each base is favored in one codon position versus another. For example, position value of A (A_pos) is calculated as follows:

C_pos, G_pos and T_pos are determined in the same way. The percentage composition of each base is also determined. These eight values are then converted into probabilities (p) of coding using a lookup table provided in the original article. Each probability is multiplied by a weight (w) for the respective base, where the value of w reflects the percentage of time each parameter alone successfully predicts coding or noncoding function for the sequences of known function. Finally, the Fickett score is calculated as follows:

The Fickett score is independent of the ORF, and when the test region is ≥200 nt in length (which includes most lncRNA), this feature alone can achieve 94% sensitivity and 97% specificity, with ‘no opinion’ on 18% of the sequences (22 (link)).
The fourth feature we used was hexamer usage bias (termed ‘hexamer score’ hereafter). This may be the most discriminating feature because of the dependence between adjacent amino acids in proteins (23 (link)). The hexamer score can be computed in numerous ways; here, we used a log-likelihood ratio to measure differential hexamer usage between coding and noncoding sequences. For a given DNA sequence, we calculated the probability of the sequence under the model of coding DNA and under the model of noncoding DNA, and then we took the logarithm of the ratio of these probabilities as the score of coding potential. We used F (h_i) (i = 0, 1, … , 4095) and F′ (h_i) (i = 0, 1, … , 4095) to represent in-frame hexamer frequency, calculated from coding and noncoding training data sets (described below), respectively. For a given hexamer sequence S = H₁, H₂, … , H_m,

Hexamer score determines the relative degree of hexamer usage bias in a particular sequence. Positive values indicate a coding sequence, whereas negative values indicate a noncoding sequence.
We build a logistic regression model using these four linguistic features as predictor variables. A χ² test was used to evaluate whether our logit model with predictors fits the training data significantly better than the null model, which had only an intercept. We built a high-confidence training data set to measure the prediction performance of our logit model. This data set contained 10 000 protein-coding transcripts selected from the RefSeq database; all transcripts had high-quality protein sequences annotated by the Consensus Coding Sequence project. We also added 10 000 randomly selected noncoding transcripts from the GenCODE database. We evaluate the model with a 10-fold cross-validation and measured its sensitivity, specificity, accuracy, precision and area under the curve (AUC) characteristics. The receiver operating characteristic (ROC) curve and precision–recall (PR) curve were generated using ROCR package (24 ). We also built a nonparametric two-graph ROC curve for selecting the optimal CPAT score threshold that maximizes the sensitivity and specificity of CPAT while minimizing misclassifications.
We built an independent test data set to compare the performance of CPAT with that of CPC, PhyloCSF and PORTRAIT. This test set composed of 4000 high-quality protein-coding genes (RefSeq annotated) and 4000 lncRNAs from a human lncRNA catalog (5 (link)). None of these 8000 genes was included in the training data set for CPAT. Assuming that all 4000 lncRNAs are truly noncoding sequences, we could compute the sensitivity, specificity, accuracy and precision of the algorithms to measure their performance. PhyloCSF could not determine the coding status of 528 (13.2%) noncoding genes. Those 528 genes were equally assigned to the true-negative and false-positive categories. The abbreviations in the equations below are as follows: FN, false negative; FP, false positive; TN, true negative; TP, true positive

We uniformly collected and preprocessed a set of Genome-Wide Association Studies(GWAS) summary statistics that (1) came from European ancestry; (2) SNV heritability h²>0.01; (3) z score of h²>4; (4) sample size >10 000. We applied linkage disequilibrium score regression (LDSC) to analyse whether the heritability of these traits enriched in common SNVs around RBP-gene (window size=100 kb). We used 1000 Genome21 (link) European population as a reference panel, only the SNVs within the HapMap322 (link) project were included, and the baseline model and other parameters of LDSC were set at default. To control the bias of incorrect SNV-to-gene mapping, we additionally applied the abstract mediation model (AMM),23 (link) an extension of LDSC that also considered the k-nearest genes of each SNV. We used the default hyperparameters provided by AMM, which were estimated by the benchmark gene set of loss-of-function intolerant genes. We directly transformed the enrichment z score of AMM into p value under a normal distribution and used it for FDR adjustment, without log-transformation of the enrichment.
To analyse whether the human-specific directional impact of common SNVs on RBP profile has a phenotypic consequence, we first calculated a ‘humanisation score’ (HS) of each 1000 genome common SNVs. As described above, we first used 13 520 465 overlapping blocks (1000 bp length each) b=1, 2, …13 520 465 to cover the full length of all protein-coding genes. For each block b, we used sequence-mode Seqweaver to calculate a vector of RBP for the hg38 sequence on ${\displaystyle b\left(hg{38}_b\left(r\right),r=1,2,\dots 217\right)}$ , a vector of RBP for common primate ancestor sequence on ${\displaystyle b\left(an{c}_b\left(r\right)\right)}$ and calculated the difference of them ${\displaystyle \mathrm{\Delta }RB{P}_b\left(r\right)=an{c}_b\left(r\right)-hg{38}_b\left(r\right)}$ . We then mapped each 1000 G common variation s to a block ${\displaystyle b=f\left(s\right)}$ , whose midpoint was closest to s (note that multiple SNVs could be mapped to the same block). If s was not on a gene body, $f\left(s\right)=0.$ We generated a DNA sequence corresponded to $s$ by SAMtools24 (link) consensus option, where the reference FASTA was the sequence of block ${\displaystyle b=f\left(s\right)}$ . We applied Seqweaver on this DNA sequence to obtain vector ${\displaystyle RB{P}_{f\left(s\right)}\left(r\right)}$ , and calculated the HS of SNV $s$ as the inner product of two vectors:
${\displaystyle H{S}_s=\left(hg{38}_{f\left(s\right)}\left(r\right)-RB{P}_{f\left(s\right)}\left(r\right)\right)\times \mathrm{\Delta }RB{P}_{f\left(s\right)}\left(r\right)}$
A negative value of ${\displaystyle H{S}_s}$ indicated that SNV $s$ modified the human RBP profile in the opposite direction to that in which all HLMs on block $f\left(s\right)$ collectively modified the ancestor RBP profile, and it made the human profile closer to the ancestor profile, termed ‘de-humanisation’. Likewise, positive ${HS}_s$ indicated that $s$ modified human RBP profile further from the ancestor RBP profile, corresponding to ‘over-humanisation’. We assumed that, if human evolution on RBP profile is polygenic, then there will be a large number of blocks $b$ whose ${\Delta RBP}_b\left(r\right)$ had a small but non-zero phenotypic association. Then, SNVs on these blocks that impact $\Delta RBP$ would also slightly impact phenotypes. Thus, ${HS}_s$ would be associated with SNV-based trait heritability enrichment. We used two approaches to evaluate this association: