启动子的分析和预测

一、摘要

1. 加深对基因启动子的理解和认知；
2. 学会如何获取已知基因的启动子序列数据；
3. 熟悉不同启动子分析软件的使用及其适用范围；
4. 学会设计启动子克隆引物。
5. 熟悉EPD和TransFac数据库的使用；
6. 学会使用已知的启动子和转录因子TransFac的HMM模型，并能够独立编程，利用该HMM模型来计算鉴别未知启动子

二、材料和方法

1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB

2、系统平台

Windows 8.1、Ubuntu

3、软件平台

【1】Primer-BLAST
【2】Softberry系列工具
【3】Promoter 2.0
【4】BDGP
【5】Cister
【6】NEBcutter

4、数据库资源

NCBI数据库：https://www.ncbi.nlm.nih.gov/
UCSC数据库：http://genome.ucsc.edu/

5、研究对象

人类谷胱甘肽硫转移酶M1的promoter区域

三、结果

基因启动子序列的获取

选择基因：谷胱甘肽硫转移酶M1(GSTM1)
概况：当携带风险基因型时，对环境毒素和致癌物质的敏感性提高，易发生DNA突变和染色体畸变，患白血病的风险因而显著增加。
首先进入UCSC genome browser 查看GSTM1上游5kb范围内有无其他基因。发现该基因的上游存在同一家族的GSTM2，所以promoter大概只有3kb。

图表 1UCSC genome browser
接下来进入Genbank，搜索GSTM1，查看该基因在基因组中的定位和基因结构。

图表 2查看基因定位和结构
打开该基因的序列信息，获取该基因的启动子序列（包含exon1）

Neural Network Promoter Prediction

进入BDGP: Neural Network Promoter Prediction网站http://www.fruitfly.org/seq_tools/promoter.html，进行启动子预测

图表 3 BDGP: Neural Network Promoter Prediction网站
一共预测出来3个启动子(这个网站预测出来的promoter都是50bp)

图表 4 BDGP预测结果

Promoter 2.0 Prediction

使用Promoter 2.0 Prediction Server http://www.cbs.dtu.dk/services/Promoter/
进行启动子预测，也是一共预测出来3个启动子

图表 5Promoter 2.0预测结果

Softberry预测

TSSW、TSSP、TSSG、FPROM都是softberry提供的启动子预测工具，进入
官网(http://www.softberry.com/)，然后点击service即可，启动子预测工具网址：
http://www.softberry.com/berry.phtml?topic=index&group=programs&subgroup=promoter

TSSW

TSSW具体网址如下(http://www.softberry.com/berry.phtml?topic=tssw&group=programs&subgroup=promoter)，输入序列进行预测即可。TSSW并没有预测出来promoter区域。

图表 6TSSW预测结果

TSSP

TSSP具体网址如下(http://www.softberry.com/berry.phtml?topic=tssp&group=programs&subgroup=promoter)，输入序列进行预测即可。共计预测出来一个promoter区域。

图表 7 TSSW预测结果

TSSG

TSSG具体网址如下(http://www.softberry.com/berry.phtml?topic=tssg&group=programs&subgroup=promoter)，输入序列进行预测即可。TSSG并没有预测出来promoter区域。

图表 8TSSG预测结果

FPROM

FPROM具体网址如下(http://www.softberry.com/berry.phtml?topic=fprom&group=programs&subgroup=promoter)，输入序列进行预测即可。FPROM并没有预测出来promoter区域。

图表 9FPROM预测结果

Cister

Transcription Elements预测平台：Cis-element Cluster Finder
https://zlab.bu.edu/~mfrith/cister.shtml
由于序列只有3kb，默认参数预测出来的转录元件太少，将average distance between clusters参数由默认的3w修改为3k，最有可能的结果还是NF-1

图表 10Cister预测结果

Match

转录因子预测集合网站http://gene-regulation.com/pub/programs.html (需要注册)
具体网址http://gene-regulation.com/cgi-bin/pub/programs/match/bin/match.cgi

图表 11Match预测结果

AliBaba 2.1

转录因子预测集合网站http://gene-regulation.com/pub/programs.html (需要注册)
具体网址http://gene-regulation.com/pub/programs/alibaba2/index.html
预测出来一大堆，下面仅截取其中一部分。

图表 12AliBaba预测结果

基因结构绘图

虽然使用了6个promoter预测在线平台，但是只有3个平台预测出promoter。
利用在线平台processon绘制基因结构图

图表 13基因结构图
接下来大致将预测结果分为4个区域，将引物设计的范围同时绘制出来（箭头为引物）

PrimerBlast引物设计

图表 14引物结构
先使用默认参数进行尝试，设置好Forward primer和Reverse primer的Range，再将PCR product size最大值调成整段序列的长度3005，同时# of primers to return参数调整为1，以方便截图。

图表 15引物位置

图表 16初步设置参数
然而由于有重复序列，经过repeat filtering，不会跑出来引物，只能将参数调宽松。

图表 17默认参数结果
在Advanced parameters的Primer Parameters栏目，将Repeat filter关掉，同样可以看出来GC含量和TM值也筛选掉很多引物，在Internal hybridization oligo parameters栏目，将引物内杂交的参数调整宽松。

图表 18高级参数调整

引物设计结果：

图表 19第一段

图表 20第二段

图表 21第三段

图表 22第四段
第一段：

第二段：

第三段：

第四段：

NEBcutter酶切位点分析

使用NEBcutter分析该启动子序列，为了更加全面，查找全部的特异性位点Enzymes to use: All specificities
http://nc2.neb.com/NEBcutter2/
保存没有酶切位点“0 cutters”的核酸内切酶数据，见附录。

图表 23NEBcutter结果

pGL4.17载体

查询pGL4.17的载体数据，获得酶切信息。查询关键词：promega pGL4.17。

图表 24pGL4.17载体信息
其中SfiI、Acc65I、KpnI、SacI、NheI、XhoI、EcoRV、BglII、HindIII这九个酶都属于在promoter内部没有酶切位点的，这些都可以选用。

引物设计

从上面九个酶中随便选两个（真实情况要考虑到切割率等问题）
选择KpnI和SacI，下面是酶切位点和保护碱基对应表，KpnI选两个保护碱基。

可以看出来酶切位点序列在反向互补以后和原序列相同，直接把这段序列加在引物前面就成。
最后按照“保护碱基+酶切序列+PCR引物”的顺序，设计用于可以转到pGL4.17载体的引物。

后续实验流程

接下来，用这四组引物，把四个promoter区域PCR出来，顺带PCR出来的还有保护碱基和酶切序列，导入pGL4.17，用双荧光素酶报告系统看看哪儿个promoter活性最高，大概会出来下面这种图，后面那张图明显说明promoter3活性最高。

再接下来，还可以用TFSEARCH，TFBSs，TRED这样的转录因子预测软件（上面也做了几个预测），看看活性最高的那段区域和哪儿些转录因子相关，或者用pubmed查查看文献，ENCODE，TRANSFAC等数据库，查找这个基因启动子区域的转录因子信息。

接下来是编程练习部分

HMM模型

TransFac是转录因子数据库，但是好像需要注册才能下载模型的矩阵。
从EPD真核生物启动子数据库下载脊椎动物TATA-box的矩阵（共计12位碱基）。
网址http://epd.vital-it.ch/promoter_elements.php
利用该矩阵建立打分模型，对上面谷胱甘肽硫转移酶M1(GSTM1)的启动子序列进行分析，具体代码见附录。
打分值：每次取出12bp序列计算，依次计算每位碱基所占比例，再累乘得到分值（由于数值太小，分值皆除以最大分值）

图表 25打分值统计图
P值计算：使用bootstrap方法，将12bp序列打乱1000次，再按照上述方法计算分值，如果1000次内有n次分值高于“打乱之前的分值”，则p值为n/1000

图表 26 p值统计图
看的出来，整段promoter区域大部分分值都为0，p值为1。之前在线预测软件中的200,400,1200,2570这四个位置，此处也能预测出来，效果还可以。

附录

”0 cutters”核酸内切酶

模型矩阵

从EPD下载，保存为data.txt
http://epd.vital-it.ch/promoter_elements.php

R代码

setwd("G:/AllShare/genomicsHomework/HMMmodel")
hmmmodel<- read.table("data.txt",header = TRUE)
rownames(hmmmodel) <- hmmmodel$Position
hmmmodel <- hmmmodel[,-1]
hmmmodel <- hmmmodel/100
hmmmodel <- t(hmmmodel)
library(seqinr)
promoter<- read.fasta(file = "promoter.fa")
seqmatrix<- as.matrix(promoter$`NC_000001.11:109684796-109687800`)
# 转成大写字母
seqmatrix<- toupper(seqmatrix)
maxseq<- strsplit("TGTATAAAAGGG",split = "")[[1]]
# 计算打分值
computeScore <- function(seq){
  score <- 1
  for(i in 1:length(seq)){
    score <- score*hmmmodel[i,seq[i]]
  }
  return(score)
}
# 使用bootstrap方法,计算p值
bootstrap<- function(seq){
  flag <- 0
  for(i in 1:1000){
    tmp <- sample(seq)
    score<- computeScore(tmp)
    if(score>=computeScore(seq)){
      flag <- flag+1
    }
  }
  return(flag/1000)
}
maxscore<- computeScore(maxseq)
bootstrap(maxseq)
scorevector <- c()
pvaluevector <- c()
for(i in 1:(length(seqmatrix)-11)){
  tmp<- seqmatrix[i:(i+11)]
  score <- computeScore(tmp)
  pvalue <- bootstrap(tmp)
  scorevector <- c(scorevector,score)
  pvaluevector <- c(pvaluevector,pvalue)
}
result <- data.frame(scorevector,pvaluevector)
result$position <- 1:(length(promoter$`NC_000001.11:109684796-109687800`)-11)
colnames(result) <- c("score","pvalue","position")
result$score <- result$score/maxscore
library(ggplot2)
ggplot(result, aes(x=position, y=score)) + 
  geom_line() + 
  geom_point(size=4, shape=20) +
  labs(title="score by HMM model")+
  theme(plot.title = element_text(hjust = 0.5))
ggsave("HMM模型打分图.pdf")

ggplot(result, aes(x=position, y=pvalue)) + 
  geom_line() + 
  geom_point(size=4, shape=20) +
  geom_hline(aes(yintercept = 0.05),colour="red",linetype="dashed")+
  labs(title="p value by HMM model")+
  scale_x_continuous(breaks = c(0,200,400,1030,1200,1632,2570,2700))+
  theme(plot.title = element_text(hjust = 0.5))
ggsave("HMM模型p值图.pdf")

序列

Promoter序列
>NC_000001.11:109684796-109687800 Homo sapiens chromosome 1, GRCh38.p7 Primary Assembly
TCTGCTCTTGAACCCTGTGTTCTGTTGTTTAAGATGTTTATCAAGACAATATGTGCACCACTGAACATAG
ACCCTTATCAGGAGTTCTACTTTTGCCTTTGTCCTGTTTCCTCAGAAGCATGTGATCTTTGTTCTGCTTT
TTGCCCTTTAAAGCATGTGATCTTCGTACCTACCCTCTGTTCGTACACCACCACCCCTTTTGCAATCCTT
AATAAAAACTTGCTGGTTTTGAGGCTCGGGCAGGCATCATGGTCCTACCGATATGTGATGTCACCCCTGG
CGGCCCAGCTGTAAAATTCTTCTCTTTGTACTCTTTCTCTTTATTTCTCAGCTGGCTGACACTTATGGAA
AATAGAAAGAACCTACGTTGAAATATTGGGGGCAGGTTCCCCAATAGCCTTGCTGAGGAAATTAAATTTA
TGTTCAAGTGCTATTTCTTTATGGCACCAAGGAACAAGTATTTCAAACAATACTAATGTAACAGTACTGG
TTCTATGTGTTTCAAAATTATTATTCTCATGAGTGTTAGCTTTCTTAAAAAATCGTTTTTTTTTTCAATT
GGATCTAGACATCTTATCTTTCACAGCTCAAGACGGATTAACTCAGAATCATAAACTCTTAATGCATAAT
GAGAAATATAATGTTTCCTAGGGCCAGGCACTTGTGTCTGTGCTGGTGCTATTGCCTCAATGCAGGAAAA
TCTATGTAAGAGTTCACTGTGAGGCCAAAACTGCTTCCTAAACATGGATACCTGCCAGGTATCTGAGCTG
GGAGTACTGCCCAGGTCTGGATGGGCGGGGAGTGTTTGCAACAAGGACTGTGCCTTGCCAGCCTCAGTGA
CACAGTGTCCAAGTGCCCCAACTTAGCAGCCACCTGCTGACCACCTGATTTCTGTGGCCTAATAGGGATG
TGATGAAGTCTACCTGTTTACTCAACCCCAAACCACACATTATCCAGGTGGTTTGAAACTTTTTTGATAT
ACTGGGTTCATCCTCTGGAGTCCTAACAATGTTTTAGCTAATTTACAAAAAACAAAACAAAACAAAACAA
AACAAAACAAAACAAAAAAACTACTTTTTTTGCAGCACAACAGCCTGGTTTACATTGCAAAATGATTTCT
CATTAAAGGTCTATCATCTATTTCCATATGTCCATTATTATTTGCAATATCCTTTAAAGCAGTCAACCCC
AGGCTAATCCATTGCACAACTCTTTTGAAAGTCTTCCTTCTACCTTGAAAGAAGAAAGTTGGCAGGTTGG
ACATTGTTCTCGTGGAGGTTGTACCATGGGTCACATATCACGGTGTGACTTCAAAGGCCACTGGAGCCAC
CGTCTCATACTGAAGAACACACATGGGTCAGGAGCCAGGTCCAGGTCCGGAATGGTGGATCTGGAGAGGG
GAGGGTCCCTGCCTGTGGTCCTGTGGGGAGCCCTCAGGCTCCTCTCTGGCCACCATCCTCTGACCTCCCT
CCTCAGCAGGACAGGGTTCTGGCTTCTCTGAGGGACAGGTTCTGTGGCAGGCCAGGCGTCACTAACACAG
GCCTCCATAACAACTGTTTCAGTACTGACTGAGTGGTGAAGTTAAATATTAAAAGCTGAAAAAAGCCAGT
ACCTTTATACAGAGGCTGGATGTAACAAAAGCCCACCAAGAGTTTTGCTTAGGCCTTTCCTGGGCCTTAA
AGCATGACAAAACAATGAAGGAATTCTTAACAGGACCTATTTAGAATTAAACAAGTTTTATTGTGAGTCT
GAAGAAACTCCCCAGGCCTCCACAAACAAGTTTATTGGGCGTCTGAAGGAACTCCCCAAACCTCCGTGAT
TTAGCAGGAGACAAGATAAGGGTAATCATCCCCCGCACCTGGACCCATTTAGATTAAATAAATAGACTGA
GGCTCCAGAATAAGGTCCTCAGGACCCAGACCTCAGTTACAGATTAAAGAAGTTAATCACTTATGTCTTT
AGATGAATGCACACTTACTTGTAGACATATACCTTAGAAGGTATATATGCTCTGGAAAACTTTGTAATAT
TGAGTTGGTCTGGTGGTAATTTCTAGGCCTTCTCCCTGTTACCGGTTGCAGAAATAAAACCTCTCTTCCT
CCCCATTTGATCTGCATCTCGTTATTGGGCCTAGAGAAATAGCAGCCGGACCCTCAGTTTGGTCCGGGAA
GTTCTTCCATCCTCCCTCGCCTGCTCTCTGTGGCCACTGCACTCACTGTTGCTGTTGCTGTTCCGGTCTC
TGTGAGGTTCACCTAGTGGACTGGCTGGACATTTCTAGGGGGCACCTCAGATACCTCACCAACTTGCTGG
ATCTGATCCTTGGATTTCGATTCATAAATTGTGCCAAAATACGAAGTGGCTAATTTACACAGTACTTAGC
CAGATGACCGAAGGACTCAGTACCCGAGGGCCCCTAACAGAAAACACAGACCACATTTCCTTTACTCTGG
CCCTTTTCCTGGGGGTCCTTCCTATACCACTGACACTGTTCCTGTGTAGGCGGGGCTAGAGGGGAGACTA
AGCCCTGGGAGTAGCTTTCGGATCAGAGGAAGTCCTGCTCTTACAGTGACAGGGGCTGAATTAAATTCCC
AGGTTGGGGCCACCACTTTTTAGTCTGACCCCTGCAGCCGGAGTCTCCCAGAGCCCTTGGGAACTCGGCA
GCGGAGAGAAGGCTGAGGGACACCGCGGGCAGGGAGGAGAAGGGAGAAGAGCTTTGCTCCGTTAGGATCT
GGCTGGTGTCTCAAGCGCACAGCCAAGTCGCTGTGGACCTAGCAAGGGCTGGATGGACTCGTGGAGCCTC
AGGGCTGGGTAGGGAAGCTGGCGAGGCCGAGCCCCGCCTTGGGCTTCTGGGCGCCCTGACTTCGCTCCCG
GAACCCTCGGGCCTGGGAGGCGGGAGGAAGTCTTACTGAGTGCAGCCCCAGGCGCCCTCTCCCGGGCCTC
CAGAATGGCGCCTTTCGGGTTGTGGCGGGCCGAGGGGCGGGGTCGCAGCAAGGCCCCGCCTGTCC