看文獻時，翻到一張保守性得分圖：

image

圖片來源：Reducing the structure bias of RNA-Seq reveals a large number of non-annotated non-coding RNA

原文描述：As expected, the most conserved group of NA_RNAs are the tRNAs, tRNA fragments and pre-tRNAs, with most tRNA conserved at least as far as the opossum (Monodelphis domestica). This is to be expected as tRNAs are the most highly conserved ncRNA class across all domains of life

剛好需要計算保守性得分，搜了下，有R包支持。

OK，這就上手。

1、安裝、加載 phastCons100way.UCSC.hg19 包

BiocManager::install("phastCons100way.UCSC.hg19")
library(GenomicRanges)
library(phastCons100way.UCSC.hg19)
ls("package:.UCSC.hg19")
phast <- phastCons100way.UCSC.hg19

注意，這個包是hg19版本的哦，也就是說你輸入的基因組位置對應(yīng)的版本也是hg19。別在這種小細節(jié)上犯錯。

2、計算指定區(qū)域的平均保守性得分

比如計算7號染色體基因組區(qū)域117232380到117232384的平均保守性得分：

gscores(phast, GRanges("chr7:117232380-117232384"))

該命令也可以寫成：
gscores(phast, GRanges(seqnames="chr7", IRanges(start=117232380, width=5)))

結(jié)果是一樣的。

顯示chr7:117232380-117232384區(qū)段平均保守性得分為0.92：

image

3、計算指定區(qū)域的保守性得分

計算7號染色體基因組區(qū)域117232380到117232384的保守性得分：

gscores(phast, GRanges(seqnames="chr7", IRanges(start=117232380:117232384, width=1)))

顯示chr7:117232380-117232384區(qū)段每個堿基的保守性得分：

image

可以看到 (0.8+0.8+1+1+1)/5=0.92，與第2步計算出來的平均保守性得分是一致的。

4、計算多個區(qū)域的平均保守性得分

計算基因組區(qū)域chr7:117232380-117232384、chr2:115262390-115262395、chr3:19597000-19597005的平均保守性得分：

gscores(phast, GRanges(c("chr7:117232380-117232384","chr2:115262390-115262395","chr3:19597000-19597005")))

結(jié)果如下所示：

image

5、結(jié)果怎么看

分值越高，越保守咯~