DNA의 많은 서열 조각이 들어있는 fasta 파일 모음이 있습니다. 각 파일에서 찾을 수있는 k-mers의 총 발생 수를 계산하려고합니다. k-mers 계산의 좋은 부분은 크기가 4 ** ** k 인 단일 배열을 만들 수 있다는 것입니다. 여기서 k는 사용되는 k-mer의 크기입니다. 처리중인 시퀀스 파일은 새로운 세대 시퀀싱 머신의 짧은 읽기 시퀀스이므로 읽기가 모두 5 '-> 3'끝에서 이루어진다 고 가정하면 수행 할 수 없습니다. 이 문제를 해결하는 가장 좋은 방법은 관찰 된 모든 k-mer를 정방향 및 역방향 보완 시퀀스에 대한 단일 인덱스로 매핑하는 것입니다.파이썬에서 DNA 시퀀스의 순방향 및 역방향 보완을 붕괴시키는 알고리즘?
원하는 매핑 : 어레이에 대한 K = 2 & 시작 인덱스
0
문자열 = 'A'이고; 색인에 매핑 -> 0
문자열 = 'tt'; 색인에 매핑 -> 0
문자열 = 'at'; AA를 : - 인덱스 지도> 1
손으로는 알아낼이 앞으로의 붕괴 모든 mers의 배열 및 보완의 길이가 10이 될 것입니다 및 특정 인덱스는 다음 mers의를 나타낼 것 역전 할 수 있었다
큰 크기의 k에 대해 가능한 mers 수를 알기 위해 일반화 된 알고리즘을 생각하는 데 문제가 있습니다. count 배열에 몇 개의 셀을 할당해야합니까?
필자는 기존의 코드에서이 세 가지 함수를 사용하여 조각을 처리하고 역 보체를 생성하며 mer (또는 역방향 보수)를 인덱스에 매핑합니다.
이 첫 번째 함수는 mer 문자열을 가져 와서 4 * k 크기 배열의 mer와 관련된 인덱스를 반환합니다.
def mer_index_finder(my_string, mer_size):
# my_string = my_string.lower()
char_value = {}
char_value["a"] = 0
char_value["t"] = 1
char_value["c"] = 2
char_value["g"] = 3
i = 0
j = 0
base_four_string = ""
while(i < mer_size):
base_four_string += str(char_value[my_string[i]])
i += 1
index = int(base_four_string, 4)
return index
이 함수는 DNA 단편을 모두 처리하고 **
def get_mer_count(mer_size, file_fragments, slidingSize):
mer_counts = {}
for fragment in file_fragments:
j = 0
max_j = len(fragment) - mer_size
while(j < max_j):
mer_frag = fragment[j:j+mer_size]
mer_frag = mer_frag.lower()
if("n" not in mer_frag):
try:
mer_counts[mer_frag] += 1
except:
mer_counts[mer_frag] = 1
j += slidingSize
myNSV = [0] * (4**mer_size)
for mer in mer_counts.keys():
mer_index = mer_index_finder(mer, mer_size)
# examples showing how to collapse,
# without shrinking the array
# rev_mer = make_complment_mer(mer)
# print rev_mer
# rev_index = mer_index_finder(rev_mer, mer_size)
# min_index = min(mer_index, rev_index)
# print mer_index,"\t",rev_index,"\t",min_index
# myNSV[min_index] += mer_counts[mer]
myNSV[mer_index] = mer_counts[mer]
return myNSV[:]
마지막이 함수는 메르 걸릴 K와 역방향 보수를 생성 사이즈 (4)의 배열의 인덱스 카운트 매핑 :
def make_complment_mer(mer_string):
nu_mer = ""
compliment_map = {"a" : "t", "c" : "g", "t" : "a", "g" : "c"}
for base in mer_string:
nu_mer += compliment_map[base]
nu_mer = nu_mer[::-1]
return nu_mer[:]
항상 배열이 앞으로 붕괴 할 때이 함께 보완 반대한다 얼마나 많은 세포를 알 수있는 확실한 방법이 있어야 것 같습니다, 그리고 시험이있다 문학에 관한 내용과이를 보여주는 일부 패키지가 완성되었습니다. 그러나 소스 코드에서이 계산을 생성 할 수있는 위치를 찾지 못했습니다.
이 질문의 두 번째 부분은 머가 모두를 확인하지 않고 순방향 또는 역방향 보완 여부를 어떻게 알 수 있습니까?
예 :
(순방향)AAGATCACGG
(보수)
TTCTAGTGCC
(후진 보수)
CCGTGATCTT
위의 코드에서 두 인덱스 중 낮은 값을 취합니다. 그러나 한 번 전진하고 한 번 역방향으로 보완하여 각 mer에 대한 인덱스를 두 번 찾을 필요없이이를 파악할 수있는 방법이있는 것처럼 보입니다.
TL; DR 정방향 및 역방향 보완이 동일한 인덱스에 매핑되는 경우 배열의 크기는 어떻게됩니까?
편집 : 내가 인덱스의 크기를 만들려면 다음 행을 포함하는 get_mer_count() 수정 답변을 사용하여 배열의 크기를 확인하려면 각 k
-mer를 들어
array_size = (4 ** mer_size)/2
if mer_size % 2 == 0:
array_size += 2**(mer_size - 1)
myNSV = [0] * array_size
굉장! 고맙습니다! 해결책을 간단한 if 문으로 변경했습니다. 편집을 참조하십시오. –