python blast在线比对、每天可达数万条

背景：在碱基序列，如

TRINITY_DN931509_c0_g1_i1
GCACGTTCTTCTCGATGCTGACGGAGCGAGGCGCGCTGGCGGAGCTGATGGC
GGCGCAGCTGGAGCGAGGGCGCACGCCGCGGGTGCGCCGGCGCCGGCCGC
GCCCCAGGCCCAGGGCCAGGCCGCGCCCGAGGCCGAGGAGTAGGAAATAAT
GATCTGTGACGAAACTAACCCACCCCTCGAGCTTGCGGAAATGCTGGAATCTG
CGACGTAGGCACTACTCGAACTTGAG

需要比对（即Basic Local Alignment Search Tool，BLAST）时，采取的方式主要有如下 ：

1. 美国国立生物技术信息中心（National Center for Biotechnology Information ，NCBI）官网提供免费在线blast服务
2. 搭建本地的blast环境。
3. 使用云服务商提供的付费BLAST服务，如blast2go

当需比对碱基序列量少时，即不考虑时间的情况下，上述三种方式均可；但如果需要大量比对时，常规操作下三种方式各有利弊，第一种和第三种的共同点——慢；第二种的比对极大消耗本地服务器集群的算力。本文试图改良第一种方式提高效率

一.啰嗦一下在线BLAST的流程

1.首先进入比对网址：https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastx&PAGE_TYPE=BlastSearch&BLAST_SPEC=&LINK_LOC=blasttab&LAST_PAGE=blastn
在线blast序列可采用两种上传方式，直接粘贴到输入框，或上传本地序列文件，如下图。

2.然后根据需要配置其他参数，点击下方的”BLAST” 按钮进入漫长的等待（主要时间为等待NCBI服务器处理提交的序列）如下图

这个页面会自动刷新，当比对完成后刷新出的截面就是结果页面；请注意图中红框标出的”Request ID” 的值，后面要用到。
3.结果页面如下图

在多选框选定某特定碱基序列的选项，页面会自动跳转到该结果的页面(通过该页面定义为GetResults()的JavaScript函数)，然后再根据需要再此页面索取信息。有时BLAST会匹配不到结果（这是序列自身与数据库中的碱基对没匹配上，请跟后面会提到的比对失败区别开），上图中，选项较浅且以 * 开头的为没有匹配到结果的序列，即这些序列的匹配结果为NULL或空。

二、分析

1.F12碱基序列提交后的等待页面发现页面在不断的发送一个POST包到指定地址，以询问是否出结果，如果出结果则跳转到结果页面（恕小弟才疏学浅没有发现POST包的构成规律）。
2.发现结果页面有GET方法可以访问https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=***********；其中***********为上文提到的Request ID；也即是说，提交碱基序列过程和信息抓取过程可以分开，浏览器不必一直等待结果。可以提交序列，然后得到Request ID，另一个程序或设备用Request ID来抓取数据。
3.访问到结果页面后，通过不断选择下拉框的方式来获得同一批提交的不同的碱基序列的结果，而这些操作，使用python+selenium是方便的。
4.不断向NCBI服务器发起比对请求，并使用多线程的方式处理不同比对RID的结果，提高程序效率。

三、实际操作

1. 环境：

   windows 10+python 3.6.5（推荐搭配PyCharm）+phantomjs +Mysql

2. 用到的python第三方扩展包：

   1. selenium
   2. pymysql
   3. openpyxl

3. 预处理：

   1. Mysql中新建table2表存放碱基序列，结构为

n	id	dna

其中n为自增的主键，id为碱基序列的名称，dna为碱基序列。
2.Mysql中新建submit表存放Request ID，结构为

RID	TTime	len

其中TTime为碱基序列提交的时间，len为本次提交的总长度。

将待处理碱基序列导入表table2中。可以参考下列python代码

import pymysql
def sp(st):#分离序列的ID和序列b=st.split('\n')ID=b[0]DNAs=''for i in b:if b.index(i)==0:continueDNAs=DNAs+ireturn ID, DNAs
def DOSQLin(n,idL,DNAs):
#生成插入SQL语句id=idL.split(' ')ID=id[0]sql = "insert into table2 (n,id,dna) values (%s,'%s','%s');"%(n,ID,DNAs)return sql
f = open('C:\\*.fasta','r') #
All=f.read()
a=All.split('>')#每一条序列都是>开头的
del a[0]#删除第一个空元素
#链接数据库
conn = pymysql.connect(host='127.0.0.1', port=3306, user='root', passwd='luanchen',db='blast')
cur = conn.cursor()
n=1
for i in a:ID,DNA=sp(i)sql=DOSQLin(n, ID, DNA)try:cur.execute(sql)except:prio
cur.close()
conn.cmmit()#提交改变结果
conn.close()

四、提交碱基序列得到Request ID

过程分为下列步骤：

1. 从mysql表中得到若干条待比对的碱基序列，再将这若干条从mysql表中删除；
2. 将1中得到的序列生成一个本地文件(*.fasta) ；
3. 通过selenium的方式打开在线BLAST网页，上传2得到的文件；
4. 得到上传文件后的Request ID，并保存；
5. 重复1~4直至mysql table2表中没有序列。

python代码（尽量加了很多注释，实际写得很乱…）

import string
from time import sleep
from selenium import webdriver
from selenium.webdriver.common.by import By
from threading import Thread
import random
import pymysql
def Fasta(lid,n=50):#读取table2表中若干条数据生成文件file = 'D:\\pills\\%s.fasta' % lid#应该事先创建该文件夹，存放随机生成名字的序列文件，下一步再将其上传到NCBI。conn = pymysql.connect(host='127.0.0.1', port=3306, user='root', passwd='luanchen',db='hp')cur = conn.cursor()cur.execute('LOCK tables hp.table2 write')#给table2上表锁，防止多个线程同时读取时，删除序列时发生错误。sql='select n,id,dna from table2 order by n limit %s'% ncur.execute(sql)lf=0di=0f = open(file, 'w') #for i in cur:lf=lf+len(i[2])if lf>30000:#控制提交的总长度，最好不超过3万，否则容易被服务器拒绝print('..')breakdi=di+1f.write('>'+str(i[0])+i[1]+'\n')f.write(i[2]+'\n')cur.execute('delete from table2 order by n limit %s'% di)#取table2的若干条序列后将其在table2中删除conn.commit()cur.execute('Unlock tables;')#解开表锁，使得其他线程可访问该数据表if di==0:return 0cur.close()conn.close()return file
def isExist(driver,id):#判断页面某元素是否存在try:driver.find_element_by_id(id)return Trueexcept:return False
def submit(RIDs,Len=0):f = open('C:\\ursb.txt', 'a+')#备份RID到本地，有可能会有用f.write(RIDs+'\n')f.close()RID = RIDs.split(' ')sql = "insert into blast.submit (RID,len,TTime) values ('%s',%s,now());" % (RID[-1], Len)conn = pymysql.connect(host='127.0.0.1', port=3306, user='root', passwd='luanchen', db='hp')cur = conn.cursor()cur.execute(sql)conn.commit()cur.close()conn.close()
def do(path):driver = webdriver.PhantomJS(executable_path="D:/phantomjs-2.1.1-windows/bin/phantomjs")#无头浏览器phantomjs的位置应该指定，或设置环境变量#driver = webdriver.Firefox()#测试时，应使用可见的浏览器，方便发现问题并解决。try:driver.set_page_load_timeout(60)driver.get('https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastx&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome')except:driver.quit()return do(path)try:upload = driver.find_element_by_id('upl')#sleep(3)upload.send_keys(path)  # send_keys 上传本地文件#print(upload.get_attribute('value'))  # check valuedriver.find_element_by_id('b1').click() #模拟点击BLAST 按钮 提交文件except:driver.quit()return do(path)r = ''for i in range(10):if isExist(driver,'statInfo'):#statInfo 页面元素内部有RIDtry:table_loc = (By.ID, 'statInfo')table_tr_list = driver.find_element(*table_loc).find_elements(By.TAG_NAME, "tr")#通过元素定位方式找到RIDfor tr in table_tr_list:LEN = len(open(path, 'r').read())print(tr.text,LEN)submit(tr.text, Len=LEN)driver.quit()breakbreakexcept:sleep(1)else:continuereturn
def go(lid):#随机生产文件名,go()执行一次提交。#lid并未用到，懒得改了= =suqs=[]Str=''.join(random.sample(string.ascii_letters + string.digits, 4))file=Fasta(Str)do(file)return
def main(n):#注意线程数量和sleep的关系，按照实际情况调节。for i in range(n):try:go(1)sleep(random.randint(6, 15))#随机等待6-15秒except:continue
Thread(target=main,args=(1000,)).start()#这一行代码创建了1个线程，线程不要创建太多，否则会被NCBI拒绝。
sleep(5)
main(1000)

提一下：应该注意自己的Mysql连接方式，并且在代码运行前，在MYSQL上试一下看相关的SQL语句能否运行，比如，第一次执行SQL的批量删除语句时会报错

You are using safe update e and you tried to update a table…

这时应该在MYSQL上执行SET SQL_SAFE_UPDATES = 0;而这种类型的错误在python的错误信息是不体现的。
运行结果：

五、通过RID得到比对结果

步骤如下：
1. 通过submit表中的RID构建GET地址（上文二、分析中已提到），并在submit表中删除该条RID。
2. 通过firefox进入1中构建的网址
3. 使用s = Select(driver.find_element_by_id('queryList'))
s.select_by_index(n)方式改变下拉框的选定，以期跳转到不同碱基序列的页面。
4. 通过元素定位的方式，得到结果页面需要的信息。
5. 重复3~4直至处理完当次RID的所有不同碱基序列的结果。
6. 重复1~5直至submit为空。

在实际操作的情况下，主要有两种比对失败的情况：

1. 一直不出结果（服务器原因）：通过GET访问结果页面却出现等待结果界面，因为五、通过RID获取结果的操作是具有人工（python程序）设置的查询长度限制和20分钟延后性，所以都应该得出结果了。这应该认为是比对失败，跳过该条RID。事实上，当20分钟没有出现结果，数小时后也不会出结果。
2. 发起的查询请求碱基序列过长，被服务器拒绝，如图：
   
   这种情况视作（通过检查相关元素是否存在来判断）比对失败，跳过该条RID。

除了比对失败的情况，还有一种页面不正常的情况，即选定下拉框选项后页面不跳转到选项对应的结果，这个时候如果抓取信息的话会得到错误的信息（误以为跳转前的信息对应跳转后的信息）；通过下图中两个红框内的值是否等同，来判断是否跳转。
如果没跳转用driver.execute_script('GetReasult()')执行页面自带的跳转js代码。

python代码

from selenium import webdriver
from time import sleep
from selenium.webdriver.support.ui import Select
from selenium import webdriver
import threading
from selenium.webdriver.common.by import By
from threading import Thread
from collections import Counter
import random
import pymysql
import os
#停词表
mu = threading.Lock()#线程锁，避免同时读取同意文件发生错误
stop_words=['hypothetical','protein','predictd:']#停词表
def isExist(driver,id):#判断元素是否存在try:driver.find_element_by_id(id)return Trueexcept:return False
def getSc(DESs,Accs):#得分函数。。这应该根据自己需要来定义all=[]for i in DESs:word=i.replace('[', ' ').replace(']', ' ')words=word.split(' ')for j in words:if j.lower() in stop_words:continueall.append(j)count=Counter(all).most_common(2)if len(count)==1:return DESs[0],Accs[0]Wmax=count[0][0]#出现频率最大Wmax0=count[1][0]#第二大#高频得分W0=count[1][1]*0.1W=(count[0][1]-count[1][1])*0.2+W0n=2score=[]for i in range(len(Accs)):sc=0if 'PREDICTED'in DESs[i]:sc=sc-0.4if Wmax in DESs[i]:sc=Welif Wmax0 in DESs[i]:sc=W0else:sc=0Tsc=n+scscore.append(Tsc)n=n-0.2ind=score.index(max(score))return DESs[ind],Accs[ind]
def dele(upid): #删除取得的RIDconn = pymysql.connect(host='47.106.121.42', port=3306, user='root', passwd='luanchen', db='blast')cur = conn.cursor()sql = "delete from submit where RID='%s'"% upidtry:cur.execute(sql)except:print(sql)conn.commit()cur.close()conn.close()
def get_table_content(driver,tableId):#得到表内内容arr=[]arr1=[]DESs=[]Accs=[]table_loc=(By.ID,tableId)
# 按行查询表格的数据，取出的数据是一整行，按空格分隔每一列的数据  try:table_tr_list=driver.find_element(*table_loc).find_elements(By.TAG_NAME,"tr")except:driver.refresh()driver.switch_to_alert().accept()n=0for tr in table_tr_list:n = n + 1if n==1:continueif n>11:breakarr1=(tr.text).split(" ")#以空格拆分成若干个(个数与列的个数相同)一维列表  DES=' '.join(arr1[:-6])Acc=arr1[-1]DESs.append(DES)Accs.append(Acc)des,acc=getSc(DESs,Accs)return des,acc
def output(suc):#输出到文件mu.acquire()#线程锁file = 'D:\\fanal.txt'f = open(file, 'a+')for i in suc:if len(i)<2:continueindex=i[0]n=index[:index.index('T')]id=index[index.index('T'):]f.write('%s\t%s\t%s\t%s\n'%(n,id,i[1],i[2]))f.close()mu.release()#释放锁
def do(driver,ht):#进入结果页面后的操作while(True):if isExist(driver, 'msgFrm') or isExist(driver, 'statInfo'):#这是两种比对错误的情况，一种是一直不出结果（服务器原因），另一种是比对序列太长被服务器拒绝print('wrong!!!')driver.quit()return -1elif isExist(driver, 'descriptions')==False:driver.quit()driver = webdriver.Firefox()driver.get('https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=%s' % ht)sleep(4)continueelif isExist(driver, 'queryList'):breakelse:continues = Select(driver.find_element_by_id('queryList'))Nop=s.optionssuc=[]shitbug=[]for i in Nop:shitbug.append(i.text)n=0b=0l=len(shitbug)while(True):#存放选项nameif isExist(driver, 'msgFrm') or isExist(driver, 'statInfo'):#比对失败print('wrong!!!')driver.quit()return -1if(n>=l):breakif '*' in shitbug[n]:#没有比对上item=[]#为了对齐item.append(shitbug[n].split(' ')[1])item.append(' ')item.append(' ')#outline(item)suc.append(item)n=n+1continuesleep(2)try:s = Select(driver.find_element_by_id('queryList'))s.select_by_index(n)sleep(5)#s = Select(driver.find_element_by_id('queryList'))item = []id = driver.find_element_by_xpath("//div[@id='querysummary']/dl/dd[3]").textif id not in shitbug[n]:try:driver.execute_script('GetReasult()')continueexcept:continuedes, acc = get_table_content(driver, 'dscTable')##这是处理得到结果表格的操作，本文主要是通过结果表格每一行的情况得到最合适的describtion和accessitem.append(id)item.append(des)item.append(acc)#outline(item)suc.append(item)n=n+1b=b+1except:driver.refresh()try:driver.switch_to_alert().accept()except:continuesleep(5)driver.quit()if b<5:return 0output(suc)print('done')return suc
def read(n):#从submit表中取一条RID并在表中删除，再构建GET方法，通过打开的浏览器访问GET方法，返回浏览器驱动，再执行do()#参数n当时用来标记的，这段代码没有用到if os.path.exists('C:\\stop.txt'):#作为程序停止的标记，当C:\\stop.txt存在时，程序停止运行return 0conn = pymysql.connect(host='127.0.0.1', port=3306, user='root', passwd='luanchen', db='hp')cur = conn.cursor()cur.execute('LOCK tables hp.submit write')#上锁sql = 'select RID,len from hp.submit where now()-TTime>1200 order by TTime limit %s' % 1#SQL语句查询上传时间超过1200秒的RIDcur.execute(sql)n = 0ht=''for i in cur:ht=(i[0])if ht=='':return 0sql = "delete from submit where RID='%s'" % httry:cur.execute(sql)except:print(sql)conn.commit()cur.execute('Unlock tables;')cur.close()conn.close()driver = webdriver.Firefox()#driver = webdriver.PhantomJS(executable_path="E:/phantomjs-2.1.1-windows/bin/phantomjs")while(True):try:driver.get('https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=%s' % ht)#get方法得到结果页面sleep(10)breakexcept:driver.quit()#river = webdriver.PhantomJS(executable_path="E:/phantomjs-2.1.1-windows/bin/phantomjs")driver = webdriver.Firefox()driver.get('https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=%s' % ht)sleep(10)breaksuc=do(driver,ht)read(n)
def test():t=[]for i in range(15):#这里设置线程数，根据CPU需要设置，建议使用线程池T=Thread(target=read, args=(i,))T.start()sleep(6)
test()

六、总结

1. 第四过程的最好使用线程池管理线程。
2. 在一次操作后，会有部分比对失败的序列，应该分析待比对序列和结果，看哪些不在结果中，再将这些漏掉的重新导入mysql table2表中，这个操作不建议在mysql中操作，而用python操作，因为sql语句执行太慢。把每次的待比对序列放入一个备份表是会更方便的
3. 根据设备和网络情况设置线程数。