from skbio import DNA

seq1 = DNA.read('data/seqs.fasta', qual='data/seqs.qual')
seq2 = DNA.read('data/seqs.fastq', variant='illumina1.8')
seq1

seq1 == seq2

"""
Explanation: A Bioinformatics Library for Data Scientists, Students, and Developers
Jai Rideout and Evan Bolyen
Caporaso Lab, Northern Arizona University
What is scikit-bio?
A Python bioinformatics library for:


data scientists


students


developers



"The first step in developing a new genetic analysis algorithm is to decide how to make the input data file format different from all pre-existing analysis data file formats." - Law's First Law

<span style='line-height:2em; word-spacing:2em'>Axt BAM SAM BED bedGraph bigBed bigGenePred table bigWig Chain GenePred table GFF GTF HAL MAF Microarray Net Personal Genome SNP format PSL VCF WIG  abi ace clustal embl fasta fastq genbank ig imgt nexus phred phylip pir seqxml sff stockholm swiss tab qual uniprot-xml emboss PhyolXML NexML newick CDAO MDL bcf caf gcproj scf SBML lsmat ordination qseq BIOM ASN.1 .2bit .nib ENCODE ... </span>
<span style='line-height:2em; word-spacing:2em'>Axt BAM SAM BED bedGraph bigBed bigGenePred table bigWig Chain GenePred table GFF GTF HAL MAF Microarray Net Personal Genome SNP format PSL VCF WIG  abi ace <span class='supio'>clustal</span> embl <span class='supio'>fasta</span> <span class='supio'>fastq</span> genbank ig imgt nexus phred <span class='supio'>phylip</span> pir seqxml sff stockholm swiss tab qual uniprot-xml emboss PhyolXML NexML <span class='supio'>newick</span> CDAO MDL bcf caf gcproj scf SBML <span class='supio'>lsmat</span> <span class='supio'>ordination</span> <span class='supio'>qseq</span> BIOM ASN.1 .2bit .nib ENCODE ... </span>
I/O in bioinformatics is hard


format redundancy (many-to-many)


format ambiguity


heterogeneous sources


How can we solve this?
An I/O Registry!
Format redundancy (many-to-many)
End of explanation
"""


import skbio.io

skbio.io.sniff('data/mystery_file.gz')

"""
Explanation: Format ambiguity
End of explanation
"""


from skbio import TreeNode

tree1 = skbio.io.read('http://localhost:8888/files/data/newick.gz', 
                      into=TreeNode)
print(tree1.ascii_art())

"""
Explanation: Heterogeneous sources
Read a gzip file from a URL:
End of explanation
"""


import io 

with io.open('data/newick.bz2', mode='rb') as open_filehandle:
    tree2 = skbio.io.read(open_filehandle, into=TreeNode)

print(tree2.ascii_art())

"""
Explanation: Read a bz2 file from a file path:
End of explanation
"""


tree3 = skbio.io.read(['((a, b, c), d:15):0;'], into=TreeNode)
print(tree3.ascii_art())

"""
Explanation: Read a list of lines:
End of explanation
"""


!cat data/yasf-seq.yml

import yaml

yasf = skbio.io.create_format('yasf')

@yasf.sniffer()
def yasf_sniffer(fh):
    return fh.readline().rstrip() == "#YASF", {}

@yasf.reader(DNA)
def yasf_to_dna(fh):
    seq = yaml.load(fh.read())
    return DNA(seq['Sequence'], metadata={
        'id': seq['ID'],
        'location': seq['Location'],
        'description': seq['Description']
    })

seq = DNA.read("data/yasf-seq.yml")
seq

"""
Explanation: Let's make a format!
YASF (Yet Another Sequence Format)
End of explanation
"""


seq.write("data/not-yasf.fna", format='fasta')
!cat data/not-yasf.fna

"""
Explanation: Convert YASF to FASTA
End of explanation
"""


from skbio.util._decorator import stable

@stable(as_of='0.4.0')
def add(a, b):
    """add two numbers.
    
    Parameters
    ----------
    a, b : int
        Numbers to add.
        
    Returns
    -------
    int
        Sum of `a` and `b`.
    
    """
    return a + b

help(add)

"""
Explanation: We are in beta - should you even use our software?
YES!
API Lifecycle
End of explanation
"""


seq = DNA("AacgtGTggA", lowercase='exon')
seq

"""
Explanation: What is stable:

skbio.io 
skbio.sequence

&nbsp;
&nbsp;
What is next:

skbio.alignment
skbio.tree
skbio.diversity
skbio.stats
&lt;your awesome subpackage!&gt;

Sequence API: putting the scikit in scikit-bio
End of explanation
"""


seq.values

"""
Explanation: Made with numpy
End of explanation
"""


seq.positional_metadata

"""
Explanation: And a pinch of pandas
End of explanation
"""


seq[seq.positional_metadata['exon']]

"""
Explanation: Slicing with positional metadata:
End of explanation
"""


aligned_seqs_fp = 'data/gg_13_8_otus/rep_set_aligned/82_otus.fasta'
taxonomy_fp = 'data/gg_13_8_otus/taxonomy/82_otu_taxonomy.txt'

from skbio import DNA

fwd_primer = DNA("GTGCCAGCMGCCGCGGTAA",
                 metadata={'label':'fwd-primer'})
rev_primer = DNA("GGACTACHVGGGTWTCTAAT",
                 metadata={'label':'rev-primer'}).reverse_complement()

def seq_to_regex(seq):
    result = []
    for base in str(seq):
        if base in DNA.degenerate_chars:
            result.append('[{0}]'.format(
                ''.join(DNA.degenerate_map[base])))
        else:
            result.append(base)

    return ''.join(result)

regex = '({0}.*{1})'.format(seq_to_regex(fwd_primer),
                            seq_to_regex(rev_primer))

import numpy as np
import skbio

starts = []
stops = []
for seq in skbio.io.read(aligned_seqs_fp, format='fasta', 
                         constructor=DNA):
    for match in seq.find_with_regex(regex, ignore=seq.gaps()):
        starts.append(match.start)
        stops.append(match.stop)
        
locus = slice(int(np.median(starts)), int(np.median(stops)))
locus

kmer_counts = []
seq_ids = []
for seq in skbio.io.read(aligned_seqs_fp, format='fasta',
                         constructor=DNA):
    seq_ids.append(seq.metadata['id'])
    sliced_seq = seq[locus].degap()
    kmer_counts.append(sliced_seq.kmer_frequencies(8))

from sklearn.feature_extraction import DictVectorizer
X = DictVectorizer().fit_transform(kmer_counts)

taxonomy_level = 3 # class
id_to_taxon = {}
with open(taxonomy_fp) as f:
    for line in f:
       id_, taxon = line.strip().split('\t')
       id_to_taxon[id_] = '; '.join(taxon.split('; ')[:taxonomy_level])

y = [id_to_taxon[seq_id] for seq_id in seq_ids]

from sklearn.feature_selection import SelectPercentile

X = SelectPercentile().fit_transform(X, y)

from sklearn.cross_validation import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                    random_state=0)

from sklearn.svm import SVC

y_pred = SVC(C=10, kernel='linear', degree=3,
             gamma=0.001).fit(X_train, y_train).predict(X_test)

from sklearn.metrics import confusion_matrix, f1_score

cm = confusion_matrix(y_test, y_pred)
cm_normalized = cm / cm.sum(axis=1)[:, np.newaxis]
plot_confusion_matrix(cm_normalized, title='Normalized confusion matrix')

print("F-score: %1.3f" % f1_score(y_test, y_pred, average='micro'))

"""
Explanation: Application: building a taxonomy classifier
End of explanation
"""


