일단 처음으로 input data는 txt파일이다.
position은 covid sequence(or codon)의 위치이다. 그리고 frequency는 그 군집에서 얼마나 변이가 나타났는지에 대한 count이다. percentage는 count/total(군집의 수)이며, 얼마나 그 군집에서 변이가 잘 나타났는지에 대한 척도이다. Entropy는 데이터의 신뢰도에 대한 척도인데, 이번 프로젝트에는 다루지 않기 때문에 신경 쓸 필요는 딱히 없다.
Position Frequency Percentages Entropy
1 0 0.0 0.0
2 0 0.0 0.0
3 0 0.0 0.0
4 0 0.0 0.0
5 1 0.33333333333333337 0.03223035521176948
6 0 0.0 0.0
7 0 0.0 0.0
8 0 0.0 0.0
9 1 0.33333333333333337 0.03223035521176948
10 0 0.0 0.0
11 0 0.0 0.0
12 0 0.0 0.0
13 13 4.333333333333334 0.291143835806015
14 1 0.33333333333333337 0.03223035521176948
이 입력 파일을 parameter로 입력하게 되면, 이 중에서 해당 percentage가 넘는 position이 있을때 그 position을 저장한다. 그 position을 ccm이라고 하자.
// ccm result
idx pos i_freq mut_n freq_sum freq_avg per_sum per_avg entr_sum entr_avg es deps
12 13 13 41 43 1 14.333333333333334 0.34959349593495936 0.9493715602930303 0.023155403909586106 2.8888888888888893 39
24 25 11 49 43 0 14.333333333333334 0.2925170068027211 0.9493715602930303 0.019374929801898578 2.444444444444444 47
25 26 15 59 63 1 21.0 0.3559322033898305 1.3218240257888483 0.022403797047268616 3.3333333333333335 57
51 52 20 90 106 1 35.333333333333336 0.3925925925925926 2.1678874040153704 0.024087637822393003 4.4444444444444455 88
96 97 29 130 90 0 30.0 0.23076923076923078 1.8159464426311578 0.013968818789470445 6.444444444444444 128
125 126 17 76 56 0 18.666666666666664 0.24561403508771926 1.0766881959896575 0.014166949947232336 3.7777777777777772 74
208 209 0 2 62 31 20.666666666666668 10.333333333333334 0.849716847882154 0.424858423941077 0.0 0
이렇게 한번 search를 하면 ccm이 나올텐데, 그 ccm을 가지고 cluster를 만들어 나간다. cluster는 좌,우로 기준(100bps)이 되는 길이를 확장시켜나가면서 percentage가 일정수준이하로 내려갈때 까지 확장한다.
percentage가 일정 수준으로 내려가면 더이상 확장을 멈추고 확장된 범위까지를 cluster로 지정하고 다음 ccm으로 넘어가 이 과정을 반복한다. 이렇게 되면 ccm마다 cluster가 만들어진다.
// clusters
ccm_pos idx lpos rpos lidx ridx length
13 12 1 31 0 30 30
25 24 9 41 8 40 32
26 25 6 45 5 44 39
52 51 30 75 29 74 45
97 96 70 124 69 123 54
126 125 106 146 105 145 40
210 209 171 249 170 248 78
220 219 215 225 214 224 10
440 439 403 478 402 477 75
프로젝트를 시작할때 받아온 데이터가 전처리가 어느정도 되어 있어서 fastq나 fasta parsing을 하지 않아도 됐다.
아래 데이터 처럼 변이가 있는 경우에는 position 오른쪽에 어떤 변이가 일어났는지에 대한 정보가 나와있었기 때문에 전처리 하는 과정이 오래걸리지는 않았다.
그리고 B.1.620처럼 covid 변이 종류를 나타내는 정보도 있었기 때문에 변이 별로 데이터를 정리 할 수도 있었다.
Cusomer ID,,COV_SCO_074
LINEAGE,,B.1.620
REF_base,REF_position,CR_21_01129_CV_CV
A,1,
T,2,
T,3,
...
A,22,
G,23,
G,24,T
...
아래 코드를 통해서 위에 있는 데이터를 mutclust를 돌리기 위해 입력파일 format으로 수정하는 코드를 작성했다. 종종 변이가 들어가있는 column에 None이 들어가있는 경우가 있었기 때문에 isna함수를 사용해서 예외처리를 진행해주었다.
header = 'Position' + '\\t' + 'Frequency' + '\\t' + 'Percentages'+ '\\t' + 'Entropy'+'\\n'
for file in path:
index = file.split('.')
name = index[0].replace('-','_')
data = pd.read_csv(init+file)
filename = data[name][0]
table[filename] = table.get(filename, 0) + 1
for i in range(2,len(data)):
if pd.isna(data[name][i]) is False: # and filename == kind[10]
if file not in scoreTable.keys():
print(file)
pass
else:
frequency[scoreTable[file]][int(data['Unnamed: 1'][i])]+=1
f = open("/data3/projects/2020_MUTCLUST/data/covid/newsscore/score_type"+str(i+1)+".txt", "w")
f.write(header)
for k in range(1,len(frequency[i])):
insert = str(k)+'\\t'+str(frequency[i][k])+'\\t'+str(frequency[i][k]/(len(data)-2)*100) + '\\t'+str(0)+'\\n'
f.write(insert)
f.close()