FASTQ 파일 다루기

 FASTQ 형식 소개에서 다뤘듯이 FASTQ 파일 형식은 서열id,서열,품질을 동시에 품고 있다. 이번 포스트에서는 ROSALIND 문제들을 통해 FASTQ 파일을 어떻게 다루는지 서술할 것이다.

 

https://rosalind.info/problems/phre/

ROSALIND | Read Quality Distribution

 

이 문제에서는 FASTQ 파일에서 여러 개의 read에서 품질(quality)를 추출하여 품질의 평균이 역치를 넘은 서열의 개수를 구하는 문제이다. 우리가 SeqIO.parse()를 통해서 파일을 읽을 때 우리는 SeqRecord 객체를 얻는다. 이 SeqRecord 객체는 여러가지 정보를 attribute로 저장하고 있다. 이를 우리는 부를 수 있다.

from Bio import SeqIO

with open(r'파일경로','r') as f:
	threshold = int(f.readline().rstrip())
    
    records = SeqIO.parse(f,'fastq')
    ans=0
    for record in records:
    	if sum(record.letter_annotations['phred_quality'])/len(record) >= threshold:
        	ans+=1
	print(ans)

이와 같이 우리는 record.letter_annotations['phred_quality']로 base 당 품질이 어떻게 되는지 부를 수 있다.

 

https://rosalind.info/problems/filt/

ROSALIND | Read Filtration by Quality

 

 이 문제는 여러개의 read중 주어진 역치를 넘은 base의 비율이 주어진 비율보다 높은 read의 개수를 세는 문제이다. 이 문제를 통해 우리는 base 하나씩 품질에 접근해서 품질이 낮은 base들은 쳐낼 수 있다.

from Bio import SeqIO

with open(r"파일경로",'r') as f:
    threshold,percentage = map(int,f.readline().rstrip().split())
    
    records = SeqIO.parse(f,'fastq')
    ans=0
    for record in records:
        total = len(record.seq)
        over=sum(i>=threshold for i in record.letter_annotations['phred_quality'])
        if over/total >= percentage/100:
            ans+=1
    
    print(ans)

 

https://rosalind.info/problems/bphr/

ROSALIND | Base Quality Distribution

 

이 문제는 반대로 여러 개의 길이가 같은 read에서 유독 품질이 떨어지는 특정 base position의 개수를 세는 문제이다. 

from Bio import SeqIO

with open(r'파일경로','r') as f:
    threshold = int(f.readline().rstrip())

    records = SeqIO.parse(f,'fastq')
    mean_q = next(records).letter_annotations['phred_quality']
    cnt=1
    for record in records:
        cnt+=1
        for idx,q in enumerate(record.letter_annotations['phred_quality']):
            mean_q[idx]+=q
    
    ans = sum(i/cnt < threshold for i in mean_q)
    print(ans)

 

https://rosalind.info/problems/bfil/

ROSALIND | Base Filtration by Quality

 

 이 문제는 주어진 read 중 품질이 떨어지는 양쪽 끝을 trimming 함으로써 그나마 쓸 수 있는 정보를 건져내는 문제이다.

from Bio import SeqIO

def trim(record,cut_off):
    quality_list = record.letter_annotations['phred_quality']
    front,back = -1,-1
    for idx,q in enumerate(quality_list):
        if q>=cut_off:
            front = idx
            break
    
    for idx,q in enumerate(quality_list[::-1]):
        if q>=cut_off:
            back = idx
            break
    
    return front,back

with open(r"파일경로",'r') as f:
    cut_off = int(f.readline().rstrip())

    records = SeqIO.parse(f,'fastq')
    ans_file = open(r"파일경로",'w')
    for record in records:
        front,back=trim(record,cut_off)
        new_record = record[front:len(record)-back]

        SeqIO.write(new_record,ans_file,'fastq')
    ans_file.close()

 

우리는 이와 같이 FASTQ 파일에서 품질에 접근하여 서열을 다룰 수 있다.