PythonDNA序列中子序列出现频率

Python DNA序列在使用的时候有很多需要我们注意的东西，其实在不断的学习中有很多问题存在，下面我们就详细的看看如何进行相关的技术学校。ms是我师弟的rotation project：给定一堆Python DNA序列，即由字符A, C, G, T组成的字符串，统计所有长度为n的子序列出现的频率。

站在用户的角度思考问题，与客户深入沟通，找到蕉岭网站设计与蕉岭网站推广的解决方案，凭借多年的经验，让设计与互联网技术结合，创造个性化、用户体验好的作品，建站类型包括：网站建设、网站设计、企业官网、英文网站、手机端网站、网站推广、主机域名、虚拟主机、企业邮箱。业务覆盖蕉岭地区。

比如 ACGTACGT，子序列长度为2，于是 AC=2, CG=2, GT=2, TA=1，其余长度为2的子序列频率为0.

最先想到的就是建一个字典，key是所有可能的子序列，value是这个子序列出现的频率。Python DNA序列但是当子序列比较长的时候，比如 n=8，需要一个有65536 (4的8次方) 个key-value pair的字典，且每个key的长度是8字符。这样ms有点浪费内存。。

于是想到，所有的长度为n的子序列是有序且连续的，所以可以映射到一个长度为4的n次方的的list里。令 A=0, C=1, G=2, T=3，则把子序列 ACGT 转换成 0*4^3 + 1*4^2 + 2*4 + 3 = 27, 映射到list的第27位。如此，list的index对应子序列，而list这个index位置则储存这个子序列出现的频率。

于是我们先要建立2个字典，表示ACGT和0123一一对应的关系：

 
 
 
  
  
  i2mD = {0:'A', 1:'C', 2:'G', 3:'T'} 
  
  
  m2iD = dict(A=0,C=1,G=2,T=3) 
  
  
  # This is just another way to initialize a 
dictionary

以及下面的子序列映射成整数函数：

 
 
 
  
  
  def motif2int(motif): 
  
  
  '''convert a sub-sequence/motif to a non-negative 
integer''' 
  
  
  total = 0
  
  
  for i, letter in enumerate(motif): 
  
  
  total += m2iD[letter]*4**(len(motif)-i-1) 
  
  
  return total 
  
  
  Test: 
  
  
  >>> motif2int('ACGT')

虽然我们内部把子序列当成正整数来存储（确切地说，其实这个整数是没有存在内存里的，而是由其在list的index表示的），为了方便生物学家们看，输出时还是转换回子序列比较好。

于是有了下面的整数映射成子序列函数，其中调用了另外一个函数baseN()，来源在此，感谢作者～

 
 
 
  
  
  def baseN(n,b): 
  
  
  '''convert non-negative decimal integer n to 
  
  
  equivalent in another base b (2-36)''' 
  
  
  return ((n == 0) and '0' ) or ( baseN(n // b, b).lstrip('0') + \ 
  
  
  "0123456789abcdefghijklmnopqrstuvwxyz"[n % b]) 
  
  
  def int2motif(n, motifLen): 
  
  
  '''convert non-negative integer n to a sub-sequence/motif with length motifLen''' 
  
  
  intBase4 = baseN(n,4) 
  
  
  return ''.join(map(lambda x: i2mD[int(x)],'0'*(motifLen-len(intBase4))+intBase4)) 
  
  
  Test: 
  
  
  >>> int2motif(27,4) 
  
  
  'ACGT'

以下代码从命令行读入一个存有DNA序列的fasta文件，以及子序列长度，并输出子序列和频率。注意以下代码需要Biopython module。

 
 
 
  
  
  if __name__ == '__main__': 
  
  
  import sys 
  
  
  from Bio import SeqIO 
  
  
  # read in the fasta file name and motif length 
  
  
  # from command line parameters 
  
  
  fastafile = sys.argv[1] 
  
  
  motifLen = int(sys.argv[2]) 
  
  
  # list to store subsequence frequency 
  
  
  frequencyL = [0]*4**motifLen 
  
  
  # go over each DNA sequence in the fasta file 
  
  
  # and count the frequency of subsequences 
  
  
  it = SeqIO.parse(open(fastafile),'fasta') 
  
  
  for rec in it: 
  
  
  chrom = rec.seq.tostring() 
  
  
  for i in range(len(chrom)-motifLen+1): 
  
  
  motif = chrom[i:i+motifLen] 
  
  
  frequencyL[motif2int(motif)] += 1 
  
  
  # print frequency result to screen 
  
  
  for i, frequency in enumerate(frequencyL): 
  
  
  print int2motif(i, motifLen), frequency

以上就是创新互联对Python DNA序列的相关介绍。

分享名称：PythonDNA序列中子序列出现频率
文章URL：http://www.shufengxianlan.com/qtweb/news17/114967.html

网站建设、网络推广公司-创新互联，是专注品牌与效果的网站制作，网络营销seo公司；服务项目有等

声明：本网站发布的内容（图片、视频和文字）以用户投稿、用户转载内容为主，如果涉及侵权请尽快告知，我们将会在第一时间删除。文章观点不代表本网站立场，如需处理请联系客服。电话：028-86922220；邮箱：631063699@qq.com。内容未经允许不得转载，或转载时需注明来源：创新互联

猜你还喜欢下面的内容