# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# <h1>biopython Introduction</h1>
# <p>
# For the whole documentation, see http://biopython.org/wiki/Biopython, for the long tutorial see http://biopython.org/DIST/docs/tutorial/Tutorial.html, for the cookbook see http://biopython.org/wiki/Category%3ACookbook.
# </p>
# <p>
# biopython is a universal bioinformatic Python package, coordinated by the Open Bioinformatics Foundation (https://www.open-bio.org/wiki/Main_Page, also coordinates - amongst other projects - BioPerl http://bioperl.org/ and BioJava http://biojava.org/). If you use the Anaconda python distribution, the easiest way to install biopython is by using the Anaconda Navigator. If you use pip, you can install biopython with the command "pip install biopython".
# </p>
# <p>
# Once you have installed biopython, you can import it to a Python script like this (not that you may have to import submodules of biopython to access all of its function):
# </p>

import Bio

# <h3>1. Storing DNA, RNA and Amino Acid Sequences</h3>

from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

# An unambiguous DNA sequence
unambiguous = Seq("ACGT", IUPAC.unambiguous_dna)
print(unambiguous)
print(unambiguous.alphabet)
print(unambiguous.count("AC")) # Returns only non-overlapping occurences
print(len(unambiguous))

# An ambiguous DNA sequence
ambiguous = Seq("ACGT", IUPAC.ambiguous_dna)
ambiguous

# An RNA sequence
# An ambiguous DNA sequence
ambiguous = Seq("ACGU", IUPAC.ambiguous_rna)
ambiguous

# An amino acid sequence
protein = Seq("TRUMP", IUPAC.protein)
protein

# <h3>2. Translation and transcription</h3>

# Complementary DNA
sequence = Seq("AGCCCTYA", IUPAC.ambiguous_dna)
sequence.complement()

# Reverse-complementary DNA
sequence.reverse_complement()

# DNA->RNA
dna = Seq("AAATTTCCCGGG", IUPAC.unambiguous_dna)
rna = dna.transcribe()
rna

# RNA->Proteins
print(dna.translate()) # The DNA's or RNA's length be a multiple of 3
print(rna.translate())

# <h3>3. Access to common databases</h3>

# NCBI Entrez...
from Bio import Entrez

# ...get all databases...
handle = Entrez.einfo()
result = handle.read()
handle.close()
result

# ...look up in pubmed
Entrez.email = ""
handle = Entrez.esearch(db="pubmed", term="Vibrio")
record = Entrez.read(handle)
record

# KEGG...
from Bio.KEGG import REST
from Bio.KEGG import Enzyme

# ...get pathways for Vibrio alginolyticus...
pathways = REST.kegg_list("pathway", "vag").read()

# ...select a pathway...
pathway_line = pathways.rstrip().split("\n")[0]
entry, description = pathway_line.split("\t")
entry

# ...and analyze this pathway
pathway_file = REST.kegg_get(entry).read()

# <h3>4. Sequence Alignments</h3>

# +
# Simple Needleman-Wunsch (note: two cells below, a more complex alternative is shown)
from Bio import pairwise2
from Bio import SeqIO

seq1 = Seq("AAATTGGGCCGG", IUPAC.protein)
seq2 = Seq("AAAGGTTAAA", IUPAC.protein)

alignments = pairwise2.align.globalxx(seq1, seq2)
print(alignments[0]) # First optimal alignment
# -

# Simple Smith-Watherman (note: two cells below, a more complex alternative is shown)
alignments2 = pairwise2.align.localxx(seq1, seq2)
print(alignments2[0]) # First optimal alignment

# +
# Better global alignment alternative with penalizing options
from Bio.SubsMat.MatrixInfo import blosum62

alignments3 = pairwise2.align.globalds(seq1, seq2, blosum62, -10, -0.5)
print(alignments3[0])
# -

# Better local alignment alternative with penalizing options
alignments4 = pairwise2.align.localds(seq1, seq2, blosum62, -10, -1)
print(alignments4[0])

# <h3>5. BLAST</h3>

# The biopython option is not reliable. Use the command line instead, e.g.
blastn_output = run("blastn", "-query", "temp_in.txt", "-out", "temp_out.txt", "-db", "nt", "-outfmt", '6 sacc pident qstart qend evalue bitscore stitle sscinames sskingdoms', "-max_target_seqs", "5")

# <h3>6. Sequence Motif Analysis</h3>

from Bio import motifs

# Generate nucleotide count report
seqset = [Seq("ACGTG"),
          Seq("ATGTT"),
          Seq("TTTTT")]
motif = motifs.create(seqset)
motif.counts

# Access a line
motif.counts["T"]

# Get consensus
motif.consensus

# Get anticonsensus
motif.anticonsensus

# Get degenerate consensus
motif.degenerate_consensus

# Create a logo (using a web service)
motif.weblogo("./biopython/motif.png")

# <img src="./biopython/motif.png"></img>

# PSB 2017

StarcoderdataJupyterScriptsDedupFiltered