---
title: "gkm-SVM"
author: "Giannina Bambaren Capurro"
date: "24/6/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=FALSE)
```

### Búsqueda de secuencias 
El motivo de P1 lo buscamos en todas las cepas encontradas de la SP  feathery mottlle, generamos secuencias solo de 200 nts  en las posiciones de 1300-1500 ya que es aqu? donde se encuentra.

```{bash}
esearch -db nucleotide -query "Sweet potato feathery mottle virus [ORGN] complete genome" | efetch -format fasta -seq_start 1300 -seq_stop 1500 > feath_200t.fasta
grep ">" feath_200t.fasta
```

Para buscar el motivo de P3 realizamos el mismo procedimiento anterior pero esta vez buscamos el motivo en todas las secuencias de potyvirus la longitud de búsqueda es de 2400 nts ya que se encuentra en distintas posiciones dependiendo de la especie
```{bash}
esearch -db nucleotide -query ?Potyvirus [ORGN]? | efetch -format fasta -seq_start 2400 -seq_stop 4800 > potyvirus_p3_meme.fasta

```

#### Para buscar coincidencias en fimo

```{bash}
# Todas las cepas de genoma completo de Feathery motte
esearch -db nucleotide -query "Sweet potato feathery mottle  virus [ORGN] complete genome" | efetch -format fasta > swp.fasta 
grep ">" swp.fasta
grep ">" swp.fasta | wc -l


# Todas las secuencias/cepas de genoma completo del g?nero Potyvirus
# - efilter : restringue los par?metros de la consulta

esearch -db nucleotide -query "virus complete genome NOT (partial cds)" | efilter -query "Potyvirus[ORGN]"| efetch -format fasta > potyvirus_fimo.fasta
grep ">" potyvirus_fimo.fasta | head
grep ">" potyvirus_fimo.fasta | wc -l

#Todas las secuencias/cepas de de genoma completo de virus h?espede de Sweet Potato

esearch -db nucleotide -query "Sweet potato virus complete genome NOT (partial cds)" | efilter -query ?Potyvirus [ORGN]?| efetch -format fasta > sweetpotato_fimo.fasta
grep ">" sweetpotato_fimo.fasta.fasta | head
grep ">" sweetpotato_fimo.fasta.fasta | wc -l

# Todas las cepas de genoma completo de la familia Potyviridae
esearch -db nucleotide -query "virus complete genome NOT (partial cds)" | efilter -query ?Potyviridae[ORGN]?| efetch -format fasta > sweetpotato_fimo.fasta
grep ">" sweetpotato_fimo.fasta.fasta | head
grep ">" sweetpotato_fimo.fasta.fasta | wc -l

```

#### Obtención de las Motivos con MEME ####

```{bash}
## Instalaci?n MEME-5.1.1 
# instalar dependencias
sudo apt-get install zlib1g-dev
sudo apt-get install autoconf
sudo apt-get install automake
sudo apt-get install libtool
sudo apt-get install default-jdk
# descargar de: http://meme-suite.org/doc/download.html
tar zxf meme-5.1.1.tar.gz
cd meme-5.1.1
./configure --prefix=$HOME/meme --with-url=http://meme-suite.org/ --enable-build-libxml2 --enable-build-libxslt
make
make test
make install
export PATH=$HOME/meme/bin:$HOME/meme/libexec/meme-5.1.1:$PATH
```

#### Parámetros utilizados

 - meme (file input)
 - (-dna) tipo de secuencia
 - (-oc) directorio para la salida
 - (-nostatus) no genera resumen en la terminal
 - (-time segs) tiempo en el que dejar? de buscar motivos si ha encontrado al menos 1, si ha superado "x seg." del CPU
 - (-mod) modo, "zoops" para encontrar motivos con al menos una ocurrencia
 - (-nmotifs) n?mero de motivos que desea obtener
 - (-minw) minima longitud del motivo
 - (-maxw) m?xima longitud del motivo
 - (-objfun) "clasic" para puntuar los motivos utilizando E-value
 - (-revcomp) para obtener el complemente reverso del motivo
 - (-markow_order) para usar el modelo de m?ximo orden de markov si no se especifica nada el programa usa el modelo "fasta-get-markov"
 
#### P1
```{bash}
mkdir /home/uoc/meme_out/meme2
meme /home/uoc/feath_200t.fasta -dna -oc /home/uoc/meme_out/meme2 -nostatus -time 3000 -mod zoops -nmotifs 20 -minw 6 -maxw 16 -objfun classic -revcomp -markov_order 0
mkdir /home/uoc/meme_out/meme3
meme /home/uoc/feath_200t.fasta -dna -oc /home/uoc/meme_out/meme3 -nostatus -time 3000 -mod zoops -nmotifs 20 -minw 6 -maxw 40 -objfun classic -revcomp -markov_order 0
```

#### P3

```{bash}
mkdir /home/uoc/meme_out/meme3_p3
meme /home/uoc/potyvirus_p3_meme.fasta -dna -oc /home/uoc/meme_out/meme3_p3 -mod zoops -nmotifs 30 -minw 20 -maxw 30 -objfun classic 
```

### FIMO: 

#### obtención del motivo en formato meme

Una vez obtenido el motivo que m?s se ajusta al de deslizamiento de la polimerasa, utilizamos la documentación de MEME para crear  el archivo que contenga el motivo en formato "MEME"  http://meme-suite.org/doc/meme-format.html  generamos los archivos:

- meme4_motif.meme -> para el motivo P1
- memep3.meme -> para P3

#### P1

```{bash}
mkdir /home/uoc/meme_out/fimo2
fimo --oc /home/uoc/meme_out/fimo2 /home/uoc/meme4_motif.meme /home/uoc/swp.fasta
mkdir /home/uoc/meme_out/fimo3
fimo --oc /home/uoc/meme_out/fimo3 /home/uoc/meme4_motif.meme /home/uoc/sweetpotato_fimo.fasta
mkdir /home/uoc/meme_out/fimo4
fimo --oc /home/uoc/meme_out/fimo4 /home/uoc/meme4_motif.meme /home/uoc/potyvirus_fimo.fasta
```

#### p3

```{bash}

mkdir /home/uoc/meme_out/fimo_p3_1
fimo --oc /home/uoc/meme_out/fimo_p3_1 /home/uoc/memep3.meme /home/uoc/sweetpotato_fimo.fasta
mkdir /home/uoc/meme_out/fimo_p3_2
fimo --oc /home/uoc/meme_out/fimo_p3_2 /home/uoc/memep3.meme /home/uoc/potyvirus_fimo.fasta
mkdir /home/uoc/meme_out/fimo_p3_3
fimo --oc /home/uoc/meme_out/fimo_p3_3 /home/uoc/memep3.meme /home/uoc/potyviridae_fimo.fasta

```


### Análisis de FIMO en R ####

Los matches obtenidos por FIMO son una aproximaci?n de la b?squeda de motivo por la secuencia
Vemos que los pvalues/qvalues obtenidos por fimo son muy altos, sin embargo muchas de las secuencias  no presentan correctamente el patrón esperado si nos basamos en los scores obtenidos por fimo vemos que estos se encuentran en un rango de 40.1 y 94.4  min y máx respectivamente, si lo comparamos con los scores de las secuencias que cumplen el patron vemos que  podríamos tener algunos falsos negativos significativos si poner un umbral de discriminaci?n por score de modo que hemos optado por generar una segunda b?squeda del patr?n de deslizamiento dentro de los  patrones obtenidos con FIMO.

De forma que tenemos dos set de secuencias postivas. 
1) 2da busqueda 2) seq. FIMO

```{r}
##### leer arhivo ####
# Generamos una funci?n que nos convierte la salida de FIMO.GGF en un dataframe
fimo.gff<-function(fimo.file){
  require(rtracklayer)
  my_tags <- c("Name", "Alias", "ID", "pvalue", "qvalue","sequence")
 df<-as.data.frame(readGFF(fimo.file, tags=my_tags))
}
#
#
##### crear las positivas con el patron correcto ####
# require como formato de entrada un archivo en formato fimo.gff en forma de df
# genera una doble b?squeda del patr?n de deslizamiento de los matches de fimo
posit2search<-function(name.df) {
  library(Biostrings)
  seq_pos<-c()
  for (i in 1:nrow(name.df)) {
    data <- name.df$sequence[i]
    dna <- DNAString(data)
    #busca el patron deseado dentro de los resultados obtenidos por fimo
    results<-matchPattern('GGAAAAAA', dna)
    empty<-isEmpty(results)
    if (empty=="FALSE")
    {
      positive<-data.frame(seqid=name.df$seqid[i],
                          score=name.df$score[i],
                          ID=name.df$ID[i],
                          start=name.df$start[i],
                          stop=name.df$end[i],
                          sequence=name.df$sequence[i])
      seq_pos<-rbind.data.frame(seq_pos,positive)
      
    }
    
  }
  return(seq_pos)
}

###### postivas 70 train / 30 test #####
# con los dataframes obtenidos generaremos la partici?n de las secuencias
# mediante caret crearemos una partici?n del 70% de las secuencias para entrenamiento 
# y 30% para evalular
library(caret)
#agregamos una semilla de aletorizaci?n
set.seed(123654)
library(seqinr)
# usamos el paquete caret para crear la particion de los datos en test y train
# establecemos un 70% de los datos para entrenamiento
# require el df para la partici?n
makeTrain<-function(df,p){
  train<- createDataPartition(df$score, 
                            p = p, 
                            list = FALSE, 
                            times = 1)
  return(train)
}
#### Secuencias positivas ####
# esta funci?n crea las secuencias positivas test y train
# requiere del df_train / df_test y un nombre para el archivo de salida ".fasta"
makePositive<-function(df,outfasta){
  for (i in 1:length(df$seqid)) {
    header<-paste(df$seqid[i],"_",
                  (df$start[i]),"_",
                  (df$stop[i]),"_rep",i,sep = "")
    write.fasta(sequences = df$sequence[i], 
                names = header, 
                file.out = outfasta,
                open = "a")
  }
}
```

Negativas  a partir de las secuencias de la particion 
 a partir de las secuencias positivas creamos un archivo con las seq. negativas
 las coordenadas de estas secuencias es completamente al azar.
 
Para esta función primero tenemos que cargar el archivo de origen donde se obtuvieron los matches de FIMO

```{r}
genomecomplete<-read.fasta(file = "genoma.fasta",seqtype = "DNA",as.string = FALSE)
#
# La funci?n require:
# - df_trais/df_test
# - nombre del archivo de salida
# - lista con el genoma completo
makenegative_p1<-function(df,outfasta,genomecomplete){
  for (i in 1:length(df$seqid)){
    rf<-df$seqid[i]
    numalt<-sample(1:9000,1,replace=F)
    genome<-unlist(genome<-seqinr::getSequence(object=genomecomplete[rf]))[numalt:(numalt+15)]
    header<-paste(df$seqid[i],"_",numalt,"_",(numalt+15),"_rep",i,sep = "")
    write.fasta(sequences = genome, names = header, file.out = outfasta,open = "a")
  }
}
```

Algoritmo gkmsvm para virus h?espedes de sweet potato con doble b?squeda de patrón

Una vez obetenidas las secuencias negativas y positivas, vamos a entrenar el algoritmo.
Antes del entrenamiento tenemos que generar un matriz kernel con las secuencias negativos y las positivas  generadas.

Para generar las secuencias que usaremos en el algoritmo debemos:

```{r}

# 1.cargar el fichero fimo.gff
old<-getwd()
setwd("GKMSVM/P1")
FM<-fimo.gff("fimo2/fimo.gff")
# distribución de los score de FIMO
library(ggplot2)
png("fimo_fm_p1.png")
ggplot(FM, aes(x=as.numeric(qvalue), group=1))+geom_density(color="darkblue")
dev.off()
# Extra. generamos un datafrene con doble b?squeda de patr?n
pos2s2<-posit2search(FM)
# 2. Generar la partici?n del df en train y test
train<-makeTrain(FM,p=0.70)
FM_train <- FM[ train,]
FM_test<- FM[-train,]
# 3. Crear las secuencias positivas
#setwd("/GKMSVM/P1")
dir.create(file.path("P1_FM"))
setwd("P1_FM")
makePositive(FM_test,"fm_pos_test.fasta")
makePositive(FM_train,"fm_pos_train.fasta")
# 4. Crear las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
library(seqinr)
swp_seq_complete<-read.fasta(file = "swp.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P1/P1_FM")
makenegative(FM_train,"neg_fm_train.fasta",swp_seq_complete)
```


### Algoritmo gkmsvm
Usando trainCV entreno el algoritmo, para ello se requiere:

 - matriz kernel
 - secuencias postivas de entrenamiento
 - secuencias ngativas de entrenamiento
 - nCv: numero de capas para la CV
 
De este algoritmo se obtiene:
"outputPDFfn" = gr?fico de la curva ROC de la predicci?n
"svmfnprfx" = score de la clasificaci?n entre las secuencias neg vs pos
"outputCVpredfn" = valores predictores CV 
"outputROCfn" = archivo con el valor roc para AUC y el valor de presici?n de rellamada

```{r}
# selecciona el directorio donde se guardar?n las salidas de este algoritmo
#setwd("/GKMSVM/P1/P1_FM")
library(gkmSVM)
# Creaci?n de la matriz kernel
gkmsvm_kernel('fm_pos_train.fasta',"neg_fm_train.fasta", 'fm_kernel.out')
# Usamos la salida "sw_pos2s_kernel.out" para entrenar el algoritmo
# usamos 5-fold CV
gkmsvm_trainCV('fm_kernel.out',
               'fm_pos_train.fasta',
               "neg_fm_train.fasta",
               nCV = 5,
               showPlots = TRUE,
               outputPDFfn = 'fm_roc.pdf',
               svmfnprfx ='fm_1X', 
               outputCVpredfn ='fm_pred.out', 
               outputROCfn ='fm_roc.out')
# Usamos el pesos de entrenamiento "fm_1x" para realizar la clasificaci?n
#setwd("P1_FM") 
gkmsvm_classify("fm_pos_test.fasta",
                svmfnprfx='fm_1x', 
                'fm_weights.out')
```

 Gr?fico de las secuencias vs. score de clasificaci?n
 
# 1. Generaremos un archivo fasta en forma de sliding windows con una de las secuencias de las cuales hemos utilizado para crear el motivo en MEME

Para ello creamos una funci?n que requiere como entrada el genoma completo 
Cargamos todo el genoma pero solo usamos 1 secuencias,
Tambi?n requiere un nombre para la salida del archivo, y los par?metros de inicio y final de la secuencia que queremos generar
generaremos secuencias en todos los casos de 10 nts ya que el la amplitud que usa por defecto el algoritmo

```{r}
setwd(old)
library(seqinr)

sliding_windows<-function(genome_compl,outfasta,p,q){
  genome_compl<-unlist(getSequence(genome_compl[1]))
  c= 0
  for (i in p:q) {
    c = c+1
    header<-paste("seq_",c,sep = "")
    genome <- genome_compl[i:(i+15)]
    write.fasta(sequences = genome, names = header, file.out = outfasta,open = "a")
  }
}

# Generamos la secuencia donde los par?metros son la region aproximada donde a donde se encuentra P1
setwd("GKMSVM/P1/P1_FM") 
sliding_windows(swp_seq_complete,"spfmv_sw.fasta",p=1,q=2000)
sliding_windows(swp_seq_complete,"spfmv_out.fasta",p=7000,q=9000)
# Clasificamos usando el modelo de entrenamiendo anterior
require(gkmSVM)
gkmsvm_classify("spfmv_sw.fasta",
                svmfnprfx='fm_1X', 
                'weight.out')
gkmsvm_classify("spfmv_out.fasta",
                svmfnprfx='fm_1X', 
                'out_weights.out')
# Cargamos el fichero obtenido que contiene las score de clasificaci?n
pos<-read.delim('fm_weights.out',header = FALSE)
png("gkmsvm_fm_p1.png")
ggplot(pos, aes(x=V2, group=1))+geom_density(color="red")+xlim(c(1.75,1.85))
dev.off()
# Cargamos el fichero obtenido que contiene las score de clasificación sliding window
pos<-read.delim('weight.out',header = FALSE)
pos1<-read.delim('out_weights.out',header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library("ggplot2")
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library(gridExtra)
png("FM_P1.png")
p<-ggplot(data =pos, aes(x =c(1:2000) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Plot of SPFMV [1:2000]")
p1<-ggplot(data =pos1, aes(x = c(7000:9000) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Plot of SPFMV [7000:9000]")
grid.arrange(p,p1, ncol = 1, nrow = 2)
dev.off()
```

#### P1: fimo contra todos los h?espedes de SP 

```{r}
setwd(old)
setwd("GKMSVM/P1")
# 1. Cargamos el fichero fimo en df
P1_SP<-fimo.gff(fimo.file = "fimo3/fimo.gff")
library(ggplot2)
png("fimo_sp_p1.png")
ggplot(P1_SP, aes(x=as.numeric(qvalue), group=1))+geom_density(color="darkblue")
dev.off()

# 2. Generamos la partici?n
train<-makeTrain(P1_SP,p=0.70)
P1_SP_train <- P1_SP[ train,]
P1_SP_test<- P1_SP[-train,]
# 3. Creamos las seq. positivas
dir.create(file.path("P1_SP"))
setwd("P1_SP")
makePositive(P1_SP_train,"sp_pos_train.fasta")
makePositive(P1_SP_test,"sp_pos_test.fasta")
# 4. Creamos las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
sp_seq_complete<-read.fasta(file = "sweetpotato_fimo.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P1/P1_SP")
makenegative(P1_SP_train,"sp_neg_train.fasta",sp_seq_complete)
makenegative(P1_SP_test,"sp_neg_test.fasta",sp_seq_complete)
require(gkmSVM)
# 5. Generamos el kernel
gkmsvm_kernel("sp_pos_train.fasta","sp_neg_train.fasta", "sp_kernel")
# 6. Entrenamos el algoritmo
gkmsvm_trainCV("sp_kernel",
               "sp_pos_train.fasta",
               "sp_neg_train.fasta",
               nCV = 5,
               showPlots = TRUE,
               outputPDFfn = "ROC_sp.pdf",
               svmfnprfx = "sp_1x", 
               outputCVpredfn = "sp_pred.out" , 
               outputROCfn ="sp_rocfn.out")
# 7. Clasificamos
gkmsvm_classify("sp_pos_test.fasta",
                svmfnprfx="sp_1x", 
                "sp_score.out")

# 8. Secuencias con mayor score
# Ordenamos los score de forma descendente
# Mostramos las secuencias que tienen el score m?s alto
# en el ID de las secuencias podemos ver su posici?n 
scores<-read.table("sp_score.out",col.names = c("seq_id","weight"))
png("gkmsvm_sp_p1.png")
ggplot(scores, aes(x=weight, group=1))+geom_density(color="red")
dev.off()
head(scores<-dplyr::arrange(scores,desc(scores$weight)))
head(scores<-dplyr::arrange(scores,(scores$weight)))
sliding_windows_R(sp_seq_complete,"sw_sp_good.fasta",1,2000,ref="MH782228.1")
gkmsvm_classify("sw_sp_good.fasta",
                svmfnprfx="sp_1x", 
                "sp_good_score.out")
# Id:KF734962.1 con peor score
sliding_windows_R(sp_seq_complete,"sw_sp_bad.fasta",1,2000,ref="MG656436.1")
gkmsvm_classify("sw_sp_bad.fasta",
                svmfnprfx="sp_1x", 
                "sp_bad_score.out")
# fuera del sitio
sliding_windows_R(pv_seq_complete,"sw_sp_out.fasta",5000,7000,ref="MH782227.1")
gkmsvm_classify("sw_sp_out.fasta",
                svmfnprfx="sp_1x", 
                "sp_out_score.out")
# 
pos1<-read.delim("sp_out_score.out",header = FALSE)
pos2<-read.delim('sp_good_score.out',header = FALSE)
pos3<-read.delim('sp_bad_score.out',header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library(ggplot2)
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library("gridExtra")
png("SP_P1.png ")
p1<-ggplot(data =pos1, aes(x = c(5000:7000), y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Out")
p2<-ggplot(data =pos2, aes(x = c(1:2000) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Good")
p3<-ggplot(data =pos3, aes(x = c(1:2000) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Bad")
grid.arrange(p2,p3,p1, ncol = 1, nrow = 3)
dev.off()
```

#### P1: fimo contra todos los Potyvirus (g?nero) 

```{r}
setwd(old)
setwd("GKMSVM/P1")
# 1. Cargamos el fichero fimo en df
P1_PV<-fimo.gff(fimo.file = "fimo4/fimo.gff")
png("fimo_pv_p1.png")
ggplot(P1_PV, aes(x=as.numeric(qvalue), group=1))+geom_density(color="darkblue")
dev.off()
# 2. Generamos la partici?n
train<-makeTrain(P1_PV,p=0.70)
P1_PV_train <- P1_PV[ train,]
P1_PV_test<- P1_PV[-train,]
# 3. Creamos las seq. positivas
dir.create(file.path("P1_PV"))
setwd("P1_PV")
makePositive(P1_PV_train,"pv_pos_train.fasta")
makePositive(P1_PV_test,"pv_pos_test.fasta")
# 4. Creamos las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
pv_seq_complete<-read.fasta(file = "potyvirus_fimo.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P1/P1_PV")
makenegative(P1_PV_train,"pv_neg_train.fasta",pv_seq_complete)
#makenegative(P1_PV_test,"pv_neg_test.fasta",pv_seq_complete)
require(gkmSVM)
# 5. Generamos el kernel
gkmsvm_kernel("pv_pos_train.fasta","pv_neg_train.fasta", "pv_kernel")
# 6. Entrenamos el algoritmo
gkmsvm_trainCV("pv_kernel",
               "pv_pos_train.fasta",
               "pv_neg_train.fasta",
               nCV = 5,
               showPlots = TRUE,
               outputPDFfn = "ROC_pv.pdf",
               svmfnprfx = "pv_1x", 
               outputCVpredfn = "pv_pred.out" , 
               outputROCfn ="pv_rocfn.out")
# 7. Clasificamos
gkmsvm_classify("pv_pos_test.fasta",
                svmfnprfx="pv_1x", 
                "pv_score.out")
# 8. Secuencias con mayor score
# Ordenamos los score de forma descendente
# Mostramos las secuencias que tienen el score más alto
# y el más bajo
scores<-read.table("pv_score.out",col.names = c("seq_id","weight"))
png("gkmsvm_pv_p1.png")
ggplot(scores, aes(x=weight, group=1))+geom_density(color="red")
dev.off()
head(scores<-dplyr::arrange(scores,(scores$weight)))
# Mostramos las secuencias peor clasificadas
head(dplyr::arrange(scores,desc(scores$weight)))
# Vamos a usar la función sliding windows para graficar los secuencias vs scores
# hemos modificado la función para poder añadirle el id de referencia de la secuencia
# genome_compl el mismo archivo con el que generamos las secuencias negativas
sliding_windows_R<-function(genome_compl,outfasta,p,q,ref){
  #genome_compl<-read.fasta(genome.fasta)
  genome_compl<-unlist(genome<-seqinr::getSequence(object=genome_compl[ref]))
  c= 0
  for (i in p:q) {
    c = c+1
    header<-paste("seq_",c,sep = "")
    genome <- genome_compl[i:(i+15)]
    write.fasta(sequences = genome, names = header, file.out = outfasta,open = "a")
  }
}
# usamos primero el ID:MH782227.1 de la secuencia con mayor score
sliding_windows(pv_seq_complete,"sw_potyvirus.fasta",1,2000,ref="MH782227.1")
gkmsvm_classify("sw_potyvirus.fasta",
                svmfnprfx="pv_1x", 
                "pv_good_score.out")
# Id:KF734962.1 con peor score
sliding_windows(pv_seq_complete,"sw_bad_potyvirus.fasta",1,2000,ref="KF734962.1")
gkmsvm_classify("sw_bad_potyvirus.fasta",
                svmfnprfx="pv_1x", 
                "pv_bad_score.out")
# fuera del sitio
sliding_windows_R(pv_seq_complete,"sw_out_potyvirus.fasta",5000,7000,ref="MH782227.1")
gkmsvm_classify("sw_out_potyvirus.fasta",
                svmfnprfx="pv_1x", 
                "pv_out_score.out")
#
# Cargamos el fichero obtenido que contiene las score de clasificaci?n

pos1<-read.delim('pv_good_score.out',header = FALSE)
pos2<-read.delim('pv_bad_score.out',header = FALSE)
pos3<-read.delim("pv_out_score.out",header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library(ggplot2)
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library("gridExtra")
png("PV_P1.png ")
p1<-ggplot(data =pos1, aes(x = c(1:2001) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Good")
p2<-ggplot(data =pos2, aes(x = c(1:2001) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Bad")
p3<-ggplot(data =pos3, aes(x = c(5000:7000) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Out")
grid.arrange(p1, p2,p3, ncol = 1, nrow = 3)
dev.off()
```

#### P3 secuencias hu?spedes de Sweet potato 

```{r}
setwd(old)
setwd("GKMSVM/P3")
# 1. Cargamos el fichero fimo en df
P3_SP<-fimo.gff(fimo.file = "fimo_p3_1/fimo.gff")
png("fimo_sp_p3.png")
ggplot(P3_SP, aes(x=as.numeric(qvalue), group=1))+geom_density(color="darkblue")
dev.off()
# 2. Generamos la partici?n
set.seed(123654)
train<-makeTrain(P3_SP,p=0.70)
P3_SP_train <- P3_SP[ train,]
P3_SP_test<- P3_SP[-train,]
# 3. Creamos las seq. positivas
dir.create(file.path("P3_SP"))
setwd("P3_SP")
makePositive(P3_SP_train ,"p3_sp_pos_train.fasta")
makePositive(P3_SP_test,"p3_sp_pos_test.fasta")
# 4. Creamos las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
sp_seq_complete<-read.fasta(file = "sweetpotato_fimo.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P3/P3_SP")
makenegative<-function(df,outfasta,genomecomplete){
  for (i in 1:length(df$seqid)){
    rf<-df$seqid[i]
    numalt<-sample(1:9000,1,replace=F)
    genome<-unlist(genome<-seqinr::getSequence(object=genomecomplete[rf]))[numalt:(numalt+29)]
    header<-paste(df$seqid[i],"_",numalt,"_",(numalt+29),"_rep",i,sep = "")
    write.fasta(sequences = genome, names = header, file.out = outfasta,open = "a")
  }
}
makenegative(P3_SP_train ,"p3_sp_neg_train.fasta",sp_seq_complete)
#makenegative(P3_SP_test,"p3_sp_neg_test.fasta",sp_seq_complete)
require(gkmSVM)
# 5. Generamos el kernel
gkmsvm_kernel("p3_sp_pos_train.fasta","p3_sp_neg_train.fasta", "p3_sp_kernel")
# 6. Entrenamos el algoritmo
gkmsvm_trainCV("p3_sp_kernel",
               "p3_sp_pos_train.fasta",
               "p3_sp_neg_train.fasta",
               nCV = 5,
               showPlots = TRUE,
               outputPDFfn = "ROC_p3_sp.pdf",
               svmfnprfx = "p3_sp_1x", 
               outputCVpredfn = "p3_sp_pred.out" , 
               outputROCfn ="p3_sp_rocfn.out")
# 7. Clasificamos
gkmsvm_classify("p3_sp_pos_test.fasta",
                svmfnprfx="p3_sp_1x", 
                "p3_sp_score.out")
# 8. Secuencias con mayor score
# Ordenamos los score de forma descendente
# Mostramos las secuencias que tienen el score m?s alto
# en el ID de las secuencias podemos ver su posici?n 
scores<-read.table("p3_sp_score.out",col.names = c("seq_id","weight"))
png("gkmsvm_sp_p3.png")
ggplot(scores, aes(x=weight, group=1))+geom_density(color="red")
dev.off()
head(scores<-dplyr::arrange(scores,desc(scores$weight)))
head(scores<-dplyr::arrange(scores,(scores$weight)))
# usamos primero el ID: de la secuencia con mayor score
sliding_windows(sp_seq_complete,"p3_sw_good.fasta",2400,4800,ref="MH023311.1")
gkmsvm_classify("p3_sw_good.fasta",
                svmfnprfx="p3_sp_1x", 
                "p3_good_score.out")
# Id: con peor score
sliding_windows(sp_seq_complete,"p3_sw_bad.fasta",2400,4800,ref="MH763687.1")
gkmsvm_classify("p3_sw_bad.fasta",
                svmfnprfx="p3_sp_1x", 
                "p3_bad_score.out")
sliding_windows(sp_seq_complete,"p3_out_bad.fasta",6000,8400,ref="MH023311.1")
gkmsvm_classify("p3_out_bad.fasta",
                svmfnprfx="p3_sp_1x", 
                "p3_out_score.out")
# Cargamos el fichero obtenido que contiene las score de clasificaci?n
pos4<-read.delim('p3_good_score.out',header = FALSE)
pos5<-read.delim('p3_bad_score.out',header = FALSE)
pos6<-read.delim('p3_out_score.out',header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library(ggplot2)
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library("gridExtra")
png("SP_P3.png ")
p4<-ggplot(data =pos4, aes(x = c(2400:4800) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Good")
p5<-ggplot(data =pos5, aes(x = c(2400:4800) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Bad")
p6<-ggplot(data =pos6, aes(x = c(6000:8400) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Out")
grid.arrange(p4, p5,p6, ncol = 1, nrow = 3)
dev.off()

```


#### P3 fimo contra todo los Potyvirus (g?nero)

```{r}
setwd(old)

setwd("GKMSVM/P3")
# 1. Cargamos el fichero fimo en df
P3_PV<-fimo.gff(fimo.file = "fimo_p3_2/fimo.gff")
png("fimo_pv_p3.png")
ggplot(P3_PV, aes(x=as.numeric(qvalue), group=1))+geom_density(color="darkblue")
dev.off()
# 2. Generamos la partici?n
train<-makeTrain(P3_PV,p=0.70)
P3_PV_train <- P3_PV[ train,]
P3_PV_test<- P3_PV[-train,]
# 3. Creamos las seq. positivas
dir.create(file.path("P3_PV"))
setwd("P3_PV")
makePositive(P3_PV_train ,"p3_pv_pos_train.fasta")
makePositive(P3_PV_test,"p3_pv_pos_test.fasta")
# 4. Creamos las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
pv_seq_complete<-read.fasta(file = "potyvirus_fimo.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P3/P3_PV")
makenegative(P3_PV_train ,"p3_pv_neg_train.fasta",pv_seq_complete)
#makenegative(P3_PV_test,"p3_pv_neg_test.fasta",pv_seq_complete)
require(gkmSVM)
# 5. Generamos el kernel
gkmsvm_kernel("p3_pv_pos_train.fasta","p3_pv_neg_train.fasta", "p3_poty_kernel")
# 6. Entrenamos el algoritmo
gkmsvm_trainCV("p3_poty_kernel",
               "p3_pv_pos_train.fasta",
               "p3_pv_neg_train.fasta",
               showPlots = TRUE,
               outputPDFfn = "ROC_p3_pv.pdf",
               svmfnprfx = "p3_pv_1x", 
               outputCVpredfn = "p3_pv_pred.out" , 
               outputROCfn ="p3_pv_rocfn.out")
# 7. Clasificamos
gkmsvm_classify("p3_pv_pos_test.fasta",
                svmfnprfx="p3_pv_1x", 
                "p3_pv_score.out")
# 8. Secuencias con mayor score
# Ordenamos los score de forma descendente
# Mostramos las secuencias que tienen el score m?s alto
# en el ID de las secuencias podemos ver su posici?n 
scores<-read.table("p3_pv_score.out",col.names = c("seq_id","weight"))
png("gkmsvm_pv_p3.png")
ggplot(scores, aes(x=weight, group=1))+geom_density(color="red")
dev.off()
head(scores<-dplyr::arrange(scores,desc(scores$weight)))
head(scores<-dplyr::arrange(scores,(scores$weight)))
# usamos primero el ID: de la secuencia con mayor score
sliding_windows(pv_seq_complete,"p3_pv_sw_good.fasta",2400,4800,ref="MH023311.1")
gkmsvm_classify("p3_pv_sw_good.fasta",
                svmfnprfx="p3_pv_1x", 
                "p3_pv_good_score.out")
# Id: con peor score
sliding_windows(pv_seq_complete,"p3_pv_sw_bad.fasta",2400,4800,ref="MG602234.1")
gkmsvm_classify("p3_pv_sw_bad.fasta",
                svmfnprfx="p3_pv_1x", 
                "p3_pv_bad_score.out")
# fuera de p3
sliding_windows_R(pv_seq_complete,"p3_pv_sw_out.fasta",7000,8400,ref="MG602234.1")
gkmsvm_classify("p3_pv_sw_out.fasta",
                svmfnprfx="p3_pv_1x", 
                "p3_pv_out_score.out")
# Cargamos el fichero obtenido que contiene las score de clasificaci?n
pos1<-read.delim('p3_pv_good_score.out',header = FALSE)
pos2<-read.delim('p3_pv_bad_score.out',header = FALSE)
pos3<-read.delim('p3_pv_out_score.out',header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library(ggplot2)
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library("gridExtra")
png("PV_P3.png ")

p1<-ggplot(data =pos1, aes(x = c(2400:4800) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Good")
p2<-ggplot(data =pos2, aes(x = c(2400:2871) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Bad")
p3<-ggplot(data =pos3, aes(x = c(7000:8400), y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Out")
grid.arrange(p1, p2, p3, ncol = 1, nrow = 3)
dev.off()
```

#### P3 fimo contra todo los Potyviridae (familia)

```{r}
setwd(old)
setwd("GKMSVM/P3")
# 1. Cargamos el fichero fimo en df
P3_PD<-fimo.gff(fimo.file = "fimo_p3_3/fimo.gff")
png("fimo_pd_p3.png")
ggplot(P3_PD, aes(x=as.numeric(qvalue), group=1))+geom_density(color="red")
dev.off()
# 2. Generamos la partici?n
train<-makeTrain(P3_PD,p=0.70)
P3_PD_train <- P3_PD[ train,]
P3_PD_test<- P3_PD[-train,]
# 3. Creamos las seq. positivas
dir.create(file.path("P3_PD"))
setwd("P3_PD")
makePositive(P3_PD_train ,"p3_pd_pos_train.fasta")
makePositive(P3_PD_test,"p3_pd_pos_test.fasta")
# 4. Creamos las negativas 
setwd(old)
setwd("GKMSVM")
# cargamos el fichero donde se encuentra el genoma completo
pd_seq_complete<-read.fasta(file = "potyviridae_fimo.fasta",seqtype = "DNA",as.string = FALSE)
setwd("P3/P3_PD")
makenegative(P3_PD_train ,"p3_pd_neg_train.fasta",pd_seq_complete)
#makenegative(P3_PD_test,"p3_pd_neg_test.fasta",pd_seq_complete)
require(gkmSVM)
# 5. Generamos el kernel
gkmsvm_kernel("p3_pd_pos_train.fasta","p3_pd_neg_train.fasta", "p3_pd_kernel")
# 6. Entrenamos el algoritmo
gkmsvm_trainCV("p3_pd_kernel",
               "p3_pd_pos_train.fasta",
               "p3_pd_neg_train.fasta",
               showPlots = TRUE,
               outputPDFfn = "ROC_p3_pd.pdf",
               svmfnprfx = "p3_pd_1x", 
               outputCVpredfn = "p3_pd_pred.out", 
               outputROCfn ="p3_pd_rocfn.out")
# 7. Clasificamos
gkmsvm_classify("p3_pd_pos_test.fasta",
                svmfnprfx="p3_pd_1x", 
                "p3_pd_score.out")
# 8. Secuencias con mayor score
# Ordenamos los score de forma descendente
# Mostramos las secuencias que tienen el score m?s alto
# en el ID de las secuencias podemos ver su posici?n 
scores<-read.table("p3_pd_score.out",col.names = c("seq_id","weight"))
png("gkmsvm_sp_p3.png")
ggplot(scores, aes(x=weight, group=1))+geom_density(color="red")
dev.off()
head(scores<-dplyr::arrange(scores,desc(scores$weight)))
head(scores<-dplyr::arrange(scores,(scores$weight)))
sliding_windows_R(pd_seq_complete,"p3_pd_sw_good.fasta",2400,4800,ref="KX279878.1")
gkmsvm_classify("p3_pd_sw_good.fasta",
                svmfnprfx="p3_pd_1x", 
                "p3_pd_good_score.out")
# Id: con peor score
sliding_windows_R(pd_seq_complete,"p3_pd_sw_bad.fasta",2400,4800,ref="LC228573.1")
gkmsvm_classify("p3_pd_sw_bad.fasta",
                svmfnprfx="p3_pd_1x", 
                "p3_pd_bad_score.out")
# fuera de p3
sliding_windows_R(pd_seq_complete,"p3_pd_sw_out.fasta",7000,8400,ref="LC228573.1")
gkmsvm_classify("p3_pd_sw_out.fasta",
                svmfnprfx="p3_pd_1x", 
                "p3_pd_out_score.out")
# Cargamos el fichero obtenido que contiene las score de clasificaci?n
pos1<-read.delim('p3_pd_good_score.out',header = FALSE)
pos2<-read.delim('p3_pd_bad_score.out',header = FALSE)
pos3<-read.delim('p3_pd_out_score.out',header = FALSE)
# Graficamos la secuencia vs. su score obtenido usando ggplot2
library(ggplot2)
# obtenemos una gr?fica donde se obtiene un pico, esto es debido a que es la zona donde se obtiene el mayor puntaje,
# es decir la secuencia que mejor representa el motivo
library("gridExtra")
png("PV_P3.png ")
p1<-ggplot(data =pos1, aes(x = c(2400:4800) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Good")
p2<-ggplot(data =pos2, aes(x = c(2400:4800) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Bad")
p3<-ggplot(data =pos3, aes(x = c(7000:8400) , y = V2,group = 1)) + xlab("Sequence")+ylab("scores")+
  geom_line() +theme_minimal()+ggtitle("Out")
grid.arrange(p1, p2, p3, ncol = 1, nrow = 3)
dev.off()
```

Utilizamos un algoritmo de doble busqueda de patrón para para comprobar las secuencias que realmente presentan G2A6 o GA6

```{r}
posit2search<-function(name.df) {
  library(Biostrings)
  seq_pos<-c()
  for (i in 1:nrow(name.df)) {
    data <- name.df$sequence[i]
    dna <- DNAString(data)
    #busca el patron deseado dentro de los resultados obtenidos por fimo
    results<-matchPattern('GGAAAAAA', dna)
    empty<-isEmpty(results)
    if (empty=="FALSE")
    {
      positive<-data.frame(seqid=name.df$seqid[i],
                           score=name.df$score[i],
                           ID=name.df$ID[i],
                           start=name.df$start[i],
                           stop=name.df$end[i],
                           sequence=name.df$sequence[i])
      seq_pos<-rbind.data.frame(seq_pos,positive)
      
    }
    
  }
  return(seq_pos)
}
p1_sp2<-posit2search(P1_SP)
p1_pv2<-posit2search(P1_PV)
p1_fm<-posit2search(FM)
posit2search<-function(name.df) {
  library(Biostrings)
  seq_pos<-c()
  for (i in 1:nrow(name.df)) {
    data <- name.df$sequence[i]
    dna <- DNAString(data)
    #busca el patron deseado dentro de los resultados obtenidos por fimo
    results<-matchPattern('GAAAAAA', dna)
    empty<-isEmpty(results)
    if (empty=="FALSE")
    {
      positive<-data.frame(seqid=name.df$seqid[i],
                           score=name.df$score[i],
                           ID=name.df$ID[i],
                           start=name.df$start[i],
                           stop=name.df$end[i],
                           sequence=name.df$sequence[i])
      seq_pos<-rbind.data.frame(seq_pos,positive)
      
    }
    
  }
  return(seq_pos)
}
p3_sp<-posit2search(P3_SP)
p3_pv2<-posit2search(P3_PV)
p3_pd2<-posit2search(P3_PD)

```

