Python总括如何开展DNA连串整理

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

Python统计在我们的使用中有很多的障碍,其中在DNA序列上的相关问题就需要我们不断的去学习。下面我们就向大家介绍有关的问题,希望在以后Python统计的使用过程中有所收获。

某随机序列[12, 4, 6, 45, 12, 78,
…..]中,找出出现次数最多的三个元素
以及他们出现的次数

DNA序列比对可能出现的问题

xiongrongchuan@126.com

威尼斯人注册官网,如果比对出来的序列,保守位点较少(在mega软件中查看,序列顶端的*较少,如下图所示)

威尼斯人注册官网 1

说明序列存在一定的问题,需要进一步核对序列,找到原因后方能进行后续操作。

威尼斯娱乐,第一种原因,可能存在重复序列,一些较长的序列不同的区域可能和自测序列都有一定的相似性,从而在“通过blast截取同源序列”的过程中,这些长序列被截取成了不同的片段。

威尼斯人注册官网 2

遇到这种原因,则删除多余的重复序列(通常是碱基数小的那一条)。

第二种原因,可能是某些序列的方向和其它序列不一致。通常这种序列在mega软件中显得特别“突兀”。

威尼斯人注册官网 3

可以把这条“可疑”序列和其它的“正常序列”中的一条上传到Genbank中相互Blast后,验证方向是否一致。

威尼斯人注册官网 4

Blast后,核对两条序列起止数字是否都是终点大于起点,如果不都是终点大于起点,说明方向反了。

威尼斯人注册官网 5

对于这样的情况,需要将可以序列进行“反向互补操作”纠正其方向,才能进行后续分析。

操作方法是,在mega中选中“可疑序列”,在“Data”菜单栏中选择“Reverse
complement”。

威尼斯人注册官网 6

核心提示:NCBI的“electronic
PCR”工具是UniSTS资源库的一部分,可以用来寻找一段目的DNA片段中的STS标记物。UniSTS
能提供所有有关STS标记物的资料,包括引物序列、产物大小、作图信息和别名。与之相链接的其他NCBI资源如Entrez、LocusLink
和MapViewer
也同样提供这些信息。e-PCR通过搜寻具有正确的方向和间距的序列且这个序列能代表用于扩增STSs的PCR引物,来寻找一段DNA序列中潜在的STSs。

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

给定一堆DNA序列,即由字符A, C, G,
T组成的字符串,统计所有长度为n的子序列出现的频率。比如
ACGTACGT,子序列长度为2,于是 AC=2, CG=2, GT=2,
TA=1,其余长度为2的子序列频率为0.

  1. 利用默认值字典,先以键值对形式存储每个元素及出现的次数,然后将键值反转存储为元组,利用排序输出后三位

先在NCBI主页上(
)找到e-PCR的主页,然后在右手栏点击“Electronic
PCR”链接。再在e-PCR主页的上端大的文本框内粘贴上目的基因序列或键入登陆号(accession
number)。例如某个序列的登录号是AF288398,结果显示该序列只包含一个STS:
stSG47693,位于此序列的2102和2232核苷之间。当点击“Marker”下标记物的名称时,从UniSTS中出现STS的详细资料。引物的信息、PCR产物大小以及标记物的替代名称也出现在主页的上端。在不同的图谱中,STSs常有不同的名称。

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

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

在“Cross-references”栏目下的 LocusLink、UniGene 和 the Genebridge
4中,将显示这个STS的定位图。在“mapping information”部分包含能链接到NCBI
的“MapViewer”浏览器。在主页的下端是“Electronic PCR
results”,显示了其他序列,包括contigs、mRNAs和包含这个STS标记物的ESTs。

于是想到,所有的长度为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位置则储存这个子序列出现的频率。

于是想到,所有的长度为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位置则储存这个子序列出现的频率。

>>> from random import randint
>>> from collections import defaultdict
>>> data = [randint(0, 10) for _ in range(0, 30)]
>>> d = defaultdict(int)
>>> for x in data:
...     d[x] += 1
...
>>> d
defaultdict(<class 'int'>, {0: 5, 1: 1, 2: 3, 3: 6, 4: 5, 5: 1, 6: 3, 7: 1, 8: 2, 9: 1, 10: 2})
>>> new_d = sorted([(v, k) for k, v in d.items()])
>>> new_d[-3:]
[(5, 0), (5, 4), (6, 3)]

为了在所有图谱中看到STS标记物及其基因组的状况,则在“Mapping
Information”部分的上端点击链接标志 “MapViewer”
,这个图谱浏览器会出现两张图谱。请注意,在这个视窗里,STS
stSG47693被称为RH92759。99�CGenebridge 4 基因图谱上有46000个STS
标记被国际放射杂交协会定位到GB4杂交面板上。

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

于是我们先要建立2个字典,Python统计表示ACGT和0123一一对应的关系:

  1. 但是对于这个问题用Counter是最简单而且高效的解决办法

STS图谱显示了如何使用e-PCR将STSs序列放置到基因组序列组装。灰色线将两个图谱的标记物连接起来,而红色线条显示STS
RH92759在两张图谱中的位置。在这个区域,STS图谱中共有211个STSs,但在这个视窗里只标记了20个。在STS图谱的右边,点击绿色和黄色圆圈会出现STS标记物的图谱。通过左边工具条的缩放工具,可以放大或缩小这个视窗。

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 
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 

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

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

>>> from random import randint
>>> from collections import Counter
>>> data = [randint(0, 10) for _ in range(0, 30)]
>>> c = Counter(data)
>>> c.most_common(3)
[(7, 5), (0, 4), (1, 4)]
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') 
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')  27 

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

以上就是对Python统计的相关介绍。虽然我们内部把子序列当成正整数来存储确切地说,其实这个整数是没有存在内存里的,而是由其在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序列的相关介绍。

DNA序列在使用的时候有很多需要我们注意的东西,其实在不断的学习中有很多问题存在,下面我们就详细的看看如何进行相关的技术学…

相关文章