生信实操丨序列进化保守性分析和可视化

1.PhastCons,Phylop介绍和小工具下载

获得保守性打分的第一步是对多物种的全基因组序列进行多序列比对，然后根据比对结果获得分值。这是一个非常繁琐的过程，好消息是UCSC基因组浏览器使用两个软件PhastCons和PhyloP把由MULTIZ产生的60个脊椎动物的多序列比对结果转换成两种单一的保守性计分和显示，这两个方法都基于已知的种系树结构和利用一个称作phylo-HMM的隐马尔可夫模型种系分析方法。其中，PhyloP能够识别加速进化和保守进化的位点，其分值有正有负；而PhastCons的计分总是一个0至1之间的值，值越大保守性越高。我们通过UCSC小工具bigWigAverageOverBed，结合PhastCons和PhyloP文件来计算保守性打分。

小鼠PhastCons和PhyloP多物种比对文件地址如下：ftp://hgdownload.soe.ucsc.edu/goldenPath/mm10/（其他物种的文件在此文件夹上层目录）

UCSC小工具bigWigAverageOverBed的下载链接如下：http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/。此外，UCSC也提供了多种其他实用的小工具，包括文件格式之间的转换，计算bigwig文件之间的相关性等等。

# 创建存储数据的文件夹

mkdir -p data_db/phastcons_phylop&&cd data_db/phastcons_phylop

# phyloP下载

wget -bc ftp://hgdownload.soe.ucsc.edu/goldenPath/mm10/phyloP60way/mm10.60way.phyloP60way.bw &

# phastCons下载

wget -bc ftp://hgdownload.soe.ucsc.edu/goldenPath/mm10/phastCons60way/mm10.60way.phastCons.bw &

# 小工具下载

wget -c http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigAverageOverBed

# 更改文件属性

chmod 744 bigWigAverageOverBed

# 检查工具是否可以运行

./bigWigAverageOverBed

可以正常运行。它的使用很简单，只是输入文件in.bw、in.bed和输出out.tab。

2.保守性分析

文件准备：

1) mm10.60way.phyloP60way.bw

2) mm10.60way.phastCons.bw

3) in.bed: 输入的bed文件

bed文件格式可参考：

http://www.genome.ucsc.edu/FAQ/FAQformat.html

输入的trans.bed文件如下：

phy=/xx/data_db/phastcons_phylop/mm10.60way.phyloP60way.bw

pha=/xx/data_db/phastcons_phylop/mm10.60way.phastCons.bw

# phyloP

nohup ./bigWigAverageOverBed $phy trans.bed trans.phylop.tab &

# phastCons

nohup ./bigWigAverageOverBed $pha trans.bed trans.phastcons.tab &

输出的结果文件：

1） trans.phylop.tab

2） trans.phastcons.tab

结果文件trans.phylop.tab和trans.phastcons.tab格式相同，这里以trans.phylop.tab文件进行展示。

trans.phylop.tab结果文件

使用awk提取name和mean列，作为下游数据分析的原始文件。需要注意的是，结果文件中的name列对应bed文件的第四列（name列）。

awk ‘{print $1,$6}’ trans.phastcons.tab > pha.res

awk ‘{print $1,$6}’ trans.phylop.tab > phy.res

3.可视化

得到打分结果之后，将数据导入R中进行处理，其中phy和pha数据框都是两列，一列是转录本ID，另一列是保守性打分。本次分析中我们认为PhastCons > 0.3,PhyloP > 0.75的转录本更加保守，这些转录本可能在特定的生物过程中扮演重要角色。

library(tidyverse)

phy <- read_delim("phy.res",delim = " ",col_names = F)

colnames(phy) <- c("trans", "phy_score")

pha <- read_delim("pha.res", delim = " ", col_names = F)

colnames(pha) <- c("trans", "pha_score")

ph <- left_join(phy, pha, by = "trans")

x <- 0.75

y <- 0.3

ph$class <- NA

ph$class[ph$phy_score > x & ph$pha_score > y] <- "phy&&pha"

ph$class[ph$phy_score > x & ph$pha_score <= y] <- "only phy"

ph$class[ph$phy_score <= x & ph$pha_score > y] <- "only pha"

ph$class[ph$phy_score <= x & ph$pha_score <= y] <- "none"

整理数据之后得到的ph数据框如下：

首先使用boxplot简单看一下数据的分布：

ph[,2:3] %>% boxplot(col = c("red", "steelblue"))

接下来，使用散点图和密度图对数据可视化：

ph_plot <- function(data2){

  library(ggplot2)

  library(ggpubr)

  library(ggsci)

  colnames(data2) <- c("trans", 'x', 'y', 'class')

  # 主图

  scatter <- ggplot() + 

    geom_point(data=data2,aes(x= x, y= y, fill= class, color = class),shape=21,size=1.5)+

    scale_fill_aaas(alpha = .6)+

    scale_color_aaas(alpha = .6) +

    geom_vline(xintercept = x,linetype = 2) +

    geom_hline(yintercept = y,linetype = 2) +

    theme_test() +

    labs(x = "phyloP", y = "phastCons") +

    theme(legend.position = c(0.9,0.2))

  # 顶部图

  hist_top <- ggplot()+ geom_density(data=data2,aes(x),colour='black',alpha=0.7,size = 1.2)+

    theme_classic()+

    theme(legend.position="none",

          axis.title.x = element_blank(),

          axis.text.x = element_blank(),

          axis.ticks.x = element_blank())

  # 右部图

hist_right <- ggplot()+  geom_density(data=data2,aes(y),colour='black',alpha=0.7,size = 1.2)+

    theme_classic()+

    coord_flip()+

    theme(legend.position="none",

          axis.title.y = element_blank(),

          axis.text.y = element_blank(),

          axis.ticks.y = element_blank())

  # 空白

  empty <- ggplot() +

    theme(panel.background=element_blank(),

          axis.title.x=element_blank(),

          axis.title.y=element_blank(),

          axis.text.x=element_blank(),

          axis.text.y=element_blank(),

          axis.ticks=element_blank())

  # 合并

  ggarrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, widths=c(4,1), heights=c(1,4),align = "hv")

}

ph_plot(ph)