In currocam/FascinRSCA: R and Python code for SCA with fascin protein
knitr::opts_chunk$set(echo = TRUE)
En este cuaderno se calculan una serie de descriptores numéricos y no numéricos (los cuales serán transformados en dummy variables para poder trabajar con ello) con un dataset de secuencias relativas a fascina (anotadas como fascina 1, fascina 2 o fascina) junto a un control negativo (formado por secuencias funcionalmente similares y secuencias aleatorias obtenidas de UniProt).
Descriptores
En primer lugar, cargamos los archivos a utilizar.
library(protr)
library(tidyverse)
library(magrittr)
library(FascinRSCA)
# load FASTA files into dfs
NonRelated <- protr::readFASTA("../inst/extdata/interim/negative_control/NonRelated.fasta")
df.Nonrelated <- tibble(fasta_header =names(NonRelated),
fasta_seq= as.character(NonRelated), annotation = "NonRelated")
Related <- protr::readFASTA("../inst/extdata/interim/negative_control/related_fascin.fasta")
df.related <- tibble(fasta_header =names(Related),
fasta_seq= as.character(Related), annotation = "Related")
df <- FascinRSCA::freezing_cluster %>%
filter(annotation != "fascin2a") %>%
filter(fasta_header != "ENSG00000186765.11|ENSP00000334665.7")%>% #Quitamos fscn2b humana
separate(fasta_header, c("Entry_from_header", "EntryName_from_header"),
sep = "\\|", remove = FALSE)%>%
mutate(Entry = is.na(Entry) %>%
if_else(Entry_from_header, Entry),
`Entry name` = is.na(`Entry name`) %>%
if_else(EntryName_from_header, `Entry name`) %>%
make.unique(),
database = file_label %>%
fct_recode(en = "Ensembl", rp = "RefProteomes", db = "SwissProt"))%>%
mutate(fasta_header = glue::glue("{database}|{Entry}|{`Entry name`}"))%>%
bind_rows(df.Nonrelated, df.related)
df$annotation %<>%
fct_collapse(fascin2 = c("fascin2"),
fascin = c("fascin", "singed", "uncharacterized"),
NegativeControl = c("Related", "NonRelated"))
df$annotation %>% fct_count()
Filtramos aquellas secuencias con aminoácidos no estándares (nos darán error al usar herramientas posteriores) y creamos archivos fasta y de identificadores con los que aplicar EMBOSS.
df %<>% filter(purrr::map_lgl(.$fasta_seq, ~protcheck(.x)))
protR_desc<- df %>%
distinct(fasta_seq, .keep_all = TRUE)%>%
select(fasta_header, fasta_seq, annotation)%>%
separate(fasta_header, c("database", "Entry", "Entry.Name"), "\\|", remove = FALSE)
ann <- protR_desc %>% pull(annotation)
fct_count(ann)
fasta_seq <- protR_desc %>%
dplyr::pull(fasta_seq) %>%
Biostrings::AAStringSet()
names(fasta_seq) <- protR_desc$fasta_header
A continuación, escribimos los archivos que vamos a usar con EMBOSS:
fasta_seq[1:500] %>%
Biostrings::writeXStringSet('fascinClassifier_1_500.fasta')
fasta_seq[501:length(fasta_seq)] %>%
Biostrings::writeXStringSet('fascinClassifier_501_568.fasta')
fasta_seq %>%
Biostrings::writeXStringSet('fascinClassifier_001_568.fasta')
write(protR_desc$Entry.Name, "fascinClassifier.txt")
EMBOSS
En esta sección se calculan una serie de descriptores usando la suite EMBOSS.
Pepstats
Pepstats calculates statistics for your protein such as molecular weight, isoelectric point etc.
pepstats.csv <- readr::read_csv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_pepstats.csv",
show_col_types = FALSE)
all(pepstats.csv$names==protR_desc$Entry.Name)
emboss.pepstats <- protR_desc %>%
left_join(pepstats.csv, by = c("Entry.Name" = "names"))%>%
select(-c(fasta_seq, Entry, fasta_header, database, annotation)) %>%
rename_with(~paste0("EMBOSS_Pepstats_", .), -c("Entry.Name"))
head(emboss.pepstats)
Pepdigest
Cálculo de los sitios de corte de enzimas proteolíticas (tripsina) usando la herramienta pepdigest.
pepdigest.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_pepdigest.tsv",
show_col_types = FALSE) %>%
dplyr::filter(SeqName != 'SeqName') ## repeated first rows
emboss.pepdigest <- pepdigest.tsv %>%
tidyr::nest(data = c(!SeqName))%>%
right_join(protR_desc %>%
select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>%
tidyr::unnest(-c("SeqName", "fasta_seq")) %>%
mutate(seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>%
select(SeqName, seq)%>%
extractEMBOSS_tsv("Trypsin.Cleavage")%>%
select_if(~sum(is.na(.))< 100)%>%
dplyr::rename("Entry.Name"= "SeqName")
dim(emboss.pepdigest)
Antigenicity
A continuación, procesamos los epítopos encontrados por EMBOSS:
#Cargamos el archivo tsv
antigenic.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_antigenic.tsv",
show_col_types = FALSE) %>%
dplyr::filter(SeqName != 'SeqName') ## repeated first rows
#Calculamos estadísticos para decidir si aplicar un filtro de score
summary(as.numeric(antigenic.tsv$Score))
emboss.antigenic <- antigenic.tsv %>%
#filter(as.numeric(Score) >1)%>%
tidyr::nest(data = c(!SeqName))%>%
right_join(protR_desc %>%
select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>%
tidyr::unnest(-c("SeqName", "fasta_seq")) %>%
mutate(seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>%
select(SeqName, seq)%>%
extractEMBOSS_tsv("Epitope")%>%
select_if(~sum(is.na(.))< 100)%>%
dplyr::rename("Entry.Name"= "SeqName")
dim(emboss.antigenic)
Sigcleave
A continuación, empleamos la información del programa sigcleave que determina el punto de escisión del péptido maduro y la secuencia señal con el método von Heijne. Incluimos una serie de variables con la secuencia del péptido maduro, del péptido escindido y una columna que nos indica si tenía un punto de escisión o no.
#Cargamos el archivo tsv
sigcleave.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_sigcleave.tsv",
show_col_types = FALSE) %>%
dplyr::filter(SeqName != 'SeqName') ## repeated first rows
emboss.sigcleave <- sigcleave.tsv %>%
tidyr::nest(data = c(!SeqName))%>%
right_join(protR_desc %>%
select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>%
tidyr::unnest(-c("SeqName", "fasta_seq")) %>%
mutate(signal_seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>%
select(SeqName, mature_peptide, signal_seq)%>%
extractEMBOSS_tsv("sigcleave")%>%
select_if(~sum(is.na(.))< 100)%>%
dplyr::rename("Entry.Name"= "SeqName")
emboss.isCleave <- tibble(Entry.Name = protR_desc$Entry.Name) %>%
mutate(isCleave = ifelse(Entry.Name %in% emboss.sigcleave$Entry.Name, TRUE, FALSE))
sum(emboss.isCleave$isCleave)
dim(emboss.sigcleave)
Generating various numerical representation schemes of protein sequences
En este cuaderno se utiliza la librería protr para calcular una serie de descriptores numéricos para las secuencias homólogas a fascina. El significado físico-químico de estos descriptores se explica de forma resumida en este cuaderno pero se puede profundizar en la documentación y en este artículo publicado el año pasado donde emplean estos descriptores.
Amino acid composition
En primer lugar, nos centraremos en la composición de aminoácidos. Usaremos las ecuaciones:
extractDC() - Dipeptide compositionextractTC() - Tripeptide composition
DC.composition <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractDC(.x))%>%
rename_with(~paste0("protR_DC_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
#glimpse(DC.composition, width = 100)
TC.composition <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractTC(.x))%>%
rename_with(~paste0("protR_TC_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
#glimpse(TC.composition)
Autocorrelation Descriptors
A continuación, vamos a emplear protr cara calcular una serie de descriptores autocorrelacionados. Estos son descriptores que nos refieren a la distribución de las propiedades de los aminoácidos a lo largo de la secuencia. En concreto, usaremos los descriptores autocorrelacionados de Moreau Broto normalizados, al ser el más sencillo.
Para una propiedad en concreto, el descriptor se definiría como
$$
AC(d) = \sum^{N-d}{i=1}{P_i P{i+d}}
$$
donde $i$ es una posición en la secuencia, $d$ es un número comprendido entre 1 y nlag (por defecto 30), al que se denomina lag y $P_r$ es la propiedad del aminoácido en la posición $r$ después de haber sido estandarizado y centralizado. Este descriptor es posteriormente normalizado siguiendo la siguiente ecuación:
$$
ATS (d) = \frac{AC(d)}{N-d}
$$
Por tanto, tendremos para cada secuencia una nueva variable cuyo nombre vendrá compuesto por la propiedad físico-química empleada y el lag, es decir, el retraso (es decir, espacio) entre los aminoácidos.
Moreau <- protR_desc %$%
purrr::map_df(fasta_seq, ~ extractMoreauBroto(.)) %>%
rename_with(~paste0("protR_Moreau_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
names(Moreau) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>%
str_replace("BHAR880101", "AverageFlexibilityIndices")%>%
str_replace("CHAM820101", "PolarizabilityParameter")%>%
str_replace("CHAM820102", " FreeEnergySolutionInWater")%>%
str_replace("CHOC760101", "ResidueAccessibleSurfaceAreaTripeptide")%>%
str_replace("BIGC670101", "ResidueVolume")%>%
str_replace("CHAM810101", "StericParameter")%>%
str_replace("DAYM780201", "RelativeMutability")
#glimpse(Moreau, width = 100)
head(Moreau)
Physicochemical properties of Composition, Transition, Distribution
Los descriptores CTD fueron desarrollados por Dubchak et al. (1995) and Dubchak et al. (1999) y representan la composición de las secuencias, la transición y la distribución en base a tres clases (para reducir dimensionalidad) en la que clasifico a los aminoácidos para cada propiedad.
| | Group 1 | Group 2 | Group 3 | |
|:---:|:---:|:---:|:---:|:---:|
| Hydrophobicity | Polar | Neutral | Hydrophobicity | |
| | R, K, E, D, Q, N | G, A, S, T, P, H, Y | C, L, V, I, M, F, W | |
| Normalized van der Waals Volume | 0-2.78 | 2.95-4.0 | 4.03-8.08 | |
| | G, A, S, T, P, D, C | N, V, E, Q, I, L | M, H, K, F, R, Y, W | |
| Polarity | 4.9-6.2 | 8.0-9.2 | 10.4-13.0 | |
| | L, I, F, W, C, M, V, Y | P, A, T, G, S | H, Q, R, K, N, E, D | |
| Polarizability | 0-1.08 | 0.128-0.186 | 0.219-0.409 | |
| | G, A, S, D, T | C, P, N, V, E, Q, I, L | K, M, H, F, R, Y, W | |
| Charge | Positive | Neutral | Negative | |
| | K, R | A, N, C, Q, G, H, I, L, M, F, P, S, T, W, Y, V | D, E | |
| Secondary Structure | Helix | Strand | Coil | |
| | E, A, L, M, Q, K, R, H | V, I, Y, C, W, F, T | G, N, P, S, D | |
| Solvent Accessibility | Buried | Exposed | Intermediate | |
| | A, L, F, C, G, I, V, W | R, K, Q, E, N, D | M, S, P, T, H, Y | |
Composition
La composición es el porcentaje de aminoácidos para cada clase y para cada proteína.
CTD.composition <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractCTDC(.x))%>%
rename_with(~paste0("protR_CTDC_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
Transition
La transición de una clase $j$ a $k$ para una propiedad es la frecuencia con la que se encuentra un dipéptido tal que $a_ja_k$ o $a_ka_j$ donde $a_i$ es un aminoácido que pertenece a la clase $i$.
CTD.transition <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractCTDT(.x))%>%
rename_with(~paste0("protR_CTDT_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
names(CTD.transition) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>%
str_replace("prop1", "hydrophobicity")%>%
str_replace("prop2", "normwaalsvolume")%>%
str_replace("prop3", " polarity")%>%
str_replace("prop4", "polarizability")%>%
str_replace("prop5", "charge")%>%
str_replace("prop6", "secondarystruct")%>%
str_replace("prop7", "solventaccess")
Distribution
La distribución de los atributos a lo largo de la secuencia. Para cada atributo hay 5 descriptores que son la posición en la secuencia entre la longitud de esta (% de la posición) del primer residuo de dicha clase para ese atributo, para el 25% de los residuos, 50%, 75% y 100%.
CTD.distribution <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractCTDD(.x)) %>%
rename_with(~paste0("protR_CTDD_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
names(CTD.distribution) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>%
str_replace("prop1", "hydrophobicity")%>%
str_replace("prop2", "normwaalsvolume")%>%
str_replace("prop3", " polarity")%>%
str_replace("prop4", "polarizability")%>%
str_replace("prop5", "charge")%>%
str_replace("prop6", "secondarystruct")%>%
str_replace("prop7", "solventaccess")
Conjoint Triad Descriptors
Se clasifican los aminoácidos según los dipolos y volúmenes de los aminoácidos.
| Nº | Dipole scale (Debye) | Volume Scale (A∘) | Class |
|:---:|:---:|:---:|:---:|
| 1 | <1.0 | < 50 | Ala, Gly, Val |
| 2 | <1.0 | > 50 | Ile, Leu, Phe, Pro |
| 3 | 1.0 < Dipole < 2.0 | > 50 | Tyr, Met, Thr, Ser |
| 4 | .0 < Dipole < 3.0; | > 50 | His, Asn, Gln, Tpr |
| 5 | > 3.0 | > 50 | Arg, Lys |
| 6 | Dipole > 3.0 with opposite orientation | > 50 | Asp, Glu |
| 7 | | > 50 | Cys |
La secuencia de proteínas se representan como un vector en un espacio multidimensional. Para ello, se trata cada tríada de aminoácidos como una única unidad y cada aminoácido como equivalente al resto de aminoácidos de su clase. Se obtiene a partir de una secuencia (en un espacio bidimensional) un vector de frecuencia $F$ alojado en un espacio multidimensional $V$. Las dimensiones de $V$ son todas las posibles combinaciones de triadas, $v_i$ y el valor del vector $F$ en la i-ésima dimensión, $f_i$ es la frecuencia de un tipo tríada en la secuencia.
ConjointTriadCalculation <- protR_desc%$%
purrr::map_df(fasta_seq, ~ extractCTriad(.x)) %>%
rename_with(~paste0("protR_Triad_", .))%>%
mutate(Entry.Name =protR_desc$Entry.Name)
Quasi-sequence-order Descriptors
A continuación, calcularemos el Sequence-order-coupling number. Este parámetro se define para un tipo de distancia físico química y para un valor de lag, $d$, de la siguiente manera:
$$
\tau_d = \sum^{N-d}{i=1}{d{i, i+d}}
$$
donde $d_{i, i+d}$ es la distancia fisicoquímica entre los aminoácidos de las posiciones $i$ e $i+d$. Por tanto, lo que se calcula es una matriz de distancias fisicoquímicas (se emplean Grantham y Schneider) y luego se calcula el anterior sumatorio.
SOCN <- protR_desc%>%
filter(nchar(fasta_seq)>60)%$%
purrr::map2_df(fasta_seq, Entry.Name,
~ extractSOCN(.x) %>%t() %>%
as_tibble() %>% mutate(Entry.Name =.y)) %>%
rename_with(~paste0("protR_SOCN_", .))%>%
rename("Entry.Name" = "protR_SOCN_Entry.Name")
Pseudo-Amino Acid Composition (PseAAC)
A continuación, calculamos la composición de pseudo aminoácidos (PAAC). Esta refleja tanto la frecuencia con la que aparece cada aminoácido como el efecto del orden y tiene una fórmula más compleja.
PAAC <- protR_desc %>%
filter(nchar(fasta_seq)>50)%$%
purrr::map2_df(fasta_seq, Entry.Name,
~ extractPAAC(.x) %>%t() %>%
as_tibble() %>% mutate(Entry.Name =.y)) %>%
rename_with(~paste0("protR_PAAC_", .))%>%
rename("Entry.Name" = "protR_PAAC_Entry.Name")
Preprocesamiento de los datos
Construimos un DataFrame que contenga todas las variables numéricas de protR y pepstats y las dummy variables de EMBOSS.
protR.df <- list(TC.composition, CTD.composition,CTD.distribution,
CTD.transition, DC.composition, TC.composition, Moreau,
ConjointTriadCalculation, emboss.pepstats, emboss.isCleave) %>%
reduce(left_join, by = "Entry.Name")
dummy.df<- emboss.sigcleave %>%
full_join(emboss.antigenic, by = "Entry.Name") %>%
full_join(emboss.pepdigest, by = "Entry.Name")
dummy.df.entry <- dummy.df$Entry.Name
dummy.df <- dummy.df %>%
select(-Entry.Name) %>%
fastDummies::dummy_cols(remove_first_dummy = TRUE,
remove_selected_columns = TRUE, ignore_na = TRUE)%>%
mutate(across(everything(), .fns = ~replace_na(.,0)))
Ahora, juntamos ambos dataframes y eliminamos de memoria las variables anteriores.
df.caret <- protR.df %>%
left_join(dummy.df %>% cbind(Entry.Name = dummy.df.entry), by= "Entry.Name")
gdata::keep(df.caret, ann, sure = TRUE)
#write_csv(df.caret.numeric, "df.caret.numeric.csv")
Caret
A continuación, vamos a preprocesar los datos utilizados para el predictor. Comenzamos cambiando de posición aquellas variables que, por bibliografía y el clustering, sabemos que tienen un significado biológico. De otro modo, la variable STYNQW, por ejemplo, sería eliminada al aplicar el combo Lineal, puesto que es combinación lineal de los aminoácidos que la componen.
spec.col <- c("EMBOSS_Pepstats_basic","EMBOSS_Pepstats_residue_n",
"EMBOSS_Pepstats_STYNQW","EMBOSS_Pepstats_charge",
"EMBOSS_Pepstats_nonPolar")
df.caret <- df.caret %>%
relocate(spec.col, .before = 1) %>%
select(where(is.numeric))
dim(df.caret)
En primer lugar, comenzaremos eliminando aquellas variables de valor único (varianza cero). Además, y para evitar problemas de capacidad de cómputo en mi portatil, también eliminaremos aquellas variables cuya varianza sea cercana al cero. De esta forma, eliminamos 25.000 variables.
library(caret)
nzv <- df.caret %>%
nearZeroVar(saveMetrics= TRUE)%>%
rownames_to_column(var = "var")
#nzv %>% filter(zeroVar | nzv)
df.caret %<>%
select(-all_of(nzv %>%
filter(zeroVar| nzv) %>% pull(var)))
dim(df.caret)
Correlación
A continuación, vamos a eliminar aquellas variables que se encuentren correlacionadas en más de un .99. Manualmente, evitaremos eliminar las variabels que habíamos determinado como importantes previamente (setdiff(highlyCorrelated, spec.col)).
`%!in%` = Negate(`%in%`)
highlyCorrelated <- findCorrelation(cor(df.caret, use = "pairwise.complete.obs", ), cutoff=0.99, names = TRUE)
df.caret <- df.caret %>%
select(-all_of(setdiff(highlyCorrelated, spec.col)))
dim(df.caret)
Combo lineal
A continuación, determinamos qué combinaciones lineales existen en nuestros datos y eliminamos una serie de variables para evitarlas. Es necesario, para este paso, eliminar una serie de secuencias que tienen NAs. Corresponden a secuencias del control negativo cuya longitud es inferior al lag utilizado en los descriptores de la librería protR.
comboInfo <-df.caret %>%
drop_na %>%
findLinearCombos()
df.caret<- df.caret[, -comboInfo$remove]
dim(df.caret)
glimpse(df.caret)
Por último, vamos a construir los datasets necesarios para entrenar los siguientes clasificadores.
df.caret <- df.caret %>% mutate(ann = ann)
df.caret %>% pull(ann) %>% fct_count()
set.seed(96)
library(rsample)
Clasificador fascina 1 Vs FSCN2 y fascina (sin control negativo)
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina 1 o, si por el contrario, se trata de una secuencia de fascina o fascina 2. No se incluyen controles negativos. Destacar que existen una serie de variables relacionadas con la composición de tripéptidos de las secuencias que tienen varianza cero si no se tienen en cuenta las secuencias de control negativo. Es decir, que ninguna de las secuencias de fascina tiene en su composición los tripéptidos KEA, KKA, QLR, KLR, VQE o AKE. Por ello, además de centrar y escalar, vamos a eliminar también dichas columnas.
FSCN1_Vs_others.df <- df.caret %>%
filter(ann %in% c("fascin", "fascin2")) %>%
group_by(ann) %>%
slice_sample(n = 41) %>%
ungroup %>%
bind_rows(df.caret %>% filter(ann %in% "fascin1"))%>%
mutate(ann = fct_collapse(ann,
not_fascina1 = c("fascin", "fascin2")))
FSCN1_Vs_others.nzv <- FSCN1_Vs_others.df %>%
nearZeroVar(saveMetrics= TRUE)%>%
filter(zeroVar) %>%
rownames_to_column(var = "var")
head(FSCN1_Vs_others.nzv)
FSCN1_Vs_others.df %<>%
select(-all_of(FSCN1_Vs_others.nzv %>% pull(var)))
FSCN1_Vs_others.df %>% pull(ann) %>% fct_count()
FSCN1_Vs_others.preProcValues <- preProcess(FSCN1_Vs_others.df, method = c("center", "scale"))
FSCN1_Vs_others.Transformed <- predict(FSCN1_Vs_others.preProcValues,
FSCN1_Vs_others.df)
write_csv(FSCN1_Vs_others.Transformed, "FSCN1_Vs_others.csv")
Clasificador fascina 1 Vs FSCN2 y fascina (con control negativo)
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina 1 o, si por el contrario, se trata de una secuencia de fascina o fascina 2. Sí se incluyen controles negativos.
FSCN1_Vs_others_and_neg.df <- df.caret %>%
filter(ann %in% c("fascin", "fascin2")) %>%
group_by(ann) %>%
slice_sample(n = 20) %>%
ungroup %>%
bind_rows(df.caret %>%
filter(ann %in% "NegativeControl")%>%
slice_sample(n = 40))%>%
bind_rows(df.caret %>% filter(ann %in% "fascin1"))%>%
mutate(ann = fct_collapse(ann,
not_fascina1 = c("fascin", "fascin2", "NegativeControl")))
FSCN1_Vs_others_and_neg.df %>% pull(ann) %>% fct_count()
FSCN1_Vs_others_and_neg.preProcValues <- preProcess(FSCN1_Vs_others_and_neg.df,
method = c("center", "scale"))
FSCN1_Vs_others_and_neg <- predict(FSCN1_Vs_others_and_neg.preProcValues,
FSCN1_Vs_others_and_neg.df)
write_csv(FSCN1_Vs_others_and_neg, "FSCN1_Vs_others_and_neg.csv")
Clasificador fascin-related Vs Control negativo
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina (en sentido amplio, incluyendo fascina, FSCN1 y FSCN2) o si por el contrario, pertenece al control negativo.
FSCN_Vs_NonFSCN.df <- df.caret %>%
filter(ann %in% "NegativeControl") %>%
group_by(ann) %>%
slice_sample(n = 271) %>%
ungroup %>%
bind_rows(df.caret %>% filter(!ann %in% "NegativeControl"))%>%
mutate(ann = fct_collapse(ann,
fascin_rel = c("fascin", "fascin2", "fascin1")))
FSCN_Vs_NonFSCN.df %>% pull(ann) %>% fct_count()
FSCN_Vs_NonFSCN.preProcValues <- preProcess(FSCN_Vs_NonFSCN.df,
method = c("center", "scale", "nzv"))
FSCN_Vs_NonFSCN.Transformed <- predict(FSCN_Vs_NonFSCN.preProcValues,
FSCN_Vs_NonFSCN.df, )
write_csv(FSCN_Vs_NonFSCN.Transformed, "FSCN_rel_Vs_NonFSCN.csv")
Multiclasificador
A continuación, se prepararan los datos para entrenar y testar un multiclasificador que prediga si una secuencia es una secuencia de fascina, fascina1, fascina2 y control negativo.
multicls.df <- df.caret %>%
filter(ann %in% "NegativeControl") %>%
slice_sample(n = 100) %>%
bind_rows(df.caret %>%
filter(!ann %in% "NegativeControl"))%>%
mutate(ann = as.factor(as.character(ann)))
multicls.df %>% pull(ann) %>% fct_count()
multicls.preProcValues <- preProcess(multicls.df, method = c("center", "scale"))
multicls.Transformed <- predict(multicls.preProcValues,
multicls.df)
write_csv(multicls.Transformed, "multiclasificador.csv")
currocam/FascinRSCA documentation built on March 21, 2022, 6:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
En este cuaderno se calculan una serie de descriptores numéricos y no numéricos (los cuales serán transformados en dummy variables para poder trabajar con ello) con un dataset de secuencias relativas a fascina (anotadas como fascina 1, fascina 2 o fascina) junto a un control negativo (formado por secuencias funcionalmente similares y secuencias aleatorias obtenidas de UniProt).
Descriptores
En primer lugar, cargamos los archivos a utilizar.
library(protr) library(tidyverse) library(magrittr) library(FascinRSCA) # load FASTA files into dfs NonRelated <- protr::readFASTA("../inst/extdata/interim/negative_control/NonRelated.fasta") df.Nonrelated <- tibble(fasta_header =names(NonRelated), fasta_seq= as.character(NonRelated), annotation = "NonRelated") Related <- protr::readFASTA("../inst/extdata/interim/negative_control/related_fascin.fasta") df.related <- tibble(fasta_header =names(Related), fasta_seq= as.character(Related), annotation = "Related") df <- FascinRSCA::freezing_cluster %>% filter(annotation != "fascin2a") %>% filter(fasta_header != "ENSG00000186765.11|ENSP00000334665.7")%>% #Quitamos fscn2b humana separate(fasta_header, c("Entry_from_header", "EntryName_from_header"), sep = "\\|", remove = FALSE)%>% mutate(Entry = is.na(Entry) %>% if_else(Entry_from_header, Entry), `Entry name` = is.na(`Entry name`) %>% if_else(EntryName_from_header, `Entry name`) %>% make.unique(), database = file_label %>% fct_recode(en = "Ensembl", rp = "RefProteomes", db = "SwissProt"))%>% mutate(fasta_header = glue::glue("{database}|{Entry}|{`Entry name`}"))%>% bind_rows(df.Nonrelated, df.related) df$annotation %<>% fct_collapse(fascin2 = c("fascin2"), fascin = c("fascin", "singed", "uncharacterized"), NegativeControl = c("Related", "NonRelated")) df$annotation %>% fct_count()
Filtramos aquellas secuencias con aminoácidos no estándares (nos darán error al usar herramientas posteriores) y creamos archivos fasta y de identificadores con los que aplicar EMBOSS.
df %<>% filter(purrr::map_lgl(.$fasta_seq, ~protcheck(.x))) protR_desc<- df %>% distinct(fasta_seq, .keep_all = TRUE)%>% select(fasta_header, fasta_seq, annotation)%>% separate(fasta_header, c("database", "Entry", "Entry.Name"), "\\|", remove = FALSE) ann <- protR_desc %>% pull(annotation) fct_count(ann) fasta_seq <- protR_desc %>% dplyr::pull(fasta_seq) %>% Biostrings::AAStringSet() names(fasta_seq) <- protR_desc$fasta_header
A continuación, escribimos los archivos que vamos a usar con EMBOSS:
fasta_seq[1:500] %>% Biostrings::writeXStringSet('fascinClassifier_1_500.fasta') fasta_seq[501:length(fasta_seq)] %>% Biostrings::writeXStringSet('fascinClassifier_501_568.fasta') fasta_seq %>% Biostrings::writeXStringSet('fascinClassifier_001_568.fasta') write(protR_desc$Entry.Name, "fascinClassifier.txt")
EMBOSS
En esta sección se calculan una serie de descriptores usando la suite EMBOSS.
Pepstats
Pepstats calculates statistics for your protein such as molecular weight, isoelectric point etc.
pepstats.csv <- readr::read_csv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_pepstats.csv", show_col_types = FALSE) all(pepstats.csv$names==protR_desc$Entry.Name) emboss.pepstats <- protR_desc %>% left_join(pepstats.csv, by = c("Entry.Name" = "names"))%>% select(-c(fasta_seq, Entry, fasta_header, database, annotation)) %>% rename_with(~paste0("EMBOSS_Pepstats_", .), -c("Entry.Name")) head(emboss.pepstats)
Pepdigest
Cálculo de los sitios de corte de enzimas proteolíticas (tripsina) usando la herramienta pepdigest.
pepdigest.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_pepdigest.tsv", show_col_types = FALSE) %>% dplyr::filter(SeqName != 'SeqName') ## repeated first rows emboss.pepdigest <- pepdigest.tsv %>% tidyr::nest(data = c(!SeqName))%>% right_join(protR_desc %>% select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>% tidyr::unnest(-c("SeqName", "fasta_seq")) %>% mutate(seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>% select(SeqName, seq)%>% extractEMBOSS_tsv("Trypsin.Cleavage")%>% select_if(~sum(is.na(.))< 100)%>% dplyr::rename("Entry.Name"= "SeqName") dim(emboss.pepdigest)
Antigenicity
A continuación, procesamos los epítopos encontrados por EMBOSS:
#Cargamos el archivo tsv antigenic.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_antigenic.tsv", show_col_types = FALSE) %>% dplyr::filter(SeqName != 'SeqName') ## repeated first rows #Calculamos estadísticos para decidir si aplicar un filtro de score summary(as.numeric(antigenic.tsv$Score)) emboss.antigenic <- antigenic.tsv %>% #filter(as.numeric(Score) >1)%>% tidyr::nest(data = c(!SeqName))%>% right_join(protR_desc %>% select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>% tidyr::unnest(-c("SeqName", "fasta_seq")) %>% mutate(seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>% select(SeqName, seq)%>% extractEMBOSS_tsv("Epitope")%>% select_if(~sum(is.na(.))< 100)%>% dplyr::rename("Entry.Name"= "SeqName") dim(emboss.antigenic)
Sigcleave
A continuación, empleamos la información del programa sigcleave que determina el punto de escisión del péptido maduro y la secuencia señal con el método von Heijne. Incluimos una serie de variables con la secuencia del péptido maduro, del péptido escindido y una columna que nos indica si tenía un punto de escisión o no.
#Cargamos el archivo tsv sigcleave.tsv <- readr::read_tsv("../inst/extdata/interim/fascin_ann_classifier_desc/fascinClassifier_001_568_sigcleave.tsv", show_col_types = FALSE) %>% dplyr::filter(SeqName != 'SeqName') ## repeated first rows emboss.sigcleave <- sigcleave.tsv %>% tidyr::nest(data = c(!SeqName))%>% right_join(protR_desc %>% select(Entry.Name, fasta_seq), by = c("SeqName" = "Entry.Name"))%>% tidyr::unnest(-c("SeqName", "fasta_seq")) %>% mutate(signal_seq = extractEMBOSS_StartEnd(Start, End, fasta_seq))%>% select(SeqName, mature_peptide, signal_seq)%>% extractEMBOSS_tsv("sigcleave")%>% select_if(~sum(is.na(.))< 100)%>% dplyr::rename("Entry.Name"= "SeqName") emboss.isCleave <- tibble(Entry.Name = protR_desc$Entry.Name) %>% mutate(isCleave = ifelse(Entry.Name %in% emboss.sigcleave$Entry.Name, TRUE, FALSE)) sum(emboss.isCleave$isCleave) dim(emboss.sigcleave)
Generating various numerical representation schemes of protein sequences
En este cuaderno se utiliza la librería protr para calcular una serie de descriptores numéricos para las secuencias homólogas a fascina. El significado físico-químico de estos descriptores se explica de forma resumida en este cuaderno pero se puede profundizar en la documentación y en este artículo publicado el año pasado donde emplean estos descriptores.
Amino acid composition
En primer lugar, nos centraremos en la composición de aminoácidos. Usaremos las ecuaciones:
extractDC()- Dipeptide compositionextractTC()- Tripeptide composition
DC.composition <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractDC(.x))%>% rename_with(~paste0("protR_DC_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name) #glimpse(DC.composition, width = 100)
TC.composition <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractTC(.x))%>% rename_with(~paste0("protR_TC_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name) #glimpse(TC.composition)
Autocorrelation Descriptors
A continuación, vamos a emplear protr cara calcular una serie de descriptores autocorrelacionados. Estos son descriptores que nos refieren a la distribución de las propiedades de los aminoácidos a lo largo de la secuencia. En concreto, usaremos los descriptores autocorrelacionados de Moreau Broto normalizados, al ser el más sencillo.
Para una propiedad en concreto, el descriptor se definiría como $$ AC(d) = \sum^{N-d}{i=1}{P_i P{i+d}} $$ donde $i$ es una posición en la secuencia, $d$ es un número comprendido entre 1 y nlag (por defecto 30), al que se denomina lag y $P_r$ es la propiedad del aminoácido en la posición $r$ después de haber sido estandarizado y centralizado. Este descriptor es posteriormente normalizado siguiendo la siguiente ecuación:
$$ ATS (d) = \frac{AC(d)}{N-d} $$ Por tanto, tendremos para cada secuencia una nueva variable cuyo nombre vendrá compuesto por la propiedad físico-química empleada y el lag, es decir, el retraso (es decir, espacio) entre los aminoácidos.
Moreau <- protR_desc %$% purrr::map_df(fasta_seq, ~ extractMoreauBroto(.)) %>% rename_with(~paste0("protR_Moreau_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name) names(Moreau) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>% str_replace("BHAR880101", "AverageFlexibilityIndices")%>% str_replace("CHAM820101", "PolarizabilityParameter")%>% str_replace("CHAM820102", " FreeEnergySolutionInWater")%>% str_replace("CHOC760101", "ResidueAccessibleSurfaceAreaTripeptide")%>% str_replace("BIGC670101", "ResidueVolume")%>% str_replace("CHAM810101", "StericParameter")%>% str_replace("DAYM780201", "RelativeMutability") #glimpse(Moreau, width = 100) head(Moreau)
Physicochemical properties of Composition, Transition, Distribution
Los descriptores CTD fueron desarrollados por Dubchak et al. (1995) and Dubchak et al. (1999) y representan la composición de las secuencias, la transición y la distribución en base a tres clases (para reducir dimensionalidad) en la que clasifico a los aminoácidos para cada propiedad.
| | Group 1 | Group 2 | Group 3 | | |:---:|:---:|:---:|:---:|:---:| | Hydrophobicity | Polar | Neutral | Hydrophobicity | | | | R, K, E, D, Q, N | G, A, S, T, P, H, Y | C, L, V, I, M, F, W | | | Normalized van der Waals Volume | 0-2.78 | 2.95-4.0 | 4.03-8.08 | | | | G, A, S, T, P, D, C | N, V, E, Q, I, L | M, H, K, F, R, Y, W | | | Polarity | 4.9-6.2 | 8.0-9.2 | 10.4-13.0 | | | | L, I, F, W, C, M, V, Y | P, A, T, G, S | H, Q, R, K, N, E, D | | | Polarizability | 0-1.08 | 0.128-0.186 | 0.219-0.409 | | | | G, A, S, D, T | C, P, N, V, E, Q, I, L | K, M, H, F, R, Y, W | | | Charge | Positive | Neutral | Negative | | | | K, R | A, N, C, Q, G, H, I, L, M, F, P, S, T, W, Y, V | D, E | | | Secondary Structure | Helix | Strand | Coil | | | | E, A, L, M, Q, K, R, H | V, I, Y, C, W, F, T | G, N, P, S, D | | | Solvent Accessibility | Buried | Exposed | Intermediate | | | | A, L, F, C, G, I, V, W | R, K, Q, E, N, D | M, S, P, T, H, Y | |
Composition
La composición es el porcentaje de aminoácidos para cada clase y para cada proteína.
CTD.composition <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractCTDC(.x))%>% rename_with(~paste0("protR_CTDC_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name)
Transition
La transición de una clase $j$ a $k$ para una propiedad es la frecuencia con la que se encuentra un dipéptido tal que $a_ja_k$ o $a_ka_j$ donde $a_i$ es un aminoácido que pertenece a la clase $i$.
CTD.transition <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractCTDT(.x))%>% rename_with(~paste0("protR_CTDT_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name) names(CTD.transition) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>% str_replace("prop1", "hydrophobicity")%>% str_replace("prop2", "normwaalsvolume")%>% str_replace("prop3", " polarity")%>% str_replace("prop4", "polarizability")%>% str_replace("prop5", "charge")%>% str_replace("prop6", "secondarystruct")%>% str_replace("prop7", "solventaccess")
Distribution
La distribución de los atributos a lo largo de la secuencia. Para cada atributo hay 5 descriptores que son la posición en la secuencia entre la longitud de esta (% de la posición) del primer residuo de dicha clase para ese atributo, para el 25% de los residuos, 50%, 75% y 100%.
CTD.distribution <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractCTDD(.x)) %>% rename_with(~paste0("protR_CTDD_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name) names(CTD.distribution) %<>% str_replace("CIDH920105", "NormalizedAverageHydrophobicityScales")%>% str_replace("prop1", "hydrophobicity")%>% str_replace("prop2", "normwaalsvolume")%>% str_replace("prop3", " polarity")%>% str_replace("prop4", "polarizability")%>% str_replace("prop5", "charge")%>% str_replace("prop6", "secondarystruct")%>% str_replace("prop7", "solventaccess")
Conjoint Triad Descriptors
Se clasifican los aminoácidos según los dipolos y volúmenes de los aminoácidos.
| Nº | Dipole scale (Debye) | Volume Scale (A∘) | Class | |:---:|:---:|:---:|:---:| | 1 | <1.0 | < 50 | Ala, Gly, Val | | 2 | <1.0 | > 50 | Ile, Leu, Phe, Pro | | 3 | 1.0 < Dipole < 2.0 | > 50 | Tyr, Met, Thr, Ser | | 4 | .0 < Dipole < 3.0; | > 50 | His, Asn, Gln, Tpr | | 5 | > 3.0 | > 50 | Arg, Lys | | 6 | Dipole > 3.0 with opposite orientation | > 50 | Asp, Glu | | 7 | | > 50 | Cys |
La secuencia de proteínas se representan como un vector en un espacio multidimensional. Para ello, se trata cada tríada de aminoácidos como una única unidad y cada aminoácido como equivalente al resto de aminoácidos de su clase. Se obtiene a partir de una secuencia (en un espacio bidimensional) un vector de frecuencia $F$ alojado en un espacio multidimensional $V$. Las dimensiones de $V$ son todas las posibles combinaciones de triadas, $v_i$ y el valor del vector $F$ en la i-ésima dimensión, $f_i$ es la frecuencia de un tipo tríada en la secuencia.
ConjointTriadCalculation <- protR_desc%$% purrr::map_df(fasta_seq, ~ extractCTriad(.x)) %>% rename_with(~paste0("protR_Triad_", .))%>% mutate(Entry.Name =protR_desc$Entry.Name)
Quasi-sequence-order Descriptors
A continuación, calcularemos el Sequence-order-coupling number. Este parámetro se define para un tipo de distancia físico química y para un valor de lag, $d$, de la siguiente manera:
$$ \tau_d = \sum^{N-d}{i=1}{d{i, i+d}} $$ donde $d_{i, i+d}$ es la distancia fisicoquímica entre los aminoácidos de las posiciones $i$ e $i+d$. Por tanto, lo que se calcula es una matriz de distancias fisicoquímicas (se emplean Grantham y Schneider) y luego se calcula el anterior sumatorio.
SOCN <- protR_desc%>% filter(nchar(fasta_seq)>60)%$% purrr::map2_df(fasta_seq, Entry.Name, ~ extractSOCN(.x) %>%t() %>% as_tibble() %>% mutate(Entry.Name =.y)) %>% rename_with(~paste0("protR_SOCN_", .))%>% rename("Entry.Name" = "protR_SOCN_Entry.Name")
Pseudo-Amino Acid Composition (PseAAC)
A continuación, calculamos la composición de pseudo aminoácidos (PAAC). Esta refleja tanto la frecuencia con la que aparece cada aminoácido como el efecto del orden y tiene una fórmula más compleja.
PAAC <- protR_desc %>% filter(nchar(fasta_seq)>50)%$% purrr::map2_df(fasta_seq, Entry.Name, ~ extractPAAC(.x) %>%t() %>% as_tibble() %>% mutate(Entry.Name =.y)) %>% rename_with(~paste0("protR_PAAC_", .))%>% rename("Entry.Name" = "protR_PAAC_Entry.Name")
Preprocesamiento de los datos
Construimos un DataFrame que contenga todas las variables numéricas de protR y pepstats y las dummy variables de EMBOSS.
protR.df <- list(TC.composition, CTD.composition,CTD.distribution, CTD.transition, DC.composition, TC.composition, Moreau, ConjointTriadCalculation, emboss.pepstats, emboss.isCleave) %>% reduce(left_join, by = "Entry.Name") dummy.df<- emboss.sigcleave %>% full_join(emboss.antigenic, by = "Entry.Name") %>% full_join(emboss.pepdigest, by = "Entry.Name") dummy.df.entry <- dummy.df$Entry.Name dummy.df <- dummy.df %>% select(-Entry.Name) %>% fastDummies::dummy_cols(remove_first_dummy = TRUE, remove_selected_columns = TRUE, ignore_na = TRUE)%>% mutate(across(everything(), .fns = ~replace_na(.,0)))
Ahora, juntamos ambos dataframes y eliminamos de memoria las variables anteriores.
df.caret <- protR.df %>% left_join(dummy.df %>% cbind(Entry.Name = dummy.df.entry), by= "Entry.Name") gdata::keep(df.caret, ann, sure = TRUE) #write_csv(df.caret.numeric, "df.caret.numeric.csv")
Caret
A continuación, vamos a preprocesar los datos utilizados para el predictor. Comenzamos cambiando de posición aquellas variables que, por bibliografía y el clustering, sabemos que tienen un significado biológico. De otro modo, la variable STYNQW, por ejemplo, sería eliminada al aplicar el combo Lineal, puesto que es combinación lineal de los aminoácidos que la componen.
spec.col <- c("EMBOSS_Pepstats_basic","EMBOSS_Pepstats_residue_n", "EMBOSS_Pepstats_STYNQW","EMBOSS_Pepstats_charge", "EMBOSS_Pepstats_nonPolar") df.caret <- df.caret %>% relocate(spec.col, .before = 1) %>% select(where(is.numeric)) dim(df.caret)
En primer lugar, comenzaremos eliminando aquellas variables de valor único (varianza cero). Además, y para evitar problemas de capacidad de cómputo en mi portatil, también eliminaremos aquellas variables cuya varianza sea cercana al cero. De esta forma, eliminamos 25.000 variables.
library(caret) nzv <- df.caret %>% nearZeroVar(saveMetrics= TRUE)%>% rownames_to_column(var = "var") #nzv %>% filter(zeroVar | nzv) df.caret %<>% select(-all_of(nzv %>% filter(zeroVar| nzv) %>% pull(var))) dim(df.caret)
Correlación
A continuación, vamos a eliminar aquellas variables que se encuentren correlacionadas en más de un .99. Manualmente, evitaremos eliminar las variabels que habíamos determinado como importantes previamente (setdiff(highlyCorrelated, spec.col)).
`%!in%` = Negate(`%in%`) highlyCorrelated <- findCorrelation(cor(df.caret, use = "pairwise.complete.obs", ), cutoff=0.99, names = TRUE) df.caret <- df.caret %>% select(-all_of(setdiff(highlyCorrelated, spec.col))) dim(df.caret)
Combo lineal
A continuación, determinamos qué combinaciones lineales existen en nuestros datos y eliminamos una serie de variables para evitarlas. Es necesario, para este paso, eliminar una serie de secuencias que tienen NAs. Corresponden a secuencias del control negativo cuya longitud es inferior al lag utilizado en los descriptores de la librería protR.
comboInfo <-df.caret %>% drop_na %>% findLinearCombos() df.caret<- df.caret[, -comboInfo$remove] dim(df.caret) glimpse(df.caret)
Por último, vamos a construir los datasets necesarios para entrenar los siguientes clasificadores.
df.caret <- df.caret %>% mutate(ann = ann) df.caret %>% pull(ann) %>% fct_count() set.seed(96) library(rsample)
Clasificador fascina 1 Vs FSCN2 y fascina (sin control negativo)
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina 1 o, si por el contrario, se trata de una secuencia de fascina o fascina 2. No se incluyen controles negativos. Destacar que existen una serie de variables relacionadas con la composición de tripéptidos de las secuencias que tienen varianza cero si no se tienen en cuenta las secuencias de control negativo. Es decir, que ninguna de las secuencias de fascina tiene en su composición los tripéptidos KEA, KKA, QLR, KLR, VQE o AKE. Por ello, además de centrar y escalar, vamos a eliminar también dichas columnas.
FSCN1_Vs_others.df <- df.caret %>% filter(ann %in% c("fascin", "fascin2")) %>% group_by(ann) %>% slice_sample(n = 41) %>% ungroup %>% bind_rows(df.caret %>% filter(ann %in% "fascin1"))%>% mutate(ann = fct_collapse(ann, not_fascina1 = c("fascin", "fascin2"))) FSCN1_Vs_others.nzv <- FSCN1_Vs_others.df %>% nearZeroVar(saveMetrics= TRUE)%>% filter(zeroVar) %>% rownames_to_column(var = "var") head(FSCN1_Vs_others.nzv) FSCN1_Vs_others.df %<>% select(-all_of(FSCN1_Vs_others.nzv %>% pull(var))) FSCN1_Vs_others.df %>% pull(ann) %>% fct_count() FSCN1_Vs_others.preProcValues <- preProcess(FSCN1_Vs_others.df, method = c("center", "scale")) FSCN1_Vs_others.Transformed <- predict(FSCN1_Vs_others.preProcValues, FSCN1_Vs_others.df) write_csv(FSCN1_Vs_others.Transformed, "FSCN1_Vs_others.csv")
Clasificador fascina 1 Vs FSCN2 y fascina (con control negativo)
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina 1 o, si por el contrario, se trata de una secuencia de fascina o fascina 2. Sí se incluyen controles negativos.
FSCN1_Vs_others_and_neg.df <- df.caret %>% filter(ann %in% c("fascin", "fascin2")) %>% group_by(ann) %>% slice_sample(n = 20) %>% ungroup %>% bind_rows(df.caret %>% filter(ann %in% "NegativeControl")%>% slice_sample(n = 40))%>% bind_rows(df.caret %>% filter(ann %in% "fascin1"))%>% mutate(ann = fct_collapse(ann, not_fascina1 = c("fascin", "fascin2", "NegativeControl"))) FSCN1_Vs_others_and_neg.df %>% pull(ann) %>% fct_count() FSCN1_Vs_others_and_neg.preProcValues <- preProcess(FSCN1_Vs_others_and_neg.df, method = c("center", "scale")) FSCN1_Vs_others_and_neg <- predict(FSCN1_Vs_others_and_neg.preProcValues, FSCN1_Vs_others_and_neg.df) write_csv(FSCN1_Vs_others_and_neg, "FSCN1_Vs_others_and_neg.csv")
Clasificador fascin-related Vs Control negativo
A continuación, se prepararan los datos para entrenar y testar un clasificador que prediga si una secuencia es una secuencia de fascina (en sentido amplio, incluyendo fascina, FSCN1 y FSCN2) o si por el contrario, pertenece al control negativo.
FSCN_Vs_NonFSCN.df <- df.caret %>% filter(ann %in% "NegativeControl") %>% group_by(ann) %>% slice_sample(n = 271) %>% ungroup %>% bind_rows(df.caret %>% filter(!ann %in% "NegativeControl"))%>% mutate(ann = fct_collapse(ann, fascin_rel = c("fascin", "fascin2", "fascin1"))) FSCN_Vs_NonFSCN.df %>% pull(ann) %>% fct_count() FSCN_Vs_NonFSCN.preProcValues <- preProcess(FSCN_Vs_NonFSCN.df, method = c("center", "scale", "nzv")) FSCN_Vs_NonFSCN.Transformed <- predict(FSCN_Vs_NonFSCN.preProcValues, FSCN_Vs_NonFSCN.df, ) write_csv(FSCN_Vs_NonFSCN.Transformed, "FSCN_rel_Vs_NonFSCN.csv")
Multiclasificador
A continuación, se prepararan los datos para entrenar y testar un multiclasificador que prediga si una secuencia es una secuencia de fascina, fascina1, fascina2 y control negativo.
multicls.df <- df.caret %>% filter(ann %in% "NegativeControl") %>% slice_sample(n = 100) %>% bind_rows(df.caret %>% filter(!ann %in% "NegativeControl"))%>% mutate(ann = as.factor(as.character(ann))) multicls.df %>% pull(ann) %>% fct_count() multicls.preProcValues <- preProcess(multicls.df, method = c("center", "scale")) multicls.Transformed <- predict(multicls.preProcValues, multicls.df) write_csv(multicls.Transformed, "multiclasificador.csv")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.