从序列查找下载（NCBI）处理到多序列比对----记录一次不太成功的尝试

Posted 2020-11-25 s-qw

一、问题呈现

找到Streptomyces属里hrdb基因的启动子（hrdbp）的保守序列，希望以此推断出-10区和-35区。

二、过程

1、下载15-20条hrdb基因的启动子序列，并处理形成一个fasta文件

1.1、以coelicolor A3(2)的hrdb基因为源头，通过blast找到得分最高的前50条序列。Download下载Hit Table(txt)格式的文件，这个文件表头会告诉你每一列显示的是什么。

接下来用excel打开这个文件，首先把 alignment length小于1500bp的全部删掉，找到subject acc.ver、s.start和s.end这三列，等一下要用这三列来生成url。

url示例:https://www.ncbi.nlm.nih.gov/nuccore/LT629768.1?report=fasta&from=6444177&to=6445864，生成url的代码如下：

 1 #读入数据
 2 fo = open(‘D:\temporary\hrdb_related\ZECB16FT01N-Alignment.csv‘,‘r‘)
 3 ls=[]
 4 for line in fo:
 5     line = line.replace(‘
‘,‘‘)
 6     ls.append(line.split(‘,‘))
 7 fo.close()
 8 #写出成url
 9 fo1 = open(‘D:\temporary\hrdb_related\output.csv‘,‘w‘)
10 for line in ls:
11     if line[-1] ==‘0‘ :    #0表示为正向序列，可以用s.start-400和s.stsrt+5来做起始和结束的位置
12         fo1.write(‘https://www.ncbi.nlm.nih.gov/nuccore/‘+line[1]+13               ‘?report=fasta&from=‘+line[9]+‘&to=‘+line[11]+‘
‘)
14 fo1.close()

1.2、本来想直接用url爬虫获得相应的二十条序列，但发现相应序列是加密过的，在网上搜了一下，好像跟异步加载有关。就直接手动打开每个url，下载相应的序列，这样我就得到了20个fasta文件，每个文件含有一条序列。

1.3、接下来我想把这20个文件合并成一个fasta文件，用于之后的多序列比对。听说在linux下一个cat命令就可以解决，所以我就想把这些个文件发送到我的linux虚拟机，试了WinSCP3，但没有连接成功，NAT模式可以联网，但换到桥接模式就连不上网络，所以WinSCP3也没有成功连接上我的linux虚拟机。后来又发现WIN10的DOS下也可以通过copy的命令实现同样的功能。只是需要把所有要合并的文件名都用加号连起来，略有些麻烦。也是用代码实现。

1 import os
2 for (dirname,subdir,subfile) in os.walk(r‘D:	emporaryhrdb_related‘) :
3     for f in subfile:
4         print(f+‘+‘,end=‘‘)

到这里就得到可以用来做多序列比对的fasta文件了。

2、meme分析和mega分析

在mega里多序列比对则发现起始密码子上游200bp都非常保守；用meme分析出来3个保守区，与clustax分析结果相比漏掉了上游100-150bp，但这一段实际上也非常保守，看来meme有一些局限性。

到这里这次不太成功的尝试就结束了，只得到起始密码子上游200bp非常保守的结论，而这个也早就被报道了。并没有得到-10区和-35区的相应序列。

3、在以上分析之外，我也试了直接预测启动子的在线软件。可能由于Streptomyces的GC含量过高，用专门的细菌启动子预测软件预测不出来或者极其不准，也没有找到合适的在线软件。