# -*- coding: utf-8 -*-
"""
Calculando el contenido en GC en oriC y otros

@author: rodri
"""


oric="atcaatgatcaacgtaagcttctaagcatgatcaaggtgctcacacagtttatccacaacctgagtggatgacatcaagataggtcgttgtatctccttcctctcgtactctcatgaccacggaaagatgatcaagagaggatgatttcttggccatatcgcaatgaatacttgtgacttgtgcttccaattgacatcttcagcgccatattgcgctggccaaggtgacggagcgggattacgaaagcatgatcatggctgttgttctgtttatcttgttttgactgagacttgttaggatagacggtttttcatcactgactagccaaagccttactctgcctgacatcgaccgtaaattgataatgaatttacatgcttccgcgacgatttacctcttgatcatcgatccgattgaagatcttcaattgttaattctcttgcctcgactcatagccatgatgagctcttgatcatgtttccttaaccctctattttttacggaagaatgatcaagctgctgctcttgatcatcgtttc"

########### 1) Algunos ejercicios simples de control de flujo y acceso a cadenas
if oric[0]=='g':
    print oric[1]
elif len(oric)<600 and len(oric)%8==0:
    print oric[len(oric)-50:len(oric)]
else:
    nuc=oric[36]
nuc


########### 2) Contar el número de adeninas
numAde=0
for nucleotido in oric:
    if nucleotido=='a':
        numAde+=1
numAde

def Acontent(seq):
    numA=0.0
    for nuc in seq:
        if nuc=='a':
            numA+=1
    return numA/len(seq)

Acontent(oric)
        

########### 3) Calcular el contenido en GC
def GCcontent(seq):
    GC=0.0
    for nuc in seq:
        if nuc=='g' or nuc=='c' or nuc=='G' or nuc=='C':
            GC+=1
    return GC/len(seq)

GCcontent(oric)

######### 4) Repetir con el genoma completo de E. coli y otras bacterias
f=open("/Users/rodri/Documents/docencia/BIE/practicas/E-coli.txt", 'r')
ec=f.read()
GCcontent(ec)

f=open("/Users/rodri/Documents/docencia/BIE/practicas/Thermotoga-petrophila.txt", 'r')
tp=f.read()
GCcontent(tp)

f=open("/Users/rodri/Documents/docencia/BIE/practicas/Vibrio_cholerae.txt", 'r')
vc=f.read()
GCcontent(vc)

f=open("/Users/rodri/Documents/docencia/BIE/practicas/Salmonella_enterica.txt", 'r')
se=f.read()
GCcontent(se)

######### 5) Contando la diferencia entre g y c
def skewTotal(seq):
    g=0
    c=0
    for i in range(len(seq)):
        if(seq[i]=='G'): 
            g=g+1
        elif(seq[i]=='C'): 
            c=c+1   
    return(g-c)
skewTotal(ec)


def skew(seq):
    count=0
    sk=[]
    sk.append(count)
    for i in range(len(seq)):
        if(seq[i]=='G'):
            count=count+1
        elif(seq[i]=='C'):
            count=count-1
        sk.append(count)    
    return(sk)
    
ecc=skew(ec)

def minimumSkew(seq):
    sk=skew(seq)
    mins=[]
    msk=min(sk)
    for i in range(len(sk)):
        if(sk[i]==msk):
            mins.append(i)
    return(mins)
minimumSkew(ec)

######### 6) Viendo la diferecia entre g y c
import matplotlib.pyplot as plt

plt.figure("Escherichia coli")
plt.plot(skew(ec))
plt.xlabel("sequence")
plt.ylabel("skew")
plt.show()

plt.plot(skew(se))
plt.title("Salmonella enterica")
plt.xlabel("sequence")
plt.ylabel("skew")
plt.show()


plt.figure("Thermotoga petrophila")
plt.plot(skew(tp))
plt.xlabel("sequence")
plt.ylabel("skew")
plt.show()