python 根据染色体起始终止点坐标来获取碱基序列

Posted 2021-01-31 Yellow_huang

在工作当中，有时候我们知到染色体编号以及染色体起始终止坐标，我们想知道这段序列是什么样的碱基。
其一，我们一般用去UCSC的genome browser里面去查询 ，其实也可以从UCSC的接口去解析网页，然后在提取序列信息
比如chr17:7676091,7676196 ，那么我只需要构造下面一个网页地址 http://genome.ucsc.edu/cgi-bin/das/hg38/dna?segment=chr17:7676091,7676196
然后 hg38可以更换成hg19，dna?segment= 后面可以按照标准格式更换，就可以返回我们想要的序列了。现在对网页返回 一个xml格式的信息，用python解析一下

  1 import requests

  2 import re
  3 from bs4 import BeautifulSoup
  4 import xlwt
  5 import xlrd
  6 from xlutils.copy import copy
  7 import os ,sys
  8 #print(sys.path)
  9 cwd=os.getcwd()
 10 
 11 
 12 
 13 def getHTMLText(url):
 14     print("111111")
 15     try:
 16         header = {‘user-agent‘: ‘Mozilla/5.0‘}
 17         r = requests.get(url,headers = header,timeout = 30 )
 18         r.raise_for_status()
 19         r.encoding =r.apparent_encoding
 20         print("get 222222222222222")
 21         return r.text
 22     except:
 23         return ""
 24 
 25 def fillDNAList(dnalist, html):
 26     # 使用正则表达式获取dna 序列的头文件
 27     match = re.search(‘SEQUENCE([sS]*?version="1.00")‘, html)
 28     print("match ok")
 29     if match:
 30         dna_header = re.search(‘SEQUENCE([sS]*?version="1.00")‘, html)
 31         #print(‘10====>‘, dna_header.group())
 32         #dna_header 存到列表
 33         dnalist.append(dna_header.group())
 34 
 35     match = re.search(‘DNA.*?length="(d*)"‘, html)
 36     if match:
 37         length_header= re.search(‘DNA.*?length="(d*)"‘, html)
 38         #print(‘11=====>‘, length_header.group())
 39         dnalist.append(length_header.group())
 40 
 41     #使用 BeautifulSoup
 42     print("BeautifulSoup tag属性 获取dna标签属性的字符串部分")
 43     soup = BeautifulSoup(html, ‘html.parser‘)
 44     tag = soup.dna
 45     seq = soup.dna.string
 46     seq = seq.replace(‘
‘,‘‘).upper()  #换行符删除掉,转换成大写
 47     # seq 存到列表
 48     dnalist.append(seq)
 49     print("final======>",dnalist)
 50     return dnalist
 51 
 52 def write_excell(dnalist,chrnum,pos):
 53     head = ‘>hg19‘ + ‘ ‘ + dnalist[0] + dnalist[1]
 54     f = xlwt.Workbook(encoding=‘utf-8‘, style_compression=0) #创建新的Excel（新的workbook）
 55     sheet = f.add_sheet(‘test8‘, cell_overwrite_ok=True) #创建新的表单
 56     # 先写第一行的头文件
 57     sheet.write(0,0,head)
 58     #再从第二行开始写，每行写入50 个字符串
 59     dna = dnalist[2]
 60     print(‘=====‘,dna,type(dna))
 61     for i in range(0,len(dna),50):
 62         sheet.write((int((i+1)/50))+1,0,dna[i:i+50])
 63 
 64     out_file = ‘chrmosome%s_%s.xls‘% (chrnum,pos)
 65     f.save(out_file)
 66     out_file_dir = os.path.join(cwd, out_file)
 67     return out_file_dir
 68 
 69 def modify_excell(out_file_dir,chrnum,pos):
 70     ‘‘‘
 71     改Excel表（xlutils模块）
 72     :return:
 73     ‘‘‘
 74     rb = xlrd.open_workbook(out_file_dir)  # 打开out_file.xls文件,创建工作簿实例对象
 75     sheet = rb.sheet_by_index(0)
 76     nrow11 = sheet.cell_value(10, 0) #修改第11行第一列，索引是10，0
 77     # 根据需要截取原单元格里面的内容与需要添加的内容进行拼接
 78     new_nrow11 = ‘[‘ + nrow11
 79     # 同理操作nrow12
 80     nrow12 = sheet.cell_value(11, 0)
 81     new_nrow12 = nrow12 + ‘]‘
 82 
 83     wb = copy(rb)
 84     ws = wb.get_sheet(0)
 85     # 往单元格中写入拼接后的新字符串内容
 86     ws.write(10,0,new_nrow11)
 87     ws.write(11, 0, new_nrow12)
 88     modify_file = ‘new_chrmosome%s_%s.xls‘ % (chrnum,pos)
 89     wb.save(modify_file)
 90 
 91 def main():
 92     hg19 = "hg19"
 93     chrnum = 17
 94     pos = 7676091
 95     start = pos - 500
 96     end = pos + 500
 97     position_DNA_list = []
 98     #url = "http://genome.ucsc.edu/cgi-bin/das/hg38/dna?segment=chr17:7676091,7676196"
 99     url = f"http://genome.ucsc.edu/cgi-bin/das/{hg19}/dna?segment=chr{chrnum}:{start},{end}"
100    
101     print(url)
102     html = getHTMLText(url)
103     dnalist = fillDNAList(position_DNA_list,html)
104     out_file_dir = write_excell(dnalist,chrnum,pos)
105     modify_excell(out_file_dir,chrnum,pos)
106 
107 main()

结果如下：

http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr17:7675591,7676591
111111
get 222222222222222
match ok
BeautifulSoup tag属性 获取dna标签属性的字符串部分
final======> [‘SEQUENCE id="chr17" start="7675591" stop="7676591" version="1.00"‘, ‘DNA length="1001"‘, ‘CCCAAGAGCCTTCAGTATACACATCAATAAAAATAATTTTAATTATTCTGATAAAAGATAAACATGAAAAGTTATGGTATGCAAAGTTGAATGACAACAACTGATACTATTTGAAATAATTGACAGAATTATATTCCGTAACAATTTATAAGCAAAGCCAAAAAAACAATGATCCCTTTGTTGAATGCACAGAACAAATCCATCTTGTCCACGGCTACTGAGCATGCCTGTGATCTCCAGGGGTCACTCAGGTTTGACTCAAAGGATCCAACAGCCTGTAGACCCTGTGCTTGAAGGCATGAGGGTCACCTCTGAGTTCACACTCACTAGTGTCCCTCCTTTCTTCAGAAAGCTAGGAACTGGGAAGACAAGGGGAAAATCAATCAAGGCCTGAGGTATGGGGCTGTAGGCTGGGAGGAAACTAACATTATTGAGAAGCTACTGATGTGAATACATTTCAATTACTACTCACATTGGTTTTTTGTTTGTTTGTTTGTTTGTTTGTTTGTTTGTTTTTTAAGACGGAGTTTTGCTCTCGTTGCCCAGGCTGGAGTGCAATGGAATGATCTAAGGTCACCACAACCTCCACCTCCCGGTTCAAGCAATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGTGTGCCACCACACCCAGCTAAGTTTGTATTTTTTTAGTAGAGACGGTGTTTCACCATGTTGGTCAGGCTGGTCTCGAACTCCTGACCTCAAGTGATCCACCCACCTCGGCCTCCCAAAGTGCTGGGATTATAGGCATGAGCCACCACACCCAGCCTCACGTTGGTTTTTGAGATGGATTTTATTGCCATTTTGTACACAAAAAGGTCAAAACTCAGTGAGGTGAATTGACATGACAGTAAGTGAAAGAACTACTATCTGATTGGGGGTCTTCTGCCGCCTGCTCTGGGACTCTTTCTGCTATGACATGAAGGACATTGGCAACCCCAGTCCTTGCAGATTTCTTTCACTGTGTGC‘]
===== CCCAAGAGCCTTCAGTATACACATCAATAAAAATAATTTTAATTATTCTGATAAAAGATAAACATGAAAAGTTATGGTATGCAAAGTTGAATGACAACAACTGATACTATTTGAAATAATTGACAGAATTATATTCCGTAACAATTTATAAGCAAAGCCAAAAAAACAATGATCCCTTTGTTGAATGCACAGAACAAATCCATCTTGTCCACGGCTACTGAGCATGCCTGTGATCTCCAGGGGTCACTCAGGTTTGACTCAAAGGATCCAACAGCCTGTAGACCCTGTGCTTGAAGGCATGAGGGTCACCTCTGAGTTCACACTCACTAGTGTCCCTCCTTTCTTCAGAAAGCTAGGAACTGGGAAGACAAGGGGAAAATCAATCAAGGCCTGAGGTATGGGGCTGTAGGCTGGGAGGAAACTAACATTATTGAGAAGCTACTGATGTGAATACATTTCAATTACTACTCACATTGGTTTTTTGTTTGTTTGTTTGTTTGTTTGTTTGTTTGTTTTTTAAGACGGAGTTTTGCTCTCGTTGCCCAGGCTGGAGTGCAATGGAATGATCTAAGGTCACCACAACCTCCACCTCCCGGTTCAAGCAATTCTCCTGCCTCAGCCTCCCAAGTAGCTGGGACTACAGGCGTGTGCCACCACACCCAGCTAAGTTTGTATTTTTTTAGTAGAGACGGTGTTTCACCATGTTGGTCAGGCTGGTCTCGAACTCCTGACCTCAAGTGATCCACCCACCTCGGCCTCCCAAAGTGCTGGGATTATAGGCATGAGCCACCACACCCAGCCTCACGTTGGTTTTTGAGATGGATTTTATTGCCATTTTGTACACAAAAAGGTCAAAACTCAGTGAGGTGAATTGACATGACAGTAAGTGAAAGAACTACTATCTGATTGGGGGTCTTCTGCCGCCTGCTCTGGGACTCTTTCTGCTATGACATGAAGGACATTGGCAACCCCAGTCCTTGCAGATTTCTTTCACTGTGTGC <class ‘str‘>