FASTQ의 경우 quality 부분과 read identifier를 제외한 순수한 서열만을 남겨둔다.
가령 uncorrected.fq -> uncorrected.read, uncorrected.qual
FASTA의 경우 read identifier를 제외한 순수한 서열만을 남겨둔다.
가령 uncorrected.fa -> uncorrected.read
위의 conversion을 수행하는 과정에서 max_read_length (가장 긴 read 길이)와 size_of_sequences (전체 서열의 길이) 저장한다.
대략적으로 64비트 머신에서 한 k-mer가 2 bit 인코딩 되었다고 가정하면 (한 read 당 리드 길이 - k-mer 크기 + 1) * 64 만큼씩 처리해야할 k-mer의 정보양이 증가한다. 대략 1000~10000개 정도의 리드를 읽어서 이 정보양의 값을 계속 더한다. 이렇게 계산된 정보양에 파일 크기를 곱한 후에 이를 지금까지 읽은 리드의 전체 길이의 합으로 나눈다. 예제 코드에서는 추가적으로 (1024^2 * 8)로 나눠서 메가 바이트 단위가 되도록 했다. 이를 다시 가용한 디스크 공간으로 나누고 1을 더한다. 그러면 n_rounds 값이 된다. 이 n_rounds 값은 최종적인 정보양이 몇 번이나 read를 다시 읽고 k-mer를 추출해야 하는지 나타낸다.
그리고 amount_of_information을 다시 n_rounds로 나누게 되면, 현재 round에서 처리해야 하는 정보양이 나온다. 이를 round_volume이라고 하자.
다음으로 이 round_volume을 사용 가능한 메모리양(이것은 물리적으로 가용한 메모리를 의미하는 것이 아니라, 지놈 크기와 k-mer당 인코딩된 비트 수를 곱한 후에 이를 바이트로 환산함)으로 나눠주고 1을 더하면 파티션의 갯수가 나온다. 파티션은 최종적으로 k-mer들이 저장되는 가장 작은 단위의 블록이다.
uint32_t kmer_nbits = 64;
uint64_t amount_of_information = 0;
uint32_t counts = 0;
uint32_t size_of_local_sequence = 0;
while (readline(inputstream, line)) {
size_of_local_sequence += line.size();
amount_of_information += (line.size() - k + 1) * kmer_nbits;
if(counts == 1000) { break; }
}
amount_of_information = (amount_of_information * file_size) / size_of_local_sequence;
amount_of_information /= 1024 / 1024 / 8;
if (0 == amount_of_information) {
amount_of_information = 1;
}
n_rounds = (amount_of_information / available_disk_space(or size_of_sequence) + 1;
round_volume = amount_of_information / n_rounds;
n_partitions = round_volume / available_memory + 1;
말로 설명하면, 다음과 같이 된다. 한 round마다 전체 read를 다 읽는다. 그리고 새로운 round가 되면 기존에 있던 파티션은 모두 비워지게 된다. 각 read를 읽어 들이면서 k-mer와 그 reverse k-mer를 추출한 후에 작은 값을 가지는 k-mer 만을 저장하려고 시도한다. 이 때, 해당 k-mer에 대한 해쉬값을 계산하고, 이 계산된 해쉬값을 n_rounds와 모듈로 연산을 했을 때, 현재 round가 나오지 않으면 해당 k-mer는 무시한다. 반면, 현재 round 값이 나왔다면, 이 해쉬값을 다시 n_rounds로 나눈 k-mer를 narrowed_kmer라는 값에 저장한다. 이제 narrowed_kmer가 어떤 파티션에 들어가야하는 지 계산하기 위해서 narrowed_kmer를 n_partitions와 모듈로 연산을 수행한다. 최종적으로는 해당 파티션에 원래의 k-mer 값을 적어두게 된다.
for (uint32_t current_round = 0; current_round < n_rounds; ++current_round) {
for(Kmer a_kmer : kmers_in_a_read) {
uint64_t hash_value = a_kmer.hash_code();
if ((hash_value % n_rounds) != current_round) {
continue;
}
Kmer narrowed_kmer = hash_value / n_rounds;
uint32_t index_of_partition = narrowed_kmer % n_partitions;
partitions[index_of_partition].add(a_kmer);
}
// counting 하는 부분
}
이제 위의 코드에서 counting 하는 부분을 살펴보자. 이 부분에서는 해당 파티션이 저장된 파일의 k-mer들을 읽어서 개별 k-mer의 frequency를 더하다가 이 값이 cutoff를 넘을 경우 solid_kmers에 저장한다. 개별 파티션에는 이미 특정 조건에 해당하는 k-mer로 분류가 되어 있기 때문에, hash_table의 크기가 매우 커지지 않는다.
댓글