Biopython SeqIO 써보기

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

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

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로 저장하는 것도 된다. (단백질은 하지말자...)

저작자표시

'Coding > Python' 카테고리의 다른 글

번외편-코딩테스트 풀이 (0)	2022.08.20
Biopython으로 MSA 해보기 (0)	2022.08.20
Biopython으로 시퀀스 레코드 생성하고 만져보기 (0)	2022.08.20
Biopython으로 시퀀스 다뤄보기 (0)	2022.08.20
Biopython으로 시퀀스 가져오기 (0)	2022.08.20

Biopython SeqIO 써보기

Iteration 활용하기

레코드 반복

레코드 목록 뽑기

데이터 추출

수정

웹에서 받아오기

파일 주세요 Genbank여...

파일 주세요 swissprot이여...

딕셔너리화하기

SeqIO.to_dict

SeqIO.index

쓰기

'Coding > Python' 카테고리의 다른 글

최근댓글

최근글

인기글

티스토리툴바

Biopython SeqIO 써보기

Iteration 활용하기

레코드 반복

레코드 목록 뽑기

데이터 추출

수정

웹에서 받아오기

파일 주세요 Genbank여...

파일 주세요 swissprot이여...

딕셔너리화하기

SeqIO.to_dict

SeqIO.index

쓰기

'Coding > Python' 카테고리의 다른 글

관련글

최근댓글

최근글

인기글

티스토리툴바