[toc]
Gerar novas análises baseados no estudo dos tomateiros, tentando associar os parâmetros baseados em 16S/ITS/WGS em resposta aos fatores ambientais, e/ou tratamento.
## Evocando bibliotecas
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.4-6library(sciplot)
library(labdsv)## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.8-23. For overview type 'help("mgcv-package")'.## This is labdsv 2.0-1
## convert existing ordinations with as.dsvord()##
## Attaching package: 'labdsv'## The following object is masked from 'package:stats':
##
## density## Diretório
setwd("/work/rcsilva/projects/saccharome-wgs/")
## Leitura das matrizes de entrada
#### Taxonomias
table_kraken <- read.delim(file = "./analises/2019-04-29_kraken_all_official/tabletax.noplants.txt")
colnames(table_kraken)[1] <- "OTU_ID"
table_16S <- read.delim(file = "../saccharome-16s-rizo/analises/outputs/P01_A01_zotus_work_in_R/output/L1_OTUs.txt",
stringsAsFactors = F)
table_ITS <- read.delim(file = "../saccharome-its-rizo/analises/E14_work_in_R/P01_A01_work_in_R/output/L1_OTUs.txt",
stringsAsFactors = F)
#### Parâmetros de solo (dgpinheiro)
load("/data/rcsilva/projects/saccharome-wgs/dados/RData/soil.dgpinheiro.RData")
#### Solo: somente os dados da amostragem
soil.df.rhizo <- soil.df[,-c(3,4,5,9,10,11)]
#### Nomes de colunas idênticos às amostras
colnames(soil.df.rhizo)[3:8] <- c(
"CRZ2", "CRZ3", "CRZ4",
"ORZ2", "ORZ3", "ORZ4"
)
#### Variáveis de amostras
samps <- c("CRZ2", "CRZ3", "CRZ4",
"ORZ2", "ORZ3", "ORZ4")
#### Orgânico e convencional
samps.c <- c("CRZ2", "CRZ3", "CRZ4")
samps.o <- c("ORZ2", "ORZ3", "ORZ4")
#### Backup da matriz de rizo original
sr_bkp <- soil.df.rhizoFunções pré-feitas pelo professor para aplicação em batelada.
# Funções de teste de hipótese (dgpinheiro)
## dgpinheiro:
# x - parâmetro obrigatória para a função que será passada para
# o método apply receberá a linha neste caso
# os demais parâmetros servem para indicar os nomes das amostras
# (orgânicas [a] e convencionais [b]), respectivamente, neste caso
## Shapiro test
## o p-valor do teste de Shapiro
s.test.pvalue <- function (x,a) {
if (sd(x[a])==0) {
x[a] <- sapply(x[a], function(w) {
return( w + runif(1,min=0,max=1e-10) ) } )
}
return (shapiro.test(as.numeric(x[a]))$p.value)
}
## Teste t
## Retorna o p-valor do teste t.
t.test.pvalue <- function (x,a,b) {
var.equal=TRUE
if (var.test(as.numeric(x[a]),
as.numeric(x[b]),
alternative="two.sided")$p.value < 0.05) {
var.equal=FALSE
}
return( t.test(as.numeric(x[a]),as.numeric(x[b]),var.equal=var.equal)$p.value );
}
## Teste de wilcox
## Retorna o p-valor do teste de Wilcoxon
w.test.pvalue <- function (x,a,b) {
return( wilcox.test(as.numeric(x[a]),as.numeric(x[b]))$p.value );
}
## Teste de Kolmogorov-Smirnov
## Testa e retorna o p-valor do teste KS
k.test.pvalue <- function (x,a) {
return ( ks.test(as.numeric(x[a]), "pnorm", mean(as.numeric(x[a])), sd(as.numeric(x[a])))$p.value )
}
## Seleção de teste
## A partir do modelo, seleciona um dos testes acima
test.selection <- function(x,p = 0.05) {
if ((x[1] >= p) && (x[2] >= p) && (x[3] >= p) && x[4] >= p) {
return('t.test')
} else {
return('w.test')
}
}
# SEM: standard error
sem <- function(x) {
return( sd(x)/sqrt(length(x)) )
}Com as matrizes prontas, a meta agora é:
Primeiro, vamos testar a normalidade de cada um dos descritores.
Vamos adicionar alguns valores que serão úteis depois no cálculo em batelada das informações interessantes ao estudo.
## Adição das médias
soil.df.rhizo$mean.o <- apply(soil.df.rhizo[, samps.o], 1, mean)
soil.df.rhizo$mean.c <- apply(soil.df.rhizo[, samps.c], 1, mean)
## Desvio padrão
soil.df.rhizo$sd.o <- apply(soil.df.rhizo[,samps.o], 1, sd)
soil.df.rhizo$sd.c <- apply(soil.df.rhizo[,samps.c], 1, sd)
## Erro padrão
soil.df.rhizo$se.o <- apply(soil.df.rhizo[,samps.o],1,sem)
soil.df.rhizo$se.c <- apply(soil.df.rhizo[,samps.c],1,sem)
## TBC: Intervalos de confiança de 95%
soil.df.rhizo$ci.low.o <- apply(soil.df.rhizo[,c("mean.o","se.o")],1, function(x) {return(x[[1]]-qnorm(0.975)*x[[2]])} )
soil.df.rhizo$ci.high.o <- apply(soil.df.rhizo[,c("mean.o","se.o")],1, function(x) {return(x[[1]]+qnorm(0.975)*x[[2]])} )
soil.df.rhizo$ci.low.c <- apply(soil.df.rhizo[,c("mean.c","se.c")],1, function(x) {return(x[[1]]-qnorm(0.975)*x[[2]])} )
soil.df.rhizo$ci.high.c <- apply(soil.df.rhizo[,c("mean.c","se.c")],1, function(x) {return(x[[1]]+qnorm(0.975)*x[[2]])} )## Para 01 variável (pH, todas as amostras)
shapiro.test(as.numeric(soil.df.rhizo[1,-c(1,2)]))##
## Shapiro-Wilk normality test
##
## data: as.numeric(soil.df.rhizo[1, -c(1, 2)])
## W = 0.61817, p-value = 2.415e-05## Testes de médias
### Teste T e W para todos
soil.df.rhizo$t.test <- apply(soil.df.rhizo[, samps], 1, t.test.pvalue, samps.o, samps.c)
soil.df.rhizo$w.test <- apply(soil.df.rhizo[, samps], 1, w.test.pvalue, samps.o, samps.c)## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(as.numeric(x[a]), as.numeric(x[b])): cannot
## compute exact p-value with ties### Teste de Shapiro para todos
soil.df.rhizo$s.test.o <- apply(soil.df.rhizo[, samps], 1, s.test.pvalue, samps.o)
soil.df.rhizo$s.test.c <- apply(soil.df.rhizo[, samps], 1, s.test.pvalue, samps.c)
### Teste de KS para todos
soil.df.rhizo$k.test.o <- apply(soil.df.rhizo[, samps], 1, k.test.pvalue, samps.o)## Warning in ks.test(as.numeric(x[a]), "pnorm", mean(as.numeric(x[a])),
## sd(as.numeric(x[a]))): ties should not be present for the Kolmogorov-
## Smirnov test## Warning in ks.test(as.numeric(x[a]), "pnorm", mean(as.numeric(x[a])),
## sd(as.numeric(x[a]))): ties should not be present for the Kolmogorov-
## Smirnov test
## Warning in ks.test(as.numeric(x[a]), "pnorm", mean(as.numeric(x[a])),
## sd(as.numeric(x[a]))): ties should not be present for the Kolmogorov-
## Smirnov testsoil.df.rhizo$k.test.c <- apply(soil.df.rhizo[, samps], 1, k.test.pvalue, samps.c)
### A partir dos resultados, escolhe o teste de médias
soil.df.rhizo$test.selected <-
apply(soil.df.rhizo[,c('s.test.o','s.test.c','k.test.o','k.test.c')], 1, test.selection )
soil.df.rhizo[['p.value.selected']] <-
apply(soil.df.rhizo[,c('w.test','t.test','test.selected')], 1, function(x) {
return(x[x['test.selected']]) }
)Assim, para os pontos 02, 03 e 04, os resultados ficam piores que considerando todos os pontos, então vamos manter a análise de todos os pontos como diferença de solo:
## Exemplo da matriz
head(soil.df.test)## param description CX C0 C1
## pH.CaCl2 pH.CaCl2 pH (CaCl[2]) 4.990000 5.480000 5.570000
## M.O. M.O. M.O. (g/dm^{3}) 21.361000 17.851400 13.057400
## P P P (mg/dm^{3}) 8.592920 36.779850 15.942920
## S S S (mg/dm^{3}) 5.400750 7.070880 7.676600
## Ca Ca Ca (mmol[c]/dm^{3}) 27.572626 36.298344 35.268909
## Mg Mg Mg (mmol[c]/dm^{3}) 7.967943 7.753899 7.779318
## C2 C3 C4 OX O0 O1 O2
## pH.CaCl2 5.38000 5.420000 5.470000 5.69000 5.820000 5.54000 5.99000
## M.O. 20.38740 16.371401 19.985800 36.13220 18.072201 20.11280 18.37940
## P 11.63536 15.976840 13.522960 32.09588 76.867706 15.94292 4.88776
## S 7.67660 19.252340 6.273480 13.15353 4.812330 8.02738 5.74731
## Ca 24.17763 31.973684 28.289474 56.65084 29.226519 36.61345 33.48684
## Mg 10.17376 7.361019 9.231083 15.35810 6.792742 10.53507 11.31516
## O3 O4 mean.o mean.c sd.o sd.c
## pH.CaCl2 5.53000 5.59000 5.694000 5.464000 0.2030516 0.07162403
## M.O. 19.18260 20.58820 19.267040 17.530680 1.0820639 2.98519205
## P 19.94080 13.52296 26.232429 18.771586 28.8379759 10.22965460
## S 5.57192 5.04575 5.840938 9.589980 1.2798461 5.43204130
## Ca 32.30263 30.72368 32.470625 31.201609 2.8195704 5.02325906
## Mg 10.69096 10.15338 9.897460 8.459815 1.7854145 1.19385928
## se.o se.c ci.low.o ci.high.o ci.low.c ci.high.c
## pH.CaCl2 0.09080744 0.03203124 5.5160207 5.871979 5.401220 5.52678
## M.O. 0.48391368 1.33501847 18.3185869 20.215494 14.914092 20.14727
## P 12.89673487 4.57484062 0.9552933 51.509565 9.805063 27.73811
## S 0.57236457 2.42928272 4.7191240 6.962752 4.828673 14.35129
## Ca 1.26095022 2.24646975 29.9992075 34.942042 26.798609 35.60461
## Mg 0.79846164 0.53391010 8.3325043 11.462416 7.413371 9.50626
## t.test w.test s.test.o s.test.c k.test.o k.test.c
## pH.CaCl2 0.04395254 0.03174603 0.17934196 0.846374439 0.6807280 0.9415448
## M.O. 0.25621213 0.30952381 0.61840095 0.567082986 0.9718428 0.9709307
## P 0.60046875 1.00000000 0.02069425 0.013769743 0.3471513 0.2891019
## S 0.20049834 0.09369262 0.10063416 0.002812027 0.5490627 0.2936029
## Ca 0.63552985 0.69047619 0.92323750 0.665934144 0.9971918 0.9761136
## Mg 0.17283304 0.22222222 0.04608693 0.262880682 0.4443207 0.6021168
## test.selected p.value.selected anova kruskal
## pH.CaCl2 t.test 0.043952541 0.04395254 0.02828012
## M.O. t.test 0.256212127 0.25621213 0.25059205
## P w.test 1.000000000 0.60046875 0.91630803
## S w.test 0.093692619 0.17146131 0.07491290
## Ca t.test 0.635529846 0.63552985 0.60150813
## Mg w.test 0.222222222 0.17283304 0.17452534O valor de espécie indicadora (IndVal) é calculado a partir do produto da especificidade de uma espécie a um grupo, comparando com a média através de grupos.
Citando Dave Roberts, o autor do pacote labdsv também usado para o cálculo do Indval:
### Matriz de exemplo
head(table_ITS)## OTU_ID CRZ2 CRZ3 CRZ4 ORZ2 ORZ3 ORZ4
## 1 ITS_Zotu57 17 16 13 45 58 52
## 2 ITS_Zotu476 4 6 0 8 2 7
## 3 ITS_Zotu322 6 2 6 10 3 3
## 4 ITS_Zotu59 51 10 16 1 12 86
## 5 ITS_Zotu20 167 49 107 12 31 75
## 6 ITS_Zotu1 717 0 1459 0 2 2### Grupos para amostras
sample_groups <- c(1,1,1,
2,2,2)
### Puxa as espécies indicadoras por grupos
## Para essa matriz, os nomes das espécies devem ser os nomes das linhas
# Backup
table_ITS_bkp <- table_ITS
# Espécies como nomes de linhas
row.names(table_ITS) <- table_ITS[,1]
# Remove a coluna de espécies
table_ITS <- table_ITS[,-1]
# Calcula o IndVal
# Quebra a máquina!
# indval_ITS <- indval(x = table_ITS,
# clustering = sample_groups)
## Estilo Marcos
ttITS <- t(table_ITS)
## Gera as indicadoras
ITS_IndVal <- indval(ttITS, sample_groups, numitr = 1000)
## Corrige os p-valores
ITS_IndVal$padj <- p.adjust(ITS_IndVal$pval)
# Quantos são menors que 0.05?
sum(ITS_IndVal$padj < 0.05)## [1] 0sum(ITS_IndVal$pval < 0.05)## [1] 0E então para 16S…
### Matriz de exemplo
head(table_16S)## OTU_ID CRZ2 CRZ3 CRZ4 ORZ2 ORZ3 ORZ4
## 1 16S_Zotu5 161 349 108 63 36 58
## 2 16S_Zotu775 1 4 4 4 1 6
## 3 16S_Zotu48 36 11 21 10 38 15
## 4 16S_Zotu522 11 7 14 3 16 11
## 5 16S_Zotu1 735 528 554 535 499 556
## 6 16S_Zotu152 12 19 24 17 30 35### Puxa as espécies indicadoras por grupos
## Para essa matriz, os nomes das espécies devem ser os nomes das linhas
# Backup
table_16S_bkp <- table_16S
# Espécies como nomes de linhas
row.names(table_16S) <- table_16S[,1]
# Remove a coluna de espécies
table_16S <- table_16S[,-1]
## Estilo Marcos
tt16S <- t(table_16S)
## Gera as indicadoras
M16S_IndVal <- indval(tt16S, sample_groups, numitr = 1000)
## Corrige os p-valores
M16S_IndVal$padj <- p.adjust(M16S_IndVal$pval,
method = "fdr")
# Quantos são menors que 0.05?
sum(M16S_IndVal$padj < 0.05)## [1] 0sum(M16S_IndVal$pval < 0.05)## [1] 0Também, nenhum valor de espécie indicadora. Vamos tentar para o WGS!
# Primeiro, ajustar os nomes de colunas
colnames(table_kraken)[1] <- "Taxa"
# Remover o taxonomy
table_kraken_bkp <- table_kraken
table_kraken <- table_kraken[,-8]
# Colocar ID como nome de linha
row.names(table_kraken) <- table_kraken[,1]
# Remove a coluna 1
table_kraken <- table_kraken[,-1]
## Estilo Marcos
ttKRK <- t(table_kraken)
## Gera as indicadoras
WGS_IndVal <- indval(ttKRK, sample_groups, numitr = 1000)
## Corrige os p-valores
WGS_IndVal$padj <- p.adjust(WGS_IndVal$pval)
# Quantos são menors que 0.05?
sum(WGS_IndVal$padj < 0.05)## [1] 0sum(WGS_IndVal$pval < 0.05)## [1] 0