Biopython FastQ #

Find similar titles

총 5회 업데이트 됨.

Edit

최초 작성자
yeye 2016-10-18 04:03:03
최근 업데이트
cloudgo 2017-01-06 11:14:24

Structured data

Category: Programming

바이오파이썬과 FastQ #

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

Biopython에서 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)

Structured data

바이오파이썬과 FastQ #

Suggested Pages #