保守序列一般预示其具有潜在的功能，或在细胞发育及调控方面可能发挥重要作用。一般来说编码区序列是高度保守的，尤其是启动子及转录起始位点（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等功能性分析（请查看历史消息），都是妥妥的套路，这些套路咱们以后慢慢揭露，尽情关注并接收更多套路文章。
网友评论