[생물정보학] FASTA gzip 파일 읽기 - 염색체 별로 염기 숫자 세기 - UCSC FASTA 다운로드 - 자바 JAVA , 파이썬 Python, 파일 읽기 속도 비교

 

안녕하세요 한주현입니다.

오늘은 자바, 파이썬으로 FASTA gzip 파일을 읽어 염기서열을 세는 프로그래밍을 해보겠습니다.

목차

1. UCSC FASTA 파일 다운로드

2. 염색체 별로 읽기

3. Python 스크립트 구현

4. JAVA 코드 구현

5. 결과

6. 파일 읽기 속도 비교 - JAVA, Python

1. UCSC FASTA 파일 다운로드 - chromosome 별

다음 경로에 들어가서 chromosome 별로 fa.gz 파일을 받습니다.

hg38
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/

hg19

http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/

사이트에 접속하시면 다음과 같은 페이지가 나오는데,

각 reference 별로

chr1.fa.gz

chr2.fa.gz

...

chr22.fa.gz

chrX.fa.gz

chrY.fa.gz 

chrM.fa.gz

의 파일들을 다운 받습니다.

리눅스 환경에서는 다음과 같이 진행해주시면 chrmosome 별로 차례대로 받으실 수 있습니다.

#!/bin/bash
 
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 X Y M
do
  wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/chr${i}.fa.gz
done

2. 염색체 별로 읽기

우리는 chromosome 별로 나뉜 파일을 하나씩 읽어 각 chromosome에 들어있는 염기의 개수를 셀 것입니다. 

염기는 A, C, G, T 그리고 N 으로 구성되어있는데, 우리는 ACGTN 구분없이 모두 읽어보겠습니다.

FASTA 파일의 형식은 > 기호로 시작하는 헤더와 그 아래 ACGTN 의 염기서열 입니다.

FASTA 파일 형식

> 헤더

ACACACGGCCNNNNNTTTTCC

CCACGGTTNNNNCCCAAAAAA

...

원하는 결과는 다음과 같이 chromosome과 염기개수로 구성된 두 개의 열이지요.

chr1   249250621

chr2   243199373

...

chrM  16571

3. Python 스크립트 구현

한 줄씩 읽으며 > 기호는 거르고 count 합니다.

#!/usr/bin/python
 
import sys
import gzip
 
class FastaReader:
    def __init__(self):
        self.fastaFile = ""
        self.count = 0
 
    def readFasta(self):
        with gzip.open(self.fastaFile,'rb') as fr:
            for line in fr:
                s = line.decode("utf-8").strip()
                if not s.startswith(">"):
                    self.count += len(s)
 
if __name__ == "__main__":
    if len(sys.argv) != 2:
        print("#usage: python %s <fasta.gz>" %sys.argv[0])
        sys.exit()
    fastaFile = sys.argv[1]
    fastaReader = FastaReader()
    fastaReader.fastaFile = fastaFile
    fastaReader.readFasta()
    print(fastaReader.count)

4. JAVA 코드 구현

한 줄씩 읽으며 > 기호는 거르고 count 합니다.

저는 static으로 해서 호출하는 방법으로 만들었습니다.

//FastaReader.java
import java.io.*;
import java.util.zip.*;
class FastaReader {
    private static long baseCount = 0;
    private static String s;
    public static long readFasta(String fastaFile) {
        try {
            BufferedReader in = new BufferedReader(new InputStreamReader(new GZIPInputStream(new FileInputStream(fastaFile))));
            while((s = in.readLine())!=null) {
                if (!(s.startsWith(">"))) {
                    baseCount += s.trim().length();
                }
            }
        } catch (IOException ex) {
            System.out.println(ex);
        }
        return baseCount;
    }
}

//FastaReaderTester.java
 
import java.io.*;
 
class FastaReaderTester {
    public static void main(String[] args) {
        if (args.length != 1) {
            System.out.println("java FastaReaderTester <fasta.gz>");
            System.exit(0);
        }
        String fastaFile = args[0];
        long baseCnt = FastaReader.readFasta(fastaFile);
        System.out.println(baseCnt);
    }
}
 

5. 결과

$ cat hg19.count.result.txt 
chr1    249250621
chr2    243199373
chr3    198022430
chr4    191154276
chr5    180915260
chr6    171115067
chr7    159138663
chr8    146364022
chr9    141213431
chr10    135534747
chr11    135006516
chr12    133851895
chr13    115169878
chr14    107349540
chr15    102531392
chr16    90354753
chr17    81195210
chr18    78077248
chr19    59128983
chr20    63025520
chr21    48129895
chr22    51304566
chrX    155270560
chrY    59373566
chrM    16571

$ cat hg38.count.result.txt 
chr1    248956422
chr2    242193529
chr3    198295559
chr4    190214555
chr5    181538259
chr6    170805979
chr7    159345973
chr8    145138636
chr9    138394717
chr10    133797422
chr11    135086622
chr12    133275309
chr13    114364328
chr14    107043718
chr15    101991189
chr16    90338345
chr17    83257441
chr18    80373285
chr19    58617616
chr20    64444167
chr21    46709983
chr22    50818468
chrX    156040895
chrY    57227415
chrM    16569

6. 파일 읽기 속도 비교 - JAVA, Python

chr21.fa.gz 파일로 각각 프로그램을 실행해보았습니다.

JAVA: 2.8초

Python3: 17.2초

가 나왔습니다.

plot 그릴 때 사용한 python 스크립트

import pandas as pd
import matplotlib.pyplot as plt
 
plt.style.use("ggplot")
 
%matplotlib inline 
 
data = {'JAVA': [2.8], 'Python3': [17.2]}
df = pd.DataFrame.from_dict(data)
 
ax = df.plot.bar(title='Read fasta.gz file')
ax.set(xlabel='Program', ylabel='Seconds')
plt.show()

JAVA, Python 소스코드는 다음 github 페이지에서 받으실 수 있습니다.

https://github.com/KennethJHan/FastaReader

오늘은 UCSC에서 chromosome 별 FASTA 파일을 받고 염색체 별로 염기서열을 세는 방법에 대해 알아보았습니다.

그럼 다음에 만나요~~ 

기부 버튼을 만들었습니다

단지 $1 의 작은 정성도 저에게는 큰 힘이 됩니다

기부해주신 분들을 기억하며

더 좋은 내용으로 보답해 드리겠습니다 :)

Donate 버튼은 paypal 결제로 paypal 계정이 없으시더라도

카드로도 기부 가능하십니다 :)

Use your credit card or bank account (where available). 옆의 continue 를 누르시면 됩니다

한주현 드림