서울여자대학교 파이썬교육 #
Find similar titles

개요
교육 내용

개요 #

주제: 생물정보 분석을 위한 파이썬 기초
일시: 2017-02-14 ~ 2017-02-16
장소: 서울여자대학교 인문사회관 (211호)
교육 특징
- 강의와 실습을 병행함. 각 단위 강의 후에 연습문제 풀이
- 문제 풀이, 참고자료 등 교육 내용의 실시간 인터넷 공유 (본 웹 페이지)
교육 진행방식
- 실습시 개인 PC에 Anaconda 설치 (Anaconda_jupyter-설치방법.pdf)
- 문서편집기는 IDLE 사용

교육 내용 #

0일차 #

강의 교재

1일차 #

강의 교재

명령행 인터페이스 작성하기 예제

#!python
import argparse
import sys
from Bio import SeqIO
from Bio.Alphabet import IUPAC

parser = argparse.ArgumentParser()
parser.add_argument('-t', '--type', type=str,
                    choices=('DNA', 'RNA', 'Protein'),
                    help="input sequence type")
parser.add_argument('-c', '--command', type=str,
                    choices=('reverse-complement', 'transcribe'),
                    help="FASTA manipulation command")
parser.add_argument('-f', '--fields', type=int,
                    default=60,
                    help='FASTA column size')
parser.add_argument('-i', '--input', type=argparse.FileType('r'),
                    default=sys.stdin,
                    help="Input file")
parser.add_argument('-o', '--output', type=argparse.FileType('o'),
                    default=sys.stdout,
                    help="Output filename")

args = parser.parse_args()

for record in SeqIO.parse(args.input, 'fasta'):
    seq = record.seq
    if args.type == 'DNA':
        seq.alphabet = IUPAC.unambiguous_dna
    elif args.type == 'RNA':
        seq.alphabet = IUPAC.unambiguous_rna
    elif args.type == 'Protein':
        seq.alphabet = IUPAC.protein

    if args.command == 'transcribe':
        record.seq = seq.transcribe()
    elif args.command == 'reverse-complement':
        record.seq = record.seq.reverse_complement()
    print(record.format('fasta'))

2일차 #

강의 내용

지난시간 복습 연습문제
BLAST
Clustal Omega
Entrez with Biopython
연습문제

강의 교재

유용한 코드 #

BLAST 결과 evalue로 필터링

#!python
import sys
from Bio.Blast import NCBIStandalone

user_evalue = float(sys.argv[1])

blast_parser = NCBIStandalone.BlastParser()
for record in NCBIStandalone.Iterator(sys.stdin, blast_parser):
    print('query {}'.format(record.query))
    for alignment in record.alignments:
        if  alignment.hsps[0].expect < float(user_evalue):
            print('   matched sequence {} [expect: {}, hsp: {}]'.format(
            alignment.title, alignment.hsps[0].expect, len(alignment.hsps)))

Entrez로 두개의 nucleotide 데이터를 받아 GC 함량 비교하기

#!python
from Bio import Entrez, SeqIO
from Bio.SeqUtils import GC
import sys
Entrez.email = "okna0215@naver.com"

def get_gc_by_id(id):
    handle = Entrez.efetch(db='nucleotide', id=id, rettype='gb',
                             retmode='text')
    record = SeqIO.read(handle, 'genbank')
    return GC(record.seq)

id1 = sys.argv[1]
id2 = sys.argv[2]

gc1 = get_gc_by_id(id1)
gc2 = get_gc_by_id(id2)

print("{}: {:.2f}".format(id1, gc1))
print("{}: {:.2f}".format(id2, gc2))

if gc1 > gc2:
    print("Result: {} > {}".format(id1, id2))
elif gc1 < gc2:
    print("Result: {} < {}".format(id1, id2))
else:
    print("Result: {} = {}".format(id1, id2))

3일차 #

실습준비

Anaconda설치
conda install bcbio-gff

원하는 위치의 유전자 찾기

#!python
import pandas as pd
import sys


def get_note(raw):
    return raw.attrs.split(";")[2]

def find_gene(chrno, pos):  
    df = pd.read_table("rice_locus.gff", sep="\t",
                  names =['chr','source','type','start',
                          'end','score','strand','phase',
                          'attrs'])
    df["note"] = df.apply(get_note, axis=1)

    aa = df[(df.chr == chrno) & (df.start <= pos) & (df.end >= pos)]
    return aa.note


chrno = sys.argv[1]
pos = int(sys.argv[2])
print(find_gene(chrno, pos))

subplots()을 이용해서 그래프 그리기

#!python
import matplotlib.pyplot as plt
%matplotlib inline

data = pd.read_csv('weight.csv')

# figure, subplot 생성
fig, axes = plt.subplots(3,1, figsize=(5,20))

# 국가별 체충의 boxplot 그리기
data.boxplot(column='Weight', by='Country', ax=axes[0])

#  체중의 histogram 그리기
data['Weight'].hist(ax=axes[1], color='gray', bins=20, range=(0,120))

# 성별로 pie차트 그리기
data['Sex'].value_counts().plot(ax=axes[2], kind='pie')

plt.subplots_adjust(hspace=0.5)

#ex5.png에 저장
plt.savefig(‘ex5.png’)

모티프 서열 위치 찾기

#!python
import sys
import urllib.request
import re

filename = []
# for line in sys.stdin:
for line in open("data/prot.txt"):
    name = line.strip()
    print(name)

    data = urllib.request.urlopen('http://www.uniprot.org/uniprot/'+name+'.fasta')
    content = data.read().decode("utf-8")
    name, *seqs = content.split('\n')
    seq = ''.join(seqs)

    pattern = re.compile(r'N[^P][ST][^P]')
    matchs = re.finditer(pattern, seq)
    for m in matchs:
        print(m.start()+1, end=" ")
    print("\n")

예제데이터

data.zip

강의교재

Table of Contents

개요 #

교육 내용 #

0일차 #

1일차 #

2일차 #

유용한 코드 #

3일차 #