保守序列一般預(yù)示其具有潛在的功能，或在細(xì)胞發(fā)育及調(diào)控方面可能發(fā)揮重要作用。一般來說編碼區(qū)序列是高度保守的，尤其是啟動子及轉(zhuǎn)錄起始位點（TSS）具有極高的保守性。
??但如何研究我們自己的序列的保守性并量化保守性的值呢，一般的方法是通過物種多序列比對，不過我們并不需要自己進(jìn)行比對，老外已經(jīng)幫我們做好了，在UCSC里已經(jīng)展示了使用PhastCons和Phylop (PHAST包)兩種方法的多物種序列（100個物種和46個物種）的進(jìn)化保守性， 多序列比對會產(chǎn)生MAF格式文件，這個比對結(jié)果文件還可以用于編碼潛能的預(yù)測，因為編碼區(qū)是高度保守嘛，比較有名的工具如phyloCSF，但不建議本地操作MAF文件，因為MAF文件實在是太TM的大了足足有上百G，建議使用galaxy在線分析序列的編碼潛能，順便安利其他幾個能預(yù)測序列編碼潛能的工具：pfamscan、CNCI和CPC，一般是預(yù)測完編碼潛能后這幾個工具取交集作為結(jié)果。
??下面開始分析序列保守性并進(jìn)行
實際操作，首先，如果我有條序列及其位置信息，如何觀察這條序列的保守性呢？很簡單，打開
UCSC，左上角點擊
Genomes,選擇基因組版本，在configure界面找到
“Comparative Genomics”欄目，把
“conservation”track調(diào)為
” full”就行，如果已經(jīng)是” full”就不用調(diào)了，點擊”submit”，然后在基因組瀏覽器界面輸入你的序列的位置信息，比如我想看看FBLIM1基因啟動子位置的保守性，我直接輸入FBLIM1或” chr1:16,075,390-16,103,219”，得到結(jié)果如下圖，可以看出FBLIM1的轉(zhuǎn)錄起始位點具有高的保守性（黑色凸起的柱子就是）。
??但是，如果我有一堆具有
類似特征的序列，比如我預(yù)測的新的轉(zhuǎn)錄本，我想看看這些新轉(zhuǎn)錄本的保守性或這些轉(zhuǎn)錄本的啟動子的保守性，那該如何表示呢？首先從UCSC獲得保守性分值的數(shù)據(jù)文件，100個物種和46個物種的都可以，我一般看到文獻(xiàn)里都是用的
PhastCons表示保守性分值，所以下載
PhastCons conservation的
wig格式的文件而不是
PhyloP conservation的wig格式的文件（當(dāng)然你想用這個也可以），下載頁面在這里：
http://hgdownload.cse./goldenPath/hg19/phyloP100way， 然后將其轉(zhuǎn)為bigwig文件，不知道怎么轉(zhuǎn)以及為什么轉(zhuǎn)的請看以前的推送。下載下一個壓縮包，解壓出來是染色體命名的
.wigfix格式的文件（
wig格式有可變步長var和固定步長fix格式兩種），然后先將每個染色體的.wigfix文件批量轉(zhuǎn)為.bw文件，bigWigMerge工具從
http://hgdownload.cse./admin/exe/linux.x86_64/ 下載。
wigFix_a=(`echo *.wigFix`)
for i in${wigFix_a[@]}
do
wigFix_out=`echo $i|cut -d . -f 1|awk '{print $1".bw"}' -`
/home/ding/tool/wigToBigWig $i/home/ding/tool/hg19.chrom.sizes $wigFix_out
done
/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
/home/ding/tool/wigToBigWig phastCons46.bedGraph /home/ding/tool/hg19.chrom.sizes phastCons46.bw
??現(xiàn)在我們已經(jīng)得到了保守性分值的.bw格式文件，接下來就是套路了，用bwtool提取要分析的序列位置的PhastCons分值的均值，然后作圖，shell代碼及注釋如下：
randomBed -l 200 -n 2000-seed 2 -g ../hg19.chrom_24.sizes|sortBed -i - >./enh_statistics/random_all_enh.bed
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轉(zhuǎn)錄起始位點的保守性最高，這是理所當(dāng)然的，隨機區(qū)域的保守性最低且沒有起伏，其他幾個預(yù)測的轉(zhuǎn)錄本集合的轉(zhuǎn)錄起始位點的保守性Seq-p2(黃線)是最高的，而且可以看出蛋白編碼基因的phastCons分值達(dá)到了0.45左右，所以我個人認(rèn)為非要給保守性高低一個閾值的化，我覺得應(yīng)該是0.4，phastCons高于0.4是具有極高的保守性的，高于0.3具有高的保守性，當(dāng)然還需要進(jìn)一步與其他功能元件進(jìn)行比較。
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()
??寫在最后，保守性分析完了可以進(jìn)行motif預(yù)測分析，然后motif又可以進(jìn)行GO等功能性分析（請查看歷史消息），都是妥妥的套路，這些套路咱們以后慢慢揭露，盡情關(guān)注并接收更多套路文章。
