2016-12-03 3 views
2

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 

답변

4

을,이 두 가능성 : 그것은 정확히 하나의 역 보완 물을 갖거나, 자신의 역의 칭찬 ("회문 변증")이다. 따라서 p 회문색 k -mer가 있다면 배열 크기는 p + (4**k - p)/2이어야합니다. 중간 염기가 자신의 칭찬이 될 수 없기 때문에

  • 이상한 k를 들어, 어떤 상동 mers의이 없습니다. 배열의 크기는 4**k/2이어야합니다.

  • k의 경우에도 j의 경우 k = 2*j입니다. 메르는 상반기가 후반기의 역 칭찬 인 경우에만 회문색을 나타낸다. 가능한 첫 번째 반쪽은 4**j이며, 따라서 p = 4**j = 2**k 회문색 k - mer가 있습니다. 따라서 배열 위의 공식을 사용하려면 크기가 p + (4**k - p)/2 = 2**k + (4**k - 2**k)/2이어야합니다.

+0

굉장! 고맙습니다! 해결책을 간단한 if 문으로 변경했습니다. 편집을 참조하십시오. –