MSA: multiple sequence alignment
여기에 관한 이론적인 설명은 나중에 또 입 털어드림 ㅇㅇ 

아 참고로 MSA 관련해서 다른건 다 결과가 제대로 나왔는데 툴 관련해서 결과가 안나왔어요 
이게 암만 찾아도 답이 없어서 좀 더 알아보고 수정할 가능성이 있음 

시범조교 앞으로

오늘은 시범조교 가짓수가 좀 많은데 그 중에서도 FASTA 파일들을 좀 올리고자 한다. 이거 말고도 pfam에서 두갠가 받았는데 그건 걍 가서 암거나 받으면 된다. 

해당 파일은 박테리아의 16s rRNA 시퀀스가 들어있는 파일이다. FASTA 파일이라 메모장 있으면 일단 열 수는 있다. 리눅스에서는 gedit으로 만들고 편집하고 다 했다. (vim 안씀) rRNA 시퀀스는 Silva에서 가져왔다. 고마워요 실바! 

Agrobacterium

A. radiobacter를 필두로 하는 뿌리혹세균들(밑에 두놈도 뿌리에 혹 만드는지는 모름)

-Agrobacterium
Agrobacterium radiobacter (구 Agrobacterium tumefaciens)
Agrobacterium agile
Agrobacterium pusense
Agrobacterium salinitolerans

-Rhizobium
Rhizobium tropici (일반적으로 알고 있는 뿌리혹박테리아)
Rhizobium hainanense
Rhizobium gallicum
Rhizobium fabae

-Hoeflea
Hoeflea alexandrii
Hoeflea halophila 웬지 짠 거 좋아하게 생겼는데? 할로박테리움같은건가 
Hoeflea trophica

-Ciceribacter
Ciceribacter lividus
Ciceribacter azotifigens
Ciceribacter thiooxidans

Enterobacter

E.coli를 필두로 하는 장내 세균들

-Enterococcus
Enterococcus faecalis (어디서 많이 봤음)
Enterococcus hirae
Enterococcus avium
Enterococcus caccae

-Escherichia
Escherichia coli (그냥 어디서나 볼 수 있는 대장균)
Escherichia coli O157:H7 (감염되면 X되는 대장균)
Escherichia albertii
Escherichia fergusonii

-Shigella
Shigella boydii
Shigella sonnei
Shigella dysenteriae
Shigella flexneri

Lactobacillus

님들 많이 드시는 그 유산균 맞습니다. 

-Lactobacillus
Lactobacillus acidophilus
Lactobacillus helveticus
Lactobacillus plantarum (이분도 김치에서 발견된다)
Lactobacillus thailandensis DSM 22698 = JCM 13996 (아마 뒤에껀 strain 이름인 듯)

-Leuconostoc
Leuconostoc carnosum
Leuconostoc mesenteroides
Leuconostoc kimchii (이건 있다)
Leuconostoc miyukkimchii (저자 나와봐요 미역김치는 무슨 저세상 학명이야)

-Bifidobacterium (얘네는 방선균인데 일단 유산균으로 섭취 하기는 함... 비피더스 뭐 이런거)
Bifidobacterium bifidum
Bifidobacterium actinocoloniiforme DSM 22766 (아마 strain 이름...)
Bifidobacterium catenulatum

아니 근데 미역김치 학명 실화냐 대체 어디서 발견하면 학명이 저렇게 되는건데요 아니 저기 김치 하나 있는거 뭔데 김치요 참고로 C.kimchii(구 L.kimchii)는 데이터 없어서 못 넣었음... 김치 파티 각 나왔다

시퀀스 가져오기

MSA를 할래도 시퀀스를 가져와야 한다. 그것도 하나 말고 뭉텅이로. 여기서도 read()와 parse()로 나뉘긴 한데, 이거는 기존에 SeqIO로 읽을 때처럼 시퀀스 갯수가 여러 개 있어도 single alignment면 read()로 읽을 수 있다. 

참고로 쿡북 예제로 스톡홀름 파일이 나오긴 했는데 FASTA도 불러올수는 있다. 

from Bio import AlignIO
alignment = AlignIO.read("/home/koreanraichu/agrobacterium.fasta", "fasta")
print(alignment)

MSA를 할 때는 SeqIO가 아니라 AlignIO로 불러오면 된다. 
​

Alignment with 14 rows and 1561 columns
AUUCUCAACUUGAGAGUUUGAUCCUGGCUCAGAACGAACGCUGG...AAG AB102732.157.1651
................................AACGAACGCUGG...... AB680818.1.1406
................................AUUGAACGCUGG...... AB680959.1.1462
................................AACGAACGCUGG...... AB682466.1.1406
................................AACGAACGCUGG...... AB682468.1.1406
......................................CGCUGG...... AB969785.1.1417
............AGAGUUUGAUCCUGGCUCAGAACGAACGCUGG...... AJ786600.1.1449
.................................................. DQ835303.1.1353
.................................................. GU564401.1.1419
................................AACGAACGCUGG...... JF957616.1.1410
.................................................. JQ230000.1.1392
.................................................. KP142169.1.1341
............AGAGUUUGAUCAUGGCUCAGAACGAACGCUGG...... KU975391.1.1451
...............................A-ACGAACGCUGG...... KX510117.1.1407

여기서 말하는 컬럼은 글자 하나. PDB파일처럼 얘도 알파벳 하나가 컬럼이다. 

for record in alignment:
    print(record.id, record.description) # 예제 코드는 Sequence와 ID를 출력하라고 했는데 기니까 다른걸로 바꿔보자.

AB102732.157.1651 AB102732.157.1651 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Agrobacterium radiobacter
AB680818.1.1406 AB680818.1.1406 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Rhizobium tropici
AB680959.1.1462 AB680959.1.1462 Bacteria;Proteobacteria;Gammaproteobacteria;Pseudomonadales;Pseudomonadaceae;Pseudomonas;Agrobacterium agile
AB682466.1.1406 AB682466.1.1406 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Rhizobium hainanense
AB682468.1.1406 AB682468.1.1406 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Rhizobium gallicum
AB969785.1.1417 AB969785.1.1417 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Agrobacterium pusense
AJ786600.1.1449 AJ786600.1.1449 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Hoeflea;Hoeflea alexandrii
DQ835303.1.1353 DQ835303.1.1353 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Rhizobium fabae
GU564401.1.1419 GU564401.1.1419 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Hoeflea;Hoeflea halophila
JF957616.1.1410 JF957616.1.1410 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Hoeflea;Hoeflea phototrophica
JQ230000.1.1392 JQ230000.1.1392 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Ciceribacter;Ciceribacter lividus
KP142169.1.1341 KP142169.1.1341 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Agrobacterium salinitolerans
KU975391.1.1451 KU975391.1.1451 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Ciceribacter;Ciceribacter thiooxidans
KX510117.1.1407 KX510117.1.1407 Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Rhizobiaceae;Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium;Ciceribacter azotifigens

이런 식으로 레코드 만들어도 된다. 예제에는 ID랑 시퀀스로 되어 있던거 ID랑 description으로 바꿨다... 시퀀스가 너무 길어... 출력 format을 지정하거나 쌩으로 레코드를 가져올 수도 있는데, 길이가... 길이가...... 정말 장난없다... 

from Bio import AlignIO
alignments = AlignIO.parse("/home/koreanraichu/lactobacillus.fasta", "fasta")
for alignment in alignments:
    print(alignment)

Alignment with 12 rows and 1606 columns
...............UUUGAUCAUGGCUCAGGACGAACGCUGGC...... AB008203.1.1553
................UUGAUCAUGGCUCAGGACGAACGCUGGC...... AB008210.1.1552
.........................................CGC...... AB022925.1.1450
..........................................GC...... AB023242.1.1446
............GAGUUUGAUCCUGGCUCAGGACGAACGCUGGC...... AB112083.1.1557
...........GGGUUUCGAUUCUGGCUCAGGAUGAACGCUGGC...... AB437354.1.1520
.............GUUUCGAUUCUGGCUCAGGAUGAACGCUGGC...... AB437356.1.1519
...............................GAUGAACGCUGGC...... AF173986.1.1505
...........AGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGC...... AYZK01000017.137.1688
UUUUUUUGUGGAGGGUUUGAUUCUGGCUCAGGAUGAACGCUGGC...AGA CP011786.1420993.1422533
UUUUUUUGUGGAGGGUUUGAUUCUGGCUCAGGAUGAACGCUGGC...AGA CP011786.1514900.1516440
...............................GAUGAACGCUGGC...... KX232108.1.1480

Parse로 가져와도 똑같긴 한데 SeqIO나 얘나 파싱으로 가져오면 for문 있어야 제대로 뜬다. 

from Bio import AlignIO
handle="/home/koreanraichu/lactobacillus.fasta"
for alignments in AlignIO.parse(handle, "fasta",seq_count=2):
    print("Alignment lentgh %i" % alignments.get_alignment_length())
    for record in alignments:
        print(record.description)

Alignment lentgh 1606
AB008203.1.1553 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Lactobacillus;Lactobacillus acidophilus
AB008210.1.1552 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Lactobacillus;Lactobacillus helveticus
Alignment lentgh 1606
AB022925.1.1450 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Leuconostoc;Leuconostoc carnosum
AB023242.1.1446 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Leuconostoc;Leuconostoc mesenteroides
Alignment lentgh 1606
AB112083.1.1557 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Lactiplantibacillus;Lactobacillus plantarum
AB437354.1.1520 Bacteria;Actinobacteriota;Actinobacteria;Bifidobacteriales;Bifidobacteriaceae;Bifidobacterium;Bifidobacterium adolescentis
Alignment lentgh 1606
AB437356.1.1519 Bacteria;Actinobacteriota;Actinobacteria;Bifidobacteriales;Bifidobacteriaceae;Bifidobacterium;Bifidobacterium bifidum
AF173986.1.1505 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Leuconostoc;Leuconostoc kimchii
Alignment lentgh 1606
AYZK01000017.137.1688 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Lacticaseibacillus;Lactobacillus thailandensis DSM 22698 = JCM 13996
CP011786.1420993.1422533 Bacteria;Actinobacteriota;Actinobacteria;Bifidobacteriales;Bifidobacteriaceae;Bifidobacterium;Bifidobacterium actinocoloniiforme DSM 22766
Alignment lentgh 1606
CP011786.1514900.1516440 Bacteria;Actinobacteriota;Actinobacteria;Bifidobacteriales;Bifidobacteriaceae;Bifidobacterium;Bifidobacterium actinocoloniiforme DSM 22766
KX232108.1.1480 Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Leuconostoc;Leuconostoc miyukkimchii

끊어뽑기도 된다. 코드에 있는 숫자 수정하면 3개 4개 묶는 것도 된다. 

쓰기

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Align import MultipleSeqAlignment
from Bio import AlignIO
align1=MultipleSeqAlignment([
    SeqRecord(Seq("GGCC"),id="HaeIII"),
    SeqRecord(Seq("G-CC"),id="id2"),
    SeqRecord(Seq("G--C"),id="id3")
])
align2=MultipleSeqAlignment([
    SeqRecord(Seq("GAATTC"),id="EcoRi"),
    SeqRecord(Seq("G-ATTC"),id="id5"),
    SeqRecord(Seq("G--TTC"),id="id6")
])
my_alignment=[align1,align2]
# 여기까지 시퀀스를 만들고

AlignIO.write(my_alignment, "/home/koreanraichu/sequence.fasta", "fasta")
# 여기서 저(별)장

나 저거 IL-16 전에 썼던거 안된 이유 알아냈음... ㅋㅋㅋㅋㅋㅋ Seq=("")가 아니라 Seq("")였다. 아무튼 이런 식으로 쓸 수 있고 주석 있어서 대충 아시겠지만 순서대로 Alignment 요소를 만들고->시퀀스 만들고->파일로 저장하는 코드이다. 

이런 식으로 저장된다. (그야 내가 FASTA로 저장했으니까...)

참고로 미리 말해두겠지만 본인 리눅스에서 FASTA, PHYLIP, 스톡홀름, clustalW 다 지에딧으로 연다. 그렇게 열리고… 

이게 필립 파일인데 얘는 변환하려면 시퀀스 길이가 다 같아야 하고, 공백이 -이어야 한다. 공백에 . 있으면 지원 안 한다고 오류 토한다. 

선생님 이거 변환 됩니까? 

변환이요? 뭐 클러스탈이나 스톡뭐시기 말하는거면 된다. 

count = AlignIO.convert("/home/koreanraichu/agrobacterium.fasta", "fasta", "/home/koreanraichu/agrobacterium.sth",
                        "stockholm")
print("Converted %i alignments" % count)

alignments = AlignIO.parse("/home/koreanraichu/agrobacterium.fasta", "fasta")
count = AlignIO.write(alignments, "/home/koreanraichu/agrobacterium.aln","clustal")
print("Converted %i alignments" % count)
# ClustalW

위는 convert를 사용해 스톡홀름 파일로 바꾼거고, 아래는 parse와 write를 사용해 클러스탈로 바꾼 것이다. 

위가 clustalw, 아래가 스톡홀름 포맷. 

alignment = AlignIO.read("/home/koreanraichu/PF00096_seed.txt", "stockholm")
AlignIO.write([alignment], "PF05371_seed.aln", "clustal")
print("Converted %i alignments" % count)
# read 후 리스트화해서 변환

read로 읽은 다음 리스트화해서 변환하는 것도 된다. 

count = AlignIO.convert("/home/koreanraichu/PF08449_seed.txt", "stockholm", "/home/koreanraichu/PF08449_seed.phy",
                        "phylip")
print("Converted %i alignments" % count)
# 이거라면 필립 될거같은데?

공백이 -이고 시퀀스 bp가 전부 동일하다는 전제하에 Phylip 포맷으로 만드는 것도 된다. 

alignment2 = AlignIO.read("/home/koreanraichu/PF08449_seed.txt", "stockholm")
name_mapping = {}
for i, record in enumerate(alignment):
    name_mapping[i] = record.id
    record.id = "seq%i" % i
AlignIO.write([alignment], "PF08449_seed_ID.phy", "phylip")
# 오 뭔진 모르겠지만 ID가 숫자가 된 건가

딕셔너리화 한 다음 아이디 새로 지정해주는 것도 된다. 위 코드를 실행하면 

이렇게 된다. 잘 보면 ID 영역에 seq(숫자)가 들어가 있다. 

출력 형식만 바꾸기

이건 뭔 소리냐... 여는 파일은 FASTA파일인데 파일 수정 없이 clustalw나 스톡홀름 형식으로 출력할 수 있다. (phylip은 아마 글자수 같고 공백 -여야 먹힐듯) 참고로 결과는 길어서 생략. 

from Bio import AlignIO
alignment = AlignIO.read("/home/koreanraichu/lactobacillus.fasta", "fasta")
print(format(alignment, "clustal"))

해당 코드는 FASTA파일을 클러스탈 형식으로 보여달라는 얘기고, 스톡홀름으로 볼 수도 있다. 

슬라이싱

넘파이 배열 자르는것처럼 얘도 슬라이싱이 된다. (넘파이 2차원 배열 잘라먹는 거 생각하면 된다)

from Bio import AlignIO
alignment = AlignIO.read("/home/koreanraichu/PF08449_seed.txt", "stockholm")
print(alignment[0:5])

ISLIPISMIMVGCCSNVISLELIMKQSQ-SH-A---------IL...TSG S35B4_DICDI/7-328
SFVLILSLVFGGCCSNVISFEHMVQGSN-INLG---------NI...YGS YEA4_KLULA/2-321
NSLKAFALVFGGCCSNVITFETLMSNET-GSIN---------NL...LGS YEA4_YEAST/3-327
MIASALSFIFGGCCSNAYALEALVREFP-SS-G---------IL...SAR YEA4_SCHPO/1-304
--GVMLSLIFGGCCSNVFALESIIKVEP-GA-G---------TL...L-- Q1K947_NEUCR/81-406

이렇게 하면 pandas의 head()처럼 위에서 다섯개만 보여준다. 

print(alignment[0,1]) # 0번째 인덱스의 두번째 컬럼
S

이런 식으로 2차원 인덱싱과 슬라이싱도 되고

print(alignment[:,6]) # 무슨 컬럼을 가져온거냐
SSASSGTCLSSASLASL

이거는 전체 인덱스에서 일곱번째 글자들을 가져온 것. (실화다) 

print(alignment[0:5,0:5]) # index+column 범위를 지정할 수 있다
ISLIP S35B4_DICDI/7-328
SFVLI YEA4_KLULA/2-321
NSLKA YEA4_YEAST/3-327
MIASA YEA4_SCHPO/1-304
--GVM Q1K947_NEUCR/81-406

이런 식으로 범위 지정해서 잘라오는 것도 되고 특정 인덱스, 컬럼부터 다 가져오는 것도 된다. 저거 저대로 잘라먹는 것도 되고, Numpy 데려와서 배열화하는 것도 가능하다. 근데 예제대로 했더니 오류남... (안되는 건 아닌데 오류뜬다) 

본격적인 정보 얻기

아까까지 대체 뭐 한겨... 

from Bio import AlignIO
alignment = AlignIO.read("/home/koreanraichu/lactobacillus.fasta", "fasta") # 아이고 유산균씨 오랜만입니다
substitutions = alignment.substitutions
print(substitutions)

       .       A       C       G       U
. 1365.0   321.0   310.5   461.5   406.0
A  321.0 19985.0   791.0  2218.5   978.5
C  310.5   791.0 18856.0   957.5  2048.0
G  461.5  2218.5   957.5 26642.0  1015.0
U  406.0   978.5  2048.0  1015.0 16007.0

이거는 rRNA라 티민 대신 우라실이 나왔다. 저건 아마도 해당 염기가 다른걸로 바뀐 비율? 횟수? 이런 걸 봐주는 것 같다. 예를 들자면 아데닌이 시토신, 구아닌, 우라실로 바뀌는 뭐 그런거. (공백 공백은 뭐여...) 저기서 옵션으로 .select() 주면 출력 순서가 바뀐다. (예: .select(“AUGC.”)으로 주면 행렬 인덱스랑 컬럼 순서가 AUCG.로 바뀐다)

Tool

이 부분이 문제가 많아요... OTL 둘 다 깔고 해야 하나... 

ClustalW

from Bio.Align.Applications import ClustalwCommandline
cline = ClustalwCommandline("clustalw2", infile="/home/koreanraichu/lactobacillus.fasta")
print(cline)

이런 식으로 불러서 

clustalw2 -infile=/home/koreanraichu/lactobacillus.fasta

니가 왜 거기서 나옴???

이거 뭐 트리 만드는 파일 만들어준다매 왜 안해줘요. 이럴거면 그냥 메가 쓰면 안되냐

MUSCLE

from Bio.Align.Applications import MuscleCommandline
cline = MuscleCommandline(input="/home/koreanraichu/lactobacillus.aln", out="/home/koreanraichu/lactobacillus.txt")
print(cline)

얘는 이런 식으로 부르면 

muscle -in /home/koreanraichu/lactobacillus.aln -out /home/koreanraichu/lactobacillus.txt

이러시는 이유가 있으실 것 아니예요. 

from Bio.Align.Applications import MuscleCommandline
cline = MuscleCommandline(input="/home/koreanraichu/lactobacillus.fasta", out="/home/koreanraichu/lactobacillus.txt",
clw=True)
print(cline)
muscle -in /home/koreanraichu/lactobacillus.fasta -out /home/koreanraichu/lactobacillus.txt -clw

이걸 쓰면 clustalW 비슷하게 출력된다. (물론 파일은 생기지 않았습니다)

from Bio.Align.Applications import MuscleCommandline
cline = MuscleCommandline(input="/home/koreanraichu/lactobacillus.fasta", out="/home/koreanraichu/lactobacillus.txt",
clwstrict=True)
print(cline)
muscle -in /home/koreanraichu/lactobacillus.fasta -out /home/koreanraichu/lactobacillus.txt -clwstrict

이걸 쓰면 clustalW headline을 쓴다는데 파일이 없다. 착한 사람만 보이는 파일인가... 그래서 얘네는 못했다. 

Pairwise

왜 블래스트에서 뭐 검색하면 결과에 줄 그어서 나오는 그거 말하는 거 맞다. 

from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("/home/koreanraichu/alpha.faa", "fasta")
seq2 = SeqIO.read("/home/koreanraichu/beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)
print(alignments)

근데 이거 파싱으로 개별 레코드끼리 비교 못하나... 나중에 다시 함 트라이해 볼 예정. 

[Alignment(seqA='MV-LSPADKTNV---K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----T--PEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=217), Alignment(seqA='MV-LSPADKTNV--K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----T-PEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LSPADKTNV--K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----TP-EEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LS-PADKTNV--K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-TP------EEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=217), Alignment(seqA='MV-LSPADKTNV--K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHLTP------EEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LSPADKTN-V-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----T-PEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LSPADKTNV-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----TPEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=215), Alignment(seqA='MV-LS-PADKTNV-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-TP-----EEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LSPADKTNV-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHLTP-----EEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=215), Alignment(seqA='MV-LSPADKT-NV-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----TPE-EKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=216), Alignment(seqA='MV-LS-PADKTNV-K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-

뭔진 모르겠지만 제가 잘못했습니다. (저거 len으로 보면 길이 장난없다)

print(pairwise2.format_alignment(*alignments[0]))

이렇게 쳐 주면 비로소 우리가 생각하는(?)

MV-LSPADKTNV---K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-
|| |     |     | |  | ||||        |  |  |||  |  |      |    |   |  |   |  |  |   | |  |      |||||  |  |  |    ||   |  | |     | ||  |  |  ||  |||  ||   |   |   ||   | ||       ||||  |  |      |    |  |      |    ||  
MVHL-----T--PEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H
  Score=72

이 결과가 나온다. BLAST에서도 이런 식으로 나온다. 

BLOSUM62와 함께 정렬...아니고 줄을 

from Bio import pairwise2
from Bio import SeqIO
from Bio.Align import substitution_matrices
blosum62 = substitution_matrices.load("BLOSUM62")
seq1 = SeqIO.read("/home/koreanraichu/alpha.faa", "fasta")
seq2 = SeqIO.read("/home/koreanraichu/beta.faa", "fasta")
alignments = pairwise2.align.globalds(seq1.seq, seq2.seq,blosum62, -10, -0.5)
print(pairwise2.format_alignment(*alignments[0]))

이거 [0] 빼면요? 뭔진 모르겠지만 잘못했어요 소리가 절로 나온다. (길이 장난없음) 

MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF-DLS-----HGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR
|| |.|..|..|.|.||||  ...|.|.|||.|.....|.|...|..| |||     .|...||.|||||..|.....||.|........||.||..||.|||.||.||...|...||.|...||||.|.|...|..|.|...|..||.
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH
  Score=292.5

.globalds는 전역(전체 시퀀스)를 보는거고 .localds로 부분만 볼 수도 있다. 

alignments = pairwise2.align.localds("LSPADKTNVKAA", "PEEKSAV", blosum62, -10, -1)

3 PADKTNV
  |..|..|
1 PEEKSAV
  Score=16

이런 식으로 부분부분 볼 수도 있다. 엥? 근데 저렇게 해 놓으면 어딘지 어떻게 알아요? 

print(pairwise2.format_alignment(*alignments[0],full_sequences=True))

LSPADKTNVKAA
  |..|..|   
--PEEKSAV---

개발자는 다 계획이 있었다. 

Pairwise aligner

from Bio import Align
aligner = Align.PairwiseAligner()
seq1="GAATTC"
seq2="GCATTC"
score = aligner.score(seq1, seq2)
print(score)
5.0

일단 모셔보래서 모셔봤습니다. 이거는 score 뽑는거고 

from Bio import Align
aligner = Align.PairwiseAligner() # 그래서 모셔왔습니다 
seq1="GAATTC"
seq2="GCATTC"
alignments = aligner.align(seq1, seq2)
for alignment in alignments:
    print(alignment)

이렇게 치면 

G-AATTC
|-|-|||
GCA-TTC

GA-ATTC
|--||||
G-CATTC

G-AATTC
|--||||
GC-ATTC

GAATTC
|.||||
GCATTC

이렇게 나온다. 전 맨 밑에껄로 보여주세요. 픽 안되나

GAATTC
||||||
GAATTC

100% 일치하는 시퀀스는 이런 식으로 보여준다. aligner.mode="local"을 주면 위에것처럼 중구난방으로 안 뽑아주고 

GAATTC
 |||||
CAATTC

이런 식으로 뽑아준다. 

Wildcard

일반적으로 검색할 때(구글이나 파일같은 거) *이 와일드카드지만, 여기서는 와일드카드를 지정할 수 있다. 

# wildcard를 지정한 다음 쓸 수도 있다
aligner.wildcard = "?"
seq3="?GCC"
seq4="GGCC"
alignments2 = aligner.align(seq3, seq4)
for alignment in alignments2:
    print(alignment)

?G-CC
 |-||
 GGCC

?GCC
 |||
GGCC

여기서는 일단 ?로 지정해보자. 

파이참 쓰다보니 키보드 먹통돼서 재부팅 여러번 했다. 이것도 갈 때 된 듯... 

*참고로 그냥 리드나 파싱으로 긁어오는 방법은 첫 글에서 썼으므로 생략한다. 

Iteration 활용하기

근데 들어도 이건 뭔 기능인가 싶긴 함. 

from Bio import SeqIO
identifiers = [seq_record.id for seq_record in SeqIO.parse("/home/koreanraichu/sequence.gb", "genbank")]
print(identifiers)

뭐야 wrap 넣어줘요... 

레코드와 for문을 조합해서 원하는 데이터(ID나 시퀀스, description)만 가져올 수 있다. 

from Bio import SeqIO
identifiers = [seq_record.id for seq_record in SeqIO.parse("/home/koreanraichu/ls_orchid.gbk", "genbank")]
print(identifiers)

['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', 'Z78526.1', 'Z78525.1', 'Z78524.1', 'Z78523.1', 'Z78522.1', 'Z78521.1', 'Z78520.1', 'Z78519.1', 'Z78518.1', 'Z78517.1', 'Z78516.1', 'Z78515.1', 'Z78514.1', 'Z78513.1', 'Z78512.1', 'Z78511.1', 'Z78510.1', 'Z78509.1', 'Z78508.1', 'Z78507.1', 'Z78506.1', 'Z78505.1', 'Z78504.1', 'Z78503.1', 'Z78502.1', 'Z78501.1', 'Z78500.1', 'Z78499.1', 'Z78498.1', 'Z78497.1', 'Z78496.1', 'Z78495.1', 'Z78494.1', 'Z78493.1', 'Z78492.1', 'Z78491.1', 'Z78490.1', 'Z78489.1', 'Z78488.1', 'Z78487.1', 'Z78486.1', 'Z78485.1', 'Z78484.1', 'Z78483.1', 'Z78482.1', 'Z78481.1', 'Z78480.1', 'Z78479.1', 'Z78478.1', 'Z78477.1', 'Z78476.1', 'Z78475.1', 'Z78474.1', 'Z78473.1', 'Z78472.1', 'Z78471.1', 'Z78470.1', 'Z78469.1', 'Z78468.1', 'Z78467.1', 'Z78466.1', 'Z78465.1', 'Z78464.1', 'Z78463.1', 'Z78462.1', 'Z78461.1', 'Z78460.1', 'Z78459.1', 'Z78458.1', 'Z78457.1', 'Z78456.1', 'Z78455.1', 'Z78454.1', 'Z78453.1', 'Z78452.1', 'Z78451.1', 'Z78450.1', 'Z78449.1', 'Z78448.1', 'Z78447.1', 'Z78446.1', 'Z78445.1', 'Z78444.1', 'Z78443.1', 'Z78442.1', 'Z78441.1', 'Z78440.1', 'Z78439.1']

예제에 있단 파일로 돌렸더니 94개 나온다. 근데 94개가 저렇게 중구난방으로 나오니까 힘들다... 이거 넘파이 배열로 바꿔서 깔쌈하게 보여주면 안되나? 

from Bio import SeqIO
import numpy as np
identifiers = [seq_record.id for seq_record in SeqIO.parse("/home/koreanraichu/ls_orchid.gbk", "genbank")]
identifiers=np.array(identifiers)
identifiers=identifiers.reshape(2,47)
print(identifiers)

그래서 해봤다. 

[['Z78533.1' 'Z78532.1' 'Z78531.1' 'Z78530.1' 'Z78529.1' 'Z78527.1'
  'Z78526.1' 'Z78525.1' 'Z78524.1' 'Z78523.1' 'Z78522.1' 'Z78521.1'
  'Z78520.1' 'Z78519.1' 'Z78518.1' 'Z78517.1' 'Z78516.1' 'Z78515.1'
  'Z78514.1' 'Z78513.1' 'Z78512.1' 'Z78511.1' 'Z78510.1' 'Z78509.1'
  'Z78508.1' 'Z78507.1' 'Z78506.1' 'Z78505.1' 'Z78504.1' 'Z78503.1'
  'Z78502.1' 'Z78501.1' 'Z78500.1' 'Z78499.1' 'Z78498.1' 'Z78497.1'
  'Z78496.1' 'Z78495.1' 'Z78494.1' 'Z78493.1' 'Z78492.1' 'Z78491.1'
  'Z78490.1' 'Z78489.1' 'Z78488.1' 'Z78487.1' 'Z78486.1']
 ['Z78485.1' 'Z78484.1' 'Z78483.1' 'Z78482.1' 'Z78481.1' 'Z78480.1'
  'Z78479.1' 'Z78478.1' 'Z78477.1' 'Z78476.1' 'Z78475.1' 'Z78474.1'
  'Z78473.1' 'Z78472.1' 'Z78471.1' 'Z78470.1' 'Z78469.1' 'Z78468.1'
  'Z78467.1' 'Z78466.1' 'Z78465.1' 'Z78464.1' 'Z78463.1' 'Z78462.1'
  'Z78461.1' 'Z78460.1' 'Z78459.1' 'Z78458.1' 'Z78457.1' 'Z78456.1'
  'Z78455.1' 'Z78454.1' 'Z78453.1' 'Z78452.1' 'Z78451.1' 'Z78450.1'
  'Z78449.1' 'Z78448.1' 'Z78447.1' 'Z78446.1' 'Z78445.1' 'Z78444.1'
  'Z78443.1' 'Z78442.1' 'Z78441.1' 'Z78440.1' 'Z78439.1']]

참고로 numpy array는 5*6을 만들려면 진짜 데이터가 30개 있어야된다. 즉, 저게 최선이다. 

레코드 반복

이것도 Iteration을 이용한다. 

from Bio import SeqIO
record_iterator = SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta", "fasta")

first_record = next(record_iterator)
print(first_record.id)
print(first_record.description)

second_record = next(record_iterator)
print(second_record.id)
print(second_record.description)

from Bio import SeqIO
first_record_2 = next(SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta", "fasta"))
print(first_record_2)

위 코드와 아래 코드 둘 다 iteration을 활용하는 건 맞다. 

2B10_1|Chains
2B10_1|Chains A, C|Cytochrome c peroxidase, mitochondrial|Saccharomyces cerevisiae (4932)
2B10_2|Chains
2B10_2|Chains B, D|Cytochrome c iso-1|Saccharomyces cerevisiae (4932)

ID: 2B10_1|Chains
Name: 2B10_1|Chains
Description: 2B10_1|Chains A, C|Cytochrome c peroxidase, mitochondrial|Saccharomyces cerevisiae (4932)
Number of features: 0
Seq('TTPLVHVASVEKGRSYEDFQKVYNAIALKLREDDEYDNYIGYGPVLVRLAWHTS...QGL')

위 코드는 next 수가 두 개라 앞에서 두 데이터+원하는 부분만 띄운 거고 아래는 판다스의 head()마냥 맨 위에걸 띄운거다. 아니 그러면 next를 일일이 쳐야되나요? for문 안됨? 있어봐요 해봅시다. 

from Bio import SeqIO
record_iterator = SeqIO.parse("/home/koreanraichu/ls_orchid.fasta", "fasta")
for i in range(0,4):
    record = next(record_iterator)
    print(record.id)
    print(record.description)

gi|2765658|emb|Z78533.1|CIZ78533
gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
gi|2765657|emb|Z78532.1|CCZ78532
gi|2765657|emb|Z78532.1|CCZ78532 C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA
gi|2765656|emb|Z78531.1|CFZ78531
gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA
gi|2765655|emb|Z78530.1|CMZ78530
gi|2765655|emb|Z78530.1|CMZ78530 C.margaritaceum 5.8S rRNA gene and ITS1 and ITS2 DNA

되는데요? 

from Bio import SeqIO
record_iterator = SeqIO.parse("/home/koreanraichu/ls_orchid.fasta", "fasta")

i=0
while i<4
    record = next(record_iterator)
    print(record.id)
    print(record.description)
    i=i+1

참고로 While문도 된다. 

레코드 목록 뽑기

from Bio import SeqIO
records = list(SeqIO.parse("/home/koreanraichu/sequence.gb", "genbank"))
print(records)

이게 무옵션 기본 코드이다. 

print("First record:",records[0])
print("Last record:",records[-1])

인덱싱도 된다. 

first=records[0]
last=records[-1]
# 변수 안 주면 풀로 나온다...
print("First record ID:",first.id)
print("Last record ID:",last.id)

인덱싱 하면서 특정 인자만 보는 것도 된다. (위 코드는 처음과 맨 마지막 데이터의 아이디를 보여주는 코드)

데이터 추출

from Bio import SeqIO
record_iterator = SeqIO.parse("/home/koreanraichu/ls_orchid.gbk", "genbank")
first_record = next(record_iterator)
# 첫 번째 레코드를 가져온다.
print(first_record.annotations) # Annotations
print(first_record.annotations.keys()) # Key
print(first_record.annotations.values()) # Values

{'molecule_type': 'DNA', 'topology': 'linear', 'data_file_division': 'PLN', 'date': '30-NOV-2006', 'accessions': ['Z78533'], 'sequence_version': 1, 'gi': '2765658', 'keywords': ['5.8S ribosomal RNA', '5.8S rRNA gene', 'internal transcribed spacer', 'ITS1', 'ITS2'], 'source': 'Cypripedium irapeanum', 'organism': 'Cypripedium irapeanum', 'taxonomy': ['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliophyta', 'Liliopsida', 'Asparagales', 'Orchidaceae', 'Cypripedioideae', 'Cypripedium'], 'references': [Reference(title='Phylogenetics of the slipper orchids (Cypripedioideae: Orchidaceae): nuclear rDNA ITS sequences', ...), Reference(title='Direct Submission', ...)]} 
# 그냥 주석을 출력했을 때
dict_keys(['molecule_type', 'topology', 'data_file_division', 'date', 'accessions', 'sequence_version', 'gi', 'keywords', 'source', 'organism', 'taxonomy', 'references'])
# 키 주세요
dict_values(['DNA', 'linear', 'PLN', '30-NOV-2006', ['Z78533'], 1, '2765658', ['5.8S ribosomal RNA', '5.8S rRNA gene', 'internal transcribed spacer', 'ITS1', 'ITS2'], 'Cypripedium irapeanum', 'Cypripedium irapeanum', ['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliophyta', 'Liliopsida', 'Asparagales', 'Orchidaceae', 'Cypripedioideae', 'Cypripedium'], [Reference(title='Phylogenetics of the slipper orchids (Cypripedioideae: Orchidaceae): nuclear rDNA ITS sequences', ...), Reference(title='Direct Submission', ...)]])
# 값 주세요

Annotation은 저장 형태가 딕셔너리라 이렇게 나온다고... 딕셔너리 형태라는 건 Key:value로 저장되어 있기 때문에 키값이나 밸류값만 볼 수 있다는 얘기이기도 하다. 

print(first_record.annotations["source"]) # Annotation+selector(키 중 하나를 선택)
Cypripedium irapeanum # 결과

이런 식으로 키 중 하나를 선택해서 거기에 대한 밸류를 볼 수 있다. 

# 아래 코드는 리스트업 코드이다.
all_species = []
for seq_record in record_iterator:
    all_species.append(seq_record.annotations["organism"])
print(all_species)

숫자만 맞으면 어레이 만들어도 된다. 

all_species = [
    seq_record.annotations["organism"]
    for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")
]
print(all_species)

for문에 던져버려도 되고. 

수정

from Bio import SeqIO
record_iterator = SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10_mod.fasta", "fasta") # 불러와서(사본 만듦)
first_record = next(record_iterator)
first_record.id="SC_2B10_1" # ID를 바꿨다.
print(first_record.id)

원본 파일에 손대기 좀 그래서 사본 만들어서 수정했다. 수정은 이런 식으로 하면 된다. 

웹에서 받아오기

아 솔직치 다운받기 귀찮은데 아이디 알면 걍 웹에서 뜯어오면 안돼요? 됨. 

파일 주세요 Genbank여... 

from Bio import Entrez
from Bio import SeqIO
Entrez.email = "blackholekun@gmail.com"
with Entrez.efetch(
    db="nucleotide", rettype="fasta", retmode="text", id="60391722"
) as handle:
    seq_record = SeqIO.read(handle, "fasta")
print(seq_record)

순순히 파일을 내놓는다면 유혈사태는 면할 것이다! (메일 내꺼임)

ID: AB018076.3
Name: AB018076.3
Description: AB018076.3 Homo sapiens hedgehog gene, complete cds
Number of features: 0
Seq('ATGTCTCCCGCCCGGCTCCGGCCCCGACTGCACTTCTGCCTGGTCCTGTTGCTG...ACC')

그래서 받아옴

from Bio import Entrez
from Bio import SeqIO
Entrez.email = "blackholekun@gmail.com"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="60391722"
) as handle:
    seq_record = SeqIO.read(handle, "gb")
print(seq_record)

파일 또 내놔! 

ID: AB018076.3
Name: AB018076
Description: Homo sapiens hedgehog gene, complete cds
Number of features: 8
/molecule_type=DNA
/topology=linear
/data_file_division=PRI
/date=24-FEB-2005
/accessions=['AB018076', 'AB010092', 'AB018075']
/sequence_version=3
/keywords=['']
/source=Homo sapiens (human)
/organism=Homo sapiens
/taxonomy=['Eukaryota', 'Metazoa', 'Chordata', 'Craniata', 'Vertebrata', 'Euteleostomi', 'Mammalia', 'Eutheria', 'Euarchontoglires', 'Primates', 'Haplorrhini', 'Catarrhini', 'Hominidae', 'Homo']
/references=[Reference(title='Expression of Sonic hedgehog and its receptor Patched/Smoothened in human cancer cell lines and embryonic organs', ...), Reference(title='Direct Submission', ...)]
/comment=On or before Mar 1, 2005 this sequence version replaced AB010092.1,
AB018075.1, AB018076.2.
Seq('ATGTCTCCCGCCCGGCTCCGGCCCCGACTGCACTTCTGCCTGGTCCTGTTGCTG...ACC')

아, 이건 genbank 포맷이다. 파스타랑 저거랑 다 된다. 단백질 이름이요? 아니 진짜 있는건데??? 

from Bio import Entrez
from Bio import SeqIO
Entrez.email = "blackholekun@gmail.com"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="60391722,60391706,1519245148" # 여러개도 된다.
) as handle:
    for seq_record in SeqIO.parse(handle, "gb"):
        print(seq_record.description)

Homo sapiens hedgehog gene, complete cds
Homo sapiens hedgehog gene, complete cds
Homo sapiens sonic hedgehog signaling molecule (SHH), transcript variant 1, mRNA

아 물론 뭉텅이로 받는것도 된다. 

파일 주세요 swissprot이여... 

솔직히 이정도면 PDB랑 유니프롯이랑 TAIR랑 펍켐도 지원해줘야된다. 개발자 일해라. 

from Bio import ExPASy
from Bio import SeqIO
with ExPASy.get_sprot_raw("O23729") as handle:
    seq_record = SeqIO.read(handle, "swiss")
print(seq_record)

스위스프롯에서는 이렇게 가져온다. 

ID: O23729
Name: CHS3_BROFI
Description: RecName: Full=Chalcone synthase 3; EC=2.3.1.74; AltName: Full=Naringenin-chalcone synthase 3;
Database cross-references: EMBL:AF007097, SMR:O23729, PRIDE:O23729, UniPathway:UPA00154, GO:GO:0016210, GO:GO:0009813, Gene3D:3.40.47.10, InterPro:IPR012328, InterPro:IPR001099, InterPro:IPR018088, InterPro:IPR011141, InterPro:IPR016039, PANTHER:PTHR11877, Pfam:PF02797, Pfam:PF00195, PIRSF:PIRSF000451, SUPFAM:SSF53901, PROSITE:PS00441
Number of features: 2
/molecule_type=protein
/accessions=['O23729']
/protein_existence=2
/date=15-JUL-1999
/sequence_version=1
/date_last_sequence_update=01-JAN-1998
/date_last_annotation_update=10-FEB-2021
/entry_version=75
/gene_name=Name=CHS3;
/organism=Bromheadia finlaysoniana (Orchid)
/taxonomy=['Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliopsida', 'Liliopsida', 'Asparagales', 'Orchidaceae', 'Epidendroideae', 'Vandeae', 'Adrorhizinae', 'Bromheadia']
/ncbi_taxid=['41205']
/comment=FUNCTION: The primary product of this enzyme is 4,2',4',6'- tetrahydroxychalcone (also termed naringenin-chalcone or chalcone) which can under specific conditions spontaneously isomerize into naringenin.
CATALYTIC ACTIVITY: Reaction=4-coumaroyl-CoA + 2 H(+) + 3 malonyl-CoA = 2',4,4',6'-   tetrahydroxychalcone + 3 CO2 + 4 CoA; Xref=Rhea:RHEA:11128,   ChEBI:CHEBI:15378, ChEBI:CHEBI:16526, ChEBI:CHEBI:57287,   ChEBI:CHEBI:57355, ChEBI:CHEBI:57384, ChEBI:CHEBI:77645; EC=2.3.1.74;   Evidence={ECO:0000255|PROSITE-ProRule:PRU10023};
PATHWAY: Secondary metabolite biosynthesis; flavonoid biosynthesis.
SIMILARITY: Belongs to the thiolase-like superfamily. Chalcone/stilbene synthases family. {ECO:0000305}.
/references=[Reference(title='Molecular cloning and sequence analysis of chalcone synthase cDNAs of Bromheadia finlaysoniana.', ...)]
/keywords=['Acyltransferase', 'Flavonoid biosynthesis', 'Transferase']
Seq('MAPAMEEIRQAQRAEGPAAVLAIGTSTPPNALYQADYPDYYFRITKSEHLTELK...GAE')

뭔진 모르겠지만 일단 잘못했어요

참고로 저거 치니까 유니프롯(Uniprot) 나오길래 여기서 아이디 제공해주는건가 했는데 저기다가 유니프롯 아이디 넣으니까 됨. 

from Bio import ExPASy
from Bio import SeqIO
with ExPASy.get_sprot_raw("Q14005") as handle:
    seq_record = SeqIO.read(handle, "swiss")
print(seq_record.name)
print(seq_record.description)

Q14005: Interleukin-16

IL16_HUMAN
RecName: Full=Pro-interleukin-16; Contains: RecName: Full=Interleukin-16; Short=IL-16; AltName: Full=Lymphocyte chemoattractant factor; Short=LCF;

뭐야 왜 돼요 ㄷㄷ 

딕셔너리화하기

뭐 하면 뭐 좋은갑지... 모르것다... (뇌클리어)

SeqIO.to_dict

from Bio import SeqIO
orchid_dict = SeqIO.to_dict(SeqIO.parse("/home/koreanraichu/sequence.gb", "genbank"))
print(orchid_dict.keys())
print(orchid_dict.values())

dict_keys(['NM_000517.6'])
dict_values([SeqRecord(seq=Seq('ACTCTTCTGGTCCCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGC...GCA'), id='NM_000517.6', name='NM_000517', description='Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA', dbxrefs=[])])

이런 식으로 딕셔너리화할 수 있고 

print(orchid_dict['2B10_1|Chains'].description)
2B10_1|Chains A, C|Cytochrome c peroxidase, mitochondrial|Saccharomyces cerevisiae (4932)

이런 식으로 셀렉터 적용도 된다. 

for record in SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta", "fasta"):
    print(record.id, seguid(record.seq))

SEGUID checksum 함수를 적용시켜 딕셔너리화 할 수도 있다는데... 

2B10_1|Chains kSIi2w+CcP3y/k3bFRVhvvWtCtM
2B10_2|Chains SA4dH3bOSXcG43088Houu01TH2U

이건 키가 일단 너무 저세상이다. InChl key 못지 않음... 

seguid_dict = SeqIO.to_dict(SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta", "fasta"),lambda rec: seguid(rec.seq))
record = seguid_dict["kSIi2w+CcP3y/k3bFRVhvvWtCtM"]
print(record.id)
2B10_1|Chains

그러니까 저 함수 적용하고 저 키를 다 외우고 있으면 찾는 게 가능하다는 얘기. 

SeqIO.index

from Bio import SeqIO
orchid_dict = SeqIO.index("/home/koreanraichu/Octodon degus insulin.gb", "genbank")
print(orchid_dict.keys())
print(orchid_dict.values())

KeysView(SeqIO.index('/home/koreanraichu/Octodon degus insulin.gb', 'genbank', alphabet=None, key_function=None))
ValuesView(SeqIO.index('/home/koreanraichu/Octodon degus insulin.gb', 'genbank', alphabet=None, key_function=None))

예제 보면 리스트로 잘만 나오던데 나는 왜 저따구여? 

쓰기

언제까지 긁어만 올텐가... 

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

rec1 = SeqRecord(
    Seq(
        "MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD"
        "GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK"
        "NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM"
        "SSAC",
    ),
    id="gi|14150838|gb|AAK54648.1|AF376133_1",
    description="chalcone synthase [Cucumis sativus]",
)

rec2 = SeqRecord(
    Seq(
        "YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ"
        "DMVVVEIPKLGKEAAVKAIKEWGQ",
    ),
    id="gi|13919613|gb|AAK33142.1|",
    description="chalcone synthase [Fragaria vesca subsp. bracteata]",
)

rec3 = SeqRecord(
    Seq(
        "MVTVEEFRRAQCAEGPATVMAIGTATPSNCVDQSTYPDYYFRITNSEHKVELKEKFKRMC"
        "EKSMIKKRYMHLTEEILKENPNICAYMAPSLDARQDIVVVEVPKLGKEAAQKAIKEWGQP"
        "KSKITHLVFCTTSGVDMPGCDYQLTKLLGLRPSVKRFMMYQQGCFAGGTVLRMAKDLAEN"
        "NKGARVLVVCSEITAVTFRGPNDTHLDSLVGQALFGDGAAAVIIGSDPIPEVERPLFELV"
        "SAAQTLLPDSEGAIDGHLREVGLTFHLLKDVPGLISKNIEKSLVEAFQPLGISDWNSLFW"
        "IAHPGGPAILDQVELKLGLKQEKLKATRKVLSNYGNMSSACVLFILDEMRKASAKEGLGT"
        "TGEGLEWGVLFGFGPGLTVETVVLHSVAT",
    ),
    id="gi|13925890|gb|AAK49457.1|",
    description="chalcone synthase [Nicotiana tabacum]",
)

my_records = [rec1, rec2, rec3]

from Bio import SeqIO
SeqIO.write(my_records, "/home/koreanraichu/example.fasta", "fasta")

아니 이거 내가 IL로 트라이했는데 에러났음... 예제껀 잘됐는데 뭔 타입에러가 반기던데? 

아무튼 이런 식으로 저장된다. (메시지 출력 따로 안했다)

from Bio import SeqIO

records = SeqIO.parse("/home/koreanraichu/sequence.gb", "genbank")
count = SeqIO.write(records, "/home/koreanraichu/my_example.fasta", "fasta")
print("Converted %i records" % count)

count2 = SeqIO.convert("/home/koreanraichu/sequence.gb", "genbank", "/home/koreanraichu/my_example2.fasta", "fasta")
print("Converted %i records" % count2)

긁어오는 방법은 두 가지고 둘 다 된다. 

from Bio import SeqIO
records = [rec.reverse_complement(id="rc_"+rec.id, description = "reverse complement")
           for rec in SeqIO.parse("/home/koreanraichu/at5g40780.gb", "genbank")]
SeqIO.write(records, "/home/koreanraichu/at5g40780_rc.fasta", "fasta")

DNA의 경우 reverse complement sequence로 저장하는 것도 된다. (단백질은 하지말자...)

오늘껀 근데 하면서 나도 이해 못했음... 

SeqRecord를 만들려면 뭘 또 모셔와야 하는데 바로 

from Bio.SeqRecord import SeqRecord

이놈이다. 

이정도면 복잡한 거 돌리려면 아주 모셔오는데만 열댓줄 쓰게 생겼는데??? 텐서플로우나 저기 그 넘파이 판다스처럼 합쳐버리지 그러냐... (걔네는 일단 한 번 데려오면 다 써먹을 수 있음) 

레코드 생성하기

이모지 넣고싶은데 이거 리눅스라 이모지 넣는 법을 모름... 아무튼 레코드는 이렇게 생성하면 된다. 

from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
EcoRI=Seq("GAATTC") # ECoRI sequence
ECoRI_r=SeqRecord(EcoRI) #로 레코드를 만들어보겠습니다 
print(EcoRI)
print(ECoRI_r)

그러면 어떻게 나오느냐 

GAATTC # 시퀀스
ID: <unknown id> # 여기서부터 레코드
Name: <unknown name>
Description: <unknown description>
Number of features: 0
Seq('GAATTC')

이렇게 나온다. 뭐가 많이 빠졌지만 일단 레코드 맞다. 여기다 넣고싶으면 일단 다음 문단으로 가시고... 아니 그럼 생성할 때 넣을 수 있나요? 예. 됩니다. 

HaeIII_r=SeqRecord(HaeIII,id="RES_HAEIII",name="HaeIII restriction site",description="Restriction of HaeIII, cuts GG/CC")

노션도 되는 wrap이 안되는 코드블록... 네이버 니네 점검하면 대체 뭐하냐... 아무튼 이렇게 만들 때 넣으면 

ID: RES_HAEIII
Name: HaeIII restriction site
Description: Restriction of HaeIII, cuts GG/CC
Number of features: 0
Seq('GGCC')

이렇게 생성된다. 

레코드에 뭔가 넣기

아니 아까 만든거에 뭐 좀 어떻게 넣어주면 안됩니까? 아뇨 됩니다. 

# 레코드에 정보를 추가해보자
ECoRI_r.id="RES_ECORI"
ECoRI_r.name="ECoRI restriction site"
ECoRI_r.description="Restriction site of ECoRI, cuts G/AATTC"
# ID랑 이름, description 정보를 추가했다.

이렇게 필요한 정보를 넣어줄 수 있다. (참고: EcoRI이 맞는 표기다)

ID: RES_ECORI
Name: ECoRI restriction site
Description: Restriction site of ECoRI, cuts G/AATTC
Number of features: 0
Seq('GAATTC')

제가 된다고 했잖아여... 

HaeIII_r.annotations["notes"]="갑자기 해삼이 먹고싶다. " # annotation
HaeIII_r.letter_annotations['letter annotations']=['G','G/','C','C'] # letter annotations

/notes=갑자기 해삼이 먹고싶다. 
Per letter annotation for: letter annotations

주석도 넣을 수 있다. 참고로 HaeIII은 Haemophilus influenzae biogroup aegyptius에서 유래한 제한효소. 본인은 해삼 혹은 해쓰리로 읽는다. (자매품 XbaI...)

FASTA에서 긁어오기

FASTA와 Genbank 파일에서 레코드를 긁어올 수도 있다. 후자쪽이 정보는 많음. 

from Bio import SeqIO
record=SeqIO.read('/home/koreanraichu/rcsb_pdb_7PNN.fasta','fasta')
print(record)

record=SeqIO.parse('/home/koreanraichu/rcsb_pdb_2B10.fasta','fasta')

FASTA 파일에 하나 들어있으면 read, 두 개 이상이면 parse로 긁어야 한다. (두 개 이상인데 read로 긁으면 ValueError: More than one record found in handle가 당신을 반길 것이다)

print(record.id)
7PNN_1|Chain

for record2 in SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta","fasta"):
    print(record2.id)
2B10_1|Chains
2B10_2|Chains

긁는 난이도치고 접근은 되게 쉽다. 근데 parse로 긁어오면 for문 써야 나온다... 

Genbank에서 긁어오기

record=SeqIO.read("/home/koreanraichu/sequence.gb",'genbank')
print(record)

print(record.description)
Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA

얘는 리드로 파스고 없이 걍 긁어오면 접근이 바로 된다. 편한데? 

Feature

시퀀스 레코드의 정보를 담고 있다고 하는데 뭔지 모르겠는 건 기분탓이 아니다. 여기서 볼 것은 .type와 .location, .qualifiers이다. 

.type

for feature in record.features:
    print(feature.type)
source
gene
CDS
sig_peptide
mat_peptide
mat_peptide
mat_peptide
regulatory
polyA_site

얘도 for문 돌려야 나온다. 타입은 말 그대로 타입. 

​

.location

for feature in record.features:
    print(feature.location)
[0:432](+)
[0:432](+)
[41:371](+)
[41:113](+)
[113:200](+)
[206:293](+)
[299:368](+)
[413:419](+)
[431:432](+)

Feature의 위치 정보. +는 뭐임? 

.qualifiers

for feature in record.features:
    print(feature.qualifiers)
OrderedDict([('organism', ['Octodon degus']), ('mol_type', ['mRNA']), ('db_xref', ['taxon:10160']), ('tissue_type', ['pancreas'])])
OrderedDict([('gene', ['insulin'])])
OrderedDict([('gene', ['insulin']), ('codon_start', ['1']), ('product', ['insulin']), ('protein_id', ['AAA40590.1']), ('translation', ['MAPWMHLLTVLALLALWGPNSVQAYSSQHLCGSNLVEALYMTCGRSGFYRPHDRRELEDLQVEQAELGLEAGGLQPSALEMILQKRGIVDQCCNNICTFNQLQNYCNVP'])])
OrderedDict([('gene', ['insulin'])])
OrderedDict([('gene', ['insulin']), ('product', ['insulin B-chain'])])
OrderedDict([('gene', ['insulin']), ('product', ['insulin C-peptide'])])
OrderedDict([('gene', ['insulin']), ('product', ['insulin A-chain'])])
OrderedDict([('regulatory_class', ['polyA_signal_sequence']), ('gene', ['insulin'])])
OrderedDict([('gene', ['insulin'])])

Feature를 딕셔너리 형태로 저장해둔 것. Octodon degus는 데덴네의 모티브가 되기도 한 데구의 학명이다. 

Position과 location

from Bio import SeqFeature
start_pos=SeqFeature.ExactPosition(15)
end_pos=SeqFeature.ExactPosition(30)
location=SeqFeature.FeatureLocation(start_pos, end_pos)
print(start_pos,end_pos,location)
15 30 [15:30]

start_pos2 = SeqFeature.AfterPosition(1)
end_pos2 = SeqFeature.BetweenPosition(9, left=8, right=9)
my_location = SeqFeature.FeatureLocation(start_pos2, end_pos2)
print(start_pos2,end_pos2,my_location)
>1 (8^9) [>1:(8^9)]

30쓰니까 에러뜸... 

포지션이 명확할 때도 있고 명확하지 않을 때도 있다. 무슨 양자역학이냐... 슈뢰딩거의 포지션

start_pos2 = SeqFeature.AfterPosition(1)
end_pos2 = SeqFeature.BeforePosition(8)
my_location = SeqFeature.FeatureLocation(start_pos2, end_pos2)
print(start_pos2,end_pos2,my_location)
>1 <8 [>1:<8]

참고로 애프터가 있다는 것은 비포가 있다는 얘기이기도 하다. (대충 펀쿨섹좌 짤)

실제 파일에 적용해보기

ID: NM_000517.6
Name: NM_000517
Description: Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA
Number of features: 44
/molecule_type=mRNA
/topology=linear
/data_file_division=PRI
/date=06-SEP-2021
/accessions=['NM_000517']
/sequence_version=6
/keywords=['RefSeq', 'MANE Select']
/source=Homo sapiens (human)
/organism=Homo sapiens
/taxonomy=['Eukaryota', 'Metazoa', 'Chordata', 'Craniata', 'Vertebrata', 'Euteleostomi', 'Mammalia', 'Eutheria', 'Euarchontoglires', 'Primates', 'Haplorrhini', 'Catarrhini', 'Hominidae', 'Homo']
/references=[Reference(title='Molecular and Hematological Analysis of Alpha- and Beta-Thalassemia in a Cohort of Mexican Patients', ...), Reference(title='Molecular Characterization and Hematological Aspects of Hb E-Myanmar [beta26(B8)Glu-->Lys and beta65(E9)Lys-->Asn, HBB: c.[79G>A;198G>C]): A Novel beta-Thalassemic Hemoglobin Variant', ...), Reference(title='delta-Globin Chain Variants Associated with Decreased Hb A2 Levels: A National Reference Laboratory Experience', ...), Reference(title='Heterozygosity for the Novel HBA2: c.*91_*92delTA Polyadenylation Site Variant on the alpha2-Globin Gene Expanding the Genetic Spectrum of alpha-Thalassemia in Iran', ...), Reference(title='Identification and characterization of two novel and differentially expressed isoforms of human alpha2- and alpha1-globin genes', ...), Reference(title='Alpha-Thalassemia', ...), Reference(title="Cloning and complete nucleotide sequence of human 5'-alpha-globin gene", ...), Reference(title='A new abnormal human hemoglobin: Hb Prato (alpha 2 31 (B12) Arg leads to Ser beta 2)', ...), Reference(title='Hemoglobin Suan-Dok (alpha 2 109 (G16) Leu replaced by Arg beta 2): an unstable variant associated with alpha-thalassemia', ...), Reference(title='Hemoglobin Tarrant: alpha126(H9) Asp leads to Asn. A new hemoglobin variant in the alpha1beta1 contact region showing high oxygen affinity and reduced cooperativity', ...)]
/comment=REVIEWED REFSEQ: This record has been curated by NCBI staff. The
reference sequence was derived from V00493.1, Z84721.1 and
AI016696.1.
This sequence is a reference standard in the RefSeqGene project.
On Aug 3, 2018 this sequence version replaced NM_000517.5.
Summary: The human alpha globin gene cluster located on chromosome
16 spans about 30 kb and includes seven loci: 5'- zeta - pseudozeta
- mu - pseudoalpha-1 - alpha-2 - alpha-1 - theta - 3'. The alpha-2
(HBA2) and alpha-1 (HBA1) coding sequences are identical. These
genes differ slightly over the 5' untranslated regions and the
introns, but they differ significantly over the 3' untranslated
regions. Two alpha chains plus two beta chains constitute HbA,
which in normal adult life comprises about 97% of the total
hemoglobin; alpha chains combine with delta chains to constitute
HbA-2, which with HbF (fetal hemoglobin) makes up the remaining 3%
of adult hemoglobin. Alpha thalassemias result from deletions of
each of the alpha genes as well as deletions of both HBA2 and HBA1;
some nondeletion alpha thalassemias have also been reported.
[provided by RefSeq, Jul 2008].
Publication Note:  This RefSeq record includes a subset of the
publications that are available for this gene. Please see the Gene
record to access additional publications.
COMPLETENESS: full length.
/structured_comment=OrderedDict([('Evidence-Data', OrderedDict([('Transcript exon combination', 'AI816040.1, AI815806.1 [ECO:0000332]'), ('RNAseq introns', 'single sample supports all introns SAMEA1965299, SAMEA1966682 [ECO:0000348]')])), ('RefSeq-Attributes', OrderedDict([('MANE Ensembl match', 'ENST00000251595.11/ ENSP00000251595.6'), ('RefSeq Select criteria', 'based on single protein-coding transcript')]))])
Seq('ACTCTTCTGGTCCCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGC...GCA')

오늘의 시범조교. 뭐야 wrap 돌려줘요.

record=SeqIO.read("/home/koreanraichu/sequence.gb",'genbank') # 여기서 뭘 뽑아봅시다
index=575
for feature in record.features:
    if index in feature:
        print("%s, %s, %s" %(feature.type,feature.location,feature.qualifiers.get("db_xref")))

여기서 575가 들어가는 feature를 찾아보자. 

source, [0:576](+), ['taxon:9606']
gene, [0:576](+), ['GeneID:3040', 'HGNC:HGNC:4824', 'MIM:141850']
exon, [337:576](+), None

아 그래서 이렇게 나온건가... (아직도 이해 못함) 

SeqFeature로 시퀀스의 특정 부분만 잘라서 보기

from Bio.Seq import Seq
from Bio.SeqFeature import SeqFeature, FeatureLocation
at5g40780=Seq("MVAQAPHDDHQDDEKLAAARQKEIEDWLPITSSRNAKWWYSAFHNVTAMVGAGVLGLPYAMSQLGWGPGIAVLVLSWVITLYTLWQMVEMHEMVPGKRFDRYHELGQHAFGEKLGLYIVVPQQLIVEIGVCIVYMVTGGKSLKKFHELVCDDCKPIKLTYFIMIFASVHFVLSHLPNFNSISGVSLAAAVMSLSYSTIAWASSASKGVQEDVQYGYKAKTTAGTVFNFFSGLGDVAFAYAGHNVVLEIQATIPSTPEKPSKGPMWRGVIVAYIVVALCYFPVALVGYYIFGNGVEDNILMSLKKPAWLIATANIFVVIHVIGSYQIYAMPVFDMMETLLVKKLNFRPTTTLRFFVRNFYVAATMFVGMTFPFFGGLLAFFGGFAFAPTTYFLPCVIWLAIYKPKKYSLSWWANWVCIVFGLFLMVLSPIGGLRTIVIQAKGYKFYS")
# 단백질 시퀀스 왜 이렇게 길어 이거
feature = SeqFeature(FeatureLocation(0, 30), type="protein", strand=1)
feature_seq = at5g40780[feature.location.start:feature.location.end]
print(feature_seq)

TAIR에서 돌아온 시범조교, LHT1의 단백질 시퀀스다. 뭐요. 파이참 wrap 일일이 켜기도 귀찮은데 걍 혀... 아무튼 feature라는 변수를 통해 30번째까지 볼 거라는 얘기다. strand가 -1이면 상보적인 시퀀스를 추출한다는데 저거 단백질이라 의미 없음. 

feature_seq2 = feature.extract(at5g40780)
print(feature_seq2)

Feature 변수 주고 extract해도 된다. 

레코드 비교하기

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
record1=SeqRecord(Seq("GAATTC"), id="ECoR1",description="sticky")
record2=SeqRecord(Seq("GGCC"), id="HaeIII",description="blunt")
record3=SeqRecord(Seq("CTCGAG"), id="XhoI",description="sticky")
record4=SeqRecord(Seq("GATC"), id="DpnII",description="sticky")

이번 시범조교는 얘네들이다. 

레코드 비교? 그럼 다이렉트로 비교 되나요? 그러면 NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest. 에러가 당신을 반겨줄 것이다. 레코드 비교는 안에 든 걸 비교하는거지 레코드를 쌩으로 비교하라는 얘기가 아니다. 

print(record1.description==record3.description)
True

EcoRI과 XhoI의 description이 같으므로 참이다.

print(len(record1.seq)==len(record2.seq))
False

EcoRI과 HaeIII의 시퀀스 길이는 다르므로 false. 

포맷 만들기

아니 컴퓨터 포맷 아님. 

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
record=SeqRecord(
    Seq(
        "ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTACTCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCCATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACACTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGATCGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTGCCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAATCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTATTTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCCATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGGGCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACAACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCGGGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCAAAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTCCCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATGTCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTCATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTCAAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTTGCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTTGGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATCTACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGTCTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAAGGATACAAGTTTTACTCATAA"),
id="at5g40780|LHT1",
description="Lysine-histidine transporter 1|Arabidopsis thaliana")
print(record.format("fasta"))

이런 식으로 입력하면 

>at5g40780|LHT1 Lysine-histidine transporter 1|Arabidopsis thaliana
ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGA
CAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTAC
TCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCC
ATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACA
CTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGAT
CGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTG
CCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAA
TCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTAT
TTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCC
ATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGG
GCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACA
ACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCG
GGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCA
AAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTC
CCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATG
TCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTC
ATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTC
AAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTT
GCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTT
GGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATC
TACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGT
CTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAA
GGATACAAGTTTTACTCATAA

이렇게 FASTA 포맷이 된다. 

슬라이싱

feature도 slicing이 된다. 어째서인지는 모르겠으나... 

from Bio import SeqIO
record = SeqIO.read("/home/koreanraichu/sequence.gb", "genbank")
print(record.features[20])

이건 근데 인덱싱 아니냐. 

sub_record = record[0:5]
print(sub_record)
ID: NM_000517.6
Name: NM_000517
Description: Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA
Number of features: 0
/molecule_type=mRNA
Seq('ACTCT')

아무튼 레코드에서 0부터 5까지 뽑아보았다. 염기 6개가 나온다. 

record = SeqIO.read("/home/koreanraichu/sequence.gb", "genbank")
sub_record = record[0:576]
print(sub_record.features[0])

범위를 풀로 바꾸고 0번째 feature를 봤다. 

type: source
location: [0:576](+)
qualifiers:
    Key: chromosome, Value: ['16']
    Key: db_xref, Value: ['taxon:9606']
    Key: map, Value: ['16p13.3']
    Key: mol_type, Value: ['mRNA']
    Key: organism, Value: ['Homo sapiens']

그러면 이렇게 나옵니다.

print(sub_record.format("genbank"))
print(sub_record.format("fasta"))

위 코드를 이용해 출력 포맷을 바꿀 수도 있다. 

>NM_000517.6 Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA
ACTCTTCTGGTCCCCACAGA
# FASTA

LOCUS       NM_000517                 20 bp    mRNA             UNK 01-JAN-1980
DEFINITION  Homo sapiens hemoglobin subunit alpha 2 (HBA2), mRNA.
ACCESSION   NM_000517
VERSION     NM_000517.6
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
ORIGIN
        1 actcttctgg tccccacaga
//
# Genbank

위는 FASTA, 아래는 Genbank 포맷. 

전에요? 아니 그건 가져온거고. 

시퀀스도 텍스트나 리스트처럼 인덱싱과 슬라이싱이 된다. 방법도 똑같다. 그래서 그건 생략할 예정임... 그리고 그거 아니어도 분량이 많아요 또. 

.count()

특정 문자열이 몇 개인지 세 준다. 뭐 시퀀스에서 아데닌이나 구아닌 갯수 세주고 그런거다. 이걸 이용하면 특정 단백질 시퀀스에서 특정 아미노산(카테고리면 겁나 세야된다...)이 차지하는 비율도 볼 수 있다. 

쉬어가는 코너-Primer에서 GC함량 구하기

아 이거 중요합니다. Primer 만들 때 중요한 요소 중 하나가 GC 함량이다. 참고로 여기서는 두 가지 방법으로 구해볼건데 첫번째는 .count()를 이용해 G와 C의 수를 세서 전체 DNA 염기 수로 나누는 거고, 두번째는 바이오파이썬의 모듈을 이용하는 방법이다. 

GC_cont=(primer.count('C')+primer.count('G'))/len(primer)
print(GC_cont*100,"%")

from Bio.SeqUtils import GC
primer=Seq('CAGCAAGCAAAGGTGTTCAA')
print(GC(primer),"%")

위는 직접 계산하는 방법이고 아래는 모듈을 사용했다. 모듈은 쓰려면 또 모셔와야 한다. 

시퀀스에 시퀀스를 더하면 1+1인가요 

시퀀스 그냥 +로 갖다 붙이거나 리스트업된 거 for문으로 붙이거나 join()쓰거나... 일단 for와 join 보고 갑시다. 

seq_list=[Seq('ATCC'),Seq('ATAT'),Seq('TNNN')]
add = Seq("")
for i in seq_list:
    add += i
print(add)

seq_list=[Seq('ATCC'),Seq('ATAT'),Seq('TNNN')]
spacer=Seq('N'*5)
print(spacer.join(seq_list))

돌연변이 만들기

시퀀스 데이터는 튜플마냥 안에 있는 내용물을 수정할 수 없다. 수정하려면 또 다른 모듈을 모셔와야 한다. 

from Bio.Seq import Seq
from Bio.Seq import MutableSeq #이걸 불러와서 
my_seq=Seq("GAATTC")
mutable_seq=MutableSeq(my_seq) #적용해줘야 한다
mutable_seq[0]="A"
print(mutable_seq)

Mutableseq에 던져둔 것은 타입도 

<class 'Bio.Seq.MutableSeq'>

로 바뀐다. 일반 시퀀스는 Seq. 

전사번역

입력한 시퀀스를 mRNA로 전사도 해 주고, 번역도 된다. 물론 다이렉트로 된다. 

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import GC
at5g40780=Seq('ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTACTCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCCATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACACTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGATCGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTGCCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAATCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTATTTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCCATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGGGCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACAACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCGGGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCAAAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTCCCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATGTCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTCATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTCAAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTTGCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTTGGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATCTACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGTCTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAAGGATACAAGTTTTACTCATAA')
at5g40780_rc=at5g40780.reverse_complement()
# DNA와 상보적인 sequence를 만들었다.

at5g40780_mrna=at5g40780.transcribe()
at5g40780_mrna_rc=at5g40780_mrna.reverse_complement()
# 여기 RNA 추가요!
print(at5g40780_mrna)
AUGGUAGCUCAAGCUCCUCAUGAUGAUCAUCAGGAUGAUGAGAAAUUAGCAGCAGCGAGACAAAAAGAGAUCGAAGAUUGGUUACCAAUUACUUCAUCAAGAAAUGCAAAGUGGUGGUACUCUGCUUUUCACAAUGUCACCGCCAUGGUCGGUGCCGGAGUUCUUGGACUCCCUUACGCCAUGUCUCAGCUCGGAUGGGGACCGGGAAUUGCAGUGUUGGUUUUGUCAUGGGUCAUAACACUAUACACAUUAUGGCAAAUGGUGGAAAUGCAUGAAAUGGUUCCUGGAAAGCGUUUUGAUCGUUACCAUGAGCUCGGACAACACGCGUUUGGAGAAAAACUCGGUCUUUAUAUCGUUGUGCCGCAACAAUUGAUCGUUGAAAUCGGUGUUUGCAUCGUUUAUAUGGUCACUGGAGGCAAAUCUUUAAAGAAAUUUCAUGAGCUUGUUUGUGAUGAUUGUAAACCAAUCAAGCUUACUUAUUUCAUCAUGAUCUUUGCUUCGGUUCACUUCGUCCUCUCUCAUCUUCCUAAUUUCAAUUCCAUCUCCGGCGUUUCUCUUGCUGCUGCCGUUAUGUCUCUCAGCUACUCAACAAUCGCAUGGGCAUCAUCAGCAAGCAAAGGUGUUCAAGAAGACGUUCAAUACGGUUACAAAGCGAAAACAACAGCCGGUACGGUUUUCAAUUUCUUCAGCGGUUUAGGUGAUGUGGCAUUUGCUUACGCGGGUCAUAAUGUGGUCCUUGAGAUCCAAGCAACUAUCCCUUCAACGCCUGAGAAACCCUCAAAAGGUCCCAUGUGGAGAGGAGUCAUCGUUGCUUACAUCGUCGUAGCGCUCUGUUAUUUCCCCGUGGCUCUCGUUGGAUACUACAUUUUCGGGAACGGAGUCGAAGAUAAUAUUCUCAUGUCACUUAAGAAACCGGCGUGGUUAAUCGCCACGGCGAACAUCUUCGUCGUGAUCCAUGUCAUUGGUAGUUACCAGAUAUAUGCAAUGCCGGUAUUUGAUAUGAUGGAGACUUUAUUGGUCAAGAAGCUAAAUUUCAGACCAACCACAACUCUUCGGUUCUUUGUCCGUAAUUUCUACGUUGCUGCAACUAUGUUUGUUGGUAUGACGUUUCCGUUCUUCGGUGGGCUUUUGGCGUUCUUUGGUGGAUUCGCGUUUGCCCCAACCACAUACUUCCUCCCUUGCGUCAUUUGGUUAGCCAUCUACAAACCCAAGAAAUACAGCUUGUCUUGGUGGGCCAACUGGGUAUGCAUCGUGUUUGGUCUUUUCUUGAUGGUCCUAUCGCCAAUUGGAGGGCUAAGGACAAUCGUUAUUCAAGCAAAAGGAUACAAGUUUUACUCAUAA

.transcribe()를 이용하면 전사도 해 준다. 예시에 있는 시퀀스는 at5g40780(의 CDS). 

from Bio.Seq import Seq
at5g40780=Seq('ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTACTCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCCATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACACTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGATCGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTGCCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAATCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTATTTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCCATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGGGCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACAACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCGGGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCAAAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTCCCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATGTCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTCATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTCAAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTTGCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTTGGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATCTACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGTCTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAAGGATACAAGTTTTACTCATAA')
# DNA sequence

at5g40780_mrna=at5g40780.transcribe()
# Transcription
print(at5g40780_mrna)
AUGGUAGCUCAAGCUCCUCAUGAUGAUCAUCAGGAUGAUGAGAAAUUAGCAGCAGCGAGACAAAAAGAGAUCGAAGAUUGGUUACCAAUUACUUCAUCAAGAAAUGCAAAGUGGUGGUACUCUGCUUUUCACAAUGUCACCGCCAUGGUCGGUGCCGGAGUUCUUGGACUCCCUUACGCCAUGUCUCAGCUCGGAUGGGGACCGGGAAUUGCAGUGUUGGUUUUGUCAUGGGUCAUAACACUAUACACAUUAUGGCAAAUGGUGGAAAUGCAUGAAAUGGUUCCUGGAAAGCGUUUUGAUCGUUACCAUGAGCUCGGACAACACGCGUUUGGAGAAAAACUCGGUCUUUAUAUCGUUGUGCCGCAACAAUUGAUCGUUGAAAUCGGUGUUUGCAUCGUUUAUAUGGUCACUGGAGGCAAAUCUUUAAAGAAAUUUCAUGAGCUUGUUUGUGAUGAUUGUAAACCAAUCAAGCUUACUUAUUUCAUCAUGAUCUUUGCUUCGGUUCACUUCGUCCUCUCUCAUCUUCCUAAUUUCAAUUCCAUCUCCGGCGUUUCUCUUGCUGCUGCCGUUAUGUCUCUCAGCUACUCAACAAUCGCAUGGGCAUCAUCAGCAAGCAAAGGUGUUCAAGAAGACGUUCAAUACGGUUACAAAGCGAAAACAACAGCCGGUACGGUUUUCAAUUUCUUCAGCGGUUUAGGUGAUGUGGCAUUUGCUUACGCGGGUCAUAAUGUGGUCCUUGAGAUCCAAGCAACUAUCCCUUCAACGCCUGAGAAACCCUCAAAAGGUCCCAUGUGGAGAGGAGUCAUCGUUGCUUACAUCGUCGUAGCGCUCUGUUAUUUCCCCGUGGCUCUCGUUGGAUACUACAUUUUCGGGAACGGAGUCGAAGAUAAUAUUCUCAUGUCACUUAAGAAACCGGCGUGGUUAAUCGCCACGGCGAACAUCUUCGUCGUGAUCCAUGUCAUUGGUAGUUACCAGAUAUAUGCAAUGCCGGUAUUUGAUAUGAUGGAGACUUUAUUGGUCAAGAAGCUAAAUUUCAGACCAACCACAACUCUUCGGUUCUUUGUCCGUAAUUUCUACGUUGCUGCAACUAUGUUUGUUGGUAUGACGUUUCCGUUCUUCGGUGGGCUUUUGGCGUUCUUUGGUGGAUUCGCGUUUGCCCCAACCACAUACUUCCUCCCUUGCGUCAUUUGGUUAGCCAUCUACAAACCCAAGAAAUACAGCUUGUCUUGGUGGGCCAACUGGGUAUGCAUCGUGUUUGGUCUUUUCUUGAUGGUCCUAUCGCCAAUUGGAGGGCUAAGGACAAUCGUUAUUCAAGCAAAAGGAUACAAGUUUUACUCAUAA

.translate()는 번역이다. (*은 종결코돈)

전사는 그냥 뭐 전사했구나 하고 넘어가도 되는데, 번역은 옵션이 좀 있으니 일단 보고 가자. 

to_stop=True: 처음 만나는 종결코돈에서 멈춘다

from Bio.Seq import Seq
at5g40780_mrna=Seq('AUGGUAGCUCAAGCUCCUCAUGAUGAUCAUCAGGAUGAUGAGAAAUUAGCAGCAGCGAGACAAAAAGAGAUCGAAGAUUGGUUACCAAUUACUUCAUCAAGAAAUGCAAAGUGGUGGUACUCUGCUUUUCACAAUGUCACCGCCAUGGUCGGUGCCGGAGUUCUUGGACUCCCUUACGCCAUGUCUCAGCUCGGAUGGGGACCGGGAAUUGCAGUGUUGGUUUUGUCAUGGGUCAUAACACUAUACACAUUAUGGCAAAUGGUGGAAAUGCAUGAAAUGGUUCCUGGAAAGCGUUUUGAUCGUUACCAUGAGCUCGGACAACACGCGUUUGGAGAAAAACUCGGUCUUUAUAUCGUUGUGCCGCAACAAUUGAUCGUUGAAAUCGGUGUUUGCAUCGUUUAUAUGGUCACUGGAGGCAAAUCUUUAAAGAAAUUUCAUGAGCUUGUUUGUGAUGAUUGUAAACCAAUCAAGCUUACUUAUUUCAUCAUGAUCUUUGCUUCGGUUCACUUCGUCCUCUCUCAUCUUCCUAAUUUCAAUUCCAUCUCCGGCGUUUCUCUUGCUGCUGCCGUUAUGUCUCUCAGCUACUCAACAAUCGCAUGGGCAUCAUCAGCAAGCAAAGGUGUUCAAGAAGACGUUCAAUACGGUUACAAAGCGAAAACAACAGCCGGUACGGUUUUCAAUUUCUUCAGCGGUUUAGGUGAUGUGGCAUUUGCUUACGCGGGUCAUAAUGUGGUCCUUGAGAUCCAAGCAACUAUCCCUUCAACGCCUGAGAAACCCUCAAAAGGUCCCAUGUGGAGAGGAGUCAUCGUUGCUUACAUCGUCGUAGCGCUCUGUUAUUUCCCCGUGGCUCUCGUUGGAUACUACAUUUUCGGGAACGGAGUCGAAGAUAAUAUUCUCAUGUCACUUAAGAAACCGGCGUGGUUAAUCGCCACGGCGAACAUCUUCGUCGUGAUCCAUGUCAUUGGUAGUUACCAGAUAUAUGCAAUGCCGGUAUUUGAUAUGAUGGAGACUUUAUUGGUCAAGAAGCUAAAUUUCAGACCAACCACAACUCUUCGGUUCUUUGUCCGUAAUUUCUACGUUGCUGCAACUAUGUUUGUUGGUAUGACGUUUCCGUUCUUCGGUGGGCUUUUGGCGUUCUUUGGUGGAUUCGCGUUUGCCCCAACCACAUACUUCCUCCCUUGCGUCAUUUGGUUAGCCAUCUACAAACCCAAGAAAUACAGCUUGUCUUGGUGGGCCAACUGGGUAUGCAUCGUGUUUGGUCUUUUCUUGAUGGUCCUAUCGCCAAUUGGAGGGCUAAGGACAAUCGUUAUUCAAGCAAAAGGAUACAAGUUUUACUCAUAA')
# 이건 mRNA 시퀀스고
at5g40780_protein=at5g40780_mrna.translate(to_stop=True)
#번역했다.
print(at5g40780_protein)
MVAQAPHDDHQDDEKLAAARQKEIEDWLPITSSRNAKWWYSAFHNVTAMVGAGVLGLPYAMSQLGWGPGIAVLVLSWVITLYTLWQMVEMHEMVPGKRFDRYHELGQHAFGEKLGLYIVVPQQLIVEIGVCIVYMVTGGKSLKKFHELVCDDCKPIKLTYFIMIFASVHFVLSHLPNFNSISGVSLAAAVMSLSYSTIAWASSASKGVQEDVQYGYKAKTTAGTVFNFFSGLGDVAFAYAGHNVVLEIQATIPSTPEKPSKGPMWRGVIVAYIVVALCYFPVALVGYYIFGNGVEDNILMSLKKPAWLIATANIFVVIHVIGSYQIYAMPVFDMMETLLVKKLNFRPTTTLRFFVRNFYVAATMFVGMTFPFFGGLLAFFGGFAFAPTTYFLPCVIWLAIYKPKKYSLSWWANWVCIVFGLFLMVLSPIGGLRTIVIQAKGYKFYS

table=2: 아래 사이트에서 제공하는 코돈 테이블로 번역할 수 있다. (기본은 1번 스탠다드)
https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi

at5g40780_protein=at5g40780_mrna.translate(table=2)
MVAQAPHDDHQDDEKLAAA*QKEIEDWLPITSS*NAKWWYSAFHNVTAMVGAGVLGLPYAMSQLGWGPGIAVLVLSWVMTLYTLWQMVEMHEMVPGKRFDRYHELGQHAFGEKLGLYIVVPQQLIVEIGVCIVYMVTGGKSLKKFHELVCDDCKPIKLTYFIMIFASVHFVLSHLPNFNSISGVSLAAAVMSLSYSTIAWASSASKGVQEDVQYGYKAKTTAGTVFNFFSGLGDVAFAYAGHNVVLEIQATIPSTPEKPSKGPMW*GVIVAYIVVALCYFPVALVGYYIFGNGVEDNILMSLKKPAWLIATANIFVVIHVIGSYQMYAMPVFDMMETLLVKKLNF*PTTTLRFFVRNFYVAATMFVGMTFPFFGGLLAFFGGFAFAPTTYFLPCVIWLAIYKPKKYSLSWWANWVCIVFGLFLMVLSPIGGL*TIVIQAKGYKFYS*

위 CDS를 table 2(Vertebrate mitochondria)로 설정하고 번역한 결과. *은 종결 코돈이다. 

stop_symbol="@": 종결코돈의 기호를 @로 바꾼다. 

at5g40780_protein=at5g40780_mrna.translate(table=2,stop_symbol="-")
MVAQAPHDDHQDDEKLAAA-QKEIEDWLPITSS-NAKWWYSAFHNVTAMVGAGVLGLPYAMSQLGWGPGIAVLVLSWVMTLYTLWQMVEMHEMVPGKRFDRYHELGQHAFGEKLGLYIVVPQQLIVEIGVCIVYMVTGGKSLKKFHELVCDDCKPIKLTYFIMIFASVHFVLSHLPNFNSISGVSLAAAVMSLSYSTIAWASSASKGVQEDVQYGYKAKTTAGTVFNFFSGLGDVAFAYAGHNVVLEIQATIPSTPEKPSKGPMW-GVIVAYIVVALCYFPVALVGYYIFGNGVEDNILMSLKKPAWLIATANIFVVIHVIGSYQMYAMPVFDMMETLLVKKLNF-PTTTLRFFVRNFYVAATMFVGMTFPFFGGLLAFFGGFAFAPTTYFLPCVIWLAIYKPKKYSLSWWANWVCIVFGLFLMVLSPIGGL-TIVIQAKGYKFYS-

위 예시에서는 -로 바꿨다. 

글자 다이렉트로 전사번역하기

시퀀스도 글자긴 한데 str과 다른 점이 있다면 상보적인 시퀀스를 만들 수 있다는 것이다. 근데 그렇다고 헐 시퀀스 안만들었다! 망했다! 이건 아니고... 

from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate #질문있는데 이걸 꼭 이렇게 불러야겠니
insulin="GCATTCTGAGGCATTCTCTAACAGGTTCTCGACCCTCCGCCATGGCCCCGTGGATGCATCTCCTCACCGTGCTGGCCCTGCTGGCCCTCTGGGGACCCAACTCTGTTCAGGCCTATTCCAGCCAGCACCTGTGCGGCTCCAACCTAGTGGAGGCACTGTACATGACATGTGGACGGAGTGGCTTCTATAGACCCCACGACCGCCGAGAGCTGGAGGACCTCCAGGTGGAGCAGGCAGAACTGGGTCTGGAGGCAGGCGGCCTGCAGCCTTCGGCCCTGGAGATGATTCTGCAGAAGCGCGGCATTGTGGATCAGTGCTGTAATAACATTTGCACATTTAACCAGCTGCAGAACTACTGCAATGTCCCTTAGACACCTGCCTTGGGCCTGGCCTGCTGCTCTGCCCTGGCAACCAATAAACCCCTTGAATGAG"
#Wrap 알아서 켜주면 안되냐고...
insulin_rc=reverse_complement(insulin)
insulin_mrna=transcribe(insulin)
insulin_prot=translate(insulin)
# 순서대로 reverse complement/전사/번역 
print(insulin_rc)
print(insulin_mrna)
print(insulin_prot)

시퀀스때랑 불러오는 건 다르지만 어쨌든 된다. 

코돈 테이블 보기

내 블로그 아미노산 배경화면에도 있긴 한데 그걸 누가 외웁니까... 아무튼 그렇다. 

from Bio.Data import CodonTable
standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
print(standard_table)

from Bio.Data import CodonTable
standard_table = CodonTable.unambiguous_dna_by_id[1]
print(standard_table)

둘 다 같은 테이블을 호출하는건데, 위에껀 이름으로 부르고 밑에껀 ID로 부른다. 이름이나 ID는 위에 올린 NCBI 주소로 가면 나오니까 그걸로 부르면 된다. 

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
--+---------+---------+---------+---------+--

아무튼 소환했다. 

색인

저거 검색도 된다. 실화다. 

print(standard_table.stop_codons)
['TAA', 'TAG', 'TGA']

해당 테이블의 Stop codon을 검색하거나 

print(mito_table.start_codons)
['ATT', 'ATC', 'ATA', 'ATG', 'GTG']

개시코돈을 검색하거나(저거 2번 테이블이다)

print(standard_table.forward_table['AAA'])
K

특정 코돈이 지정하는 아미노산을 찾을 수 있다. 

참고로 다른 코딩글은 여기다가 안 올린다. 코드블럭도 없고 카테고리도 없어서 어지간한 코딩 이야기는 노션과 워드프레스에 올리는 중... 아니면 미디움이나. 근데 이건 우째도 올렸네? 아 생물학 카테고리가 있잖아요. 

Biopython은 생물정보학에서 쓰는건데 마침 심심하던 차에 잘됐다 써봐야징 해서 어제 깔았다. 리눅스의 경우 pip install bio(안되면 biopython)으로 걍 터미널에서 깔면 된다(우분투 20.02 LTS 기준). 윈도우마냥 pip 있는 데 안 찾아가도 됨 ㄹㅇ 편함. 근데 파이참 바로가기 없는 건 너무한 거 아니냐고... 내가 터미널 열고 직접 모시러 가야것냐... 

*Notes: Biopython을 사용하려면 사용할 기능에 맞는 모듈을 모셔와야 한다. 코드에 그게 생략되는 경우가 있는데(부분삽입의 경우) from ~ import ~가 생략된 것이다. 

일단 오늘의 시범조교를 모셔보자. 

ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTACTCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCCATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACACTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGATCGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTGCCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAATCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTATTTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCCATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGGGCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACAACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCGGGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCAAAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTCCCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATGTCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTCATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTCAAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTTGCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTTGGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATCTACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGTCTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAAGGATACAAGTTTTACTCATAA

At5g40780의 CDS 시퀀스이다. (코드블럭은 있는데 wrap이 안돼서 걍 생으로 붙였음)

이걸 시퀀스 오브젝트로 생성하는 방법도 간단하다. 그냥 모듈 모셔온 다음 Seq=("")로 감싸버리면 된다. 

at5g40780=Seq("ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGTGGTACTCTGCTTTTCACAATGTCACCGCCATGGTCGGTGCCGGAGTTCTTGGACTCCCTTACGCCATGTCTCAGCTCGGATGGGGACCGGGAATTGCAGTGTTGGTTTTGTCATGGGTCATAACACTATACACATTATGGCAAATGGTGGAAATGCATGAAATGGTTCCTGGAAAGCGTTTTGATCGTTACCATGAGCTCGGACAACACGCGTTTGGAGAAAAACTCGGTCTTTATATCGTTGTGCCGCAACAATTGATCGTTGAAATCGGTGTTTGCATCGTTTATATGGTCACTGGAGGCAAATCTTTAAAGAAATTTCATGAGCTTGTTTGTGATGATTGTAAACCAATCAAGCTTACTTATTTCATCATGATCTTTGCTTCGGTTCACTTCGTCCTCTCTCATCTTCCTAATTTCAATTCCATCTCCGGCGTTTCTCTTGCTGCTGCCGTTATGTCTCTCAGCTACTCAACAATCGCATGGGCATCATCAGCAAGCAAAGGTGTTCAAGAAGACGTTCAATACGGTTACAAAGCGAAAACAACAGCCGGTACGGTTTTCAATTTCTTCAGCGGTTTAGGTGATGTGGCATTTGCTTACGCGGGTCATAATGTGGTCCTTGAGATCCAAGCAACTATCCCTTCAACGCCTGAGAAACCCTCAAAAGGTCCCATGTGGAGAGGAGTCATCGTTGCTTACATCGTCGTAGCGCTCTGTTATTTCCCCGTGGCTCTCGTTGGATACTACATTTTCGGGAACGGAGTCGAAGATAATATTCTCATGTCACTTAAGAAACCGGCGTGGTTAATCGCCACGGCGAACATCTTCGTCGTGATCCATGTCATTGGTAGTTACCAGATATATGCAATGCCGGTATTTGATATGATGGAGACTTTATTGGTCAAGAAGCTAAATTTCAGACCAACCACAACTCTTCGGTTCTTTGTCCGTAATTTCTACGTTGCTGCAACTATGTTTGTTGGTATGACGTTTCCGTTCTTCGGTGGGCTTTTGGCGTTCTTTGGTGGATTCGCGTTTGCCCCAACCACATACTTCCTCCCTTGCGTCATTTGGTTAGCCATCTACAAACCCAAGAAATACAGCTTGTCTTGGTGGGCCAACTGGGTATGCATCGTGTTTGGTCTTTTCTTGATGGTCCTATCGCCAATTGGAGGGCTAAGGACAATCGTTATTCAAGCAAAAGGATACAAGTTTTACTCATAA")
print(at5g40780)

이런 식으로. 

print(type(at5g40780))
<class 'Bio.Seq.Seq'>

누가봐도 시퀀스다. 

엥 근데 저거 걍 텍스트랑 뭔 차이임? Biopython을 통해 시퀀스로 생성한 데이터는 complement sequence와 reverse complement sequence를 만들 수 있다. 즉, 시퀀스로 불러온 CDS에 상보적인 시퀀스를 명령어 한 방에 만들 수 있는 것이다. 텍스트는 저런거 안된다. 

from Bio.Seq import Seq
at5g40780=Seq("ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGT")
print("Sequence",at5g40780)
print("Complement sequence",at5g40780.complement())
print("Reverse complement sequence",at5g40780.reverse_complement())

이렇게 입력하면 

Sequence ATGGTAGCTCAAGCTCCTCATGATGATCATCAGGATGATGAGAAATTAGCAGCAGCGAGACAAAAAGAGATCGAAGATTGGTTACCAATTACTTCATCAAGAAATGCAAAGTGGT
Complement sequence TACCATCGAGTTCGAGGAGTACTACTAGTAGTCCTACTACTCTTTAATCGTCGTCGCTCTGTTTTTCTCTAGCTTCTAACCAATGGTTAATGAAGTAGTTCTTTACGTTTCACCA
Reverse complement sequence ACCACTTTGCATTTCTTGATGAAGTAATTGGTAACCAATCTTCGATCTCTTTTTGTCTCGCTGCTGCTAATTTCTCATCATCCTGATGATCATCATGAGGAGCTTGAGCTACCAT

짜자잔 

complement는 CDS 시퀀스에 그냥 상보적인거고, reverse complement는 그 상보적인 시퀀스 순서를 뒤집은 모양이다. 

아니 근데 DNA가 겁나 바이트가 길어요... 저 시범조교 CDS도 1300자가 넘는데 어느 세월에 그거 복붙할거임... 그래서 FASTA와 Genbank의 데이터를 직접 열어볼 것이다. 

FASTA 시범조교는 PDB ID 2B10, cytochrome c peroxiase이다. 

>2B10_1|Chains A, C|Cytochrome c peroxidase, mitochondrial|Saccharomyces cerevisiae (4932)
TTPLVHVASVEKGRSYEDFQKVYNAIALKLREDDEYDNYIGYGPVLVRLAWHTSGTWDKHDNTGGSYGGTYRFKKEFNDPSNAGLQNGFKFLEPIHKEFPWISSGDLFSLGGVTAVQEMQGPKIPWRCGRVDTPEDTTPDNGRLPDADKDADYVRTFFQRLNMNDREVVALMGAHALGKTHLKNSGYEGPWGAANNVFTNEFYLNLLNEDWKLEKNDANNEQWDSKSGYMMLPTDYSLIQDPKYLSIVKEYANDQDKFFKDFSKAFEKLLENGITFPKDAPSPFIFKTLEEQGL
>2B10_2|Chains B, D|Cytochrome c iso-1|Saccharomyces cerevisiae (4932)
TEFKAGSAKKGATLFKTRCLQCHTVEKGGPHKVGPNLHGIFGRHSGQAEGYSYTDANIKKNVLWDENNMSEYLTNPKKYIPGTKMASGGLKKEKDRNDLITYLKKACE

FASTA 파일 형식은 이렇게 생겼다. 일단 저 파일을 어떻게 가져올거냐... PDB 가서 2B10으로 검색하고 FASTA 파일을 받아보자. 

for seq_record in SeqIO.parse("/home/koreanraichu/rcsb_pdb_2B10.fasta","fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))
2B10_1|Chains
Seq('TTPLVHVASVEKGRSYEDFQKVYNAIALKLREDDEYDNYIGYGPVLVRLAWHTS...QGL')
294
2B10_2|Chains
Seq('TEFKAGSAKKGATLFKTRCLQCHTVEKGGPHKVGPNLHGIFGRHSGQAEGYSYT...ACE')
108

FASTA 파일은 이런 식으로 가져올 수 있다. (seqIO 데려오면 된다) 예? 경로요? 아니 이거 리눅슨디? 

사실 분량때문에 넣지는 못하지만 

print(seq_record)

만 치면 전체 파일을 볼 수 있다. 처음 하시는 분들은 전체를 보고 저 예제를 따라해보는 걸 추천한다. 

for seq_record in SeqIO.parse("/home/koreanraichu/sequence.gb","genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

Genbank에서 뭐 검색하고 그거 gk(gbk)파일로 받으면 그게 genbank 파일이다. 이것도 일단 다운받아서 불러와야 한다. 웹에 있는거요? 그거 바로 못갖고온다. 

아무튼 저 코드를 실행하면 

NM_000517.6
Seq('ACTCTTCTGGTCCCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGC...GCA')
576

이렇게 된다. 그리고 얘도 FASTA와 마찬가지로 

print(seq_record)

를 주면 전체 파일을 볼 수는 있는데 분량이 어유...

문제

https://www.acmicpc.net/problem/10989

카운팅 정렬로 숫자 정렬하는 문제. 그러나 이 문제에는 함정카드가 하나 있다. (대충 함정좌 짤)

메모리 제한이 8MB밖에 안된다. 자바랑 코틀린만 많이 준다.

Reference

https://8iggy.tistory.com/123

https://seongonion.tistory.com/130

https://www.acmicpc.net/board/view/26132

https://scarlett-choi.tistory.com/9

풀이

문제를 보면 아시겠지만 이전의 정렬 문제와 달리 숫자가 중복이 된다. 우리야 숫자 같으니까 똑같이 두지 뭐 하지만 컴퓨터는 이게 뭐여 하면서 밥상뒤집기를 시전할 수도 있다. 그래서 이럴 때 필요한 게 카운팅 정렬임… 함정 얘기는 이따가 해드릴게여.

a = [1, 2, 6, 5, 6, 6, 4, 3, 3, 4, 2]

def count_sort(arr):
    m = max(arr)
    # 최댓값 도출
    C = [0] * (m + 1)
    # Counting Array 생성
    for a in arr:
        C[a] += 1
        print("Counting Array: {}".format(C))
        # counting Array에 채워넣는다
    for i in range(1, m + 1):
        C[i] += C[i - 1]
        print("Array's Sum: {}".format(C))
        # 누적합으로 바꾼다
    result = [0] * len(arr)
    # 정렬을 위한 배열
    for a in arr:
        result[C[a] - 1] = a
        C[a] -= 1
        print("Sorting... :{}".format(result))
        # 정렬 콜
    return result

print(count_sort(a))

이게 카운팅 정렬의 코드이다. 뭐야 최댓값이 왜 나와요 저 배열 뭔데요…

카운팅 정렬이 진행되는 과정에서는 배열이 두 개 필요하다. 첫번째가 Counting Array, 두번째는 결과를 만들기 위한 배열이다. Counting Array는 말 그대로 이 배열에 어떤 원소가 몇 개 있느냐를 표시하기 위한 배열로, 배열 전체를 돌면서 원소를 세고 있을때마다 카운트를 하나씩 올리는 방식이다. 이렇게 배열 내 원소를 전부 셌다면, 그 다음은 누적합(개수의 누적합)을 구해 배열화하게 된다. 근데 이것만 하면 정렬이 되나요? 된다.

누적합 배열을 통해 개수를 유추할 수 있으니, 어디에서 어디까지 어떤 원소가 들어가게 되는 지 알 수 있기 때문에 알아서 채우면 된다.

그럼 저거 내면 되나요?

메모리 초과: YOU JUST ACTIVATED MY TRAP CARD (브금은 셀프)

위에서 자바랑 코틀린 말고 메모리가 8메가라고 했는데, 저거 일일이 받아서 배열에 저장하면 초과 뜬다.

import sys
N = int(sys.stdin.readline())
array = [0] * 10001

for i in range(N):
    input_num = int(sys.stdin.readline())
    array[input_num] = array[input_num] + 1

for i in range(10001): 
    if array[i] != 0: 
        for j in range(array[i]): 
            print(i)

그래서 이거 내야됨. 아닛 이 짤막한 코드로 된다고? ㅇㅇ 됨 내가 봤음.

위에서 카운팅 정렬은 입력받은 배열에서 세기 위한 배열과 결과를 저장할 배열 두 개를 만든다고 했는데, 요 작고 소듕한 코드는 그 만드는 과정을 생략해버렸다. 

그래서 만드는 배열은 0이 10001개인 배열 하나뿐이고, 숫자가 입력될때마다 배열의 입력받은 값에 해당하는 인덱스의 숫자를 하나 추가한다. 이게 뭔 소리냐면 배열이 다 0이고 입력값이 7이면 배열[7]을 0에서 1로 바꾼다는 얘기. 그리고 배열을 싹 돌아서 안에 0이 아닌 숫자가 있으면 그만큼 해당 숫자를 출력한다.

문제

https://www.acmicpc.net/problem/2750

주어진 수를 오름차순으로 정렬하기(2751번도 같은 문제인데 버블/선택/삽입으로는 시간초과 뜬다)

풀이

정렬 알고리즘과 관련된 이론적인 설명은 아래를 참고할 것.

https://koreanraichu.sfuhost.com/2022/6650/

위 글에서 다섯가지 대표적인 정렬 알고리즘을 봤는데, 문제가 하나 있었다. ‘버블 알고리즘은 배열의 길이가 커질수록 시간복잡도가 무지막지하게 증가한다’는 것. 근데 시간제한이 1초여… 그래서 선택 알고리즘으로 했다. 

import sys
N = int(sys.stdin.readline())
array = []

for i in range(N):
    i = int(sys.stdin.readline())
    array.append(i)

입력은 들어오는 만큼 받아서 리스트에 넣으면 된다.

import sys
N = int(sys.stdin.readline())
array_in = []

for i in range(N):
    i = int(sys.stdin.readline())
    array_in.append(i)

def selection_sort(array):
	n = len(array)
	for i in range(n):
		min_index = i
		for j in range(i + 1, n):
			if array[j] < array[min_index]:
				min_index = j
		array[i], array[min_index] =  array[min_index], array[i]
	return array

array_out = selection_sort(array_in)

for i in array_out:
    print(i)

선택 정렬을 함수로 정의하고 돌려서 출력을 하면 되는’데’… 출력이 한줄에 하나씩 나와야 하잖음? 그래서 배열을 통으로 출력하면 안되고 for문으로 출력해야 한다.

2751번도 같은 문제인데, 버블/선택/삽입 정렬로 하면 시간초과 뜬다. 병합이나 퀵을 쓰자.

문제

https://www.acmicpc.net/problem/1436

Reference

https://hongcoding.tistory.com/108

https://pearlluck.tistory.com/523

풀이

실로 골때리는데 해결방안이 정말 의외인 문제. 아니 농담 아니고 진짜 그렇다니까?

참고로 666 다음이 1666 2666 3666 4666 5666이니까 이 다음은 6666이 올 것 같은가? 놉. 6660 6661 6662… 중간에 순서가 꼬인다. 근데 그때까지 영화가 안 망할 수 있을까 예제에 187 500 나오던데 여기까지 찍기 전에 감독님 죽어요… 

아무튼 처음에 시도했던 건 이거였다.

1. 1부터 6660000(10000이 제일 큰 수)까지 중 666이 들어가는 수를
2. 배열에 넣고
3. 배열의 N-1번째 요소를 반환한다

그래서 이케 했지.

import sys
N = int(sys.stdin.readline())
terminate_list =[]
terminate_end = (666 * 10000) + 1

for i in range(1,terminate_end):
    if str(i).find('666') != -1:
        terminate_list.append(i)

print(terminate_list[N-1])

근데 IndexError가 날 반길 줄은 몰랐지. 그래서 방법을 찾아봤는데 어떻게 하느냐…

그냥 1부터 쭉 더해가면서 666이 들어가는 수를 찾으면 됨.

import sys
N = int(sys.stdin.readline())
count = 0
num = 1
while True: 
    if '666' in str(num):
        count += 1
    if count == N:
        print(num)
        break;
    num += 1

농담 아니고 진짜다.

문제

https://www.acmicpc.net/problem/1018

보드에서 최소한으로 재도색하면서 체스판을 만들 수 있는 가짓수는?

Reference

https://bambbang00.tistory.com/43

풀이

체스판은 검흰검흰이거나 흰검흰검으로 반복된다. …앞에앞에 문제에 기물 주운 애랑 합치면 체스 세트 완성인데? 아무튼. 근데 이제 브루트 포스맛이면 8*8을 다 훑어보면서 체크하게 된다.

import sys
N, M = map(int, sys.stdin.readline().split())
board = []
cut_board = []
for _ in range(N):
    row = sys.stdin.readline().strip()
    board.append(row)
# 전체 보드 만들기

for m in range(N-7):
    for n in range(M-7):
        board_index_1 = 0
        board_index_2 = 0
        for o in range(m, m+8):
            for p in range(n, n+8):
                if (o + p) % 2 == 0:
                    if board[o][p] != 'W':
                        board_index_1 += 1
                    if board[o][p] != 'B': 
                        board_index_2 += 1
                else: 
                    if board[o][p] != 'B':
                        board_index_1 += 1
                    if board[o][p] != 'W':
                        board_index_2 += 1
        cut_board.append(min(board_index_1,board_index_2))

print(min(cut_board))

이게 전체 코드. 보드는 원래 보드(줍줍한 보드), 컷보드는 8*8로 자른 상태이다. 근데 뭔 for문이 이렇게 많냐고? 일단 위에 있는거 말고 요 두놈을 보자.

1. for m in range/for n in range
2. for o in range/for p in range

두 for문은 각각

1. 체스보드의 시작점(8*8로 자르는거지 주어진 게 8*8이 아님)
2. 시작점~시작점+8(range는 a 이상 b 미만)까지를 range로 해서 확인

이렇게 되어 있다. 위에도 썼지만 체스판이 8*8이니까 그걸로 잘라서 만드는거지 주어지는 판이 8*8이라는 법은 없다. 그럼 밑에 If가 네갈래인 이유는요?

if (o + p) % 2 == 0:
    if board[o][p] != 'W':
    	board_index_1 += 1
	if board[o][p] != 'B':
    	board_index_2 += 1
else :
    if board[o][p] != 'B':
    	board_index_1 += 1
	if board[o][p] != 'W':
    	board_index_2 += 1

체스판의 시작점이 검정색일때와 흰색일 때, 각각 인접해야 하는 색이 다르다. 즉 흰색 옆에는 검정, 검정 옆에는 흰색이 와야 해서 if가 저렇게 갈라졌다고 보면 된다. 

문제

https://www.acmicpc.net/problem/25304

구매한 각 물건의 가격과 총합이 일치하는지 확인해야 한다.

풀이

야 우리집에서 이마트 털러 가도 26만원어치는 못사… 게임팩을 털어왔나 아무튼 그래요… 가랑비에 옷 젖는 줄 모른다는 말이 이럴 때 쓰는 말임…

아무튼 로직이 2단인데

1. 다 더한다
2. 비교한다

이게 다다.

import sys
total_price = int(sys.stdin.readline())
what_buy = int(sys.stdin.readline())
# 총합과 물건 종류 수
if_total = 0
for i in range(what_buy):
    price, amount = map(int, sys.stdin.readline().split())
    if_total += price * amount
# 더함
if if_total - total_price == 0:
    print("Yes")
else: 
    print("No")
# 끝

그래서 이게 다다. 엥? if에 왜 if_total – total_price == 0이 들어갔어요? 두개 똑같으면 뺐을때 0이잖음.

import sys
total_price = int(sys.stdin.readline())
what_buy = int(sys.stdin.readline())
# 총합과 물건 종류 수
if_total = 0
i = 1
while i <= what_buy:
    price, amount = map(int, sys.stdin.readline().split())
    if_total += price * amount
    i += 1
# 더함
if if_total - total_price == 0:
    print("Yes")
else: 
    print("No")
# 끝

While은 이거.

import sys
total_price = int(sys.stdin.readline())
what_buy = int(sys.stdin.readline())
# 총합과 물건 종류 수
if_total = 0
while True:
    try:
        price, amount = map(int, sys.stdin.readline().split())
        if_total += price * amount
    except:
        break
# 더함
if if_total - total_price == 0:
    print("Yes")
else: 
    print("No")
# 끝

While True는 트라이 익셉 주면 되는데, 이건 틀린다. 왜냐하면 입력할 개수를 지정해줬기 때문에 While True가 의미가 없기 때문.

문제

https://www.acmicpc.net/problem/3003

체스 기물의 수를 입력하면 몇 개가 모자라거나 남는지 출력하기.

풀이

일단 체스는 흑백의 킹, 퀸, 룩, 비숍, 나이트, 폰으로 이루어져 있다. 킹 하나, 퀸 하나, 룩/비숍/나이트 둘에 폰 여덟이라 16개. 그나저나 이 문제 if로 가야 하는 거 아님? 왜 여기 있음? 아니 if 안가도 됨…

일단 각 케이스를 보자.

0 1 2 2 2 7 -> 1 0 0 0 0 1
2 1 2 1 2 1 -> -1 0 0 1 0 7

1번 케이스: 킹과 폰이 하나씩 모자람

2번 케이스: 킹이 하나 남고 비숍, 폰이 각각 1개, 7개 모자람

그러니까 입력값과 출력값을 합해서 1 1 2 2 2 8이 나와야 한다.

import sys
king, queen, rook, bishop, knight, pawn = map(int, sys.stdin.readline().split())
king = 1 - king
queen = 1 - queen
rook = 2 - rook
bishop = 2 - bishop
knight = 2 - knight
pawn = 8 - pawn
print(king, queen, rook, bishop, knight, pawn)

그러므로 입력받은 기물의 수를 1 1 2 2 2 8에서 빼면 된다.

문제

https://www.acmicpc.net/problem/7568

키와 몸무게를 이용해 덩치 등수를 산출한다.

Reference

https://bgspro.tistory.com/61

풀이

지금은 출석번호를 가나다순으로 매기지만(근데 대학 학번은 또 가나다순이 아님) 예에에에에에에에에전에는 키순으로 매겼었다. 나도 1학년때는 키순이었던걸로 기억함… 그래서 지금은 출석번호가 뒤쪽입니다 그러면 이름이 최씨 하씨 뭐 이런… 자음 순서에서도 뒤쪽이구나 하지만(보통 황씨, 황보씨가 맨 뒤) 예전같았으면 뭘 먹고 키가 그렇게 크셨어요? 했다.

당시에는 키가 비슷해보이는 애들끼리 비교했지만 이게 브루트 포스의 특성(개 노가다)과 맞물리면 이제 대단위 비교가 들어가게 된다. 

일단 입력 인자가 키랑 몸무게 두 개이다. 그러면 어캄? 아 어레이 불러야죠.

import sys
a = int(sys.stdin.readline().strip())
array = [[0 for i in range(a)] for j in range (a)]
for i in array:
    for j in i:
        print(j,end="")
    print()

Array는 이런 식으로 불러오는 게 맞다. 근데 키와 몸무게 두 개만 입력받는거고, 그 리스트를 a(입력값)만큼 층층이 쌓는… 그런 2차원 어레이를 만들거기때문에 컬럼값이 고정된다. 즉

import sys
a = int(sys.stdin.readline().strip())
array = [[0 for i in range(2)] for j in range (a)]
for i in array:
    for j in i:
        print(j,end="")
    print()

입력값이 n일 때, 우리는 2 * n 배열을 만들 것이다. (컬럼/로우)

import sys
a = int(sys.stdin.readline().strip())
array = [[list(map(int, sys.stdin.readline().split()))] for j in range (a)]
for i in array:
    for j in i:
        print(j,end="")
    print()

키랑 몸무게를 입력받을거니까 이렇게 해 주자. 다음으로는 덩치 등수를 매길건데… 일단 예시를 보자.

어? 왜 3등 아니고 5등임? 그건 문제에 나와있다. ‘덩치 등수는 나보다 덩치 큰 사람 k명 +1’이고, E는 5명 중 덩치가 제일 작으니까 자기보다 큰 사람이 4명 있어서 5등이다. 2등 세명은 씁 애매한데? 라고 보면 된다. 문제에서는 키와 몸무게 둘 다 클 때 덩치가 크다고 정의하고 있는데, 한 쪽은 크나 한 쪽이 작은 경우에는 비교하기 애매하기 때문. 체적을 비교하죠 어떻게요 아르키메데스좌 출동 

import sys
a = int(sys.stdin.readline().strip())
array = [list(map(int, sys.stdin.readline().split())) for j in range (a)]
ans = []
for i in range(a):
    count = 0
    for j in range(a):
        if array[i][0] < array[j][0] and array[i][1] < array[j][1]:
            count += 1
    ans.append(count+1)

for k in ans:
    print(k,end=" ")

그니까 나보다 큰 사람이 몇 명인가’만’ 세면 된다. 그 count를 배열에 넣으면 되고, 그 배열을 end=” “로 뽑으면 된다. 엥? 그거 공백 빼면 안돼요? end=””로 하면 숫자 다 붙어서 나온다.

11
1.2 55.5 #텅비드
2.4 333.6 #매시붕
1.8 25.0 #페로코체
3.8 100.0 #전수목
9.2 999.9 #철화구야
0.3 0.1 #종이신도
5.5 888.0 #악식킹
0.6 1.8 #베베놈(베베벱베)
3.6 150.0 #아고용
5.5 820.0 #차곡차곡(이름임)
1.8 13.0 #두파팡(대왕츕팝츕스)
7 4 7 4 1 11 2 10 4 2 7

float로 바꾸고 울비들 덩치 등수 매긴 결과. 순서는 도감 번호다. 울비들 덩치는 철화구야 > 악식킹/차곡차곡 > 매시붕/전수목/아고용 > 텅비드/페로코체/두파팡 > 베베놈 > 종이신도.

문제

일단 생성자가 뭔지를 알아야 이 문제를 풀 수 있다. 분해합은 어떤 수와 그 자릿수를 더하는 것으로, 예를 들자면

255 + 2 + 5 + 5 = 267

이런 식으로 더하는 것. 이 때, 255는 267의 생성자이다. 이런 식으로 어떤 수와 그 수의 자릿수를 다 더했을 때 입력한 수가 되면 그게 생성자이다. 그럼 ‘제일 작은 생성자’가 뭔지는 알 것 같은데, 생성자가 없는 수도 있나요? 1이랑 255 없었음.

import sys
a = int(sys.stdin.readline().strip())
sum = 0
for i in range(1,a+1):
    print(i)

sum은 나중가면 다른걸로 바꾼다. 배열 합이 그래서… 아무튼 일단 입력받은 수까지 for문 뺑뺑이를 돌려보자.

import sys
a = int(sys.stdin.readline().strip())
result = 0
for i in range(1,a+1):
    A = list(map(int, str(i)))
    result = sum(A) + i
    if result == a:
        print(i)
        break
    elif i == a:
        print(0)

이게 답인데… 숫자를 배열하해서 각각의 자리수를 리스트에 넣고 그 리스트의 합(sum)을 구해 원래 수와 더하게 된다. 그걸 한땀한땀 정성스럽게 계산하다가 생성자가 나오면 출력하고 땡이고, 계산 다 해봤는데 i가 입력한 수까지 가게 되면 씁 이건 에반데가 된다. 그러면 여기서 여러분들은 한가지 궁금증이 생기게 된다.

그런데 break가 굳이 있어야 할 이유가 있나요? 사족 아니예요? 

문제에 이런 얘기가 있었다.

물론, 어떤 자연수의 경우에는 생성자가 없을 수도 있다. 반대로, 생성자가 여러 개인 자연수도 있을 수 있다.

그리고 생성자가 ‘없는’ 경우를 elif에서 처리했고, ‘가장 작은 생성자’에 대한 처리를 break로 했다고 보면 된다.

1024 # 입력한 값
998
1016

1024는 생성자가 두 개거든. 이런 경우 응애 애기출력초과! 가 여러분을 반길 수 있다.

import sys
a = int(sys.stdin.readline().strip())
result = 0
i = 1

while i <= a:
    A = list(map(int, str(i)))
    result = sum(A) + i
    if result == a:
        print(i)
        break
    elif i == a:
        print(0)
    i += 1

이건 While버전. for문은 1부터 입력값+1이고, while은 i가 입력값과 같아질때까지 계속 한다. if문에 break가 들어가는 건 똑같다.

문제

https://www.acmicpc.net/problem/2798

블랙잭… 진짜로 그냥 블랙잭이다. 카드 장 수와 마지노선, 그리고 카드가 주어질 때 카드 세 장의 합이 마지노선을 넘지 않으면서 제일 큰 수를 구하는 문제다.

Reference

https://go-coding.tistory.com/67

https://duwjdtn11.tistory.com/297

풀이

일단 마지노선 얘기가 왜 나왔냐면… 블랙잭의 룰이 그렇다. 블랙잭은 카드의 합이 21이 넘으면 지는 게임이기 때문. 즉, 카드 세 장의 숫자를 합쳤을 때 마지노선으로 주어진 수보다 크면 진다. 예시 입력에서 5 21이 주어졌을 때 답은 21이지만 실제로 다 계산해보면 24까지 나온다. 근데 24는 21보다 크니까 져요… 

import sys
card,maginot = map(int,sys.stdin.readline().split())
card_list = list(map(int,sys.stdin.readline().split()))
result = 0

for i in range(card):
    for j in range(i+1,card):
        for k in range(j+1,card):
            if card_list[i]+card_list[j]+card_list[k] > maginot:
                continue
            else: 
                result = max(result,card_list[i]+card_list[j]+card_list[k])
#이거 3중for로 되는거 실화? 

print(result)

이걸로 한방컷 실화…

문제

https://www.acmicpc.net/problem/11729

하노이의 탑에서 이동 경로와 최종 이동 횟수를 출력하시오. (입력: 원판 개수)

Reference

https://ko.wikipedia.org/wiki/하노이의_탑

https://shoark7.github.io/programming/algorithm/tower-of-hanoi

https://8iggy.tistory.com/100

https://soohyun6879.tistory.com/m/190

풀이

정말 의외로 심플하게 풀렸다. (실화)

하노이의 탑? 

퍼즐 게임이다. 막대기 세 개와 크기가 다른 원판이 있고 이걸 맨 끝에 있는 막대로 옮기면 되는데 규칙이 있다.

1. 원판은 한번에 하나씩 옮길 수 있다.
2. 원판은 맨 위에 것만 옮길 수 있다.
3. 큰 원판이 작은 원판 위에 오면 안된다.

참 쉽죠?

def hanoi(number_of_disks_to_move, from_, to_, via_):
    if number_of_disks_to_move == 1:
        print(from_, "->", to_)
    else:
        hanoi(number_of_disks_to_move-1, from_, via_, to_)
        print(from_, "->", to_)
        hanoi(number_of_disks_to_move-1, via_, to_, from_)

하노이의 탑 코드를 보면 이렇게 입력 인자가 4개이다. (원판 개수, 어디서, 어디로, 어디를 통해) 개수 말고 세 개는 시작 기둥, 끝 기둥, 보조 기둥이다. 그럼 입력인자 다이어트가 되나요?

import sys
a = int(sys.stdin.readline())
def hanoi(n,start,to,via):
    if n == 1:
        print('{} {}'.format(start,to))
    else:
        hanoi(n-1,start,via,to)
        print('{} {}'.format(start,to))
        hanoi(n-1,via,to,start)
print(2 ** a - 1)
hanoi(a,1,3,2)

아뇨 그냥 입력인자 중 기둥을 고정값으로 두시면 됩니다. 이동 횟수도 재귀하면서 셀 필요 없이 2^n-1 하면 된다. (n=원판 개수)