#基因组干货#之保守性分析及作图

保守序列一般预示其具有潜在的功能，或在细胞发育及调控方面可能发挥重要作用。一般来说编码区序列是高度保守的，尤其是启动子及转录起始位点（TSS）具有极高的保守性。
  但如何研究我们自己的序列的保守性并量化保守性的值呢，一般的方法是通过物种多序列比对，不过我们并不需要自己进行比对，老外已经帮我们做好了，在UCSC里已经展示了使用PhastCons和Phylop (PHAST包)两种方法的多物种序列（100个物种和46个物种）的进化保守性， 多序列比对会产生MAF格式文件，这个比对结果文件还可以用于编码潜能的预测，因为编码区是高度保守嘛，比较有名的工具如phyloCSF，但不建议本地操作MAF文件，因为MAF文件实在是太TM的大了足足有上百G，建议使用galaxy在线分析序列的编码潜能，顺便安利其他几个能预测序列编码潜能的工具：pfamscan、CNCI和CPC，一般是预测完编码潜能后这几个工具取交集作为结果。

  下面开始分析序列保守性并进行实际操作，首先，如果我有条序列及其位置信息，如何观察这条序列的保守性呢？很简单，打开UCSC，左上角点击Genomes,选择基因组版本，在configure界面找到“Comparative Genomics”栏目，把“conservation”track调为” full”就行，如果已经是” full”就不用调了，点击”submit”，然后在基因组浏览器界面输入你的序列的位置信息，比如我想看看FBLIM1基因启动子位置的保守性，我直接输入FBLIM1或” chr1:16,075,390-16,103,219”，得到结果如下图，可以看出FBLIM1的转录起始位点具有高的保守性（黑色凸起的柱子就是）。
  但是，如果我有一堆具有类似特征的序列，比如我预测的新的转录本，我想看看这些新转录本的保守性或这些转录本的启动子的保守性，那该如何表示呢？首先从UCSC获得保守性分值的数据文件，100个物种和46个物种的都可以，我一般看到文献里都是用的PhastCons表示保守性分值，所以下载PhastCons conservation的wig格式的文件而不是PhyloP conservation的wig格式的文件（当然你想用这个也可以），下载页面在这里：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phyloP100way， 然后将其转为bigwig文件，不知道怎么转以及为什么转的请看以前的推送。下载下一个压缩包，解压出来是染色体命名的.wigfix格式的文件（wig格式有可变步长var和固定步长fix格式两种），然后先将每个染色体的.wigfix文件批量转为.bw文件，bigWigMerge工具从http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ 下载。

#获得所有wigFix的文件列表并存入数组
wigFix_a=(`echo *.wigFix`)
#写个循环批量实现
for i in${wigFix_a[@]}
do
  #用染色体名定义要输出的文件名。
  wigFix_out=`echo $i|cut -d . -f 1|awk '{print $1".bw"}' -`
  #将每个.wigFix文件转为.bw文件
  /home/ding/tool/wigToBigWig $i/home/ding/tool/hg19.chrom.sizes $wigFix_out
done
#使用UCSC的bigWigMerge工具合并各个染色体的bw文件为.bedGraph文件，因为我没找到能合并为.bw的工具。。
/home/ding/tool/bigWigMergechr10.bw  chr11.bw  chr12.bw chr13.bw  chr14.bw  chr15.bw chr16.bw  chr17.bw  chr18.bw chr19.bw  chr1.bw  chr20.bw chr21.bw  chr22.bw  chr2.bw chr3.bw  chr4.bw  chr5.bw chr6.bw  chr7.bw  chr8.bw chr9.bw  chrX.bw  chrY.bw phastCons46.bedGraph
#将.bedgraph文件转为.bw文件。
/home/ding/tool/wigToBigWig  phastCons46.bedGraph  /home/ding/tool/hg19.chrom.sizes  phastCons46.bw

  现在我们已经得到了保守性分值的.bw格式文件，接下来就是套路了，用bwtool提取要分析的序列位置的PhastCons分值的均值，然后作图，shell代码及注释如下：

#使用bedtools获得随机序列2000条
randomBed -l 200 -n 2000-seed 2 -g ../hg19.chrom_24.sizes|sortBed -i -  >./enh_statistics/random_all_enh.bed
#使用bwtool获得各个bed文件中心附近2000bp的保守性分值。
bwtool agg 2000:2000  ./enh_statistics/random_all_enh.bed,./enh_statistics/my_sequence.bed,./enh_statistics/uni_cluster.bed,../gencode/protein_gene2.ss,./enh_statistics/other_cluster.bed../phastCons46/phastCons46.bw  /dev/stdout >./enh_statistics/conservation2_result.mean_signal

  使用R作图及代码如下,从图中可以看出，蛋白编码Gene转录起始位点的保守性最高，这是理所当然的，随机区域的保守性最低且没有起伏，其他几个预测的转录本集合的转录起始位点的保守性Seq-p2(黄线)是最高的，而且可以看出蛋白编码基因的phastCons分值达到了0.45左右，所以我个人认为非要给保守性高低一个阈值的化，我觉得应该是0.4，phastCons高于0.4是具有极高的保守性的，高于0.3具有高的保守性，当然还需要进一步与其他功能元件进行比较。

rm(list=ls())
setwd("/home")
library(ggplot2)
library(reshape)
library(Cairo)
mean_signal<-as.data.frame(read.table("./enh_statistics/conservation2_result.mean_signal",header=F,sep="\t",stringsAsFactors= F))
apply(mean_signal,2,as.numeric)
 
 
#画图:
png_path="./figure/conservation_isspe.png"
CairoPNG(png_path, width =6.2, height =6, units='in', dpi=600)
 
ggplot(data=mean_signal_melt,aes(x=loc,y=value,colour=type))+
  geom_line(size=0.6,alpha=1)+
  xlab("Distance to transcriptionstart sites (bp)")+
  ylab("PhastConsScore")+
  xlim(c(-1000,1000))+
  scale_colour_manual(limits=c("seq1","seq2","seq2","gene TSS","random"),labels=c(expression(seq1-p1),expression(seq-p2),expression(seq-p3),"Gene TSS","random"),values =c("red","gold","green","blue","black"))+
  theme_bw()+
  theme(
    axis.text=element_text(size = rel(1.3)),
    axis.title=element_text(size=rel(1.3)),
    # axis.title.x=element_blank(),
    legend.text=element_text(size=rel(1.3)),
    legend.title=element_blank(),
    legend.background=element_blank(),
    legend.key = element_blank(),
    legend.margin=margin(0,0,0,0,"mm"),
    legend.box.spacing =unit(1,"mm"),
    #legend.position=c(1,1),legend.justification=c(1,1),
    legend.position="top"
  )
dev.off()

---

  写在最后，保守性分析完了可以进行motif预测分析，然后motif又可以进行GO等功能性分析（请查看历史消息），都是妥妥的套路，这些套路咱们以后慢慢揭露，尽情关注并接收更多套路文章。