一、問題呈現
找到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('\n','') 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]+'\n') 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:\temporary\hrdb_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含量過高,用專門的細菌啟動子預測軟件預測不出來或者極其不准,也沒有找到合適的在線軟件。