这里承诺的是我写的脚本。作为参考,我试图用Burrows-Wheeler变换在DNA串中进行重复匹配。基本上这个想法是采取一段长度为M的DNA链,并在该串中找到所有的重复序列。所以,举个例子,如果我有奇怪的acaacg并且搜索了大小为2的所有重复的子字符串,我会得到1的计数和0,3的起始位置。然后你可以键入字符串[0:2]和字符串[3:5]来验证它们确实匹配,并且它们的结果是“ac”。
如果有人有兴趣了解的挖洞轮车,就可以了维基百科搜索产生非常有用的结果。这是斯坦福大学的另一个来源,也很好地解释了这一点。 http://www.stanford.edu/class/cs262/notes/lecture5.pdf
现在,我还没有解决一些问题。首先,我使用n^2空间来创建BW变换。另外,我创建了一个后缀数组,对其进行排序,然后用数字替换,以创建可能占用一定空间的数字。然而,最后我只是真的存储了occ矩阵,最后一列和单词本身。
尽管内存问题的字符串大于4^7(得到这个工作的字符串大小为40,000,但没有更大的...),我会称这是一个成功,看到像以前的星期一,我唯一的新的在python中如何做是为了让它打印我的名字和你好世界。
# generate random string of DNA
def get_string(length):
string=""
for i in range(length):
string += random.choice("ATGC")
return string
# Make the BW transform from the generated string
def make_bwt(word):
word = word + '$'
return ''.join([x[-1] for x in
sorted([word[i:] + word[:i] for i in range(len(word))])])
# Make the occurrence matrix from the transform
def make_occ(bwt):
letters=set(bwt)
occ={}
for letter in letters:
c=0
occ[letter]=[]
for i in range(len(bwt)):
if bwt[i]==letter:
c+=1
occ[letter].append(c)
return occ
# Get the initial starting locations for the Pos(x) values
def get_starts(word):
list={}
word=word+"$"
for letter in set(word):
list[letter]=len([i for i in word if i < letter])
return list
# Single range finder for the BWT. This produces a first and last position for one read.
def get_range(read,occ,pos):
read=read[::-1]
firstletter=read[0]
newread=read[1:len(read)]
readL=len(read)
F0=pos[firstletter]
L0=pos[firstletter]+occ[firstletter][-1]-1
F1=F0
L1=L0
for letter in newread:
F1=pos[letter]+occ[letter][F1-1]
L1=pos[letter]+occ[letter][L1] -1
return F1,L1
# Iterate the single read finder over the entire string to search for duplicates
def get_range_large(readlength,occ,pos,bwt):
output=[]
for i in range(0,len(bwt)-readlength):
output.append(get_range(word[i:(i+readlength)],occ,pos))
return output
# Create suffix array to use later
def get_suf_array(word):
suffix_names=[word[i:] for i in range(len(word))]
suffix_position=range(0,len(word))
output=zip(suffix_names,suffix_position)
output.sort()
output2=[]
for i in range(len(output)):
output2.append(output[i][1])
return output2
# Remove single hits that were a result of using the substrings to scan the large string
def keep_dupes(bwtrange):
mylist=[]
for i in range(0,len(bwtrange)):
if bwtrange[i][1]!=bwtrange[i][0]:
mylist.append(tuple(bwtrange[i]))
newset=set(mylist)
newlist=list(newset)
newlist.sort()
return newlist
# Count the duplicate entries
def count_dupes(hits):
c=0
for i in range(0,len(hits)):
sum=hits[i][1]-hits[i][0]
if sum > 0:
c=c+sum
else:
c
return c
# Get the coordinates from BWT and use the suffix array to map them back to their original indices
def get_coord(hits):
mylist=[]
for element in hits:
mylist.append(sa[element[0]-1:element[1]])
return mylist
# Use the coordinates to get the actual strings that are duplicated
def get_dupstrings(coord,readlength):
output=[]
for element in coord:
temp=[]
for i in range(0,len(element)):
string=word[element[i]:(element[i]+readlength)]
temp.append(string)
output.append(temp)
return output
# Merge the strings and the coordinates together for one big list.
def together(dupstrings,coord):
output=[]
for i in range(0,len(coord)):
merge=dupstrings[i]+coord[i]
output.append(merge)
return output
Now run the commands as follows
import random # This is needed to generate a random string
readlength=12 # pick read length
word=get_string(4**7) # make random word
bwt=make_bwt(word) # make bwt transform from word
occ=make_occ(bwt) # make occurrence matrix
pos=get_starts(word) # gets start positions of sorted first row
bwtrange=get_range_large(readlength,occ,pos,bwt) # Runs the get_range function over all substrings in a string.
sa=get_suf_array(word) # This function builds a suffix array and numbers it.
hits=keep_dupes(bwtrange) # Pulls out the number of entries in the bwt results that have more than one hit.
dupes=count_dupes(hits) # counts hits
coord=get_coord(hits) # This part attempts to pull out the coordinates of the hits.
dupstrings=get_dupstrings(coord,readlength) # pulls out all the duplicated strings
strings_coord=together(dupstrings,coord) # puts coordinates and strings in one file for ease of viewing.
print dupes
print strings_coord
感谢您的帮助。我能够按照我想要的方式继续下去。我不太熟悉.append命令,这正是我所需要的。此外,空对象的分配也有所帮助。完成后我会发布我的脚本。 – doggysaywhat 2012-04-09 10:57:12
看看你可以用列表[这里]做什么(http://docs.python.org/tutorial/datastructures.html#more-on-lists)。 – 2012-04-09 11:00:14