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

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

引言

我们在做基因组相关的研究时,经常需要通过高通量筛选的方法得到所需的蛋白编码基因或长非编码RNA(LncRNA)。我们在挑选基因或者lncRNA时,可以从进化保守性角度考虑。一般来说,保守性强的基因或lncRNA可能具有更加重要的功能。本篇文章中,以小鼠为例,对已经通过CPAT预测得到的非编码RNA(ncRNA),筛选长度大于200bp,外显子数目大于2的ncRNA作为lncRNA进行保守性分析,进而筛选可能会在小鼠特定的发育过程中发挥重要作用的lncRNA。

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)

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

注:此推文未经许可禁止转载!
阅读推荐:
  • 生信实操|一个生信菜鸟的上道经验分享-转录组测序(富集分析篇)

  • 生信实操|一个生信素人的上道经验分享-转录组测序(可变剪接篇)

  • 生信实操|一个生信素人的上道经验分享-转录组测序(绘图篇)

  • 生信实操|一个生信素人的上道经验分享-转录组测序(新转录本鉴定篇)

  • 生信实操丨一个生信菜鸟的上道经验分享-转录组测序(富集分析绘图篇)

生信实操丨序列进化保守性分析和可视化
生物医学科研方法

刚刚!2021年度国家自然科学基金项目指南公布!

2021-1-16 23:14:59

生物医学科研方法

分享|如何“复制”顶刊文章漂亮的配色?

2021-1-17 0:04:39