2 Direct comparison of phenotypes
2.1 Phenotypes associated to the gene candidate
Let’s use our prior knowledge to find disorders, and eventually human phenotypes, associated to the gene candidate CC2D2A:
genDis <- traitDef[
match(
geneByTrait[which(geneByTrait$entrez==genId), "id"],
traitDef$id
),
]
genDis## db id name
## 10479 OMIM 249000 Meckel syndrome type 1
## 10484 OMIM 612284 Meckel syndrome type 6
## 8809 OMIM 216360 COACH syndromegenHpDef <- hpDef[
match(
geneByHp[which(geneByHp$entrez==genId), "hp"],
hpDef$id
),
]
genHp <- genHpDef$id
dim(genHpDef)## [1] 119 2head(genHpDef)## id name
## 1765 HP:0002342 Intellectual disability, moderate
## 1974 HP:0002650 Scoliosis
## 136 HP:0000154 Wide mouth
## 1696 HP:0002240 Hepatomegaly
## 2157 HP:0002896 Neoplasm of the liver
## 1716 HP:0002269 Abnormality of neuronal migrationHow many of these 119 HP terms are shared with the 8 HP related to the syndrome under focus?
genHpDef[which(genHpDef$id %in% hpOfInterest),]## id name
## 506 HP:0000639 Nystagmus
## 1817 HP:0002419 Molar tooth sign on MRI
## 972 HP:0001250 SeizuresHowever some of the 116 CC2D2A associated HPs which are not among the 8 HPs of interest present strong similarity with the last ones. For example the Dandy-Walker malformation (HP:0001305) phenotype, associated to CC2D2A, is an indirect but closely related descendant of Ventriculomegaly (HP:0002119) phenotype of interest.
The following steps describe a way to measure similarity between different HP terms.
2.2 Information content and semantic similariy
Here we measure semantic similarity between HPs using gene information content (IC). The formula below shows how information content is computed for each HP term \(p\):
\[ IC_{p} = -ln\left(\frac{|p|}{|root|}\right) \]
Where \(|p|\) is the number of gene associated to the HP term \(p\) and all its descendants. \(root\), in our case, is Phenotypic abnormality (HP:0000118). By definition: \(IC_{root} = -ln\left(\frac{|root|}{|root|}\right) = 0\).
Let’s use the computeHpIC function to compute IC for all HP terms descendants of Phenotypic abnormality in the human phenotype ontology. This function needs to know the genes associated to each HP and the descendants of each HP term.
info <- unstack(geneByHp, entrez~hp)
ic <- computeHpIC(info, hp_descendants)Let’s have a look at the distribution of IC:
IC is an measure of the specificity of genes associated to HPs. The higher IC, the more specific.
Semantic similarity (\(SS_{p_{1}p_{2}}\)) between two HP terms is then defined as the IC of the most informative common ancestor (MICA) (i.e. showing the higher IC).
Let’s use the clacHpSim function to compute the semantic similarity between different couples of HP terms:
hp1 <- "HP:0000518"
hp2 <- "HP:0030084"
hpDef[which(hpDef$id %in% c(hp1, hp2)), 1:2]## id name
## 401 HP:0000518 Cataract
## 9645 HP:0030084 ClinodactylycalcHpSim(hp1, hp2, IC=ic, ancestors=hp_ancestors)## [1] 0hp1 <- "HP:0002119"
hp2 <- "HP:0001305"
hpDef[which(hpDef$id %in% c(hp1, hp2)), 1:2]## id name
## 1012 HP:0001305 Dandy-Walker malformation
## 1601 HP:0002119 VentriculomegalycalcHpSim(hp1, hp2, IC=ic, ancestors=hp_ancestors)## [1] 2.8500152.3 Comparing two sets of phenotype
Now, we can compute semantic similarity between all HP of interest and CC2D2A associated HPs using the compMat function:
compMat <- compareHPSets(
hpSet1=genHp, hpSet2=hpOfInterest,
IC=ic,
ancestors=hp_ancestors,
method="Resnik",
BPPARAM=SerialParam()
)
dim(compMat)## [1] 119 8head(compMat)## HP:0000483 HP:0000510 HP:0000518 HP:0000639 HP:0001249
## HP:0002342 0 0 0 0 1.4502253
## HP:0002650 0 0 0 0 0.0000000
## HP:0000154 0 0 0 0 0.0000000
## HP:0002240 0 0 0 0 0.0000000
## HP:0002896 0 0 0 0 0.0000000
## HP:0002269 0 0 0 0 0.4262711
## HP:0001250 HP:0002119 HP:0002419
## HP:0002342 0.5169219 0.4262711 0.4262711
## HP:0002650 0.0000000 0.0000000 0.0000000
## HP:0000154 0.0000000 0.0000000 0.0000000
## HP:0002240 0.0000000 0.0000000 0.0000000
## HP:0002896 0.0000000 0.0000000 0.0000000
## HP:0002269 0.4262711 0.9002305 0.9002305Then we compute the symmetric semantic similarity score of the matrix to get single value corresponding to similarity between the two sets of HP terms: the HP terms of interest and CC2D2A associated HPs.
hpSetCompSummary(compMat, method="bma", direction="symSim")## [1] 1.357548Unfortunately it is not easy to interpret such a score and to assess it’s significance. To do it we need to compare the score of the candidate gene (CC2D2A) to the score of all the other genes for which we can compute it. Let’s compute the score for all the genes:
## Compute semantic similarity between HP of interest and all HP terms
## This step is time consumming and can be parallelized.
## Use the BPPARAM parameter to specify your own
## back-end with appropriate number of workers.
hpGeneResnik <- compareHPSets(
hpSet1=names(ic), hpSet2=hpOfInterest,
IC=ic,
ancestors=hp_ancestors,
method="Resnik",
BPPARAM=SerialParam()
)
## Group the results by gene
hpByGene <- unstack(geneByHp, hp~entrez)
hpMatByGene <- lapply(
hpByGene,
function(x){
hpGeneResnik[x, , drop=FALSE]
}
)
## Compute the corresponding scores
resnSss <- unlist(lapply(
hpMatByGene,
hpSetCompSummary,
method="bma", direction="symSim"
))
## Get the score of the gene candidate
candScore <- resnSss[genId]
candScore## 57545
## 1.357548And now, we can compare the score of the candidate to all the others:
candRank <- sum(resnSss >= candScore)
candRank## [1] 130candRelRank <- candRank/length(resnSss)
candRelRank## [1] 0.04086765
According to a direct comparison, the candidate gene CC2D2A is in the top 4.1% genes the most relevant for the set of HPs of interest. This result can be used for candidate prioritization.
