바이오파이썬 예제코드

Taeyoon Kim

2019-06-06 08:21

아래 내용은 책 한주현, 바이오파이썬으로 만나는 생물정보학, 비제이퍼블릭, 2019을 읽고 실습한 코드를 정리한 것입니다.

In [1]:

import Bio

print(Bio.__doc__)

Collection of modules for dealing with biological data in Python.

The Biopython Project is an international association of developers
of freely available Python tools for computational molecular biology.

http://biopython.org

DNA 서열 다루기¶

In [2]:

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
with Entrez.efetch(
    db="nucleotide", rettype="fasta", retmode="text", id="1490011893"
) as handle:
    seq_record = SeqIO.read(handle, "fasta")

seq_record

Out[2]:

SeqRecord(seq=Seq('ATGGAGGAGATGCTGCCCCTCTTTGAGCCCAAGGGCCGGGTCCTCCTGGTGGAC...GAG', SingleLetterAlphabet()), id='MG727867.1', name='MG727867.1', description='MG727867.1 Synthetic construct Taq DNA polymerase gene, partial cds', dbxrefs=[])

서열 개수 세기¶

In [3]:

taq = seq_record.seq
len(taq)

Out[3]:

GC_content 계산하기¶

In [4]:

from Bio.SeqUtils import GC

GC(taq)

Out[4]:

67.70833333333333

대소문자 변환하기¶

In [5]:

taq.upper()

Out[5]:

Seq('ATGGAGGAGATGCTGCCCCTCTTTGAGCCCAAGGGCCGGGTCCTCCTGGTGGAC...GAG', SingleLetterAlphabet())

In [6]:

taq.lower()

Out[6]:

Seq('atggaggagatgctgcccctctttgagcccaagggccgggtcctcctggtggac...gag', SingleLetterAlphabet())

DNA 서열 Translation, Transcription¶

In [7]:

taq.transcribe()

Out[7]:

Seq('AUGGAGGAGAUGCUGCCCCUCUUUGAGCCCAAGGGCCGGGUCCUCCUGGUGGAC...GAG', RNAAlphabet())

In [8]:

taq.translate()

Out[8]:

Seq('MEEMLPLFEPKGRVLLVDGHHLAYRTFHALKGLTTSRGEPVQAVYGFAKSLLKA...AKE', ExtendedIUPACProtein())

Molecular weight 구하기¶

In [9]:

from Bio.SeqUtils import molecular_weight

molecular_weight(taq)

Out[9]:

771589.062500014

In [10]:

molecular_weight(taq.translate())

Out[10]:

94387.32270000054

코돈 테이블¶

In [11]:

from Bio.Data import CodonTable

print(CodonTable.unambiguous_dna_by_name["Standard"])

Table 1 Standard, SGC0

  |  T      |  C      |  A      |  G      |
--+---------+---------+---------+---------+--
T | TTT F   | TCT S   | TAT Y   | TGT C   | T
T | TTC F   | TCC S   | TAC Y   | TGC C   | C
T | TTA L   | TCA S   | TAA Stop| TGA Stop| A
T | TTG L(s)| TCG S   | TAG Stop| TGG W   | G
--+---------+---------+---------+---------+--
C | CTT L   | CCT P   | CAT H   | CGT R   | T
C | CTC L   | CCC P   | CAC H   | CGC R   | C
C | CTA L   | CCA P   | CAA Q   | CGA R   | A
C | CTG L(s)| CCG P   | CAG Q   | CGG R   | G
--+---------+---------+---------+---------+--
A | ATT I   | ACT T   | AAT N   | AGT S   | T
A | ATC I   | ACC T   | AAC N   | AGC S   | C
A | ATA I   | ACA T   | AAA K   | AGA R   | A
A | ATG M(s)| ACG T   | AAG K   | AGG R   | G
--+---------+---------+---------+---------+--
G | GTT V   | GCT A   | GAT D   | GGT G   | T
G | GTC V   | GCC A   | GAC D   | GGC G   | C
G | GTA V   | GCA A   | GAA E   | GGA G   | A
G | GTG V   | GCG A   | GAG E   | GGG G   | G
--+---------+---------+---------+---------+--

Tm 값 계산하기¶

In [13]:

from Bio.SeqUtils import MeltingTemp
from Bio.Seq import Seq

primer_DNA = Seq("ATGGAGGAGATGCTGCCCCTCT")
MeltingTemp.Tm_Wallace(primer_DNA)

Out[13]:

70.0

아미노산 서열 약자 변환하기¶

In [14]:

from Bio.SeqUtils import seq1, seq3

seq3("MEEMLPLFEPKGRVLLVDGHHLAYRTFHALKGLTTSRGEPVQAVYGFAKSLLKA")

Out[14]:

'MetGluGluMetLeuProLeuPheGluProLysGlyArgValLeuLeuValAspGlyHisHisLeuAlaTyrArgThrPheHisAlaLeuLysGlyLeuThrThrSerArgGlyGluProValGlnAlaValTyrGlyPheAlaLysSerLeuLeuLysAla'

In [15]:

seq1(
    "MetGluGluMetLeuProLeuPheGluProLysGlyArgValLeuLeuValAspGlyHisHisLeuAlaTyrArgThrPheHisAlaLeuLysGlyLeuThrThrSerArgGlyGluProValGlnAlaValTyrGlyPheAlaLysSerLeuLeuLysAla"
)

Out[15]:

'MEEMLPLFEPKGRVLLVDGHHLAYRTFHALKGLTTSRGEPVQAVYGFAKSLLKA'

Weblogo 그리기¶

In [16]:

from Bio import AlignIO
from Bio.motifs import Motif
from Bio import motifs
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
from IPython.display import Image

alignment = AlignIO.read("./data/HBA.aln", "clustal")
instance = []

for record in alignment:
    s = Seq(str(record.seq), IUPAC.protein)
    instance.append(s)

m = motifs.create(instance)
Motif.weblogo(m, "HBA_logo.png")
Image("HBA_logo.png")

Out[16]:

No description has been provided for this image

계통수(phylogenetic tree) 그리기¶

In [17]:

%matplotlib inline
from Bio import Phylo

tree = Phylo.read("./data/HBA.newick", "newick")
Phylo.draw(tree)

Entrez 데이터 베이스 검색¶

In [3]:

from Bio import Entrez

Entrez.email = "your@email.com"
handle = Entrez.esearch(db="pubmed", term="machine learning")
record = Entrez.read(handle)
print("Pubmed에 machine learning를 검색하면 총 {}개의 결과".format(record["Count"]))

Pubmed에 machine learning를 검색하면 총 29978개의 결과

KEGG API 사용¶

In [4]:

from Bio.KEGG import REST

human_pathways = REST.kegg_list("pathway", "hsa").read()

pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description.lower():
        pathways.append(entry)
        print(entry, description)
print(pathways)

path:hsa03410 Base excision repair - Homo sapiens (human)
path:hsa03420 Nucleotide excision repair - Homo sapiens (human)
path:hsa03430 Mismatch repair - Homo sapiens (human)
['path:hsa03410', 'path:hsa03420', 'path:hsa03430']

In [8]:

genes = []
for pathway in pathways:
    pathway_file = REST.kegg_get(pathway).read()
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()
        if not section == "":
            current_section = section
            if current_section == "GENE":
                gene_identifiers, gene_description = line[12:].split("; ")
                gene_id, gene_symbol = gene_identifiers.split()
                if gene_symbol not in genes:
                    genes.append(gene_symbol)
print("pathway에 연관된 유전자는 {} 이다".format(",".join(genes)))

pathway에 연관된 유전자는 OGG1,RBX1,SSBP1 이다