<jupyter_start><jupyter_text>Labs - Biopython and data formats Outline- Managing dependencies in Python with environments- Biopython     - Sequences (parsing, representation, manipulation)    - Structures (parsing, representation, manipulation) 1. Python environments- handles issues with dependencies versions- ensures reproducibility- does not clutter users' global site-packages directory`python3 -m venv venv/        Creates an environment called venv/``source venv/bin/activate``pip install biopython``pip freeze > requirements.txt``(venv) % deactivate`On a different machine, the environment can be replicated by creating a new environment and running`pip install -r requirements.txt` 2. BiopythonBiopython is a library consisting of tools for both sequence and structure bioinformatics. Among other things it enables parsing, handling and storing molecular data present in common formats such as FASTA, PDB or mmCIF.Install biopython using `pip install biopython`Functionality divided into packages list of which is available in the [docs](https://biopython.org/docs/1.75/api/Bio.html). Main sequence and structure packages: - [Bio.Seq](https://biopython.org/docs/latest/api/Bio.Seq.html) - [Bio.Align](https://biopython.org/docs/latest/api/Bio.Align.html)  - [Bio.SeqIO](https://biopython.org/docs/latest/api/Bio.SeqIO.html) - [Bio.PDB](https://biopython.org/docs/latest/api/Bio.PDB.html)     Sequences   Loading a sequence from a string:<jupyter_code>from Bio.Seq import Seq
seq = Seq("AGTACACTG")
print(seq)<jupyter_output>AGTACACTG<jupyter_text>This creates a [sequence object](https://biopython.org/docs/latest/api/Bio.Seq.html) with a couple of fancy methods, especially when it comes to nuclotide sequences such as `reverse_complement` or `translate`.<jupyter_code>print(seq.translate())
print(seq.reverse_complement().transcribe())
print(seq.reverse_complement().transcribe().translate())
coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
print(coding_dna.translate())
print(coding_dna.translate(to_stop=True))
print(coding_dna.translate(table=2))
print(coding_dna.translate(table=2, to_stop=True))<jupyter_output>MAIVMGR*KGAR*
MAIVMGR
MAIVMGRWKGAR*
MAIVMGRWKGAR<jupyter_text>Notice, in the example above we used different genetic tables. Check [NCBI genetic codes](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi) for details.To list all the methods, run, e.g., one of the following:<jupyter_code>print(dir(seq))
print(help(seq))<jupyter_output>['__abstractmethods__', '__add__', '__bytes__', '__class__', '__contains__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__imul__', '__init__', '__init_subclass__', '__le__', '__len__', '__lt__', '__module__', '__mul__', '__ne__', '__new__', '__radd__', '__reduce__', '__reduce_ex__', '__repr__', '__rmul__', '__setattr__', '__sizeof__', '__slots__', '__str__', '__subclasshook__', '__weakref__', '_abc_impl', '_data', 'back_transcribe', 'complement', 'complement_rna', 'count', 'count_overlap', 'encode', 'endswith', 'find', 'index', 'join', 'lower', 'lstrip', 'replace', 'reverse_complement', 'reverse_complement_rna', 'rfind', 'rindex', 'rsplit', 'rstrip', 'split', 'startswith', 'strip', 'tomutable', 'transcribe', 'translate', 'ungap', 'upper']
b'AGTACACTG'<jupyter_text>Methods for accessing by position  are available as well.<jupyter_code>print(seq[3])
print(seq[3:5])
print(seq[::-1])<jupyter_output>A
AC
GTCACATGA<jupyter_text>If needed, the Seq object can be converted into a string.<jupyter_code>print(str(seq))
print(str(seq).translate({65: 88}))<jupyter_output>AGTACACTG
XGTXCXCTG<jupyter_text>To parse sequence from a file, you can use [Bio.SeqIO](https://biopython.org/docs/latest/api/Bio.SeqIO.html). [Here](https://biopython.org/wiki/SeqIOfile-formats) is the list of supported formats. The format name is passed into the `parse` method.<jupyter_code>from Bio import SeqIO

sars2_it = SeqIO.parse("R1A-B_SARS2.fasta", "fasta")
for seq_record in sars2_it:    
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))
sars2_seq_recs = list(sars2_it)<jupyter_output>sp|P0DTD1|R1AB_SARS2
Seq('MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLV...VNN')
7096
sp|P0DTC1|R1A_SARS2
Seq('MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLV...FAV')
4405<jupyter_text>The result is an iterator of [SeqRecord](https://biopython.org/docs/latest/api/Bio.SeqRecord.html)s. Other attributes of `SeqRecord` such as features or annotations are more relevant for other formats, such as genbank. The underlying gene for the two isoforms (R1A_SARS2/P0DTC1 and R1AB_SARS2/P0DTD1) is ORF1ab and the two isoforms are caused by ribosomal slippage during translation (see, e.g., [here](https://www.science.org/doi/full/10.1126/science.abf3546)). Both reading frames R1A_SARS2 and R1AB_SARS2 are polyproteins and are encoded by the same [gene](https://www.ncbi.nlm.nih.gov/gene/43740578). Let's explore this.<jupyter_code>gb_rec = list(SeqIO.parse("NC_045512.gb", "genbank"))[0]
print(gb_rec.)
print(gb_rec.annotations)
print(gb_rec.features)
gb_rec.features<jupyter_output><empty_output><jupyter_text>Let's obtain all CDS (coding sequence) features.<jupyter_code>print(gb_rec.features)
cds = [seq_feature for seq_feature in gb_rec.features if seq_feature.type == 'CDS']
cds[0].extract(gb_rec.seq).translate()<jupyter_output><empty_output><jupyter_text>Now, let's get the DNA sequence for the the polyprotein 1ab.<jupyter_code>aa_seq = cds[0].extract(gb_rec.seq).translate()
print(aa_seq[-10:])
print(gb_rec.seq.translate()[-10:])<jupyter_output>ISSDVLVNN*
KKKKKKKKKK<jupyter_text>To write a sequence into a file, use `SeqIO.write`.<jupyter_code>SeqIO.write([gb_rec, SeqIO.SeqRecord(aa_seq, id="id", description="aa")], "fasta_from_gb.fasta", "fasta")<jupyter_output><empty_output><jupyter_text>Now, carry out the following tasks by yourselfs:- Obtain the protein sequnece for polyprotein 1ab and check with UniProt that it matches (just by eyeballing).- Obtain the protein sequence for the polyprotien 1a.- Obtain protein sequences for all the proteins and list them together with their names<jupyter_code>print(cds[1].extract(gb_rec.seq).translate())
#print(cds[1].extract(gb_rec.seq).translate())

peptides = [print("{}: {}".format(ft.qualifiers['protein_id'], ft.extract(gb_rec.seq).translate())) for ft in gb_rec.features if ft.type == 'mat_peptide']<jupyter_output>['YP_009725297.1']: MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLVEVEKGVLPQLEQPYVFIKRSDARTAPHGHVMVELVAELEGIQYGRSGETLGVLVPHVGEIPVAYRKVLLRKNGNKGAGGHSYGADLKSFDLGDELGTDPYEDFQENWNTKHSSGVTRELMRELNGG
['YP_009725298.1']: AYTRYVDNNFCGPDGYPLECIKDLLARAGKASCTLSEQLDFIDTKRGVYCCREHEHEIAWYTERSEKSYELQTPFEIKLAKKFDTFNGECPNFVFPLNSIIKTIQPRVEKKKLDGFMGRIRSVYPVASPNECNQMCLSTLMKCDHCGETSWQTGDFVKATCEFCGTENLTKEGATTCGYLPQNAVVKIYCPACHNSEVGPEHSLAEYHNESGLKTILRKGGRTIAFGGCVFSYVGCHNKCAYWVPRASANIGCNHTGVVGEGSEGLNDNLLEILQKEKVNINIVGDFKLNEEIAIILASFSASTSAFVETVKGLDYKAFKQIVESCGNFKVTKGKAKKGAWNIGEQKSILSPLYAFASEAARVVRSIFSRTLETAQNSVRVLQKAAITILDGISQYSLRLIDAMMFTSDLATNNLVVMAYITGGVVQLTSQWLTNIFGTVYEKLKPVLDWLEEKFKEGVEFLRDGWEIVKFISTCACEIVGGQIVTCAKEIKESVQTFFKLVNKFLALCADSIIIGGAKLKALNLGETFVTHSKGLYRKCVKSREETGLLMPLKAPKEIIFLEGETLPTEVLTEEVVLKTGDLQPLEQPTSEAVEAPLVGTPVCINGLMLLEIKDTEKYCALAPNMMVTNNTFTLKGG
['YP_009725299.1']: APTKVTFGDDTVIEVQGYKSVNITFELDERIDKVLNEKCSAYTVELGTEVNEFACVVADAVIKTLQPVSELLTPLGIDLDEWSMATYYLFDESGEFKLASHMYCSFYPPDEDEEEGDCEE[...]<jupyter_text>StructuresStructure processing is managed by the [Bio.PDB](https://biopython.org/docs/latest/api/Bio.PDB.html) package.To read a structure from a PDB file, use the `PDBParser`. We will be using the 3C-like proteinase protein, which is one of the processed proteins present in the ORF1a discussed above. One of it's structures is [7ALH](https://www.ebi.ac.uk/pdbe/entry/pdb/7alh). To see all the proteins, I suggest checking out the PDBe-KB page for [P0DTD1](https://www.ebi.ac.uk/pdbe/pdbe-kb/protein/P0DTD1).<jupyter_code>from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure = parser.get_structure("7alh", "7alh.ent")<jupyter_output><empty_output><jupyter_text>As the PDB format is considered deprecated, one should use the mmCIF file instead. This is done the same way as in case of PDB files.<jupyter_code>from Bio.PDB.MMCIFParser import MMCIFParser
parser = MMCIFParser()
structure = parser.get_structure("7alh", "7alh.cif")<jupyter_output><empty_output><jupyter_text>To retrieve the individual CIF dictionary fields, one can use the `MMCIF2Dict` module.<jupyter_code>from Bio.PDB.MMCIFParser import MMCIF2Dict
mmcif_dict = MMCIF2Dict("7alh.cif")
print(mmcif_dict["_citation.title"])<jupyter_output>['Crystal structure of the main protease (3CLpro/Mpro) of SARS-CoV-2 at 1.65A resolution (spacegroup C2).']<jupyter_text>The structure record has the structure->model->chain->residue architecture.![SMRCA](http://biopython.org/DIST/docs/tutorial/images/smcra.png) Each of the levels in the hierarchy is represented by a submodule in Bio.PDB, namely [Bio.Structure](https://biopython.org/docs/latest/api/Bio.PDB.Modul.html), [Bio.Module](https://biopython.org/docs/latest/api/Bio.PDB.Module.html),[Bio.Chain](https://biopython.org/docs/latest/api/Bio.PDB.Chain.html),[Bio.Residue](https://biopython.org/docs/latest/api/Bio.PDB.Residue.html) and [Bio.Atom](https://biopython.org/docs/latest/api/Bio.PDB.Atom.html). For details regarding IDs, check the [section on ID](https://biopython.org/docs/1.75/api/Bio.PDB.Entity.htmlBio.PDB.Entity.Entity.get_full_id) of the Entity class which is the superclass of the Module/Chain/Residue/Atom classes.<jupyter_code>print(structure.get_list())

print('---------- MODEL INFO ----------')

model = structure[0]
print("Full ID: {}\nID: {}".format(model.get_full_id(), model.get_id()))
print(model.get_list())

print('---------- CHAIN INFO ----------')
chain = model['A']
print("Full ID: {}\nID: {}".format(chain.get_full_id(), chain.get_id()))
print(chain.get_list())

print('---------- RESIDUE INFO ----------')
res = chain[(' ',1,' ')]
print("Full ID: {}\nID: {}".format(res.get_full_id(), res.get_id()))
print(res.get_resname())
res = chain[1]
print(res.get_resname())

print(res.get_list())
print('---------- ATOM INFO ----------')
atom=res['CA']
print("Full ID: {}\nID: {}".format(atom.get_full_id(), atom.get_id()))
print("{}\n{}\n{}\n{}".format(atom.get_name(), atom.get_id(), atom.get_coord(), atom.get_fullname()))
print(atom.get_vector())<jupyter_output>[<Model id=0>]
---------- MODEL INFO ----------
Full ID: ('7alh', 0)
ID: 0
[<Chain id=A>]
---------- CHAIN INFO ----------
Full ID: ('7alh', 0, 'A')
ID: A
[<Residue SER het=  resseq=1 icode= >, <Residue GLY het=  resseq=2 icode= >, <Residue PHE het=  resseq=3 icode= >, <Residue ARG het=  resseq=4 icode= >, <Residue LYS het=  resseq=5 icode= >, <Residue MET het=  resseq=6 icode= >, <Residue ALA het=  resseq=7 icode= >, <Residue PHE het=  resseq=8 icode= >, <Residue PRO het=  resseq=9 icode= >, <Residue SER het=  resseq=10 icode= >, <Residue GLY het=  resseq=11 icode= >, <Residue LYS het=  resseq=12 icode= >, <Residue VAL het=  resseq=13 icode= >, <Residue GLU het=  resseq=14 icode= >, <Residue GLY het=  resseq=15 icode= >, <Residue CYS het=  resseq=16 icode= >, <Residue MET het=  resseq=17 icode= >, <Residue VAL het=  resseq=18 icode= >, <Residue GLN het=  resseq=19 icode= >, <Residue VAL het=  resseq=20 icode= >, <Residue THR het=  resseq=21 icode= >, <Residue CYS het=  resseq=22 icode[...]<jupyter_text>To download a file from PDB, one can use the PDBList module.<jupyter_code>from Bio.PDB.PDBList import PDBList
pdbl = PDBList()
pbl_7lkr=pdbl.retrieve_pdb_file("7LKR", file_format="mmCif", pdir=".")
from Bio.PDB.MMCIFParser import MMCIFParser
parser = MMCIFParser()
structure = parser.get_structure("7lkr", "7lkr.cif")<jupyter_output>c:\python39\lib\site-packages\Bio\PDB\StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line 4467.
  warnings.warn(
c:\python39\lib\site-packages\Bio\PDB\StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain B is discontinuous at line 4531.
  warnings.warn(
c:\python39\lib\site-packages\Bio\PDB\StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line 4608.
  warnings.warn(
c:\python39\lib\site-packages\Bio\PDB\StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain B is discontinuous at line 4774.
  warnings.warn(
{'fraction_non_alphanumeric': 0.14104824491104342, 'fraction_numerical': 0.017871453758615163, 'mean_word_length': 6.508423586040915, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 81, 'https://': 21, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '9520', 'license': 'MIT', 'n_tokens_mistral': 5000, 'n_tokens_neox': 4691, 'n_words': 948, 'path': '01-data-formats/lab.ipynb', 'repo_name': 'davidhoksza/bioinformatics-algo-labs', 'size': '12478'}