{ Les snippets de bioinfo-fr.net }

#!/usr/bin/env python
# coding: utf-8

# # Installing PyEnsembl
# `PyEnsembl`is a Python interface to Ensembl reference genome metadata such as exons and transcripts. `PyEnsembl` downloads GTF and FASTA files from the Ensembl FTP server and loads them into a local database. PyEnsembl can also work with custom reference data specified using user-supplied GTF and FASTA files.

# In[52]:

import pyensembl


# In[53]:

from pyensembl import EnsemblRelease


# In[54]:

data = EnsemblRelease(78)


# In[55]:

import os, sys


# In[56]:

os.path.isfile("TruSeq.txt")


# In[57]:

f = open("TruSeq.txt",'r',encoding='utf8')


# In[58]:

L = f.readlines()


# In[59]:

f.close()


# In[60]:

genes = [str(L[i]).rstrip("\n") for i in range(1,len(L))]


# In[61]:

liste = [str(genes[i]).split('\t') for i in range(len(genes))]


# In[62]:

liste_simple = []
for elt in range(0,len(liste)):
    liste_simple = liste_simple + liste[elt]


# In[63]:

print(liste_simple)


# ## get all Genes ID from TrueSeq Amplicon cancer panel

# In[64]:

f = open("TruSeq_geneIDs.txt",'w',encoding = 'utf8')
tTruSeq = []

for i in range(0,len(liste_simple)):
    gene_courant = str(liste_simple[i]).rstrip()
    gene_ID = pyensembl.ensembl_grch38.genes_by_name(gene_courant)[0].gene_id
    f.write("gene: {} , gene_ID : {} \n".format(str(liste_simple[i]),gene_ID))
    tTruSeq = [gene_ID, pyensembl.ensembl_grch38.genes_by_name(gene_courant)[0].transcripts] + tTruSeq
f.close()
    

# In[65]:

f = open("genotype_transcripts.txt",'w',encoding = 'utf8')
f.write(str(tTruSeq))
f.close()


# ## compare genotype Transcript IDs vs NGS phenotypic transcripts (VariantsNGS.tsv)

# In[66]:

import os, sys


# In[67]:

import pyensembl


# In[68]:

os.path.isfile("VariantsNGS.tsv")


# In[69]:

f = open("VariantsNGS.tsv",'r',encoding='utf8')


# In[70]:

L2 = f.readlines() 


# In[71]:

f.close()


# In[72]:

genes = [str(L2[i]).rstrip("\n") for i in range(1,len(L2))]


# In[73]:

liste2 = [str(genes[i]).split('\t') for i in range(len(genes))]


# In[74]:

liste_simple2 = []
for elt in range(0,len(liste2)):
    liste_simple2 = liste_simple2 + liste2[elt]

# In[75]:

print(liste_simple2)

# ## manual cross-check of each transcript IDs

# In[76]:

VHL = pyensembl.ensembl_grch38.genes_by_name("VHL")


# In[77]:

tVHL = VHL[0].transcripts


# In[78]:

print(tVHL[0].transcript_id)


# same Transcript IDs found

# In[79]:

SETD2 = pyensembl.ensembl_grch38.genes_by_name("SETD2")
tSETD2 = SETD2[0].transcripts
print(tSETD2[0].transcript_id)


# different transcript IDs found : ENST00000431180 != ENST00000409792

# l'automatisation de la comparaison des deux listes simples n'est pas encore terminee

import Bio
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

#on importe la sequence mystere avec les outils en ligne

def seq_codante_de_phase(phase):
    meilleure_seq=compteur=i=j=start=stop=0
    longueur = len(phase)-1
    for compteur in range (i,longueur):
        if (phase[i] != 'M'):
            i += 1
        else:
            start = i
            j = start
            while (phase[j] != '*' and j < len(phase)-1):
                j += 1
            if (phase[j] == "*"):
                stop = j
            seq_codante = phase[i:j]
            if (len(seq_codante) > meilleure_seq):
                meilleure_seq = len(seq_codante)
                meilleur_start = start
                meilleur_stop = stop
                proteine_de_phase = phase[meilleur_start+1:meilleur_stop]
            i += 1
    return(proteine_de_phase)        

#main

mysterious_sequence = SeqIO.read(open('my_sequence.fasta'),
'fasta',
alphabet=generic_dna).seq


phase_1 = mysterious_sequence[0::]
phase_2 = mysterious_sequence[1::]
phase_3 = mysterious_sequence[2::]

#on renverse la sequence mystere pour creer les 3 dernieres phases 
complement_sequence = mysterious_sequence.complement()
reverse_sequence = complement_sequence[::-1]

#sequence à partir du dernier nucleotide en sens contraire
phase_4 = reverse_sequence[0::]
#sequence à partir de l'avant-dernier nucleotide en sens contraire
phase_5 = reverse_sequence[1::]
#sequence à partir de l'antépénultième nucleotide en sens contraire 
phase_6 = reverse_sequence[2::]

prot_1 = str(phase_1.translate())
prot_2 = str(phase_2.translate())
prot_3 = str(phase_3.translate())
prot_4 = str(phase_4.translate())
prot_5 = str(phase_5.translate())
prot_6 = str(phase_6.translate())

liste=[seq_codante_de_phase(prot_1),seq_codante_de_phase(prot_2),seq_codante_de_phase(prot_3),seq_codante_de_phase(prot_3),seq_codante_de_phase(prot_4),seq_codante_de_phase(prot_5),seq_codante_de_phase(prot_6)]
sorted(liste, key=len)
#la phase vraisemblable est la plus longue, donc la première de la liste triee par longueur
print (liste[1])

import functools
import gc
import itertools
import sys
from timeit import default_timer as _timer


def timeit(_func=None, *, repeat=3, number=1000, file=sys.stdout):
    """Decorator: prints time from best of `repeat` trials.
    Mimics `timeit.repeat()`, but avg. time is printed.
    Returns function result and prints time.
    You can decorate with or without parentheses, as in
    Python's @dataclass class decorator.
    kwargs are passed to `print()`.
    >>> @timeit
    ... def f():
    ...     return "-".join(str(n) for n in range(100))
    ...
    >>> @timeit(number=100000)
    ... def g():
    ...     return "-".join(str(n) for n in range(10))
    ...
    """

    _repeat = functools.partial(itertools.repeat, None)

    def wrap(func):
        @functools.wraps(func)
        def _timeit(*args, **kwargs):
            # Temporarily turn off garbage collection during the timing.
            # Makes independent timings more comparable.
            # If it was originally enabled, switch it back on afterwards.
            gcold = gc.isenabled()
            gc.disable()

            try:
                # Outer loop - the number of repeats.
                trials = []
                for _ in _repeat(repeat):
                    # Inner loop - the number of calls within each repeat.
                    total = 0
                    for _ in _repeat(number):
                        start = _timer()
                        result = func(*args, **kwargs)
                        end = _timer()
                        total += end - start
                    trials.append(total)

                # We want the *average time* from the *best* trial.
                # For more on this methodology, see the docs for
                # Python's `timeit` module.
                #
                # "In a typical case, the lowest value gives a lower bound
                # for how fast your machine can run the given code snippet;
                # higher values in the result vector are typically not
                # caused by variability in Python’s speed, but by other
                # processes interfering with your timing accuracy."
                best = min(trials) / number
                print(
                    "Best of {} trials with {} function"
                    " calls per trial:".format(repeat, number)
                )
                print(
                    "Function `{}` ran in average"
                    " of {:0.3f} seconds.".format(func.__name__, best),
                    end="\n\n",
                    file=file,
                )
            finally:
                if gcold:
                    gc.enable()
            # Result is returned *only once*
            return result

        return _timeit

    # Syntax trick from Python @dataclass
    if _func is None:
        return wrap
    else:
        return wrap(_func)

import numpy as np
import os


def genlouvain(B, seed=None):
    '''
    The optimal community structure is a subdivision of the network into
    nonoverlapping groups of nodes which maximizes the number of within-group
    edges and minimizes the number of between-group edges.
    This function is a fast an accurate multi-iterative generalization of the
    louvain community detection algorithm. This function subsumes and improves
    upon modularity_[louvain,finetune]_[und,dir]() and additionally allows to
    optimize other objective functions (includes built-in Potts Model i
    Hamiltonian, allows for custom objective-function matrices).
    Parameters
    ----------
    B : str | NxN np.arraylike
        string describing objective function type, or provides a custom
        objective-function matrix. builtin values 'modularity' uses Q-metric
        as objective function, or 'potts' uses Potts model Hamiltonian.
        Default value 'modularity'.
    seed : int | None
        random seed. default value=None. if None, seeds from /dev/urandom.
    Returns
    -------
    ci : Nx1 np.array
        final community structure
    q : float
        optimized q-statistic (modularity only)
    '''
    np.random.seed(seed)
    st0 = np.random.get_state()

    n = len(B)
    ci = np.arange(n) + 1
    Mb = ci.copy()

    B = np.squeeze(np.asarray(B))
    B = (B + B.T)/2.0
    Hnm = np.zeros((n, n))
    for m in range(1, n + 1):
        Hnm[:, m - 1] = np.sum(B[:, ci == m], axis=1)  # node to module degree
    H = np.sum(Hnm, axis=1)  # node degree
    Hm = np.sum(Hnm, axis=0)  # module degree

    q0 = -np.inf
    # compute modularity
    q = np.sum(B[np.tile(ci, (n, 1)) == np.tile(ci, (n, 1)).T])

    first_iteration = True

    while q - q0 > 1e-10:
        it = 0
        flag = True
        while flag:
            it += 1
            if it > 1000:
                raise ValueError('Modularity infinite loop style G. '
                                    'Please contact the developer.')
            flag = False
            for u in np.random.permutation(n):
                ma = Mb[u] - 1
                dQ = Hnm[u, :] - Hnm[u, ma] + B[u, u]  # algorithm condition
                dQ[ma] = 0

                max_dq = np.max(dQ)
                if max_dq > 1e-10:
                    flag = True
                    mb = np.argmax(dQ)
                    Hnm[:, mb] += B[:, u]
                    Hnm[:, ma] -= B[:, u]  # change node-to-module strengths

                    Hm[mb] += H[u]
                    Hm[ma] -= H[u]  # change module strengths

                    Mb[u] = mb + 1

        _, Mb = np.unique(Mb, return_inverse=True)
        Mb += 1

        M0 = ci.copy()
        if first_iteration:
            ci = Mb.copy()
            first_iteration = False
        else:
            for u in range(1, n + 1):
                ci[M0 == u] = Mb[u - 1]  # assign new modules

        n = np.max(Mb)
        b1 = np.zeros((n, n))
        for i in range(1, n + 1):
            for j in range(i, n + 1):
                # pool weights of nodes in same module
                bm = np.sum(B[np.ix_(Mb == i, Mb == j)])
                b1[i - 1, j - 1] = bm
                b1[j - 1, i - 1] = bm
        B = b1.copy()

        Mb = np.arange(1, n + 1)
        Hnm = B.copy()
        H = np.sum(B, axis=0)
        Hm = H.copy()

        q0 = q
        q = np.trace(B) # compute modularity
        dir = 'output'
        if not os.path.isdir(dir): os.makedirs(dir)
        
        output = open (os.path.join(dir, 'gen_louvain_seed.txt'),'a+')
        print >> output, "%i" % (st0[1][0])


    return ci, q