Skip to content

biopython FastQ #

Find similar titles

2회 업데이트 됨.

Edit
  • 작성자
    lca

You are seeing an old version of the page. Go to latest version

Structured data

Category
Programming

바이오파이썬과 FastQ #

대용량 NGS파일을 분석하기 위해서 R과 마찬가지로 Biopython에서는 관련 기능을 제공한다. 다만 R과 같이 다양한 기능을 제공하는 것이 아니라 최소 기능만 제공하고 있다.

바이오파이썬에서 FastQ파일은 읽기 위해서 기본적으로 SeqIO기능을 이용한다.

from Bio import SeqIO
count = 0
for rec in SeqIO.parse("SRR020192.fastq", "fastq"):
    count += 1
print("%i reads" % count)

위의 예제에서는 기본 서열의 입출력과 마찬가지로 SeqIO.parse(서열이름, 포맷)의 형태로 사용하면 된다. 위의 예제에서는 FastQ파일을 읽어들여서 서열의 갯수를 파악해보는 예제이다.

from Bio import SeqIO
good_reads = (rec for rec in \
          SeqIO.parse("SRR020192.fastq", "fastq") \
          if min(rec.letter_annotations["phred_quality"]) >= 20)
count = SeqIO.write(good_reads, "good_quality.fastq", "fastq")
print("Saved %i reads" % count)

위의 예제는 역시 서열의 정보를 읽은 후에 서열 정보 중 QC를 수행하는 예제이다. 즉 서열의 품질 정도를 파악해서 Q20 이상의 고퀄리티 정보를 추출하는 예제이다. letter_annotations["phred_quality"])에서 처럼 관련 정보 즉 phred_quality정보를 추출하여 이를 비교하고 별도의 파일로 저장하는 예제이다.

Q20은 99% 즉 1%의 에러만 허용하겠다는 뜻이다.

다음의 예제에서는 서열의 Adaptor서열을 trimming하기 위한 예제이다. Trimming을 하기 위해서 startswith("GATGACGGTGT") 즉 특정 서열로 시작하는 서열을 찾아내어 그 서열만 따로 추출해내는 예제이다.

from Bio import SeqIO
primer_reads = (rec for rec in \
            SeqIO.parse("SRR020192.fastq", "fastq") \
            if rec.seq.startswith("GATGACGGTGT"))
count = SeqIO.write(primer_reads, "with_primer.fastq", "fastq")
print("Saved %i reads" % count)

아래 예제에서는 위와 같이 startswith("GATGACGGTGT")로 시작하는 서열을 찾고 즉 서열 길이 10만큼 "GATGACGGTGT" 이후의 서열만을 추출한다. 즉 [11:]과 같이 서열 중 12번째 서열부터 끝까지 정보를 Slicing을 통해서 정보를 추출한다.그리고 이 서열들을 SeqIO.write()기능을 이용하여 별도의 파일로 저장하는 예제이다.

from Bio import SeqIO
trimmed_primer_reads = (rec[11:] for rec in \
                    SeqIO.parse("SRR020192.fastq", "fastq") \
                    if rec.seq.startswith("GATGACGGTGT"))
count = SeqIO.write(trimmed_primer_reads, "with_primer_trimmed.fastq", "fastq")
print("Saved %i reads" % count)

다음의 예제는 위에서 소개했던 내용을 함수 생성을 통해 재사용이 가능하게 만든 예제이다. def trim_primer() 함수를 생성하고 여기에 레코드 객체, 프라이머 서열을 인수로 넘겨주게 된다 그러면 프라이머 서열을 찾아서 그 서열이 실제 서열에 존재하는지 체크를 하게 되고 이를 별도의 파일로 저장하는 예제이다.

from Bio import SeqIO
def trim_primer(record, primer):
  if record.seq.startswith(primer):
    return record[len(primer):]
else:
    return record

trimmed_reads = (trim_primer(record, "GATGACGGTGT") for record in \
             SeqIO.parse("SRR020192.fastq", "fastq"))
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)

같은 기능을 하지만 조금 다르게 코딩을 해보면 다음과 같이 할 수도 있다. 여기서는 yield 기능을 사용하여 generator를 생성하고 이를 리턴하는 예제이다. 파이썬의 generator에 대해서는 별도로 다루고자 한다.

from Bio import SeqIO
def trim_primers(records, primer):
"""Removes perfect primer sequences at start of reads.

This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
    len_primer = len(primer) #cache this for later
    for record in records:
    if record.seq.startswith(primer):
        yield record[len_primer:]
    else:
        yield record

original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_primers(original_reads, "GATGACGGTGT")
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)

다음의 예제는 위와 동일하지만 조금 다른 형태로 작성된 예제이다. 특정 서열 즉 프라이머 서열을 찾기위해서 startswith() 함수를 사용하지 않고 find() 함수를 사용하였다. find() 함수는 특정 타겟 문자를 찾은 후에 그 찾은 위치를 리턴해준다. 즉 -1이라면 찾지 못한 것이 되고 양의 수를 리턴하면 그 찾은 위치가 된다. 즉 동일한 기능이지만 다른 형태로 작성된 예제이다.

from Bio import SeqIO
def trim_adaptors(records, adaptor):
"""Trims perfect adaptor sequences.

This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
    len_adaptor = len(adaptor) #cache this for later
    for record in records:
      index = record.seq.find(adaptor)
      if index == -1:
        #adaptor not found, so won't trim
        yield record
     else:
        #trim off the adaptor
        yield record[index+len_adaptor:]

original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_adaptors(original_reads, "GATGACGGTGT")
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)

아래 예제는 서열 중에 프라이머 서열을 찾고 이 중에서 서열의 길이가 100미만인 것들만 찾아서 따로 파일로 정리하는 내용이다. 특정 서열의 길이가 너무 짧다면 문제가 되기 때문에 해당 서열의 길이를 구한 후 이를 최소 길이와 비교하여 데이터를 정리한다.

from Bio import SeqIO

def trim_adaptors(records, adaptor, min_len):
"""Trims perfect adaptor sequences, checks read length.

This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
    len_adaptor = len(adaptor) #cache this for later
    for record in records:
       len_record = len(record) #cache this for later
       if len(record) < min_len:
          #Too short to keep
          continue
        index = record.seq.find(adaptor)
        if index == -1:
          #adaptor not found, so won't trim
          yield record
        elif len_record - index - len_adaptor >= min_len:
          #after trimming this will still be long enough
          yield record[index+len_adaptor:]

original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_adaptors(original_reads, "GATGACGGTGT", 100)
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)

다음은 각 장비별 데이터를 표준 FastQ로 바꿔주는 예제이다. 즉 해당 머신에서 생산된 데이터가 조금씩 다른 형태를 가지고 있기 때문에 표준 포맷으로 변환한다. SeqIO.convert(입력, 포맷, 출력, 포맷)의 형태로 사용이 되는데 입력된 파일과 포맷을 지정하고 저장될 파일과 포맷을 지정해주면 자동으로 변환한다.

from Bio import SeqIO
SeqIO.convert("solexa.fastq", "fastq-solexa", "standard.fastq", "fastq")


from Bio import SeqIO
SeqIO.convert("illumina.fastq", "fastq-illumina", "standard.fastq", "fastq")

다음 예제는 fastq파일에서 서열만 추출하여 별도의 fasta포맷으로 전환하는 예제이다. 그 아래는 qual파일 즉 퀄리티 파일만 따로 추출하는 예제이다 from Bio import SeqIO SeqIO.convert("example.fastq", "fastq", "example.fasta", "fasta")

from Bio import SeqIO
SeqIO.convert("example.fastq", "fastq", "example.qual", "qual")

간혹 fastq파일이 qual과 fasta로 구분되어 가지고 있는 경우가 있는 분석 등의 이유로 인해서 두 파일을 하나도 합칠 필요도 있다. 이를 위해 PairedFastaQualIterator를 통해서 서열 파일과 퀄리티 파일을 하나도 합쳐서 사용할 수도 있다.

from Bio.SeqIO.QualityIO import PairedFastaQualIterator
for record in PairedFastaQualIterator(open("example.fasta"), open("example.qual")):
print(record)

아래의 예제에서는 동작하는 다른 예제이다. 즉 두 파일, 퀄리티 파일과 fasta파일을 하나로 합쳐 새로운 fastq파일 포맷으로 저정하는 예제이다. 경우에 따라서 상황에 맞게끔 사용하면된다.

from Bio import SeqIO
from Bio.SeqIO.QualityIO import PairedFastaQualIterator
handle = open("temp.fastq", "w") #w=write
records = PairedFastaQualIterator(open("example.fasta"), open("example.qual"))
count = SeqIO.write(records, handle, "fastq")
handle.close()
print("Converted %i records" % count)
0.0.1_20240318_1_v95