Bioinformatic workflow for Nanopore metabarcoding analysis

install.packages("NAME")

cp -r barcode12 barcode24 barcode36 barcode48 ~/PATH-to-store-barcodes

zcat FAW16131_pass_barcode12_1dba8d5b_de32243e_0.fastq.gz | head -2

zcat * > barcode12.fastq.gz

for name in ./*.fastq.gz; do

porechop-runner.py -i "$name" -o "${name}.pc.fastq.gz"

done

fastqc -o ./FastQC_results *.pc.fastq.gz

chopper --minlength 1000 --maxlength 2000 -i $FILE > $output

$ apptainer exec ./pm2.sif parallel-meta -h

gunzip *

for name in ./*.fastq; do


parallel-meta -r ${name}.fastq -d X -e 1e-30 -o "${name}_parallel"

done

find . -name '*classification.txt' | sed 's/.txt//' > classification_list

#!/bin/bash

while read name; do

cut -f5 ${name}.txt | grep -v 'Classification' |sort | uniq -c | sort -nr | sed 's/; /;/g; s/;$//' | sed 's/^ \+//' | sed 's/ /\t/' | awk -F '\t' '{print $2 "\t" $1}' > ${name}_count_matrix.txt


echo "Finished for ${name}"

done < classification_list

find . -name '*classification_count_matrix.txt' > matrix_list.txt

#!/usr/bin/perl
use strict;
use warnings;

scalar@ARGV == 1 || die "usage: $0 <file_list.txt>
      file_list.txt      List of input files, each one is a matrix with the same format:
                         - No header
                         - Lineage is in the first column separated by ';'
                         - Abundance in integer numbers in the second column
";

my $files = $ARGV[0];
my @names;
my %each_matrix;
open (VALUES, $files) or die ("I can not open the file $files\n") ;
while (<VALUES>) {
	chomp;
	my $name = $_;
	PARSE($name);
	push(@names, $name);
}
close(VALUES);

## Print the matrix
open MATCOUNT, ">tmp" or die ("I can not create the file integrated_matrix.txt, check write permissions on current directory\n"); 
my $header = join("\t", @names);
print MATCOUNT "SampleID\t$header\n"; #<STDIN>;

my $lastname=scalar@names;
my $cont_names=1;
foreach my $tax (sort keys %each_matrix) {
	print MATCOUNT "$tax\t";
	foreach my $file (@names) {
		if ($each_matrix{$tax}{$file}) {
			print MATCOUNT "$each_matrix{$tax}{$file}\t";
		} else {
			print MATCOUNT "0\t";
		}
	}
	print MATCOUNT "\n";
}
close(MATCOUNT);

system("sed 's/\t\$//' tmp > integrated_matrix.txt && rm tmp");

### Subroutine to parse the file
sub PARSE {
	my $file=shift;
#	my @new_array=();
	open (FILE, $file) or die ("I can not open the file $file\n") ;
	while (<FILE>) {
		chomp;
		my($taxo, $abundance)=split('\t', $_);
#		my@lineage=split(/;/, $taxo);
#		foreach my$each_name (@lineage){
#			if ($each_name !~ /epsilon_subdivisions/ and $each_name !~ /[Gg]roup/  and $each_name !~ /unclassified/ and $each_name !~ /[Ii]ncertae_[Ss]edis/ ) {
#				push(@new_array, $each_name);
#			}
#		}
#		if (scalar@new_array == 6) {
#			my$sp=join("_", $new_array[-1], "sp\.");
#			push(@new_array, $sp);
#			$taxo=join(";", @new_array);
#		}else{
#			$taxo=join(";", @new_array);
#		}
		$each_matrix{$taxo}{$file}=$abundance;
#		@new_array=();
	}
	close(FILE);
}

perl matrix.pl matrix_list.txt

sed 's/classification_count_matrix.txt//g ; s/[./]//g ;'s/\'//g' ; s/_S[0-9][0-9][0-9]//g ; s/_S[0-9][0-9]//g ; s/_S[0-9]//g ; /Eukaryota/d;/Mitochondria/d;/Chloroplast/d' integrated_matrix.txt > integrated_matrix_good.txt

#!/usr/bin/perl
use strict;
use warnings;

scalar@ARGV == 2 || die "usage: $0 <matrix.txt> <counts|abun>

      Programa que separa una matrix de abundancias o conteos en niveles taxonomicos

      matrix.txt      Matriz de abundancias donde la primer columna es la asignacion
                      taxonomica y el resto son las abundancias para cada muestra

      counts|abun     Especificar si la matriz es de conteos crudos o abundancias, en 
                      el segundo caso, reporta la proporcion de abundancia asignada a 
                      niveles taxonomicos mas altos del corte dado

";

my $file=$ARGV[0];
my $tipo_mat=$ARGV[1];

unless ($tipo_mat =~ "counts"  or $tipo_mat =~ "abun") {
        die "El tipo de matriz solo puede ser counts o abun en minusculas\n";
}

## Armando el hash de la matriz de abundancia original
my %matrix=();
open (FILE, $file) or die ("No puedo abrir el archivo $file\n");
my $header=<FILE>;
while (<FILE>) {
	chomp;
	my@line=split('\t', $_);
	my$tax_name=shift@line;
	if (exists $matrix{$tax_name}){	
		my@sums=();
		for (my $i=0; $i<scalar@line; $i++){
			my$add=$matrix{$tax_name}[$i]+$line[$i];
			push(@sums, $add);
		}
		$matrix{$tax_name}=[@sums];
	}else{
		$matrix{$tax_name}=[@line];
	}
}
close(FILE);

my @samples=split("\t", $header);
my $samples_number=scalar@samples -1 ;

## Armando las matrices para cada nivel taxonomico
my %domain=SEPARA(\%matrix, 1, $samples_number, $tipo_mat);
open DOM, ">dom_tmp.txt";
print DOM "$header";
foreach my $dom (keys %domain){
	print DOM "$dom\t" ;
	foreach my $count (@{ $domain{$dom} }) {
		print DOM "$count\t";
	} 
		print DOM "\n";
}
close(DOM);
system("sed 's/\t\$//' dom_tmp.txt > domain_matrix.txt");
system("rm dom_tmp.txt");

my %phylum=SEPARA(\%matrix, 2, $samples_number, $tipo_mat);
open PHY, ">phy_tmp.txt";
print PHY "$header";
foreach my $phy (keys %phylum){
	print PHY "$phy\t" ;
	foreach my $count (@{ $phylum{$phy} }) {
		print PHY "$count\t";
	} 
		print PHY "\n";
}
close(PHY);
system("sed 's/\t\$//' phy_tmp.txt > phylum_matrix.txt");
system("rm phy_tmp.txt");

my %class=SEPARA(\%matrix, 3, $samples_number, $tipo_mat);
open CLASS, ">class_tmp.txt";
print CLASS "$header";
foreach my $cla (sort keys %class){
	print CLASS "$cla\t" ;
	foreach my $count (@{ $class{$cla} }) {
		print CLASS "$count\t";
	} 
		print CLASS "\n";
}
close(CLASS);
system("sed 's/\t\$//' class_tmp.txt > class_matrix.txt");
system("rm class_tmp.txt");

my %order=SEPARA(\%matrix, 4, $samples_number, $tipo_mat);
open ORD, ">order_tmp.txt";
print ORD "$header";
foreach my $ord (sort keys %order){
	print ORD "$ord\t" ;
	foreach my $count (@{ $order{$ord} }) {
		print ORD "$count\t";
	} 
		print ORD "\n";
}
close(ORD);
system("sed 's/\t\$//' order_tmp.txt > order_matrix.txt");
system("rm order_tmp.txt");

my %family=SEPARA(\%matrix, 5, $samples_number, $tipo_mat);
open FAM, ">family_tmp.txt";
print FAM "$header";
foreach my $fam (sort keys %family){
	print FAM "$fam\t" ;
	foreach my $count (@{ $family{$fam} }) {
		print FAM "$count\t";
	} 
		print FAM "\n";
}
close(FAM);
system("sed 's/\t\$//' family_tmp.txt > family_matrix.txt");
system("rm family_tmp.txt");

my %genus=SEPARA(\%matrix, 6, $samples_number, $tipo_mat);
open GEN, ">genus_tmp.txt";
print GEN "$header";
foreach my $gen (sort keys %genus){
	print GEN "$gen\t" ;
	foreach my $count (@{ $genus{$gen} }) {
		print GEN "$count\t";
	} 
		print GEN "\n";
}
close(GEN);
system("sed 's/\t\$//' genus_tmp.txt > genus_matrix.txt");
system("rm genus_tmp.txt");

my %species=SEPARA(\%matrix, 7, $samples_number, $tipo_mat);
open SPC, ">species_tmp.txt";
print SPC "$header";
foreach my $spc (sort keys %species){
	print SPC "$spc\t" ;
	foreach my $count (@{ $species{$spc} }) {
		print SPC "$count\t";
	} 
		print SPC "\n";
}
close(SPC);
system("sed 's/\t\$//' species_tmp.txt > species_matrix.txt");
system("rm species_tmp.txt");

my %subspecies=SEPARA(\%matrix, 8, $samples_number, $tipo_mat);
open SSPC, ">subspecies_tmp.txt";
print SSPC "$header";
foreach my $subspc (sort keys %subspecies){
	print SSPC "$subspc\t" ;
	foreach my $count (@{ $subspecies{$subspc} }) {
		print SSPC "$count\t";
	} 
		print SSPC "\n";
}
close(SPC);
system("sed 's/\t\$//' subspecies_tmp.txt > subspecies_matrix.txt");
system("rm subspecies_tmp.txt");

sub SEPARA {
	my %matriz_sub=%{$_[0]};
	my $level=$_[1];
	my $muestras=$_[2];
	my $type_matrix=$_[3];
	my %regresa=();
	my(@sumas, @taxa_levels, @diferencia);
	my($level_key, $pos, $diff,$taxa2print);
	my$suma_val=0;

	foreach my $llaves (sort keys %matriz_sub) {
		@taxa_levels=split(';', $llaves);

		my $tamano=scalar@taxa_levels;
		if ($tamano >= $level) {
			$pos=$level-1;
			
			if ($taxa_levels[$pos] =~ /no_rank/){
				$pos=$level;
			}
			
			$taxa2print=join(';', @taxa_levels[0..$pos]);	### variable nueva para imprimir

			if (!exists($regresa{$taxa2print})) {
				$regresa{$taxa2print}=$matriz_sub{$llaves};
			}else{ 
			
				for (my $i=0; $i<$muestras; $i++) { 
					$suma_val=$regresa{$taxa2print}[$i]+$matriz_sub{$llaves}[$i];
	            			push(@sumas, $suma_val);
				}
			$regresa{$taxa2print}=[@sumas];
			@sumas=();
			}
		}
	}

	if ($type_matrix =~ 'abun'){
		my @suma_per_sample=();
		my $suma_val=0;

		for (my $i=0; $i<$muestras; $i++) { 
			foreach my $key (sort keys %regresa){
				$suma_val=$regresa{$key}[$i]+$suma_val;
			}
			push(@suma_per_sample, $suma_val);
			$suma_val=0;
		}
	
		foreach my $cada_suma (@suma_per_sample){
			my $resta=100-$cada_suma;
			if ($resta > 0) {
				my $higher_level= sprintf "%.6f", $resta;
				push(@diferencia, $higher_level);
			}else{
				push(@diferencia, "0");
			}
		}
		$regresa{"Classified_at_a_higher_level"}=[@diferencia];
		return %regresa;
	}else{
		return %regresa;	
	}
}

perl taxa_levels.pl integrated_matrix_good.txt counts

#!/bin/bash

awk '{print $1}' species_matrix.txt > 7taxonomy_table.txt

cut -f1 species_matrix.txt > 7taxonomy_table.txt
sed 's/;/\t/g' 7taxonomy_table.txt | sed 's/SampleID/Domain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecie/g' > 7_taxonomy_table1.txt
awk 'BEGIN{ FS = OFS = "\t" } { print (NR==1? "TAXA" : "TAXA"NR-1), $0 }' 7_taxonomy_table1.txt > taxonomy_table.txt
cut -f2- species_matrix.txt > 7_otu_table.txt
awk 'BEGIN{ FS = OFS = "\t" } { print (NR==1? "TAXA" : "TAXA"NR-1), $0 }' 7_otu_table.txt > otu_table.txt

taxonkit name2taxid --data-dir $PATH/taxonkit/database/ -i 8 --show-rank taxonomy_table.txt > out_species.txt

awk -F'\t' '$9 == ""' out_species.txt

taxonkit reformat --data-dir /home/pardolab/Apps/taxonkit/database/ -I 9 -f "{K}|{p}|{c}|{o}|{f}|{g}|{s}" -F out_species.txt > species_corrected.txt

cut -f1,11 species_corrected.txt  | sed 's/|/\t/g' | sed "s/'//g" | sed 's/TAXA\t/TAXA    Domain  Phylum  Class   Order   Family  Genus   Specie/' > taxonomytable.txt

head -1 integrated_matrix_good.txt | tr "\t" "\n" > metadata.txt

Class_1up = filter_taxa(Class_perc, function(x) sum(x) > 0.1, TRUE)
Class_1up
#40
Class_top = aggregate_top_taxa2(Class_perc,top = 40, level = "Class")

  facet_nested(. ~ Origin + Depth + Platform , scales = "free_x",space="free_x") +

#libraries
library(phyloseq)
library(dplyr)
library(vegan)
library(tidyverse)
library(tibble)
library(cowplot)
library(ggh4x)
library(factoextra)
library(reshape2)
library(dplyr)
library(microbiome  )
library(ggplot2)
library(ranacapa)
library(microViz)
library(microbiomeutilities)

#2 change working directory
setwd("$PATH_to_files")

################################################################################
#3 importar data 

#16S
otu16s_df<-read.table("otu_table.txt",
                      header=TRUE, row.names=1, check.names=FALSE)
otu16s_mat<-as.matrix(otu16s_df)
tax16s_df<-read.table("taxonomytable.txt",
                      header=TRUE, row.names=1, check.names=FALSE, sep = '\t')
tax16s_mat<-as.matrix(tax16s_df)
metadata16s_df<-read.table("metadata.txt",
                           sep="\t", header=TRUE,  row.names=1)
OTU16s = otu_table(otu16s_mat, taxa_are_rows = TRUE)
TAX16s = tax_table(tax16s_mat)
ps_16s = phyloseq(OTU16s,TAX16s, sample_data(metadata16s_df))
ps_16s

ps_16s_genus = ps_16s %>% aggregate_taxa(level = "Genus")


##### Curves ###############
fuente_rare = theme(text = element_text(size= 12, face="bold"),
                    axis.text.x = element_text(size = 12,face="bold"),
                    axis.text.y = element_text(size = 12,face="bold"),
                    panel.background = element_rect(fill = "white",colour = "grey50"), 
                    panel.grid.major = element_blank(),panel.grid.minor = element_blank())

#genus
ggrare(ps_16s_genus, step = 209, plot = T,
       parallel = T, se = F, color = "Sample" ) + 
  fuente_rare + labs(y="Géneros observados", x="Tamaño de la muestra") 


##################alpha#######################3
alphadiv <- plot_richness(ps_16s, x = "Sample",
                          measures=c("Observed","Chao1", "Shannon", "Simpson"), color = "Sample") + 
  labs(x="Origin", y="Alpha diversity")
alphadiv

##############
#16s 
genus_comp = ps_16s_genus %>% transform(transform = "compositional")
genus_perc = transform_sample_counts(genus_comp, function(x) x*100)
genus_perc
genus_1up = filter_taxa(genus_perc, function(x) sum(x) > 1, TRUE)
genus_1up
#70
genus_5up = filter_taxa(genus_perc, function(x) sum(x) > 2, TRUE)
genus_5up
#30

genus_top = aggregate_top_taxa2(genus_perc,top = 32, level = "Genus")
genus_top
data = psmelt(genus_top) # create dataframe from phyloseq object
data$Genus <- as.character(data$Genus) #convert to character
data$Genus[data$Genus == "Other"] <- "< 2% overall rel. ab."

mycolors= c( "#F6F6F6","#68228B","#2BDEBB","#C0C0C0","#EEAD0E","#CD3333",
             "#CDC8B1","#B565A7","#fff5b5","#66CD00","#53868B",
             "#fcff5d","#8B4513","#009B77","#FF6F61","#CAFF70",
             "#98F5FF","#EFC050","#FFD29B","#00BFFF","#722F37","#BF3EFF",
             "#FF1493","#096e24","#BDFFBC","#93FF33","#FF8C00",
             "#5F9EA0","#955251","#00008B","#FFD700","#CCCCFF",
             "#8B8278","#A6CF4A","#CD661D","#FF0000","#1E90FF","#CDAA7D","#D65076",
             "#FF4040","#DFFF00","#1874CD","#CD950C","#00CDCD",
             "#458B00","#FF3030","#8B6508","#cb7cc3","#76EE00",
             "#2f2aa0","#cb7c80","#92A8D1","#8B0A50","#F7CAC9",
             "#C1FFC1","#fe9000","#2BDE8A","#FA8072","#909aa6", "black")
bars_genus_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Genus)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_genus_top

###############PHYLA####################
ps_16s_Phylum = ps_16s %>% aggregate_taxa(level = "Phylum")

Phylum_comp = ps_16s_Phylum %>% transform(transform = "compositional")
Phylum_perc = transform_sample_counts(Phylum_comp, function(x) x*100)
Phylum_perc
Phylum_1up = filter_taxa(Phylum_perc, function(x) sum(x) > 1, TRUE)
Phylum_1up
#30 y 10
Phylum_top = aggregate_top_taxa2(Phylum_perc,top = 30, level = "Phylum")
Phylum_top
data = psmelt(Phylum_top) # create dataframe from phyloseq object
data$Phylum <- as.character(data$Phylum) #convert to character
data$Phylum[data$Phylum == "Other"] <- "< 2% overall rel. ab."
bars_Phylum_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_Phylum_top

############################################################################

###############Class####################
ps_16s_Class = ps_16s %>% aggregate_taxa(level = "Class")

Class_comp = ps_16s_Class %>% transform(transform = "compositional")
Class_perc = transform_sample_counts(Class_comp, function(x) x*100)
Class_perc
Class_1up = filter_taxa(Class_perc, function(x) sum(x) > 0.1, TRUE)
Class_1up
#64 y 40
Class_top = aggregate_top_taxa2(Class_perc,top = 40, level = "Class")
Class_top
data = psmelt(Class_top) # create dataframe from phyloseq object
data$Class <- as.character(data$Class) #convert to character
data$Class[data$Class == "Other"] <- "< 0.5% overall rel. ab."
bars_Class_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Class)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_Class_top

############################################################################

###############Order####################
ps_16s_Order = ps_16s %>% aggregate_taxa(level = "Order")

Order_comp = ps_16s_Order %>% transform(transform = "compositional")
Order_perc = transform_sample_counts(Order_comp, function(x) x*100)
Order_perc
Order_1up = filter_taxa(Order_perc, function(x) sum(x) > 1, TRUE)
Order_1up
#42
Order_top = aggregate_top_taxa2(Order_perc,top = 42, level = "Order")
Order_top
data = psmelt(Order_top) # create dataframe from phyloseq object
data$Order <- as.character(data$Order) #convert to character
data$Order[data$Order == "Other"] <- "< 1% overall rel. ab."
bars_Order_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Order)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_Order_top

############################################################################

###############Family####################
ps_16s_Family = ps_16s %>% aggregate_taxa(level = "Family")

Family_comp = ps_16s_Family %>% transform(transform = "compositional")
Family_perc = transform_sample_counts(Family_comp, function(x) x*100)
Family_perc
Family_1up = filter_taxa(Family_perc, function(x) sum(x) > 2, TRUE)
Family_1up
#30 y 10
Family_top = aggregate_top_taxa2(Family_perc,top = 35, level = "Family")
Family_top
data = psmelt(Family_top) # create dataframe from phyloseq object
data$Family <- as.character(data$Family) #convert to character
data$Family[data$Family == "Other"] <- "< 2% overall rel. ab."
bars_Family_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Family)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_Family_top

############################################################################

###############Specie####################
ps_16s_Specie = ps_16s %>% aggregate_taxa(level = "Specie")

Specie_comp = ps_16s_Specie %>% transform(transform = "compositional")
Specie_perc = transform_sample_counts(Specie_comp, function(x) x*100)
Specie_perc
Specie_1up = filter_taxa(Specie_perc, function(x) sum(x) > 1, TRUE)
Specie_1up
#30 y 10
Specie_top = aggregate_top_taxa2(Specie_perc,top = 54, level = "Specie")
Specie_top
data = psmelt(Specie_top) # create dataframe from phyloseq object
data$Specie <- as.character(data$Specie) #convert to character
data$Specie[data$Specie == "Other"] <- "< 1% overall rel. ab."
bars_Specie_top <- ggplot(data, aes(x=Sample, y=Abundance, fill=Specie)) + 
  geom_bar(aes(), stat="identity", position="stack", color="black") +
  facet_nested(. ~ Origin + Platform, scales = "free_x",space="free_x") +
  labs(x = "", y = "Abundancia relativa (%)") +
  scale_fill_manual(values=mycolors) + guides(x = guide_axis(angle = 90)) + 
  theme(panel.spacing=unit(0.1,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        panel.border=element_rect(color="grey90", fill = NA),
        axis.ticks.x=element_blank()) +
  scale_y_continuous(expand = c(0.01, 0.01))
bars_Specie_top

#genus###################
ps_16s_genus
genus_df <- data.frame(phyloseq::otu_table(ps_16s_genus))
#write.table(genus_transpose,file="genus_nmds.txt", sep = "\t", row.names = TRUE, col.names = NA)
genus_transpose_m <- as.data.frame(t(as.matrix(genus_df)))
#write.table(genus_transpose,file="genus_nmds_transpose.txt", sep = "\t", row.names = TRUE, col.names = NA)

genus_nmds <- metaMDS(genus_transpose_m, k=2, trymax=100, sfgrmin = 1e-9)
genus_nmds

stressplot(genus_nmds) # Produces a Shepards diagram
plot(genus_nmds)
genus_nmds
ordiplot(genus_nmds,type="n")
orditorp(genus_nmds,display="sites",cex=1.25,air=0.01)