HackMD
Meta'omics - Ab-V5
Ferrarezi et al. (2023): https://doi.org/10.3389/fpls.2023.1172839
Background
The beneficial plant-microbes association has been widely explored with the purpose of understanding the mechanisms involved in this interaction, in order to develop bioinoculants for agriculture aiming at sustainability and reduction of chemical inputs use. Plant growth promoting bacteria (PGPBs) have been reported as an interesting alternative to mineral fertilizers for a range of crops. The genus Azospirillum has shown potential in plant growth promotion for maize. However, the lack of holistic comprehension about PGPBs-plant-native soil microbiota can lead to inconsistency of results in field conditions. Recent ecological theories revealed that plant microbiomes are organized as microbial hubs with keystone species and helpers. The objective of this work is to characterize the microbial community associated with a commercial maize genotype (30A37PW) under influence of the potential keystone PGPB A. brasilense Ab-v5 to support the original development of a synthetic minimum microbiome including this strain as keystone for plant growth promotion. Aiming to understand the PGPB-plant-microbiome system, we conducted a greenhouse experiment in which the synergistic effect of Ab-V5 and soil native microbial community were evaluated through soil and rhizosphere metataxonomics, metagenomics and plant phenomics. Our results revealed that maize seed bacterization with the PGPB Ab-V5 promotes plant growth. This effect is boosted by the interaction with specific microbial groups present in the soil native community. The strain Ab-V5 persists in the maize rhizosphere until 25 days after sowing. In order to understand which microbial groups are important during the interaction Ab-V5-plant-native soil microbial community, 16S rDNA metataxonomics and metagenomics analysis will be conducted in maize soil and rhizosphere samples. This unprecedented omics approach for the holistic view of the complex system composed by inoculant-plant-microbiota and data integration will allow the development of new strategies for optimization of plant biostimulants recommendation with enhanced efficiency for establishment in the field.
Useful R packages:
library("ExpDes.pt")
library("agricolae")
library("laercio")
library("ggplot2")
library("nortest")
library("car")
library("ggfortify")
library("doBy")
library("dplyr")
library("tibble")
library("tidyverse")
library("ggpubr")
library("rstatix")
library("tinytex")
library("factoextra")
library("wesanderson")
require(gridExtra)
require(patchwork)
library("devtools")
library("factoextra")
library(RColorBrewer)
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("mixOmics")
library(mixOmics)
Sparse Partial Least Squares regression (PLS) for discrimination analysis (sPLS-DA, mixOmics)
de<-read.csv("WinRhizo_julho2021 - Dilu-extin_AbV5.csv")
splsda_de<-splsda(de[,4:50],de$Tratamento,keepX=47)
splsda<-plotIndiv(splsda_de,ind.names=TRUE,legend=TRUE,ellipse=TRUE,
star=TRUE,title='sPLS-DA',X.label='PLS-DA 1',Y.label='PLS-DA 2')
splsda<-plotIndiv(splsda_de,pch.levels=c("NS+Ab-V5", "NS 10⁻³+Ab-V5",
"NS 10⁻⁶+Ab-V5", "NS 10⁻⁹+Ab-V5", "SS+Ab-V5"),col.per.group=c("chocolate1",
"coral3", "chartreuse4","deepskyblue4", "lightcoral"),ind.names=TRUE,legend=TRUE,
ellipse=TRUE,star=TRUE,style="ggplot2",title='sPLS-DA',X.label='PLS-DA 1',
Y.label='PLS-DA 2',legend.title="Treatments")
splsda_de$explained_variance
$X
comp 1 comp 2
0.4839414 0.1361531
$Y
comp 1 comp 2
0.2508990 0.2389356
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
Figure 1. Sparse partial least squares discriminant analysis (sPLS-DA) among five treatments based on plant phenotypic traits 25 days after sowing bacterialized seeds. NS: Natural soil (microbial community); SS: Sterilized soil; NS 10⁻³+Ab-V5: dilution 10⁻³ of natural soil; NS 10⁻⁶+Ab-V5: dilution 10⁻⁶ of natural soil; NS 10⁻⁹+Ab-V5: dilution 10⁻⁹ of natural soil.
Descriptive analyses of experimental data
shapiro.test checks the normality of the regression residuals
Ho: the regression residuals follow a normal distribution
Ha: the regression residuals do not follow a normal distribution
Significance alpha: 0.05
SDM: Shoot dry matter (g)
(shapiro.test(de$SDM_g))
Shapiro-Wilk normality test
data: de$SDM_g
W = 0.96468, p-value = 0.08933
leveneTest checks whether the homogeneity of variance between groups is even
Ho: the groups show homogeneity of variance
Ha: the groups do not show homogeneity of variance
Significance alpha: 0.05
(leveneTest(SDM_g ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 6.2562 | 0.0003*** |
| 53 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model.de.sdm<-aov(SDM_g ~ Tratamento, data=de, weights=1/SDM_g)
anova(model.de.sdm)
Analysis of Variance Table
Response: SDM_g
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamento | 4 | 12.1351 | 3.03378 | 21.509 | 1.384e-10 *** |
| Residuals | 53 | 7.4756 | 0.14105 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Test for significant mean differences among groups
LSD test is a non-parametric test
out.de.sdm<-LSD.test(model.de.sdm, c("Tratamento"), main="SDM_g ~ Tratamento, console=TRUE")
out.de.sdm$groups
| SDM_g | groups | |
|---|---|---|
| SS+Ab-V5 | 2.129042 | a |
| NS 10⁻⁹+Ab-V5 | 1.913036 | ab |
| NS 10⁻³+Ab-V5 | 1.701350 | b |
| NS 10⁻⁶+Ab-V5 | 1.604918 | b |
| NS+Ab-V5 | 0.739400 | c |
de.sdm.groups<-out.de.sdm$groups
de.sdm.groups["Tratamento"]<-c("SS+Ab-V5", "NS 10⁻⁹+Ab-V5", "NS 10⁻³+Ab-V5",
"NS 10⁻⁶+Ab-V5", "NS+Ab-V5")
de.sdm.means<-out.de.sdm$means
de$Tratamento<-factor(de$Tratamento,ordered=TRUE,levels=c("NS+Ab-V5",
"NS 10⁻³+Ab-V5","NS 10⁻⁶+Ab-V5","NS 10⁻⁹+Ab-V5","SS+Ab-V5"))
de.sdm<-ggplot(data=de, mapping=aes(x=Tratamento, y=SDM_g,color=Tratamento))+
geom_jitter(data=de,position=position_jitter(0.2))+
scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral"))+
geom_boxplot(width=0.3) + geom_text(data=de.sdm.groups, aes(label=groups),
nudge_y = 0.8, size=5, hjust=1, color="black") + theme_classic() +
theme(legend.title=element_blank())+ labs(y="Shoot dry matter (g)") + theme(legend.position = "none") +
theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10),
text = element_text(size = 10))
de.sdm
de.sdm+coord_flip()
png("de_sdm.png")
de.sdm
dev.off()
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
RDM: Root dry matter (mg)
(shapiro.test(de$RDM_mg))
Shapiro-Wilk normality test
data: de$RDM_mg
W = 0.97529, p-value = 0.2823
(leveneTest(RDM_mg ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 0.9715 | 0.4309 |
| 53 |
model.de.rdm<-aov(RDM_mg ~ Tratamento, data=de)
anova(model.de.rdm)
Analysis of Variance Table
Response: RDM_mg
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamento | 4 | 5208.2 | 1302.05 | 2.4043 | 0.06109 . |
| Residuals | 53 | 28702.3 | 541.55 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
out.de.rdm<-LSD.test(model.de.rdm, c("Tratamento"), main="RDM_mg ~ Tratamento, console=TRUE")
out.de.rdm$groups
| RDM_mg | groups | |
|---|---|---|
| NS 10⁻³+Ab-V5 | 90.82500 | a |
| NS 10⁻⁶+Ab-V5 | 81.50909 | ab |
| SS+Ab-V5 | 73.53333 | ab |
| NS+Ab-V5 | 67.70000 | b |
| NS 10⁻⁹+Ab-V5 | 64.96364 | b |
de.rdm.groups<-out.de.rdm$groups
de.rdm.groups["Tratamento"]<-c("NS 10⁻³+Ab-V5", "NS 10⁻⁶+Ab-V5", "SS+Ab-V5", "NS+Ab-V5", "NS 10⁻⁹+Ab-V5")
de.rdm.means<-out.de.rdm$means
de.rdm<-ggplot(data=de, mapping=aes(x=Tratamento, y=RDM_mg,color=Tratamento))+geom_jitter(data=de,position=position_jitter(0.2))+scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral")) + geom_boxplot(width=0.3) + geom_text(data=de.rdm.groups, aes(label=groups), nudge_y = 40, size=5, hjust=1, color="black") + theme_classic() + theme(legend.title=element_blank())+ labs(y="Root dry matter (mg)") + theme(legend.position = "none") + theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10), text = element_text(size = 10))
de.rdm+coord_flip()
Shoot to root ratio (g.mg⁻¹)
(shapiro.test(de$ShootRoot_ratio))
Shapiro-Wilk normality test
data: de$ShootRoot_ratio
W = 0.95436, p-value = 0.02895
(leveneTest(ShootRoot_ratio ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 3.231 | 0.01908 * |
| 53 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model.de.srr<-aov(ShootRoot_ratio ~ Tratamento, data=de, weights=1/ShootRoot_ratio)
anova(model.de.srr)
Analysis of Variance Table
Response: ShootRoot_ratio
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamento | 4 | 159.32 | 39.829 | 17.832 | 2.494e-09 *** |
| Residuals | 53 | 118.38 | 2.234 |
out.de.srr<-LSD.test(model.de.srr,c("Tratamento"), main="ShootRoot_ratio ~ Tratamento, console=TRUE")
out.de.srr$groups
| ShootRoot_ratio | groups | |
|---|---|---|
| NS 10⁻⁹+Ab-V5 | 32.36755 | a |
| SS+Ab-V5 | 30.06475 | b |
| NS 10⁻⁶+Ab-V5 | 20.26033 | c |
| NS 10⁻³+Ab-V5 | 19.43915 | c |
| NS+Ab-V5 | 11.64676 | d |
de.srr.groups<-out.de.srr$groups
de.srr.groups["Tratamento"]<-c("NS 10⁻⁹+Ab-V5", "SS+Ab-V5", "NS 10⁻⁶+Ab-V5", "NS 10⁻³+Ab-V5", "NS+Ab-V5")
de.srr.means<-out.de.srr$means
de.srr<-ggplot(data=de, mapping=aes(x=Tratamento, y=ShootRoot_ratio, color=Tratamento))+geom_jitter(data=de, position=position_jitter(0.2))+scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral")) + geom_boxplot(width=0.3) + geom_text(data=de.srr.groups, aes(label=groups),nudge_y=10, size=5,hjust=1.5,color="black")+theme_classic() + theme(legend.title=element_blank())+ labs(y="Shoot/root ratio") +theme(legend.position = "none") + theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10), text = element_text(size = 10))
SRL: specific root length (g.cm⁻¹)
(shapiro.test(de$SRL))
Shapiro-Wilk normality test
data: de$SRL
W = 0.93665, p-value = 0.004638
(leveneTest(SRL ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 3.6528 | 0.01061 * |
| 53 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model.de.srl<-aov(SRL ~ Tratamento, data=de, weights=1/SRL)
anova(model.de.srl)
Analysis of Variance Table
Response: SRL
| | Df | Sum Sq |Mean Sq |F value | Pr(>F)|
|Tratamento | 4 |3.1886| 0.79715 | 4.0861| 0.005843 **|
|Residuals |53 |10.3397 |0.19509|
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
out.de.srl<-LSD.test(model.de.srl,c("Tratamento"), main="SRL ~ Tratamento, console=TRUE")
out.de.srl$groups
| SRL | groups | |
|---|---|---|
| NS 10⁻⁶+Ab-V5 | 4.289778 | a |
| SS+Ab-V5 | 3.774809 | b |
| NS 10⁻⁹+Ab-V5 | 3.712263 | bc |
| NS 10⁻³+Ab-V5 | 3.369935 | c |
| NS+Ab-V5 | 2.746755 | d |
de.srl.groups<-out.de.srl$groups
de.srl.groups["Tratamento"]<-c("NS 10⁻⁶+Ab-V5", "SS+Ab-V5", "NS 10⁻⁹+Ab-V5", "NS 10⁻³+Ab-V5", "NS+Ab-V5")
de.srl.means<-out.de.srl$means
de.srl<-ggplot(data=de,mapping=aes(x=Tratamento, y=SRL,color=Tratamento))+geom_jitter(data=de, position=position_jitter(0.2))+scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral")) + geom_boxplot(width=0.3) + geom_text(data=de.srl.groups,aes(label=groups),nudge_y=1.5, size=5,hjust=1.5,color="black")+theme_classic() + theme(legend.title=element_blank())+ labs(y="Specific root length (mg.cm⁻¹)")+theme(legend.position = "none") + theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10), text = element_text(size = 10))
SRSA: specific root surface area (g.cm⁻²)
(shapiro.test(de$SRSA))
Shapiro-Wilk normality test
data: de$SRSA
W = 0.85598, p-value = 6.258e-06
(leveneTest(SRSA ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 4.0478 | 0.006157 ** |
| 53 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model.de.srsa<-aov(SRSA ~ Tratamento, data=de, weights=1/SRSA)
anova(model.de.srsa)
Analysis of Variance Table
Response: SRSA
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamento | 4 | 0.46398 | 0.115994 | 5.4424 | 0.0009547 *** |
| Residuals | 53 | 1.12958 | 0.021313 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
out.de.srsa<-LSD.test(model.de.srsa,c("Tratamento"), main="SRSA ~ Tratamento, console=TRUE")
out.de.srsa$groups
| SRSA | groups | |
|---|---|---|
| NS 10⁻⁶+Ab-V5 | 0.7738787 | a |
| NS 10⁻⁹+Ab-V5 | 0.6468011 | b |
| SS+Ab-V5 | 0.6154526 | bc |
| NS 10⁻³+Ab-V5 | 0.6064365 | bc |
| NS+Ab-V5 | 0.5040541 | c |
de.srsa.groups<-out.de.srsa$groups
de.srsa.groups["Tratamento"]<-c("NS 10⁻⁶+Ab-V5", "NS 10⁻⁹+Ab-V5", "SS+Ab-V5", "NS 10⁻³+Ab-V5", "NS+Ab-V5")
de.srsa.means<-out.de.srsa$means
de.srsa<-ggplot(data=de,mapping=aes(x=Tratamento, y=SRSA,color=Tratamento))+geom_jitter(data=de,position=position_jitter(0.2))+scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral")) + geom_boxplot(width=0.3) + geom_text(data=de.srsa.groups,aes(label=groups),nudge_y=0.5, size=5,hjust=1.5,color="black")+theme_classic() + theme(legend.title=element_blank())+ labs(y="Specific root surface area (mg.cm⁻²)")+theme(legend.position = "none") + theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10), text = element_text(size = 10))
Root length (cm)
(shapiro.test(de$Length.cm.))
Shapiro-Wilk normality test
data: de$Length.cm.
W = 0.97292, p-value = 0.2196
(leveneTest(Length.cm. ~ Tratamento, de))
Levene's Test for Homogeneity of Variance (center = median)
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 4 | 0.6285 | 0.6443 |
| 53 |
model.de.leng<-aov(Length.cm. ~ Tratamento, data=de)
anova(model.de.leng)
Analysis of Variance Table
Response: Length.cm.
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamento | 4 | 131004 | 32751 | 4.47 | 0.003468 ** |
| Residuals | 53 | 388324 | 7327 |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
out.de.leng<-LSD.test(model.de.leng,c("Tratamento"),main="Length.cm. ~ Tratamento, console=TRUE")
out.de.leng$groups
| Length cm | groups | |
|---|---|---|
| NS 10⁻⁶+Ab-V5 | 318.4564 | a |
| NS 10⁻³+Ab-V5 | 301.2228 | a |
| SS+Ab-V5 | 278.6917 | a |
| NS 10⁻⁹+Ab-V5 | 245.8799 | ab |
| NS+Ab-V5 | 184.6067 | b |
de.leng.groups<-out.de.leng$groups
de.leng.groups["Tratamento"]<-c("NS 10⁻⁶+Ab-V5", "NS 10⁻³+Ab-V5", "SS+Ab-V5", "NS 10⁻⁹+Ab-V5", "NS+Ab-V5")
de.leng.means<-out.de.leng$means
de.leng<-ggplot(data=de,mapping=aes(x=Tratamento,y=Length.cm.,color=Tratamento))+geom_jitter(data=de,position=position_jitter(0.2))+scale_colour_manual(values=c("deepskyblue4","chocolate1", "coral3", "chartreuse4","lightcoral")) + geom_boxplot(width=0.3) + geom_text(data=de.leng.groups, aes(label=groups),nudge_y=100,size=5,hjust=1.5,color="black")+theme_classic()+theme(legend.title=element_blank())+ labs(y="Root length (cm)")+theme(legend.position = "none") + theme(axis.title.x=element_blank(), axis.text.y=element_text(size=10), text = element_text(size = 10))
Joining images
part1<-ggarrange(de.sdm,de.rdm,de.srl,de.srsa,de.leng,de.vol,labels=c("A","B","C","D","E","F"),ncol=2,nrow=3)
ggarrange(splsda$graph,part1,ncol=2)
Figure 1. Effect of soil native microbial community grandient amended with Azospirillum brasilense Ab-V5 on maize phenotypic traits. A. Sparse partial least squares discriminant analysis (sPLS-DA) among five treatments based on plant phenotypic traits 25 days after sowing bacterialized seeds. B-G: Maize phenotypic traits evaluated 25 days after sowing bacterialized seeds. NS: Natural soil (microbial community); SS: Sterilized soil; NS 10⁻³+Ab-V5: dilution 10⁻³ of natural soil; NS 10⁻⁶+Ab-V5: dilution 10⁻⁶ of natural soil; NS 10⁻⁹+Ab-V5: dilution 10⁻⁹ of natural soil.
qPCR quantification of 16S rRNA gene and Ab-V5
Soil samples
Soil and rhizosphere samples
Community profiling with Divisive Amplicon Denoising Algorithm 2 - DADA2
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("dada2")
library(dada2); packageVersion("dada2")
Loading required package: Rcpp
[1] ‘1.24.0’ (em 25 de agosto de 2022)
library(ShortRead); packageVersion("ShortRead")
Loading required package: BiocGenerics
Loading required package: parallel
...
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/CN_samples/")
path_CN<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/CN_samples")
sort(list.files(path_CN))
"SA1CN_S49_L001_R1_001.fastq" "SA1CN_S49_L001_R2_001.fastq"
"SA2CN_S50_L001_R1_001.fastq" "SA2CN_S50_L001_R2_001.fastq"
"SA3CN_S51_L001_R1_001.fastq" "SA3CN_S51_L001_R2_001.fastq"
"SA4CN_S52_L001_R1_001.fastq" "SA4CN_S52_L001_R2_001.fastq"
"SE1CN_S53_L001_R1_001.fastq" "SE1CN_S53_L001_R2_001.fastq"
"SE2CN_S54_L001_R1_001.fastq" "SE2CN_S54_L001_R2_001.fastq"
"SE3CN_S55_L001_R1_001.fastq" "SE3CN_S55_L001_R2_001.fastq"
"SE4CN_S56_L001_R1_001.fastq" "SE4CN_S56_L001_R2_001.fastq"
fnFs_CN <- sort(list.files(path_CN, pattern="_R1_001.fastq"))
fnRs_CN <- sort(list.files(path_CN, pattern="_R2_001.fastq"))
sampleNames_CN <- sapply(strsplit(basename(fnFs_CN), "_"), `[`, 1)
sampleNames_CN
"SA1CN" "SA2CN" "SA3CN" "SA4CN" "SE1CN" "SE2CN" "SE3CN" "SE4CN"
fnFs_CN<-file.path(path_CN,fnFs_CN)
fnRs_CN<-file.path(path_CN,fnRs_CN)
png("Fastq_CN_F.png")
plotQualityProfile(fnFs_CN)
dev.off()
Bulk soil samples collected when the native community microbial grandient was inoculated in pots with irradiated soil - Forward reads
png("Fastq_CN_R.png")
plotQualityProfile(fnRs_CN)
dev.off()
Bulk soil samples collected when the native community microbial grandient was inoculated in pots with irradiated soil - Reverse reads
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/0DAP_samples/")
path_0D<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/0DAP_samples")
sort(list.files(path_0D))
"SA10DAP_S57_L001_R1_001.fastq" "SA10DAP_S57_L001_R2_001.fastq"
"SA20DAP_S58_L001_R1_001.fastq" "SA20DAP_S58_L001_R2_001.fastq"
"SA30DAP_S59_L001_R1_001.fastq" "SA30DAP_S59_L001_R2_001.fastq"
"SA40DAP_S60_L001_R1_001.fastq" "SA40DAP_S60_L001_R2_001.fastq"
"SB10DAP_S61_L001_R1_001.fastq" "SB10DAP_S61_L001_R2_001.fastq"
"SB20DAP_S62_L001_R1_001.fastq" "SB20DAP_S62_L001_R2_001.fastq"
"SB30DAP_S63_L001_R1_001.fastq" "SB30DAP_S63_L001_R2_001.fastq"
"SB40DAP_S64_L001_R1_001.fastq" "SB40DAP_S64_L001_R2_001.fastq"
"SC10DAP_S65_L001_R1_001.fastq" "SC10DAP_S65_L001_R2_001.fastq"
"SC20DAP_S66_L001_R1_001.fastq" "SC20DAP_S66_L001_R2_001.fastq"
"SC30DAP_S67_L001_R1_001.fastq" "SC30DAP_S67_L001_R2_001.fastq"
"SC40DAP_S68_L001_R1_001.fastq" "SC40DAP_S68_L001_R2_001.fastq"
"SD10DAP_S69_L001_R1_001.fastq" "SD10DAP_S69_L001_R2_001.fastq"
"SD20DAP_S70_L001_R1_001.fastq" "SD20DAP_S70_L001_R2_001.fastq"
"SD30DAP_S71_L001_R1_001.fastq" "SD30DAP_S71_L001_R2_001.fastq"
"SD40DAP_S72_L001_R1_001.fastq" "SD40DAP_S72_L001_R2_001.fastq"
"SE10DAP_S73_L001_R1_001.fastq" "SE10DAP_S73_L001_R2_001.fastq"
"SE20DAP_S74_L001_R1_001.fastq" "SE20DAP_S74_L001_R2_001.fastq"
"SE30DAP_S75_L001_R1_001.fastq" "SE30DAP_S75_L001_R2_001.fastq"
"SE40DAP_S76_L001_R1_001.fastq" "SE40DAP_S76_L001_R2_001.fastq"
fnFs_0D<-sort(list.files(path_0D, pattern="_R1_001.fastq"))
fnRs_0D<-sort(list.files(path_0D, pattern="_R2_001.fastq"))
sampleNames_0D<-sapply(strsplit(basename(fnFs_0D), "_"), `[`, 1)
sampleNames_0D
"SA10DAP" "SA20DAP" "SA30DAP" "SA40DAP" "SB10DAP" "SB20DAP" "SB30DAP"
"SB40DAP" "SC10DAP" "SC20DAP" "SC30DAP" "SC40DAP" "SD10DAP" "SD20DAP"
"SD30DAP" "SD40DAP" "SE10DAP" "SE20DAP" "SE30DAP" "SE40DAP"
fnFs_0D<-file.path(path_0D,fnFs_0D)
fnRs_0D<-file.path(path_0D,fnRs_0D)
png("Fastq_0D_F.png")
plotQualityProfile(fnFs_0D)
dev.off()
Bulk soil samples collected when on the same day of seed sowing - Forward reads
png("Fastq_0D_R.png")
plotQualityProfile(fnRs_0D)
dev.off()
Bulk soil samples collected when on the same day of seed sowing - Reverse reads
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/15DAP_rhizos_samples")
path_r15D<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/15DAP_rhizos_samples")
sort(list.files(path_r15D))
"RiA115DAP_S117_L001_R1_001.fastq" "RiA115DAP_S117_L001_R2_001.fastq"
"RiA215DAP_S118_L001_R1_001.fastq" "RiA215DAP_S118_L001_R2_001.fastq"
"RiA315DAP_S119_L001_R1_001.fastq" "RiA315DAP_S119_L001_R2_001.fastq"
"RiA415DAP_S120_L001_R1_001.fastq" "RiA415DAP_S120_L001_R2_001.fastq"
"RiB115DAP_S121_L001_R1_001.fastq" "RiB115DAP_S121_L001_R2_001.fastq"
"RiB215DAP_S122_L001_R1_001.fastq" "RiB215DAP_S122_L001_R2_001.fastq"
"RiB315DAP_S123_L001_R1_001.fastq" "RiB315DAP_S123_L001_R2_001.fastq"
"RiB415DAP_S124_L001_R1_001.fastq" "RiB415DAP_S124_L001_R2_001.fastq"
"RiC115DAP_S125_L001_R1_001.fastq" "RiC115DAP_S125_L001_R2_001.fastq"
"RiC215DAP_S126_L001_R1_001.fastq" "RiC215DAP_S126_L001_R2_001.fastq"
"RiC315DAP_S127_L001_R1_001.fastq" "RiC315DAP_S127_L001_R2_001.fastq"
"RiC415DAP_S128_L001_R1_001.fastq" "RiC415DAP_S128_L001_R2_001.fastq"
"RiD115DAP_S129_L001_R1_001.fastq" "RiD115DAP_S129_L001_R2_001.fastq"
"RiD215DAP_S130_L001_R1_001.fastq" "RiD215DAP_S130_L001_R2_001.fastq"
"RiD315DAP_S131_L001_R1_001.fastq" "RiD315DAP_S131_L001_R2_001.fastq"
"RiD415DAP_S132_L001_R1_001.fastq" "RiD415DAP_S132_L001_R2_001.fastq"
"RiE115DAP_S133_L001_R1_001.fastq" "RiE115DAP_S133_L001_R2_001.fastq"
"RiE215DAP_S134_L001_R1_001.fastq" "RiE215DAP_S134_L001_R2_001.fastq"
"RiE315DAP_S135_L001_R1_001.fastq" "RiE315DAP_S135_L001_R2_001.fastq"
"RiE415DAP_S136_L001_R1_001.fastq" "RiE415DAP_S136_L001_R2_001.fastq"
fnFs_r15D<-sort(list.files(path_r15D,pattern="_R1_001.fastq"))
fnRs_r15D<-sort(list.files(path_r15D,pattern="_R2_001.fastq"))
sampleNames_r15D<-sapply(strsplit(basename(fnFs_r15D),"_"), `[`, 1)
sampleNames_r15D
"RiA115DAP" "RiA215DAP" "RiA315DAP" "RiA415DAP" "RiB115DAP" "RiB215DAP"
"RiB315DAP" "RiB415DAP" "RiC115DAP" "RiC215DAP" "RiC315DAP" "RiC415DAP"
"RiD115DAP" "RiD215DAP" "RiD315DAP" "RiD415DAP" "RiE115DAP" "RiE215DAP"
"RiE315DAP" "RiE415DAP"
fnFs_r15D<-file.path(path_r15D,fnFs_r15D)
fnRs_r15D<-file.path(path_r15D,fnRs_r15D)
png("Fastq_r15D_F.png")
dev.off()
Rhizospheric soil samples collected 15 days after seed sowing - Forward reads
png("Fastq_r15D_R.png")
plotQualityProfile(fnRs_r15D)
dev.off()
Rhizospheric soil samples collected 15 days after seed sowing - Reverse reads
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/15DAP_soil_samples/")
path_s15D<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/15DAP_soil_samples")
sort(list.files(path_s15D))
"SA115DAP_S77_L001_R1_001.fastq" "SA115DAP_S77_L001_R2_001.fastq"
"SA215DAP_S78_L001_R1_001.fastq" "SA215DAP_S78_L001_R2_001.fastq"
"SA315DAP_S79_L001_R1_001.fastq" "SA315DAP_S79_L001_R2_001.fastq"
"SA415DAP_S80_L001_R1_001.fastq" "SA415DAP_S80_L001_R2_001.fastq"
"SB115DAP_S81_L001_R1_001.fastq" "SB115DAP_S81_L001_R2_001.fastq"
"SB215DAP_S82_L001_R1_001.fastq" "SB215DAP_S82_L001_R2_001.fastq"
"SB315DAP_S83_L001_R1_001.fastq" "SB315DAP_S83_L001_R2_001.fastq"
"SB415DAP_S84_L001_R1_001.fastq" "SB415DAP_S84_L001_R2_001.fastq"
"SC115DAP_S85_L001_R1_001.fastq" "SC115DAP_S85_L001_R2_001.fastq"
"SC215DAP_S86_L001_R1_001.fastq" "SC215DAP_S86_L001_R2_001.fastq"
"SC315DAP_S87_L001_R1_001.fastq" "SC315DAP_S87_L001_R2_001.fastq"
"SC415DAP_S88_L001_R1_001.fastq" "SC415DAP_S88_L001_R2_001.fastq"
"SD115DAP_S89_L001_R1_001.fastq" "SD115DAP_S89_L001_R2_001.fastq"
"SD215DAP_S90_L001_R1_001.fastq" "SD215DAP_S90_L001_R2_001.fastq"
"SD315DAP_S91_L001_R1_001.fastq" "SD315DAP_S91_L001_R2_001.fastq"
"SD415DAP_S92_L001_R1_001.fastq" "SD415DAP_S92_L001_R2_001.fastq"
"SE115DAP_S93_L001_R1_001.fastq" "SE115DAP_S93_L001_R2_001.fastq"
"SE215DAP_S94_L001_R1_001.fastq" "SE215DAP_S94_L001_R2_001.fastq"
"SE315DAP_S95_L001_R1_001.fastq" "SE315DAP_S95_L001_R2_001.fastq"
"SE415DAP_S96_L001_R1_001.fastq" "SE415DAP_S96_L001_R2_001.fastq"
fnFs_s15D<-sort(list.files(path_s15D,pattern="_R1_001.fastq"))
fnRs_s15D<-sort(list.files(path_s15D,pattern="_R2_001.fastq"))
sampleNames_s15D<-sapply(strsplit(basename(fnFs_s15D),"_"), `[`, 1)
sampleNames_s15D
"SA115DAP" "SA215DAP" "SA315DAP" "SA415DAP" "SB115DAP" "SB215DAP"
"SB315DAP" "SB415DAP" "SC115DAP" "SC215DAP" "SC315DAP" "SC415DAP"
"SD115DAP" "SD215DAP" "SD315DAP" "SD415DAP" "SE115DAP" "SE215DAP"
"SE315DAP" "SE415DAP"
fnFs_s15D<-file.path(path_s15D,fnFs_s15D)
fnRs_s15D<-file.path(path_s15D,fnRs_s15D)
png("Fastq_s15D_F.png")
plotQualityProfile(fnFs_s15D)
dev.off()
Bulk soil samples collected 15 days after seed sowing - Forward reads
png("Fastq_s15D_R.png")
plotQualityProfile(fnRs_s15D)
dev.off()
Bulk soil samples collected 15 days after seed sowing - Reverse reads
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/25DAP_rhizos_samples")
path_r25D<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/25DAP_rhizos_samples")
sort(list.files(path_r25D))
"RiA125DAP_S137_L001_R1_001.fastq" "RiA125DAP_S137_L001_R2_001.fastq"
"RiA225DAP_S138_L001_R1_001.fastq" "RiA225DAP_S138_L001_R2_001.fastq"
"RiA325DAP_S139_L001_R1_001.fastq" "RiA325DAP_S139_L001_R2_001.fastq"
"RiA425DAP_S140_L001_R1_001.fastq" "RiA425DAP_S140_L001_R2_001.fastq"
"RiB125DAP_S141_L001_R1_001.fastq" "RiB125DAP_S141_L001_R2_001.fastq"
"RiB225DAP_S142_L001_R1_001.fastq" "RiB225DAP_S142_L001_R2_001.fastq"
"RiB325DAP_S143_L001_R1_001.fastq" "RiB325DAP_S143_L001_R2_001.fastq"
"RiB425DAP_S144_L001_R1_001.fastq" "RiB425DAP_S144_L001_R2_001.fastq"
"RiC125DAP_S145_L001_R1_001.fastq" "RiC125DAP_S145_L001_R2_001.fastq"
"RiC225DAP_S146_L001_R1_001.fastq" "RiC225DAP_S146_L001_R2_001.fastq"
"RiC325DAP_S147_L001_R1_001.fastq" "RiC325DAP_S147_L001_R2_001.fastq"
"RiC425DAP_S148_L001_R1_001.fastq" "RiC425DAP_S148_L001_R2_001.fastq"
"RiD125DAP_S149_L001_R1_001.fastq" "RiD125DAP_S149_L001_R2_001.fastq"
"RiD225DAP_S150_L001_R1_001.fastq" "RiD225DAP_S150_L001_R2_001.fastq"
"RiD325DAP_S151_L001_R1_001.fastq" "RiD325DAP_S151_L001_R2_001.fastq"
"RiD425DAP_S152_L001_R1_001.fastq" "RiD425DAP_S152_L001_R2_001.fastq"
"RiE125DAP_S153_L001_R1_001.fastq" "RiE125DAP_S153_L001_R2_001.fastq"
"RiE225DAP_S154_L001_R1_001.fastq" "RiE225DAP_S154_L001_R2_001.fastq"
"RiE325DAP_S155_L001_R1_001.fastq" "RiE325DAP_S155_L001_R2_001.fastq"
"RiE425DAP_S156_L001_R1_001.fastq" "RiE425DAP_S156_L001_R2_001.fastq"
fnFs_r25D<-sort(list.files(path_r25D,pattern="_R1_001.fastq"))
fnRs_r25D<-sort(list.files(path_r25D,pattern="_R2_001.fastq"))
sampleNames_r25D<-sapply(strsplit(basename(fnFs_r25D),"_"), `[`, 1)
sampleNames_r25D
"RiA125DAP" "RiA225DAP" "RiA325DAP" "RiA425DAP" "RiB125DAP" "RiB225DAP"
"RiB325DAP" "RiB425DAP" "RiC125DAP" "RiC225DAP" "RiC325DAP" "RiC425DAP"
"RiD125DAP" "RiD225DAP" "RiD325DAP" "RiD425DAP" "RiE125DAP" "RiE225DAP"
"RiE325DAP" "RiE425DAP"
fnFs_r25D<-file.path(path_r25D,fnFs_r25D)
fnRs_r25D<-file.path(path_r25D,fnRs_r25D)
png("Fastq_r25D_F.png")
plotQualityProfile(fnFs_r25D)
dev.off()
Rhizospheric soil samples collected 25 days after seed sowing - Forward reads
png("Fastq_r25D_R.png")
plotQualityProfile(fnRs_r25D)
dev.off()
Rhizospheric soil samples collected 25 days after seed sowing - Reverse reads
setwd("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/25DAP_soil_samples")
path_s25D<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/25DAP_soil_samples")
sort(list.files(path_s25D))
"SA125DAP_S97_L001_R1_001.fastq" "SA125DAP_S97_L001_R2_001.fastq"
"SA225DAP_S98_L001_R1_001.fastq" "SA225DAP_S98_L001_R2_001.fastq"
"SA325DAP_S99_L001_R1_001.fastq" "SA325DAP_S99_L001_R2_001.fastq"
"SA425DAP_S100_L001_R1_001.fastq" "SA425DAP_S100_L001_R2_001.fastq"
"SB125DAP_S101_L001_R1_001.fastq" "SB125DAP_S101_L001_R2_001.fastq"
"SB225DAP_S102_L001_R1_001.fastq" "SB225DAP_S102_L001_R2_001.fastq"
"SB325DAP_S103_L001_R1_001.fastq" "SB325DAP_S103_L001_R2_001.fastq"
"SB425DAP_S104_L001_R1_001.fastq" "SB425DAP_S104_L001_R2_001.fastq"
"SC125DAP_S105_L001_R1_001.fastq" "SC125DAP_S105_L001_R2_001.fastq"
"SC225DAP_S106_L001_R1_001.fastq" "SC225DAP_S106_L001_R2_001.fastq"
"SC325DAP_S107_L001_R1_001.fastq" "SC325DAP_S107_L001_R2_001.fastq"
"SC425DAP_S108_L001_R1_001.fastq" "SC425DAP_S108_L001_R2_001.fastq"
"SD125DAP_S109_L001_R1_001.fastq" "SD125DAP_S109_L001_R2_001.fastq"
"SD225DAP_S110_L001_R1_001.fastq" "SD225DAP_S110_L001_R2_001.fastq"
"SD325DAP_S111_L001_R1_001.fastq" "SD325DAP_S111_L001_R2_001.fastq"
"SD425DAP_S112_L001_R1_001.fastq" "SD425DAP_S112_L001_R2_001.fastq"
"SE125DAP_S113_L001_R1_001.fastq" "SE125DAP_S113_L001_R2_001.fastq"
"SE225DAP_S114_L001_R1_001.fastq" "SE225DAP_S114_L001_R2_001.fastq"
"SE325DAP_S115_L001_R1_001.fastq" "SE325DAP_S115_L001_R2_001.fastq"
"SE425DAP_S116_L001_R1_001.fastq" "SE425DAP_S116_L001_R2_001.fastq"
fnFs_s25D<-sort(list.files(path_s25D, pattern="_R1_001.fastq"))
fnRs_s25D<-sort(list.files(path_s25D, pattern="_R2_001.fastq"))
sampleNames_s25D<-sapply(strsplit(basename(fnFs_s25D),"_"), `[`, 1)
sampleNames_s25D
"SA125DAP" "SA225DAP" "SA325DAP" "SA425DAP" "SB125DAP" "SB225DAP"
"SB325DAP" "SB425DAP" "SC125DAP" "SC225DAP" "SC325DAP" "SC425DAP"
"SD125DAP" "SD225DAP" "SD325DAP" "SD425DAP" "SE125DAP" "SE225DAP"
"SE325DAP" "SE425DAP"
fnFs_s25D<-file.path(path_s25D,fnFs_s25D)
fnRs_s25D<-file.path(path_s25D,fnRs_s25D)
png("Fastq_s25D_F.png")
plotQualityProfile(fnFs_s25D)
dev.off()
Bulk soil samples collected 25 days after seed sowing - Forward reads
png("Fastq_s25D_R.png")
plotQualityProfile(fnRs_s25D)
dev.off()
Bulk soil samples collected 25 days after seed sowing - Reverse reads
I split bulk and rhizospheric soil samples in two different directories
path<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/soil_samples")
sort(list.files(path))
fnFs <- sort(list.files(path, pattern="_R1_001.fastq"))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq"))
sampleNames <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
sampleNames
"SA10DAP" "SA115DAP" "SA125DAP" "SA1CN" "SA20DAP" "SA215DAP"
"SA225DAP" "SA2CN" "SA30DAP" "SA315DAP" "SA325DAP" "SA3CN"
"SA40DAP" "SA415DAP" "SA425DAP" "SA4CN" "SB10DAP" "SB115DAP"
"SB125DAP" "SB20DAP" "SB215DAP" "SB225DAP" "SB30DAP" "SB315DAP"
"SB325DAP" "SB40DAP" "SB415DAP" "SB425DAP" "SC10DAP" "SC115DAP"
"SC125DAP" "SC20DAP" "SC215DAP" "SC225DAP" "SC30DAP" "SC315DAP"
"SC325DAP" "SC40DAP" "SC415DAP" "SC425DAP" "SD10DAP" "SD115DAP"
"SD125DAP" "SD20DAP" "SD215DAP" "SD225DAP" "SD30DAP" "SD315DAP"
"SD325DAP" "SD40DAP" "SD415DAP" "SD425DAP" "SE10DAP" "SE115DAP"
"SE125DAP" "SE1CN" "SE20DAP" "SE215DAP" "SE225DAP" "SE2CN"
"SE30DAP" "SE315DAP" "SE325DAP" "SE3CN" "SE40DAP" "SE415DAP"
"SE425DAP" "SE4CN"
fnRs <- file.path(path, fnRs)
filt_path <- file.path(path, "filtered")
filtFs <- file.path(filt_path, paste0(sampleNames, "_F_filt.fastq"))
filtRs <- file.path(filt_path, paste0(sampleNames, "_R_filt.fastq"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(245,245), maxN=0, maxEE=c(2,2),truncQ=2, rm.phix=TRUE,compress=TRUE, multithread=TRUE)
png("soil_filtered.png", width=16, height=8, units='in', res=300)
plotQualityProfile(filtFs)
dev.off()
truncLen=c(245,245)
errF<-learnErrors(filtFs,multithread=TRUE)
101364585 total bases in 413733 reads from 13 samples will be used for learning the error rates.
errR<-learnErrors(filtRs,multithread=TRUE)
101364585 total bases in 413733 reads from 13 samples will be used for learning the error rates.
png("errors.png")
plotErrors(errF, nominalQ=TRUE)
derepFs <- derepFastq(filtFs)
derepRs <- derepFastq(filtRs)
names(derepFs) <- sampleNames
names(derepRs) <- sampleNames
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
[1] 68 40709
table(nchar(getSequences(seqtab)))
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
Identified 27783 bimeras out of 40709 input sequences.
dim(seqtab.nochim)
[1] 68 12926
sum(seqtab.nochim)/sum(seqtab)
[1] 0.6643015
table(nchar(getSequences(seqtab.nochim)))
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sampleNames
head(track)
| input | filtered | denoisedF | denoisedR | merged | nonchim | |
|---|---|---|---|---|---|---|
| SA10DAP | 25333 | 22590 | 19950 | 19153 | 12628 | 10430 |
| SA115DAP | 26252 | 23601 | 20614 | 20295 | 12855 | 10599 |
| SA125DAP | 23078 | 20474 | 17929 | 17140 | 10775 | 8549 |
| SA1CN | 26257 | 23451 | 19992 | 19277 | 11421 | 10338 |
| SA20DAP | 24037 | 21445 | 18817 | 18012 | 11393 | 10175 |
| SA215DAP | 22504 | 19993 | 17090 | 16279 | 9926 | 8934 |
write.table(track, "track.txt", sep="\t", quote=F)
taxa <- assignTaxonomy(seqtab.nochim,"/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/silva_nr99_v138.1_train_set.fa", multithread=TRUE)
write.table(taxa, "taxa.txt", sep="\t", quote=F)
taxa.
t <- taxa
rownames(taxa.print) <- NULL
head(taxa.print)
| Kingdom | Phylum | Class | Order | Family | Genus |
|---|---|---|---|---|---|
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Cupriavidus" |
| "Bacteria" | "Firmicutes" | "Bacilli" | "Bacillales" | NA | NA |
| "Bacteria" | "Proteobacteria" | "Alphaproteobacteria" | "Rhizobiales" | "Devosiaceae" | "Devosia" |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Pseudomonadales" | "Pseudomonadaceae" | "Pseudomonas" |
| "Bacteria" | "Acidobacteriota" | "Vicinamibacteria" | "Vicinamibacterales" | "Vicinamibacteraceae" | NA |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Pseudomonadales" | "Pseudomonadaceae" | "Pseudomonas" |
asv_seqs <- colnames(seqtab.nochim)
asv_headers <- vector(dim(seqtab.nochim)[2], mode="character")
for (i in 1:dim(seqtab.nochim)[2]) {
asv_headers[i] <- paste(">ASV", i, sep="_")
}
asv_fasta <- c(rbind(asv_headers, asv_seqs))
write(asv_fasta, "ASVs_solo.fa")
asv_tab <- t(seqtab.nochim)
row.names(asv_tab) <- sub(">", "", asv_headers)
write.table(asv_tab, "ASVs_solo_counts.txt", sep="\t", quote=F)
asv_tax <- taxa
row.names(asv_tax) <- sub(">", "", asv_headers)
write.table(asv_tax, "ASVs_solo_taxonomy.txt", sep="\t", quote=F)
counts <- read.table("ASVs_solo_counts.txt", header = TRUE, sep = "\t", dec = ".")
taxonomy <- read.table("ASVs_solo_taxonomy.txt", header = TRUE, sep = "\t", dec = ".")
counts_taxonomy <- as.data.frame(c(taxonomy, counts)) #concatenar
write.table(counts_taxonomy, "solo_counts_taxonomy.txt", sep="\t", quote=F)
library(vegan); packageVersion("vegan")
rarefaction <- t(counts)
S <- specnumber(rarefaction)
raremax <- min(rowSums(rarefaction))
Srare <- rarefy(rarefaction, raremax)
png("rarefaction.png")
rarecurve(rarefaction, step = 1, sample = raremax, col = "blue", cex = 0.6)
empty rows removed
dev.off()
Rarefaction curves
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("phyloseq")
library(Biostrings)
library(ggplot2)
library(phyloseq)
samples.out <- rownames(seqtab.nochim)
subject <- sapply(strsplit(samples.out, "AP"), `[`, 1)
treatment<-substr(subject,2,3)
day<-substr(subject,4,7)
samdf<-data.frame(Treatment=treatment,Day=day)
treatment[treatment=="iA"]="A"
treatment[treatment=="iB"]="B"
treatment[treatment=="iC"]="C"
treatment[treatment=="iD"]="D"
treatment[treatment=="iE"]="E"
treatment[treatment=="A1"]="A"
treatment[treatment=="A2"]="A"
treatment[treatment=="A3"]="A"
treatment[treatment=="A4"]="A"
treatment[treatment=="A5"]="A"
treatment[treatment=="B1"]="B"
treatment[treatment=="B2"]="B"
treatment[treatment=="B3"]="B"
treatment[treatment=="B4"]="B"
treatment[treatment=="B5"]="B"
treatment[treatment=="C1"]="C"
treatment[treatment=="C2"]="C"
treatment[treatment=="C3"]="C"
treatment[treatment=="C4"]="C"
treatment[treatment=="C5"]="C"
treatment[treatment=="D1"]="D"
treatment[treatment=="D2"]="D"
treatment[treatment=="D3"]="D"
treatment[treatment=="D4"]="D"
treatment[treatment=="D5"]="D"
treatment[treatment=="E1"]="E"
treatment[treatment=="E2"]="E"
treatment[treatment=="E3"]="E"
treatment[treatment=="E4"]="E"
treatment[treatment=="E5"]="E"
day[day=="115D"]<-"15D"
day[day=="215D"]<-"15D"
day[day=="315D"]<-"15D"
day[day=="415D"]<-"15D"
day[day=="425D"]<-"25D"
day[day=="325D"]<-"25D"
day[day=="225D"]<-"25D"
day[day=="125D"]<-"25D"
samdf$Treatments[samdf$Treatment=="A"]<-"Native soil + Ab-V5"
samdf$Treatments[samdf$Treatment=="B"]<-"Native soil 10⁻³ + Ab-V5"
samdf$Treatments[samdf$Treatment=="C"]<-"Native soil 10⁻⁶ + Ab-V5"
samdf$Treatments[samdf$Treatment=="D"]<-"Native soil 10⁻⁹ + Ab-V5"
samdf$Treatments[samdf$Treatment=="E"]<-"Sterilized soil + Ab-V5"
neworder<-c("Native soil + Ab-V5", "Native soil 10⁻³ + Ab-V5", "Native soil 10⁻⁶ + Ab-V5", "Native soil 10⁻⁹ + Ab-V5", "Sterilized soil + Ab-V5")
neworderday<-c("CN", "0D", "15D", "25D")
samdf$Treatments<-as.character(samdf$Treatments)
samdf$Treatments<-factor(samdf$Treatments, levels=neworder)
samdf$Day<-as.character(samdf$Day)
samdf$Day<-factor(samdf$Day,levels=neworderday)
rownames(samdf)<-samples.out
ps<-phyloseq(otu_table(seqtab.nochim,taxa_are_rows=FALSE),sample_data(samdf),tax_table(taxa))
ps
dna<-Biostrings::DNAStringSet(taxa_names(ps))
names(dna)<-taxa_names(ps)
ps<-merge_phyloseq(ps,dna)
taxa_names(ps)<-paste0("ASV",seq(ntaxa(ps)))
ps
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 12926 taxa and 68 samples ] |
| sample_data() | Sample Data: | [ 68 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 12926 taxa by 6 taxonomic ranks ] |
| refseq() | DNAStringSet: | [ 12926 reference sequences ] |
azospirillum<-subset_taxa(ps,Genus=="Azospirillum")
azospirillum
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 52 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 52 taxa by 7 taxonomic ranks ] |
| refseq() | DNAStringSet: | [ 52 reference sequences ] |
abv5<-subset_taxa(ps,Species=="Azospirillum_brasilense_AbV5")
rz2ms9<-subset_taxa(ps,Species=="Bacillus_thurigiensis_RZ2MS9")
bacillus<-subset_taxa(ps,Genus=="Bacillus")
bacillus
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 143 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 143 taxa by 7 taxonomic ranks ] |
| refseq() | DNAStringSet: | [ 143 reference sequences ] |
png("plot_richness.png", width=8, height=4,units = 'in', res = 300)
plot_richness(ps,x="Day",measures=c("Shannon", "Simpson"),color="Treatments")
dev.off()
Transform data to proportions as appropriate for Bray-Curtis distances
ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")
png("top20.png", width=8, height=4,units = 'in', res = 300)
plot_bar(ps.top20, x="Treatments", fill="Family") + facet_wrap(~Day, scales="free_x")
Warning message:
Removed 20 rows containing missing values (position_stack).
dev.off()
path<-("/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/rhizosphere_samples")
sort(list.files(path))
fnFs<-sort(list.files(path,pattern="_R1_001.fastq"))
fnRs<-sort(list.files(path,pattern="_R2_001.fastq"))
sampleNames<-sapply(strsplit(basename(fnFs),"_"), `[`, 1)
sampleNames
"RiA115DAP" "RiA125DAP" "RiA215DAP" "RiA225DAP" "RiA315DAP" "RiA325DAP"
"RiA415DAP" "RiA425DAP" "RiB115DAP" "RiB125DAP" "RiB215DAP" "RiB225DAP"
"RiB315DAP" "RiB325DAP" "RiB415DAP" "RiB425DAP" "RiC115DAP" "RiC125DAP"
"RiC215DAP" "RiC225DAP" "RiC315DAP" "RiC325DAP" "RiC415DAP" "RiC425DAP"
"RiD115DAP" "RiD125DAP" "RiD215DAP" "RiD225DAP" "RiD315DAP" "RiD325DAP"
"RiD415DAP" "RiD425DAP" "RiE115DAP" "RiE125DAP" "RiE215DAP" "RiE225DAP"
"RiE315DAP" "RiE325DAP" "RiE415DAP" "RiE425DAP"
fnFs<-file.path(path,fnFs)
fnRs<-file.path(path,fnRs)
png("Fastq_rhizos.png")
plotQualityProfile(fnFs)
filt_path<-file.path(path,"filtered")
filtFs<-file.path(filt_path,paste0(sampleNames,"_F_filt.fastq"))
filtRs<-file.path(filt_path,paste0(sampleNames,"_R_filt.fastq"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(245,240), maxN=0, maxEE=c(2,2),truncQ=2, rm.phix=TRUE,compress=TRUE, multithread=TRUE)
Creating output directory: /workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/rhizosphere_samples/filtered
png("rhizos_filtered.png",width=16,height=8,units='in', res=300)
plotQualityProfile(filtFs)
errF<-learnErrors(filtFs,multithread=TRUE)
101659320 total bases in 414936 reads from 13 samples will be used for learning the error rates.
> errR <- learnErrors(filtRs, multithread=TRUE)
105035520 total bases in 437648 reads from 14 samples will be used for learning the error rates.
png("errors_rhizo.png")
plotErrors(errF, nominalQ=TRUE)
derepFs <- derepFastq(filtFs)
derepRs <- derepFastq(filtRs)
names(derepFs) <- sampleNames
names(derepRs) <- sampleNames
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
40 36926
table(nchar(getSequences(seqtab)))
seqtab.nochim<-removeBimeraDenovo(seqtab,method="consensus",multithread=TRUE,verbose=TRUE)
Identified 26817 bimeras out of 36926 input sequences.
dim(seqtab.nochim)
40 10109
sum(seqtab.nochim)/sum(seqtab)
0.6143463
table(nchar(getSequences(seqtab.nochim)))
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sampleNames
head(track)
| input | filtered | denoisedF | denoisedR | merged | nonchim | |
|---|---|---|---|---|---|---|
| RiA115DAP | 27485 | 24492 | 21942 | 21139 | 17481 | 11889 |
| RiA125DAP | 36829 | 33314 | 30040 | 29276 | 22362 | 16139 |
| RiA215DAP | 45639 | 40605 | 36722 | 35786 | 27716 | 20101 |
| RiA225DAP | 33937 | 30835 | 28757 | 27981 | 22923 | 14823 |
| RiA315DAP | 41524 | 36336 | 32828 | 31987 | 25715 | 18324 |
| RiA325DAP | 34988 | 32045 | 30663 | 30388 | 28717 | 17201 |
write.table(track, "track.txt", sep="\t", quote=F)
taxa <- assignTaxonomy(seqtab.nochim,"/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/silva_nr99_v138.1_train_set.fa", multithread=TRUE)
is.chloro<-taxa[,"Order"] %in% "Chloroplast"
seqtab.nochloro<-seqtab.nochim[,!is.chloro]
taxa <- assignTaxonomy(seqtab.nochloro,"/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/silva_nr99_v138.1_train_set.fa", multithread=TRUE)
is.mito<-taxa[,"Family"] %in% "Mitochondria"
seqtab.noplasts<-seqtab.nochloro[,!is.mito]
taxa <- assignTaxonomy(seqtab.noplasts,"/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/silva_nr99_v138.1_train_set.fa", multithread=TRUE)
write.table(taxa, "taxa.txt", sep="\t", quote=F)
taxa.print <- taxa
rownames(taxa.print) <- NULL
head(taxa.print)
| Kingdom | Phylum | Class | Order | Family | Genus |
|---|---|---|---|---|---|
| "Bacteria" | "Cyanobacteria" | "Cyanobacteriia" | "Chloroplast" | NA | NA |
| "Bacteria" | "Cyanobacteria" | "Cyanobacteriia" | "Chloroplast" | NA | NA |
| "Bacteria" | "Actinobacteriota" | "Actinobacteria" | "Micrococcales" | "Micrococcaceae" | "Pseudarthrobacter" |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Pseudomonadales" | "Pseudomonadaceae" | "Pseudomonas" |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" |
asv_seqs <- colnames(seqtab.noplasts)
asv_headers <- vector(dim(seqtab.noplasts)[2], mode="character")
for (i in 1:dim(seqtab.noplasts)[2]) {
asv_headers[i] <- paste(">ASV", i, sep="_")
}
asv_fasta <- c(rbind(asv_headers, asv_seqs))
write(asv_fasta,"ASVs_rhizo.fa")
asv_tab<-t(seqtab.noplasts)
row.names(asv_tab) <- sub(">", "", asv_headers)
write.table(asv_tab,"ASVs_rhizo_counts.txt",sep="\t", quote=F)
asv_tax<-taxa
row.names(asv_tax) <- sub(">", "", asv_headers)
write.table(asv_tax,"ASVs_rhizo_taxonomy.txt",sep="\t", quote=F)
counts<-read.table("ASVs_rhizo_counts.txt",header = TRUE, sep = "\t", dec = ".")
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
[1] 40 36926
table(nchar(getSequences(seqtab)))
seqtab.nochim<-removeBimeraDenovo(seqtab,method="consensus",multithread=TRUE,verbose=TRUE)
Identified 26817 bimeras out of 36926 input sequences.
is.chloro<-taxa[,"Order"] %in% "Chloroplast"
seqtab.nochloro<-seqtab.nochim[,!is.chloro]
is.mito<-taxa[,"Family"] %in% "Mitochondria"
seqtab.noplasts<-seqtab.nochloro[,!is.mito]
dim(seqtab.noplasts)
40 10006
sum(seqtab.noplasts)/sum(seqtab)
0.5835919
table(nchar(getSequences(seqtab.noplasts)))
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.noplasts))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonplasts")
rownames(track) <- sampleNames
head(track)
| input | filtered | denoisedF | denoisedR | merged | nonplasts | |
|---|---|---|---|---|---|---|
| RiA115DAP | 27485 | 24492 | 21942 | 21139 | 17481 | 10402 |
| RiA125DAP | 36829 | 33314 | 30040 | 29276 | 22362 | 15680 |
| RiA215DAP | 45639 | 40605 | 36722 | 35786 | 27716 | 18940 |
| RiA225DAP | 33937 | 30835 | 28757 | 27981 | 22923 | 14585 |
| RiA315DAP | 41524 | 36336 | 32828 | 31987 | 25715 | 18082 |
| RiA325DAP | 34988 | 32045 | 30663 | 30388 | 28717 | 16523 |
write.table(track, "track.txt", sep="\t", quote=F)
taxa <- assignTaxonomy(seqtab.noplasts,"/workspace/mquecine/JessicaFerrarezi/Metataxonomics/dados/brutos/silva_nr99_v138.1_train_set.fa", multithread=TRUE)
write.table(taxa, "taxa.txt", sep="\t", quote=F)
taxa.print <- taxa
rownames(taxa.print) <- NULL
head(taxa.print)
| Kingdom | Phylum | Class | Order | Family | Genus | |
|---|---|---|---|---|---|---|
| "Bacteria" | "Actinobacteriota" | "Actinobacteria" | "Micrococcales" | "Micrococcaceae" | "Pseudarthrobacter" | |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" | |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Pseudomonadales" | "Pseudomonadaceae" | "Pseudomonas" | |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" | |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" | |
| "Bacteria" | "Proteobacteria" | "Gammaproteobacteria" | "Burkholderiales" | "Burkholderiaceae" | "Ralstonia" |
asv_seqs <- colnames(seqtab.noplasts)
asv_headers<-vector(dim(seqtab.noplasts)[2],mode="character")
for (i in 1:dim(seqtab.noplasts)[2]) {
asv_headers[i] <- paste(">ASV", i, sep="_")
}
asv_fasta <- c(rbind(asv_headers, asv_seqs))
write(asv_fasta,"ASVs_rhizo.fa")
asv_tab<-t(seqtab.noplasts)
row.names(asv_tab) <- sub(">", "", asv_headers)
write.table(asv_tab,"ASVs_rhizo_counts.txt",sep="\t", quote=F)
asv_tax<-taxa
row.names(asv_tax) <- sub(">", "", asv_headers)
write.table(asv_tax,"ASVs_rhizo_taxonomy.txt",sep="\t", quote=F)
counts<-read.table("ASVs_rhizo_counts.txt",header = TRUE, sep = "\t", dec = ".")
taxonomy<-read.table("ASVs_rhizo_taxonomy.txt",header=TRUE,sep="\t",dec=".")
counts_taxonomy<-as.data.frame(c(taxonomy,counts))
write.table(counts_taxonomy,"rhizo_counts_taxonomy.txt", sep="\t",quote=F)
library(vegan)
rarefaction<-t(counts)
S<-specnumber(rarefaction)
raremax<-min(rowSums(rarefaction))
Srare<-rarefy(rarefaction,raremax)
png("rarefaction_r.png")
rarecurve(rarefaction,step=1,sample=raremax,col="blue",cex=0.6)
library(Biostrings)
library(ggplot2)
library(phyloseq)
samples.out <- rownames(seqtab.noplasts)
subject <- sapply(strsplit(samples.out, "AP"), `[`, 1)
treatment<-substr(subject,2,3)
day<-substr(subject,4,7)
samdf<-data.frame(Treatment=treatment,Day=day)
samdf$Treatment[samdf$Treatment=="iA"]<-"A"
samdf$Treatment[samdf$Treatment=="iB"]<-"B"
samdf$Treatment[samdf$Treatment=="iC"]<-"C"
samdf$Treatment[samdf$Treatment=="iD"]<-"D"
samdf$Treatment[samdf$Treatment=="iE"]<-"E"
samdf$Treatment[samdf$Treatment=="E1"]<-"E"
samdf$Treatment[samdf$Treatment=="D1"]<-"D"
samdf$Treatment[samdf$Treatment=="C1"]<-"C"
samdf$Treatment[samdf$Treatment=="B1"]<-"B"
samdf$Treatment[samdf$Treatment=="A1"]<-"A"
samdf$Treatment[samdf$Treatment=="A2"]<-"A"
samdf$Treatment[samdf$Treatment=="A3"]<-"A"
samdf$Treatment[samdf$Treatment=="A4"]<-"A"
samdf$Treatment[samdf$Treatment=="B2"]<-"B"
samdf$Treatment[samdf$Treatment=="B3"]<-"B"
samdf$Treatment[samdf$Treatment=="B4"]<-"B"
samdf$Treatment[samdf$Treatment=="C2"]<-"C"
samdf$Treatment[samdf$Treatment=="C3"]<-"C"
samdf$Treatment[samdf$Treatment=="C4"]<-"C"
samdf$Treatment[samdf$Treatment=="D2"]<-"D"
samdf$Treatment[samdf$Treatment=="D3"]<-"D"
samdf$Treatment[samdf$Treatment=="D4"]<-"D"
samdf$Treatment[samdf$Treatment=="E2"]<-"E"
samdf$Treatment[samdf$Treatment=="E3"]<-"E"
samdf$Treatment[samdf$Treatment=="E4"]<-"E"
samdf$Day[samdf$Day=="115D"]<-"15D"
samdf$Day[samdf$Day=="215D"]<-"15D"
samdf$Day[samdf$Day=="315D"]<-"15D"
samdf$Day[samdf$Day=="415D"]<-"15D"
samdf$Day[samdf$Day=="125D"]<-"25D"
samdf$Day[samdf$Day=="225D"]<-"25D"
samdf$Day[samdf$Day=="325D"]<-"25D"
samdf$Day[samdf$Day=="425D"]<-"25D"
samdf$Niche<-c("Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Rhizosphere","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil","Bulk soil")
samdf$Treatments[samdf$Treatment=="A"]<-"Native community"
samdf$Treatments[samdf$Treatment=="B"]<-"10⁻³-Native community"
samdf$Treatments[samdf$Treatment=="C"]<-"10⁻⁶-Native community"
samdf$Treatments[samdf$Treatment=="D"]<-"10⁻⁹-Native community"
samdf$Treatments[samdf$Treatment=="E"]<-"Sterilized soil"
neworder<-c("Native community", "10⁻³-Native community", "10⁻⁶-Native community", "10⁻⁹-Native community", "Sterilized soil")
neworderday<-c("CN","0D","15D","25D")
neworderniche<-c("Bulk soil", "Rhizosphere")
samdf$Treatments<-as.character(samdf$Treatments)
samdf$Treatments<-factor(samdf$Treatments, levels=neworder)
samdf$Day<-as.character(samdf$Day)
samdf$Day<-factor(samdf$Day,levels=neworderday)
samdf$Niche<-as.character(samdf$Niche)
samdf$Niche<-factor(samdf$Niche,levels=neworderniche)
rownames(samdf)<-samples.out
ps<-phyloseq(otu_table(seqtab.noplasts,taxa_are_rows=FALSE),sample_data(samdf),tax_table(taxa))
dna<-Biostrings::DNAStringSet(taxa_names(ps))
names(dna)<-taxa_names(ps)
ps<-merge_phyloseq(ps,dna)
taxa_names(ps)<-paste0("ASV",seq(ntaxa(ps)))
png("plot_richness.png", width=8, height=4,units = 'in', res = 300)
plot_richness(ps,x="Day",measures=c("Shannon", "Simpson"),color="Treatments")
richnessdf<-data.frame(samdf,estimate_richness(ps)) ### Tabela completa com índices de riqueza
anova_richness<-aov(Shannon ~ Treatment + Day + Niche,richnessdf)
summary(anova_richness)
library(agricolae)
lsd_richness<-LSD.test(anova_richness,c("Treatment","Day","Niche"),main="Chao1 ~ Treatment+Day+Niche,console=TRUE")
$groups
| Shannon | groups | |
|---|---|---|
| D:15D:Rhizosphere | 6.344453 | a |
| D:25D:Rhizosphere | 6.177422 | ab |
| C:15D:Rhizosphere | 6.137612 | ab |
| B:15D:Rhizosphere | 6.122270 | abc |
| C:25D:Rhizosphere | 6.092669 | abc |
| E:15D:Rhizosphere | 5.902019 | abc |
| B:25D:Rhizosphere | 5.895528 | abc |
| E:25D:Rhizosphere | 5.814584 | abc |
| D:0D:Bulk soil | 5.800548 | abcd |
| E:0D:Bulk soil | 5.750761 | abcde |
| C:0D:Bulk soil | 5.729492 | abcde |
| E:15D:Bulk soil | 5.678020 | abcde |
| D:15D:Bulk soil | 5.610393 | abcdef |
| C:15D:Bulk soil | 5.575967 | abcdef |
| D:25D:Bulk soil | 5.540153 | abcdef |
| A:25D:Rhizosphere | 5.446532 | bcdef |
| B:0D:Bulk soil | 5.442790 | bcdef |
| B:15D:Bulk soil | 5.417646 | bcdef |
| B:25D:Bulk soil | 5.378118 | bcdef |
| E:25D:Bulk soil | 5.364804 | bcdef |
| A:15D:Rhizosphere | 5.310821 | cdef |
| A:15D:Bulk soil | 4.977376 | defg |
| E: CN:Bulk soil | 4.974569 | efg |
| A:0D:Bulk soil | 4.974246 | efg |
| A:25D:Bulk soil | 4.817122 | fg |
| A: CN:Bulk soil | 4.813979 | fg |
| C:25D:Bulk soil | 4.152103 | g |
lsd_richness<-LSD.test(anova_richness,c("Treatment", "Niche"),main="Shannon ~ Treatment+Day+Niche,console=TRUE")
$groups
| Shannon | groups | |
|---|---|---|
| D:Rhizosphere | 6.260937 | a |
| C:Rhizosphere | 6.115140 | ab |
| B:Rhizosphere | 6.008899 | ab |
| E:Rhizosphere | 5.858302 | abc |
| D:Bulk soil | 5.650365 | bc |
| E:Bulk soil | 5.442038 | cd |
| B:Bulk soil | 5.412851 | cd |
| A:Rhizosphere | 5.378676 | cde |
| C:Bulk soil | 5.152521 | de |
| A:Bulk soil | 4.895681 | e |
lsd_richness<-LSD.test(anova_richness,c("Niche", "Day"),main="Shannon ~ Treatment+Day+Niche,console=TRUE")
$groups
| Shannon | groups | |
|---|---|---|
| Rhizosphere:15D | 5.963435 | a |
| Rhizosphere:25D | 5.885347 | ab |
| Bulk soil:0D | 5.539567 | bc |
| Bulk soil:15D | 5.451880 | c |
| Bulk soil:25D | 5.050460 | d |
| Bulk soil:CN | 4.894274 | d |
png("richness_alltogether.png",width=8,height=4,units='in',res=300)
plot_richness(ps,color="Day",shape="Niche",x="Treatments",measures=c("Shannon","Simpson"))
ps.prop <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")
plot_ordination(ps.prop,ord.nmds.bray,color="Treatments",title="Bray NMDS")
top50<-names(sort(taxa_sums(ps),decreasing=TRUE))[1:50]
ps.top50<-transform_sample_counts(ps,function(OTU) OTU/sum(OTU))
ps.top50<-prune_taxa(top50,ps.top50)
png("top50_rhizos.png",width=16,height=8,units='in',res=300)
plot_bar(ps.top50,x="Day",fill="Order")+facet_wrap(~Treatments,scales="free_x")
library(ape)
random_tree<-rtree(ntaxa(ps),rooted=TRUE,tip.label=taxa_names(ps))
ps<-merge_phyloseq(ps,random_tree)
ps
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 10006 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 10006 taxa by 7 taxonomic ranks ] |
| phy_tree() | Phylogenetic Tree: | [ 10006 tips and 10005 internal nodes ] |
| refseq() | DNAStringSet: | [ 10006 reference sequences ] |
png("phytree.png", width=8, height=4,units = 'in', res = 300) plot_tree(ps,color="Treatments",shape="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3)
azospirillum<-subset_taxa(ps,Genus=="Azospirillum")
random_tree_abv5<-rtree(ntaxa(azospirillum),rooted=TRUE,tip.label=taxa_names(azospirillum))
png("random_tree_abv5.png",width=8, height=4,units = 'in', res = 300)
plot(random_tree_abv5)
azospirillum<-merge_phyloseq(azospirillum,random_tree_abv5)
azospirillum
png("phytree_abv5.png", width=8, height=4,units = 'in', res = 300)
plot_tree(azospirillum,color="Treatments",shape="Niche",size="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3, title="Azospirillum spp.")
bacillus<-subset_taxa(ps,Genus=="Bacillus")
bacillus
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 143 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 143 taxa by 7 taxonomic ranks ] |
| refseq() | DNAStringSet: | [ 143 reference sequences ] |
random_tree_bacillus<-rtree(ntaxa(bacillus),rooted=TRUE,tip.label=taxa_names(bacillus))
png("random_tree_bacillus.png", width=8, height=4,units = 'in', res = 300)
plot(random_tree_bacillus)
bacillus<-merge_phyloseq(bacillus, random_tree_bacillus)
bacillus
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 143 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 143 taxa by 7 taxonomic ranks ] |
| phy_tree() | Phylogenetic Tree: | [ 143 tips and 142 internal nodes ] |
| refseq() | DNAStringSet: | [ 143 reference sequences ] |
png("phytree_bacillus.png",width=8, height=4,units = 'in', res = 300)
plot_tree(bacillus, color="Treatments",shape="Niche", size="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3)
paenibacillus<-subset_taxa(ps,Genus=="Paenibacillus")
paenibacillus
| phyloseq-class experiment-level object | ||
|---|---|---|
| otu_table() | OTU Table: | [ 33 taxa and 40 samples ] |
| sample_data() | Sample Data: | [ 40 samples by 3 sample variables ] |
| tax_table() | Taxonomy Table: | [ 33 taxa by 7 taxonomic ranks ] |
| phy_tree() | Phylogenetic Tree: | [ 33 tips and 32 internal nodes ] |
| refseq() | DNAStringSet: | [ 33 reference sequences ] |
random_tree_paenibacillus<-rtree(ntaxa(paenibacillus),rooted=TRUE,tip.label=taxa_names(paenibacillus))
png("random_tree_paenibacillus.png",width=8, height=4,units = 'in', res = 300)
plot(random_tree_paenibacillus)
png("phytree_paenibacillus.png", width=8, height=4,units = 'in', res = 300)
plot_tree(paenibacillus,color="Treatments",shape="Day", label.tips="taxa_names",ladderize="right",plot.margin=0.3)
priestia<-subset_taxa(ps,Genus=="Priestia")
png("phytree_priestia.png",width=8, height=4,units = 'in', res = 300)
plot_tree(priestia,color="Treatments",shape="Day", label.tips="taxa_names",ladderize="right",plot.margin=0.3)
streptomyces<-subset_taxa(ps,Genus=="Streptomyces")
png("phytree_streptomyces.png",width=8, height=4,units = 'in', res = 300)
plot_tree(streptomyces,color="Treatments",shape="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3)
siphonobacter<-subset_taxa(ps,Genus=="Siphonobacter")
png("phytree_siphonobacter.png",width=8, height=4,units = 'in', res = 300)
plot_tree(siphonobacter,color="Treatments",shape="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3)
amycolatopsis<-subset_taxa(ps,Genus=="Amycolatopsis")
png("phytree_amycolatopsis.png",width=8, height=4,units = 'in', res = 300)
plot_tree(amycolatopsis,color="Treatments",shape="Day",label.tips="taxa_names",ladderize="right",plot.margin=0.3)
PICRUSt:
In the sheet counts, it is necessary to add "#OTU ID" in the first cell (first row/first column)
nohup picrust2_pipeline.py -s ASVs_alltogether.fa -i ASVs_alltogether_counts.txt -o picrust2_out_pipeline -p 10 &
Inkscape installation
If your Ubuntu version or derivative has not packaged Inkscape 1.2.1 for you yet, you can install the current version from our ppa for Ubuntu 22.04, 20.04 and 22.10 (Ubuntu development version):
Desktop Instructions
Start the application "Software and updates"
Select the tab "Other Software"
Click the button "Add…"
Enter the following text
ppa:inkscape.dev/stable
Click the button "Add Source"
Click on the button "Close"
When the pop up warns "The information about available software is out-of-date" click the button "Reload"
The latest stable version of Inkscape is now avaialble to be installed from the software center.
Command Line
sudo add-apt-repository ppa:inkscape.dev/stable
sudo apt update
sudo snap install inkscape
