正如这里所承诺的,这是我写的完成的脚本。作为参考,我正在尝试使用 Burrows-Wheeler 变换在 DNA 字符串中进行重复匹配。基本上,这个想法是取一段长度为 M 的 DNA 链,并在该字符串中找到所有重复。因此,例如,如果我有奇怪的 acaacg 并搜索所有重复的大小为 2 的子字符串,我会得到 1 的计数和 0,3 的起始位置。然后您可以输入 string[0:2] 和 string[3:5] 来验证它们确实匹配并且它们的结果是“ac”。
如果有人有兴趣了解 Burrows-Wheeler,在 Wikipedia 上搜索它会产生非常有用的结果。这是斯坦福大学的另一个来源,它也很好地解释了它。 http://www.stanford.edu/class/cs262/notes/lecture5.pdf
现在,有几个问题我没有解决。首先,我使用 n^2 空间来创建 BW 变换。另外,我正在创建一个后缀数组,对其进行排序,然后用数字替换它,因此创建它可能会占用一些空间。但是,最后我只真正存储了 occ 矩阵、结束列和单词本身。
尽管大于 4^7 的字符串存在 RAM 问题(使其适用于 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