Python으로 NGS 분석하기 위한 프로그래밍 실무 교육

개요 #

목표

• Python 프로그래밍을 통한 농생명 생명정보 분석 실무 역량 강화
• 기본 문법 최소화, 실무 응용 집중
• 프로그래밍 과정 시연과 직접 실습, 과제 제출을 통해 확실히 내것으로!

교육환경

• Soongsil University 벤처관 707호 실습실 Ubuntu linux
• Python 2.7.8
• Editor: IDLE or vim

필수 사전 온라인 학습 (반드시 전체 내용을 실습할 것)

1. Free Interactive Python Tutorial "Learn the Basics"

선택 사전 온라인 학습 (가능한 실습하기를 권장함)

1. Free Interactive Python Tutorial "Advanced Tutorials"
2. 왕초보를 위한 파이썬 2.7 1장 ~ 8장

교육 과정 #

0주차: 파이썬 기본 문법 이해 #

강의자료: W0-Python-Basic.pdf

• Python introduction
• Data type
  • Numeric type
  • Sequence type
  • Collection type
  • Boolean type
• Control flow
  • if statement
  • for statement
  • while statement
  • Sort
• Function
• Mudule, Package
• String formatting, File, Unicode
• Exception
• Unit test

문제풀이 #

Q1. 다음을 출력하는 파이썬 코드를 작성하시오.

$ python p1.py
*
**
***
****
*****

range 함수의 결과를 for 문으로 루프.

#!python
for each in range(1, 6):
    print '*' * each

Q2. Q1을 약간 변형함

$ python p1.py
*****
 ****
  ***
   **
    *

앞 부분 공백도 함께 변수로 처리

#!python
for each in range(5):
    print ' ' * each + '*' * (5 - each)

Q3. Euclidean distance 계산하기

emumerate를 쓰면 해당 리스트의 요소 위치(인덱스)가 함께 표시됨.

#!python
import math

def get_euclidean_distance(lista, listb):
    total = 0
    for i, each in enumerate(lista):
        total += (each - listb[i]) ** 2
    return math.sqrt(total)

Q4. 2014년이 윤년인지 아닌지 확인 (calendar 모듈 사용)

https://docs.python.org/2/library/calendar.html 문서를 통해 사용방법 익힘

#!python
>>> import calendar
>>> calendar.isleap(2014)
False
>>>

Q5. 2016년 7월 29일은 무슨요일? (calendar 모듈 사용)

#!python
>>> import calendar
>>> calendar.weekday(2016, 7, 29)
4
>>>

위 모듈 문서에 의하면 0이 월요일이다. 4는 금요일.

Q6. 오늘부터 1000일 후가 며칠인지 알려주는 함수 만들기.

#!python
>>> import datetime
>>> def date_from_today(days):
...     return datetime.date.today() + datetime.timedelta(1000)
...
>>> date_from_today(1000)
datetime.date(2017, 5, 1)
>>>

Q7. 임의의 텍스트 파일내 단어의 빈도를 조사하여 가장 많이 출현한 단어부터 정렬하여 출력하는 프로그램 (특수문자 제외, 소문자로 통일)

#!python
import sys

word_frequencies = {}
omits = ['.', ',', "'", ":"]

for line in sys.stdin:
    for word in line.strip().split():
        word = word.lower()
        for omit in omits:
            if omit in word:
                word = word.replace(omit, '')
        if word in word_frequencies:
            word_frequencies[word] += 1
        else:
            word_frequencies[word] = 1

for word, frequency in sorted(word_frequencies.items(),
        key=lambda x: x[1], reverse=True)[:10]:
    print word, frequency

Q8. 구구단 출력하기 (임의의 단을 명령행으로 입력받음)

$ python gugudan.py 4
4 x 1 = 4
4 x 2 = 8
4 x 3 = 12
4 x 4 = 16
4 x 5 = 20
4 x 6 = 24
4 x 7 = 28
4 x 8 = 32

문자열 포매팅 (.format)을 사용하여 결과 표시.

#!python
import sys

x = int(sys.argv[1])
for y in range(1,10):
    print '{} x {} = {}'.format(x, y, x * y)

1주차: 파이썬 객체지향 프로그래밍 이해 - FASTA 다루기 #

강의자료: W1-Python-OOP.pdf

• Class
• Decorator
• Iterator
• Generator

문제풀이 #

Q1. FASTA 파일내 염기서열의 GC content를 객체지향과 단위테스트로 구현하기 (Single FASTA)

#!python
import unittest
import sys

class Fasta:
    def __init__(self, title, sequence):
        self.title = title
        self.sequence = sequence

    def get_gc_content(self):
        gc = self.sequence.count('G') + self.sequence.count('C')
        return gc / float(len(self.sequence))


class FastaTest(unittest.TestCase):
    def test_fasta_object(self):
        my_fasta = Fasta('test1', 'AGTC')
        self.assertEqual(my_fasta.title, 'test1')
        self.assertEqual(my_fasta.sequence, 'AGTC')

    def test_gc_content(self):
        my_fasta = Fasta('test1', 'AGTC')
        self.assertEqual(my_fasta.get_gc_content(), 0.5)

    def test_reverse_complement(self):
        my_fasta = Fasta('test1', 'AGTC')
        self.assertEqual(my_fasta.reverse_complement(), 'GACT')

def main():
    sequence = []
    for line in sys.stdin:
        line = line.strip()
        if line.startswith('>'):
            title = line[1:]
        else:
            sequence.append(line)
    sequence = ''.join(sequence)
    my_fasta = Fasta(title, sequence)
    print my_fasta.get_gc_content()

if __name__ == '__main__':
    unittest.main()
    #main()

Q2. FASTA 입력파일의 Reverse complement 출력하기 (Multiple FASTA)

#!python
import unittest
from cStringIO import StringIO

class Fasta:
    def __init__(self, title, sequence):
        self.title = title
        self.sequence = sequence

    def get_gc_content(self):
        gc = self.sequence.count('G') + self.sequence.count('C')
        return gc / float(len(self.sequence))

    def reverse_complement(self):
        result = []
        for nucleotide in self.sequence:
            result.append(self.base_pairs[nucleotide])
        return ''.join(reversed(result))

    def get_rc_fasta(self):
        return Fasta(self.title + ' reverse complement', 
                self.reverse_complement())

    def write(self, cols=60):
        result = []
        result.append('>{}'.format(self.title))
        sequence = self.sequence
        while sequence:
            result.append(sequence[:cols])
            sequence = sequence[cols:]
        return '\n'.join(result) + '\n'

    @classmethod
    def iterator(cls, file):
        title = file.next().strip()[1:]
        sequence = []
        for line in file:
            line = line.strip()
            if line.startswith('>'):
                sequence = ''.join(sequence)
                fasta = cls(title, sequence)
                yield fasta
                title = line[1:]
                sequence = []
            sequence.append(line)
        sequence = ''.join(sequence[1:])
        fasta = cls(title, sequence)
        yield fasta


class DNA(Fasta):
    base_pairs = {
        'A': 'T',
        'T': 'A',
        'G': 'C',
        'C': 'G',
    }

class RNA(Fasta):
    base_pairs = {
        'A': 'U',
        'U': 'A',
        'G': 'C',
        'C': 'G',
    }


class FastaTest(unittest.TestCase):
    def test_fasta_object(self):
        my_fasta = Fasta('test1', 'AGTC')
        self.assertEqual(my_fasta.title, 'test1')
        self.assertEqual(my_fasta.sequence, 'AGTC')

    def test_gc_content(self):
        my_dna = DNA('test1', 'AGTC')
        self.assertEqual(my_dna.get_gc_content(), 0.5)

    def test_reverse_complement(self):
        my_dna = DNA('test1', 'AGTC')
        self.assertEqual(my_dna.reverse_complement(), 'GACT')
        my_rna = RNA('test1', 'AGUC')
        self.assertEqual(my_rna.reverse_complement(), 'GACU')

    def test_fasta_iterate(self):
        input = StringIO("""\
>test1
AGTC
AGTC
>test2
AAGG
GGTT
""")
        fi = Fasta.iterator(input)
        f1 = fi.next()
        self.assertEqual(f1.title, 'test1')
        self.assertEqual(f1.sequence, 'AGTCAGTC')
        f2 = fi.next()
        self.assertEqual(f2.title, 'test2')
        self.assertEqual(f2.sequence, 'AAGGGGTT')

    def test_fasta_write(self):
        my_fasta = Fasta('test1', 'AGTCAGTC')
        self.assertEqual(my_fasta.write(cols=4), '''\
>test1
AGTC
AGTC
''')


def main(infile, outfile):
    for dna in DNA.iterator(infile):
        outfile.write(dna.get_rc_fasta().write(cols=20))

if __name__ == '__main__':
    #unittest.main()
    import sys
    main(sys.stdin, sys.stdout)

Q3. Q2 코드를 이용하여, Multiple FASTA의 gc content를 다음처럼 구현하기

$ python gc_content.py < input.fasta
>test1 GC content
0.55
>test2 GC content
0.45

위 Q1의 코드 main 함수를 다음처럼 변경

#!python
def main(infile, outfile):
    for dna in DNA.iterator(infile):
        outfile.write('>{}\n'.format(dna.title))
        outfile.write('{}\n'.format(dna.get_gc_content()))

Q4. Q2 코드를 이용하여, Multiple FASTA의 Transcripted RNA를 다음처럼 구현하기

$ python transcribe.py < input.fasta
>test1 transcripted RNA
AGUAUAGCC
AGUUAACAG
>test2 transcripted RNA
AUAGCAGCU
ACGAGGACU

위 Q1의 DNA 클래스를 다음 처럼 변경 후, main 함수 변경

#!python
class DNA(Fasta):
    base_pairs = {
        'A': 'T',
        'T': 'A',
        'G': 'C',
        'C': 'G',
    }

    def transcribe(self):
        return RNA('{} transcribed'. format(self.title),
                self.sequence.replace('U', 'T'))


def main(infile, outfile):
    for dna in DNA.iterator(infile):
        rna = dna.transcribe()
        outfile.write(rna.write(cols=20))

2주차: Biopython을 이용한 생물정보 입출력 #

강의자료: W2-Python-Biopython.pdf

• FASTA manipulation
• Standard libraries
• 3rd party libraries
• Biopython
  • Seq
  • SeqRecord
  • SeqFeature
  • SeqIO
  • Entrez

문제풀이 #

Q1. 1주차 Q2 코드를 변경하여, 명령행에서 '-c' 옵션으로 명령(reverse-complement, gc-content, transcribe)를 넣을 수 있도록 구현. (-i 옵션으로 입력파일, -o 옵션으로 출력파일을 선택 입력할 수 있도록)

1주차 Q2 코드의 main 함수를 다음처럼 변경

#!python
def main():
    import sys
    import argparse
    parser = argparse.ArgumentParser(description='FASTA manipulator')
    parser.add_argument('-c', '--command', type=str, 
            help='Manuplation command', 
            choices=('reverse-complement', 'gc-content', 'transcribe'), 
            required=True)
    parser.add_argument('-i', '--infile', type=argparse.FileType('r'),
            default=sys.stdin,
            help='Input file, if not specified, stdin used')
    parser.add_argument('-o', '--outfile', type=argparse.FileType('w'),
            default=sys.stdout,
            help='Output file, if not specified, stdout used')
    args = parser.parse_args()

    for dna in DNA.iterator(args.infile):
        if args.command == 'reverse-complement':
            rc_dna = dna.get_rc_fasta()
            args.outfile.write(rc_dna.write(cols=20))
        elif args.command == 'gc-content':
            args.outfile.write('>{} GC content\n'.format(dna.title))
            args.outfile.write('{:f}\n'.format(dna.get_gc_content()))
        elif args.command == 'transcribe':
            rna = dna.transcribe()
            args.outfile.write(rna.write(cols=20))

Q2. Biopython을 써서 Multiple FASTA 파일의 reverse complemented 출력 프로그램 구현

#!python
import sys
from Bio import SeqIO

for record in SeqIO.parse(sys.stdin, 'fasta'):
    record.seq = record.seq.reverse_complement()
    record.description = record.description + ' reverse complemented'
    sys.stdout.write(record.format('fasta'))

Q3. 특정 키워드를 명령행 인수로 입력하면, 해당 키워드의 2014년 논문들의 초록을 pubmed_id.txt 파일로 기록하는 프로그램 구현

#!python
import sys
from Bio import Entrez
from Bio import Medline

keyword = sys.argv[1].strip()
handle = Entrez.esearch(term=keyword, db='pubmed')
record = Entrez.read(handle)
idlist = record['IdList']

handle = Entrez.efetch(id=idlist, db='pubmed', rettype='medline', retmode='text')
for record in Medline.parse(handle):
    year = int(record['DP'].split()[0])
    if year == 2014:
        abstract = record['AB']
        pubmed_id = record['PMID']

        with open('{}.txt'.format(pubmed_id), 'w') as outfile:
            otufile.write(abstract)

3주차: 외부 프로세스 다루기 - BLAST 자동화 #

강의자료: W3-Python-subprocess.pdf

• blast+
• BLAST result parsing using biopython
• Python subprocess

문제풀이 #

Q1. 임의의 BLAST 결과(pairwise)를 표준입력으로 입력하면 다음과 같은 BLAST 결과 요약정보를 출력하는 프로그램 구현

$ python blast_summary.py 0.0005 < my_query.blout.pairwise
query title1
    matched sequence title1 [e-value: 0.005, hsps:4]
    matched sequence title2 [e-value: 0.009, hsps:4]
query title2
    matched sequence title3 [e-value: 0.001, hsps:1]
query title3
    matched sequence title4 [e-value: 0.002, hsps:2]

pairwise 출력이므로, NCBIStandalone 모듈을 이용하여 BLAST 결과를 파싱해야 함.

#!python
import sys
from Bio.Blast import NCBIStandalone

blast_parser = NCBIStandalone.BlastParser()
for record in NCBIStandalone.Iterator(sys.stdin, blast_parser):
    print '{}'.format(record.query)
    for alignment in record.alignments:
        print '    {} [e-value: {:f}, hsps: {:d}]'.format(
                alignment.title, 
                alignment.hsps[0].expect, 
                len(alignment.hsps))

Q2. NCBI에서 특정 키워드로 검색한 nucleotide 정보들을 (상위 10개만) 각각 blastx로 swissprot DB에 검색했을 때 나오는 단백질의 개수는 총 몇 개인지 계산해주는 스크립트를 작성하시오. (E-value 10^-5 미만, 단백질 개수를 중복계산하지 않도록 유의할 것)

$ python search.py apoptosis
123

명령행 옵션으로 입력된 키워드로 Entrez nucleotide DB를 검색하고, 검색결과를 genbank 형식으로 받은 뒤, 이를 SeqIO 모듈로 하나하나 FASTA 로 변환하고, blastx 프로그램의 입력으로 이용한다. blastx 결과는 NCBIStandalone 모듈을 이용하여 파싱한 뒤, 매칭 단백질 갯수 (count alignments)를 합산 한 뒤 결과를 출력한다.

#!python
import sys
from Bio import Entrez
from Bio import SeqIO
from Bio.Blast import NCBIStandalone
from subprocess import Popen, PIPE
from cStringIO import StringIO

## Search esearch by keyword
Entrez.email = 'yong27@gmail.com'
keyword = sys.argv[1]
handle = Entrez.esearch(db='nucleotide', term=keyword)
result = Entrez.read(handle)
idlist = result['IdList'][:10]   ## search result top 10

blast_parser = NCBIStandalone.BlastParser()
count_matched_protein = 0

## Fetch each result
handle = Entrez.efetch(db='nucleotide', id=idlist, 
        rettype='gb', retmode='text')
for seq_record in SeqIO.parse(handle, 'gb'):
    print seq_record.name
    if len(seq_record.seq) > 100000:  ## pass big sequence 
        continue
    input_fasta = seq_record.format('fasta')

    p = Popen(['blastx', '-db', 'swissprot', '-outfmt', '0'],
            stdin=PIPE, stdout=PIPE)
    blast_output = p.communicate(input=input_fasta)[0]
    blast_record = blast_parser.parse(StringIO(blast_output))
    print '    number of hit: {}'.format(len(blast_record.alignments))
    for alignment in blast_record.alignments:
        if alignment.hsps[0] < 0.0005:
            count_matched_protein += 1

print count_matched_protein

4주차: NGS 데이터 다루기 #

강의자료: W4-Python-subprocess.pdf

• NGS data analysis
• FASTQ
• BAM
  • SAMtools
  • pysam
• VCF
  • VCFtools
  • PyVCF
• GFF3
  • bcbio-gff
  • bedtools
  • pybedtools
• Galaxy tool

후기 #

Biopython이 Python 3을 지원한다는 사실을 나중에 알게됨. 다음 기회에는 Python 3으로 시작할 수 있도록 준비해야 함.

Suggested Pages #

• 0.025 Hyungyong Kim/.vimrc
• 0.025 C
• 0.025 Design patterns
• 0.025 GTF
• 0.025 Python for Data Analysis
• 0.025 TinyDB
• 0.025 C++
• 0.025 DNA sequencing
• 0.025 Insilicogen, Inc.
• 0.025 Nucleic Acids Research
• More suggestions...

Other Posts #

• newer PyCon Korea 2014 “DNA 데이터로 가족찾기” 발표 후기
• older 개인유전체 검사 결과를 설명하는 방법