<jupyter_start><jupyter_text>basic classifier for Dna2Vec represenationBasic bacterial genome classification using Dna2Vec representation<jupyter_code>%load_ext autoreload
%autoreload 2
from faigen.genomic import sequence 
import sys
sys.path.append("..")
from faigen.data import sequence 
# from genomic import sequence
from faigen.data.sequence import regex_filter, count_filter, Dna2VecDataBunch
from functools import partial
import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn import manifold,neighbors
from scipy.cluster.hierarchy import dendrogram, linkage  
from matplotlib import pyplot as plt
import seaborn as sns; sns.set(color_codes=True)
import plotly.plotly as py
import plotly.graph_objs as go
from fastai import *
from fastai.data_block import *
from fastai.basic_train import *
from fastai.layers import *
from fastai.metrics import *
import torch.nn as nn
import torch.nn.functional as F<jupyter_output><empty_output><jupyter_text>Load Data<jupyter_code>from gensim.models import Word2Vec
# word_vectors = Word2Vec.load_word2vec_format('/data/genomes/embeddings/dna2vec-20190612-1611-k10to10-100d-10c-4870Mbp-sliding-kPR.w2v') 

word_vectors = Word2Vec.load_word2vec_format('/data/genomes/embeddings/dna2vec-20190611-1940-k8to8-100d-10c-4870Mbp-sliding-LmP.w2v') 


filters=[partial(regex_filter, rx="Streptomyces|Bacillus|Vibrio|Streptococcus|Rhizobium|Staphylococcus"),partial(regex_filter, rx="plasmid", keep=False)]
items = sequence.Dna2VecList.from_folder("/data/genomes/GenSeq_fastas/train",filters=filters,agg=partial(np.mean, axis=0),n_cpus=7,emb=word_vectors)
processors = [
    sequence.GSFileProcessor(),
    sequence.GSTokenizeProcessor(tokenizer=sequence.GSTokenizer(ngram=10, skip=0, n_cpus=7)),
    sequence.Dna2VecProcessor(emb=word_vectors)]
for p in processors:
    print("\n", type(p))
    %time p.process(items)
bunch.sanity_check()
DB="/data/genomes/GenSeq_fastas"
# DB='/home/serge/development/genomes/ncbi-genomes-2019-04-07/bacterial genomes'

filters=[partial(regex_filter, rx="Bacillus",keep=False),partial(regex_filter, rx="plasmid?\s", keep=False)]    
#        partial(count_filter,num_fastas=(1,1), keep=1)]

bunch = Dna2VecDataBunch.from_folder(DB,
             filters=filters,
             labeler=lambda x: x.split()[1],
             emb=word_vectors,ngram=8,skip=0,
             n_cpus=7,agg=partial(np.mean, axis=0), one_hot=True,bs=5)

bunch.one_batch()

class GenomeClassifier(nn.Module):
    def __init__(self):
        super(GenomeClassifier, self).__init__()
        self.l1 = nn.Linear(100, 100)
        self.linear_classifier = nn.Linear(100, bunch.train_dl.y.c)

    def forward(self, x):
        x = F.relu(self.l1(x))
        x = self.linear_classifier(x)
        return x
class BactClassifier(SequentialEx):
    pass
layers=[nn.Linear(100,100),nn.ReLU(),nn.Linear(100,17)]
bc = SequentialEx(*layers)
bc
learn = Learner(bunch, bc)
learn.lr_find()
learn.recorder.plot()
learn.fit(1)
%debug
learn.fit_one_cycle(25)
processors = [
    sequence.GSFileProcessor(),
    sequence.GSTokenizeProcessor(tokenizer=sequence.GSTokenizer(ngram=8, skip=0, n_cpus=4)),
    sequence.Dna2VecProcessor()]
%time for p in processors: p.process(data)
len(data.items)<jupyter_output><empty_output><jupyter_text>Genome vectors<jupyter_code>data.items
def log_scale(X):
    x=np.asarray(X);e=1e-6
    return np.log10(x+np.abs(x.min())+e) 


x=np.asarray(data.items)
bad_fastas = np.where(np.mean(x,axis=1) == 0.)[0]
X = np.delete(x, bad_fastas,0)
labelList=[" ".join(i.split()[1:3]) for i in data.descriptions]
labelList=np.delete(np.asarray(labelList), bad_fastas)
vocab=list(np.unique(labelList))
y=[vocab.index(x) for x in labelList]<jupyter_output><empty_output><jupyter_text>Correlation Distance in log-scaled space  tSNE<jupyter_code>tsne = manifold.TSNE(n_components=3, init='pca', perplexity=10, metric="correlation",random_state=0)
%time X3 = tsne.fit_transform(log_scale(X))
genus = [i.split()[0] for i in labelList]
genus_vocab=list(np.unique(genus))
y=[genus_vocab.index(x) for x in genus]
genus_vocab

X3_df = pd.DataFrame(data=X3, columns=["pc1",'pc2','pc3'], index=labelList)

X3_df["genus"]=genus
X3_df["y"]=y

genus_df=X3_df.groupby("genus").agg({"pc1": list, "pc2":list,"pc3":list,"y":np.mean})<jupyter_output><empty_output><jupyter_text>Correlation Distance visualisation<jupyter_code>data=[]
for g in genus_df.index:
    trace  = go.Scatter3d(
        name = str(g),
        x=genus_df.loc[g,"pc1"],
        y=genus_df.loc[g,"pc2"],
        z=genus_df.loc[g,"pc3"],
        mode='markers',
        marker=dict(
            size=8,
            color=genus_df.loc[g,"y"],                # set color to an array/list of desired values
            colorscale='Jet',           # choose a colorscale
            opacity=0.5)
    )

    data.append(trace)
    

layout = go.Layout(
    width=1000,
    height=1000,
    margin=dict(
        l=0,
        r=0,
        b=0,
        t=0
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='correlation distance ncbi-genomes-2019-04-07 Escherichia,Klebsiella,Bacillus')<jupyter_output>/home/serge/anaconda3/envs/bio/lib/python3.6/site-packages/IPython/core/display.py:689: UserWarning:

Consider using IPython.display.IFrame instead<jupyter_text>Eucleadian Distance in log-scaled space  tSNE<jupyter_code>tsne = manifold.TSNE(n_components=3, init='pca', perplexity=30,random_state=0)
%time X3 = tsne.fit_transform(log_scale(X))
genus = [i.split()[0] for i in labelList]
genus_vocab=list(np.unique(genus))
y=[genus_vocab.index(x) for x in genus]
genus_vocab

X3_df = pd.DataFrame(data=X3, columns=["pc1",'pc2','pc3'], index=labelList)

X3_df["genus"]=genus
X3_df["y"]=y

genus_df=X3_df.groupby("genus").agg({"pc1": list, "pc2":list,"pc3":list,"y":np.mean})<jupyter_output><empty_output><jupyter_text>Eucleadian Distance Visualisation<jupyter_code>data=[]
for g in genus_df.index:
    trace  = go.Scatter3d(
        name = str(g),
        x=genus_df.loc[g,"pc1"],
        y=genus_df.loc[g,"pc2"],
        z=genus_df.loc[g,"pc3"],
        mode='markers',
        marker=dict(
            size=8,
            color=genus_df.loc[g,"y"]+1,                # set color to an array/list of desired values
            colorscale='YlGnBu',           # choose a colorscale
            opacity=0.5)
    )

    data.append(trace)
    

layout = go.Layout(
    width=1000,
    height=1000,
    margin=dict(
        l=0,
        r=0,
        b=0,
        t=0
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='eucledian distance metric by genus Escherichia|Klebsiella|Bacillus')<jupyter_output>/home/serge/anaconda3/envs/bio/lib/python3.6/site-packages/IPython/core/display.py:689: UserWarning:

Consider using IPython.display.IFrame instead<jupyter_text>Eucleadian Distance in unmodified space  tSNE<jupyter_code>tsne = manifold.TSNE(n_components=3, init='pca', perplexity=30,random_state=0)
%time X3 = tsne.fit_transform(X)
genus = [i.split()[0] for i in labelList]
genus_vocab=list(np.unique(genus))
y=[genus_vocab.index(x) for x in genus]
genus_vocab

X3_df = pd.DataFrame(data=X3, columns=["pc1",'pc2','pc3'], index=labelList)

X3_df["genus"]=genus
X3_df["y"]=y

genus_df=X3_df.groupby("genus").agg({"pc1": list, "pc2":list,"pc3":list,"y":np.mean})<jupyter_output><empty_output><jupyter_text>Eucleadian Distance Visualisation<jupyter_code>data=[]
for g in genus_df.index:
    trace  = go.Scatter3d(
        name = str(g),
        x=genus_df.loc[g,"pc1"],
        y=genus_df.loc[g,"pc2"],
        z=genus_df.loc[g,"pc3"],
        mode='markers',
        marker=dict(
            size=8,
            color=genus_df.loc[g,"y"],                # set color to an array/list of desired values
            colorscale='Jet',           # choose a colorscale
            opacity=0.5)
    )

    data.append(trace)
    

layout = go.Layout(
    width=1000,
    height=1000,
    margin=dict(
        l=0,
        r=0,
        b=0,
        t=0
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='eucledian distance in native space Escherichia|Klebsiella|Bacillus')<jupyter_output>/home/serge/anaconda3/envs/bio/lib/python3.6/site-packages/IPython/core/display.py:689: UserWarning:

Consider using IPython.display.IFrame instead<jupyter_text>Genome Inventory<jupyter_code>all_fastas = sequence.Dna2VecList.from_folder(DB)
next(iter(all_fastas))
a = [x for x in all_fastas if "plasmid" not in x['description'].split()] 
a
a = pd.DataFrame(data=list(a.items), columns=["file", "description","id","name"]).groupby("file").agg({"id":'count'}).sort_values('id', ascending=False)
a.describe()
inventory = pd.DataFrame(data=[l.split()[1:3] for l in all_fastas], columns=["genus","species" ])
inventory.groupby("genus").agg({"species":"count"}).sort_values("species",ascending=False)
all_fasts
inventory.groupby(["genus", "species"]).agg({"species": "count"})
inventory.columns=["count"]
inventory
all_fastas
counts = inventory.reset_index().groupby("genus").agg({"count", sum}).drop(("species"), axis=1)
counts.columns=["n_sequences","species"]
counts.sort_values("n_sequences", ascending=False)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup