###################Importation of qiime objects using Silva 138######################################

#### metadata
metadata<- (read_q2metadata("JFFS_metadata_qiime2.tsv"))
metadata

#### feature table
ASVs<-read_qza("qiime/merged/merged-table.qza")
names(ASVs)
ASVs$data

#### Tree (rooted)
tree<-read_qza("qiime/merged/rooted-tree.qza")
tree$data 

#### Taxonomy
taxonomy138<-read_qza("qiime/merged/taxonomy_silva138.qza")
taxonomy132<-read_qza("qiime/merged/taxonomy_silva132.qza")

#tax table extraction for Silva 138
taxtable138 <- taxonomy138$data %>% 
  as_tibble() %>% 
  separate (Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))

taxtable138
# # A tibble: 595 x 9
#    Feature.ID        Kingdom   Phylum      Class        Order      Family      Genus      Species       Confidence
#    <fct>             <chr>     <chr>       <chr>        <chr>      <chr>       <chr>      <chr>              <dbl>
#  1 00b189dfeb8d57ac… d__Bacte… " p__Fusob… " c__Fusoba… " o__Fuso… " f__Fusob… " g__Fuso… " s__gut_met…      0.995
#  2 00f22f5da05122ba… d__Bacte… " p__Firmi… " c__Clostr… " o__Osci… " f__Rumin… " g__UBA1… " s__uncultu…      0.988
#  3 01ef68e4dfd8711c… d__Bacte… " p__Firmi… " c__Bacill… " o__Lact… " f__Strep… " g__Stre… " s__Strepto…      0.781
#  4 020bc3d0d8004d13… d__Bacte… " p__Firmi… " c__Clostr… " o__Pept… " f__Pepto… " g__Pept… " s__uncultu…      1.00 
#  5 036a90f816e0e135… d__Bacte… " p__Firmi… " c__Clostr… " o__Osci… " f__Oscil… " g__NK4A… "Ambiguous_t…      0.748
#  6 03f74c0ea1f06547… d__Bacte… " p__Prote… " c__Gammap… " o__Burk… " f__Sutte… " g__Sutt… " s__uncultu…      0.738
#  7 045b30d05d01cc5a… d__Bacte… " p__Prote… " c__Gammap… " o__Burk… " f__Sutte… " g__Sutt…  NA                0.971
#  8 0575f73c4fd1dcf7… d__Bacte… " p__Prote… " c__Gammap… " o__Aero… " f__Succi… " g__Anae…  NA                0.992
#  9 05b0d51e42b60014… d__Bacte… " p__Firmi… " c__Clostr… " o__Lach… " f__Lachn… " g__Marv…  NA                0.938
# 10 05b4471632cd64da… d__Bacte… " p__Bacte… " c__Bacter… " o__Bact… " f__Riken… " g__Alis… "Ambiguous_t…      0.859

str(taxtable138)
# Classes ‘tbl_df’, ‘tbl’ and 'data.frame':	595 obs. of  9 variables:
#  $ Feature.ID: Factor w/ 595 levels "00b189dfeb8d57ac51adc17b62b8bdb9",..: 1 2 3 4 5 6 7 8 9 10 ...
#  $ Kingdom   : chr  "d__Bacteria" "d__Bacteria" "d__Bacteria" "d__Bacteria" ...
#  $ Phylum    : chr  " p__Fusobacteriota" " p__Firmicutes" " p__Firmicutes" " p__Firmicutes" ...
#  $ Class     : chr  " c__Fusobacteriia" " c__Clostridia" " c__Bacilli" " c__Clostridia" ...
#  $ Order     : chr  " o__Fusobacteriales" " o__Oscillospirales" " o__Lactobacillales" " o__Peptococcales" ...
#  $ Family    : chr  " f__Fusobacteriaceae" " f__Ruminococcaceae" " f__Streptococcaceae" " f__Peptococcaceae" ...
#  $ Genus     : chr  " g__Fusobacterium" " g__UBA1819" " g__Streptococcus" " g__Peptococcus" ...
#  $ Species   : chr  " s__gut_metagenome" " s__uncultured_organism" " s__Streptococcus_pyogenes" " s__uncultured_bacterium" ...
#  $ Confidence: num  0.995 0.988 0.781 1 0.748 ...

#### constructing phyloseq object with phyloseq()
phy138 <- phyloseq(otu_table(ASVs$data, taxa_are_rows = T),
                phy_tree(tree$data),
                tax_table( taxtable138 %>%
                             select(-Confidence) %>%
                             column_to_rownames(var = "Feature.ID") %>%
                             as.matrix()),
                sample_data(metadata %>% 
                             remove_rownames() %>%
                             column_to_rownames(var = "SampleID"))
                )
phy138

phy138
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 595 taxa and 57 samples ]
# sample_data() Sample Data:       [ 57 samples by 3 sample variables ]
# tax_table()   Taxonomy Table:    [ 595 taxa by 7 taxonomic ranks ]
# phy_tree()    Phylogenetic Tree: [ 595 tips and 593 internal nodes ]              


#### direct way with qiime2R
phy2R138<-qza_to_phyloseq("qiime/merged/merged-table.qza",
 "qiime/merged/rooted-tree.qza", 
 "qiime/merged/taxonomy_silva138.qza",
 "JFFS_metadata_qiime2.tsv")

phy2R138
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 595 taxa and 57 samples ]
# sample_data() Sample Data:       [ 57 samples by 3 sample variables ]
# tax_table()   Taxonomy Table:    [ 595 taxa by 7 taxonomic ranks ]
# phy_tree()    Phylogenetic Tree: [ 595 tips and 593 internal nodes ]

#So far they look the same, however:
all.equal(phy138, phy2R138)
# [1] "Attributes: < Component “tax_table”: 'is.NA' value mismatch: 434 in current 378 in target >"

#lets look at NAs directly in taxtable138
sum(is.na(taxtable138))
sum(is.na(taxtable138))
#[1] 378