NumPy

class: center, middle, inverse, title-slide

# NumPy
### Licenciatura en Ciencias Genómicas,UNAM
### First version: 2021-10-10; Last update: 2021-10-26

---

## Contenido de la unidad

1. Ejemplo: scRNAseq
  
  2. Seeds
  
  3. Ejercicio: Busqueda sitios DNA 
  
---

## Objetivo

Familiarizarse con los contenidos aprendidos sobre NumPy con algunos ejerc

---
## Ejemplo scRNAseq

Tomamos un subset de nuestros datos de scRNA-Seq para explicar cómo funciona. Trabajaremos con datos de 6 células del tejido XX

.pull-left[
<img src="imgs/clase_5_pt2/scrna.jpeg" width="250px" style="display: block; margin: auto;" />
]
.pull-right[
<img src="imgs/clase_5_pt2/scrna3.jpeg" width="3000px" style="display: block; margin: auto;" />
]
<br>
**- ¿Qué podemos observar de nuestros datos?**

**- ¿Qué podemos obtener de los datos?**
---

## Datos de scRNAseq

Podemos leer nuestros datos con numpy con `genfromtxt` o `loadtxt`
siendo `loadtxt` más rápida pero menos flexible (ej. problemas con valores faltantes)

```python
import numpy as np
count_matrix = np.genfromtxt("files/scrnaseq.csv", dtype="int", delimiter=",") 
cm2 =np.loadtxt("files/scrnaseq.csv", delimiter=",", dtype = "int")
np.array_equal(count_matrix, cm2)
```

```
## True
```
--
Tambien podemos ingresarlos de esta manera

```python
import numpy as np
count_matrix = np.array([ [3, 3, 0],
                          [0, 0, 1],
                          [1, 1, 0],
                          [0, 0, 1],
                          [1, 0, 4],
                          [1, 2, 0]])
```

---
## ¡Nos enfocaremos en sacar una matríz de coexpresión!

Por lo que nos interesa saber en qué células se encuentran coexistiendo

Haremos una modificación a nuestra matriz, un 1 representando expresión y un 0 representando no expresión

```python
# cuando sea mayor a 0 será 1, si no será 0
count_matrix = np.where(count_matrix > 0, 1, 0)
count_matrix
```

---

## Gen 1 vs Gen 1

Para adentrarnos haremos la siguiente demostración

```python
print(count_matrix.T[0]) # checamos los valores de gen1
```

```
## [3 0 1 0 1 1]
```

Sacamos el producto punto de gen 1 contra gen 1

```python
print(np.dot(count_matrix.T[0], count_matrix.T[0]))
```

```
## 12
```

---

Ahora veremos de dónde sale ese 4 (núm de interesecciones)

```python

print(np.vstack((count_matrix.T[0],count_matrix.T[0])) )
```

```
## [[3 0 1 0 1 1]
##  [3 0 1 0 1 1]]
```

Podríamos sacar un resultado similar con operaciones booleanas

```python
print((count_matrix.T[0]) & (count_matrix.T[0]))
```

```
## [3 0 1 0 1 1]
```

Claro que tenemos que sumar dichos números (¡nos da 4!)

```python
print(sum((count_matrix.T[0]) & (count_matrix.T[0])))
```

```
## 6
```

---

### Gen 1 vs Gen 3

Hagamos la misma comparación, ahora el gen 1 y 3

```python
print(np.dot(count_matrix.T[0], count_matrix.T[2]))
```

```
## 4
```

Solo hay una posición en dónde intersectan

```python
print(np.vstack((count_matrix.T[0],count_matrix.T[2])) ) # para visualizar
```

```
## [[3 0 1 0 1 1]
##  [0 1 0 1 4 0]]
```

Sacamos lo mismo con operadores booleanos

```python
sum((count_matrix.T[0]) & (count_matrix.T[2]))
```

```
## 0
```
---

## multiplicación de matriz

Todo lo anterior se puede resumir en esta operación

```python
# np.matmul(count_matrix.T, count_matrix)
np.dot(count_matrix.T, count_matrix)
```

```
## array([[12, 12,  4],
##        [12, 14,  0],
##        [ 4,  0, 18]])
```

Esta es nuestra matríz de coexpresión. Al analizarla podremos concluir que...

---
## Seeds

Pensando en manejo de valores aleatorios, ¿Qué hacer si necesitamos obtener resultados reproducibles y estos cambian cada que corres el código?

**Usar un número determinado de seed nos permitirá obtener los mismos resultados, siempre que la llamemos antes de pedir un valor**

```python
import numpy as np
# seed = 10
np.random.seed(10)
print(np.random.rand(1), np.random.rand(1), np.random.rand(1), np.random.rand(1))
# Volvamos a poner seed en 10 
```

```
## [0.77132064] [0.02075195] [0.63364823] [0.74880388]
```

```python
np.random.seed(10) 
print(np.random.rand(1)) # mismo resultado que primer intento
```

```
## [0.77132064]
```

---
## Ejercicio

Análisis estadístico básico de búsqueda de secuencias de interés en un set dado:

- Crear función para obtener un set de secuencias de DNA aleatorias

- Buscar X patrón en dichas secuencias, obtener información(cantidad de veces encontrada, ¿posición?)

- Matriz de coexpresión (¿qué patrones coexisten en cada secuencia? ¿cuáles no?)

---
## Secuencias aleatorias DNA

Crearemos función usando NumPy y evitando ciclos for

```python
from Bio.Seq import Seq
import numpy as np

def secuencia_aleatoria(tamano = 100, p = [0.5, 0.5, 0.5, 0.5],seed = None):
    np.random.seed(seed) # posibilidad de reproducibilidad
    DNA = list("ATGC")
    #secuencia random con distribucion p
    secuencia = Seq(''.join(np.random.choice(DNA, tamano, p)))
    return(secuencia)

# probemos
secuencia = secuencia_aleatoria(25, p=[0.1,0.2,0.4,0.3])
secuencia
```

```
## Seq('GATAAATTTCCGTGACGTACGAAAC')
```
---

## Cómo sería sin numpy?

```python
def secuencia_aleatoria2(tamano=100):
  return(Seq(''.join( random.choice("AGTC") for x in range(tamano) )))
```
--
Otra manera sería:

```python
def secuencia_aleatoria3(tamano=100):
  sequence=""
  for i in range(tamano):
      sequence+=random.choice("AGTC")
  return(Seq(sequence))
```

---
## Guardemos las secuencias

```python
def seq_aleatorias(tamano = 100, p = [0.5, 0.5, 0.5, 0.5], seed = None,
                   num_seq = 20, archivo = "SeqAleat.fasta"):
    #escribamos archivo de secuencias!!
    with open(archivo, 'w') as file:
        for i in range(num_seq):
            sec = secuencia_aleatoria(tamano, p, seed)
            file.write(">Seq" + str(i) + "\n")
            file.write(str(sec))
            file.write("\n")
    file.close()
```

---

## Creando matriz! .tiny[buscando patrones]

```python
from Bio import SeqIO, SeqUtils
import csv
import re
# TFs a buscar
tf_interes_mod = { "TF_1": 'ATG[GG|TAG]', "TF_2" : 'T[TC|AA]GAAT', 
                  "TF_3" : "GTATGCGGGG", "TF_4" : "TAT[GT]CC",
                  "TF_5" : "TATATA[GT|TG]" }
              
outfile = open("ejercicioNP.csv", "w")
writer = csv.writer(outfile)

for rec in SeqIO.parse("SeqAleat.fasta", "fasta"):
    #secuencia por analizar
    seq = rec.seq
    #array para guardar numero de coindicencias
    tf_counts = np.empty((0), dtype = "int")
    for tf in tf_interes_mod.values():
        counts = len( re.findall(tf, str( seq ) ) )
        tf_counts = np.append(tf_counts, [counts])

# escribir counts encontrados (iterable)
    writer.writerows( [tf_counts] )
outfile.close()
```
---
## Leamos nuestro archivo

usaremos `loadtxt`, especificamos delimitador y tipo de dato

```python
count_matrix =np.loadtxt("files/ejercicioNP.csv", delimiter=",", dtype = "int")
#shape (20,5)
count_matrix
```

```
## array([[ 7,  1,  0,  1,  0],
##        [11,  0,  0,  0,  0],
##        [ 7,  0,  0,  0,  0],
##        [ 5,  0,  0,  0,  0],
##        [ 9,  0,  0,  1,  0],
##        [11,  1,  0,  1,  0],
##        [ 5,  0,  0,  0,  0],
##        [10,  0,  0,  0,  0],
##        [ 6,  1,  0,  0,  0],
##        [11,  2,  0,  0,  0],
##        [ 9,  0,  0,  0,  0],
##        [10,  0,  0,  1,  0],
##        [ 8,  0,  0,  0,  1],
##        [ 7,  1,  0,  0,  0],
##        [10,  1,  0,  2,  0],
##        [ 7,  0,  0,  0,  1],
##        [ 7,  1,  0,  1,  0],
##        [14,  0,  0,  0,  0],
##        [ 9,  1,  0,  0,  0],
##        [ 4,  0,  0,  0,  0]])
```

---

## Chequemos suma de columnas y filas

```python
# total de cada TF
print(np.sum(count_matrix, axis = 0))
```

```
## [167   9   0   7   2]
```

```python
# cuantos TFs se encontraron por secuencia
print(np.sum(count_matrix, axis = 1))
```

```
## [ 9 11  7  5 10 13  5 10  7 13  9 11  9  8 13  8  9 14 10  4]
```

---

## Visualicemos nuestra matriz de coexpresion de TFs

Lo que haremos será volver cambiar los valores de nuestra matriz. A 1 si se encontró y 0 si no

```python
from matplotlib.pyplot import imshow
# matriz 1s y 0s
matrix_uno = np.where(count_matrix > 0, 1, 0 )

coexpresion = np.matmul(matrix_uno.T, matrix_uno)
print(coexpresion)
```

```
## [[20  8  0  6  2]
##  [ 8  8  0  4  0]
##  [ 0  0  0  0  0]
##  [ 6  4  0  6  0]
##  [ 2  0  0  0  2]]
```

```python
#si solo te interesa ver los sitios no compartidos 
# plt.imshow(resultado.astype(bool))
plt.imshow(coexpresion)
plt.colorbar()
plt.show()
```
---
<br><br><br>

<img src="imgs/clase_5_pt2/tfMat.jpeg" width="350px" style="display: block; margin: auto;" />
---

#### Otra manera de visualizar lo anterior seria con un grafo

```python
import networkx as nx

G = nx.DiGraph(coexpresion)
nx.draw(G, node_size = 900,  with_labels=True)
```
<img src="imgs/clase_5_pt2/grafo.jpeg" width="300px" style="display: block; margin: auto;" />