[cpp] FASTA 파일을 인덱스를 사용하여 읽는 방법

 

들어가며

안녕하세요 한주현입니다

 

오늘은 cpp로 fasta 파일을 읽는데 index를 사용하여 전체 파일을 읽지 않고,

미리 인덱스 파일을 읽어 해당 위치로 이동하여 읽을 수 있는 코드에 대해 말해보겠습니다.

 

FASTA 파일 설명

먼저 FASTA 파일과 FASTA 파일의 index 파일에 대해서 설명해보겠습니다.

FASTA 파일은 아래와 같이 생겼습니다.

 

> 로 시작하는 header 줄이 있고

 

그 다음은 염기서열이 있는 줄 들이 있습니다.

각각의 >header 와 서열들을 하나의 record라고 부릅니다.

FASTA 파일에는 하나의 record가 들어있을 수도 있고 여러개의 records가 있을 수도 있습니다.

예시에서는 세 개의 records가 있고 각각의 record는 한 줄의 염기 개수가 60자로 되어있습니다.

$ cat sample.fasta
>sample1
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>sample2
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>sample3
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT                                                                                                                                                

samtools를 사용하여 FASTA 파일의 indexing을 해보겠습니다.

$ samtools faidx sample.fasta

FASTA 파일이 indexing이 되었습니다.

이 파일을 살펴보면 아래 와 같이 생겼습니다.

$ cat sample.fasta.fai
sample1 80      9       60      61
sample2 84      100     60      61
sample3 87      195     60      61

첫 번째 컬럼은 header 줄에 있는 내용.

두 번째 컬럼은 각 record에 해당하는 염기서열의 숫자.

세 번째 컬럼은 각 record의 염기서열이 시작하는 offset.

네 번째 컬럼은 한 줄에 포함된 염기서열의 개수.

마지막 컬럼은 한 줄에 포함된 염기서열의 개수와 가장 끝의 new line의 개수까지 포함한 값.

입니다.

 

파일을 전체를 읽지 않고 원하는 부분만 읽을 수 있는 전략

예를 들어 >sample1 의 offset, 즉 record에서 염기서열이 시작하는 지점은 9 입니다.

그 이유는 >sample1\n 의 각각 character 개수를 세어본다면 9 글자 이기 때문입니다.

그래서 9번째 부터 읽게 되면 sample1의 record의 염기를 읽을 수 있습니다.

 

 

cpp 코드

이번 글에서 핵심적이라고 말할 수 있는 부분은 

read_fasta 함수에서 ifs.seekg(offset, std::ios::beg) 부분입니다.

 

offset만큼을 이동하여 읽는 부분으로

이렇게 하게 되면 일반적으로 human reference genome fasta 파일이 3GB 인데 이를 전부 읽지 않고도 특정 위치만 읽어낼 수 있습니다.

즉 메모리 효율적으로 읽을 수 있다는 점이죠.

 

오늘은 offset을 사용하여 FASTA 파일을 읽는 방법에 대해 알아보았습니다.

 

여러분들께 도움이 되셨으면 합니다.

 

그럼 다음 포스팅에서 또 만나요.