---
# System prepended metadata

title: New data 16S rRNA - CRC WG3 cost
tags: [COST Action, 16S rRNA, ML]

---

# New data 16S rRNA - CRC WG3 cost 

###### tags: `16S rRNA` `COST Action` `ML`

### Qiime2 objects to parse to phyloseq 

**Metadata** 

[Metadata all studies](https://drive.google.com/file/d/1Koz4Ft7amQdfc_c8jC4HXeVtVq6oHc58/view?usp=share_link)

**Single-end data**

[Taxonomy table.](https://drive.google.com/file/d/1T_HvmChQFamt7spqmcW8nAsloDTm3mMh/view?usp=share_link)
[ASV table.](https://drive.google.com/file/d/1lFLhJXcKfMAa_yUG__22J3_phQHac2Jr/view?usp=share_link) 

[Stats.](https://drive.google.com/file/d/14hKnLxw1c5wAnnH-IkuXNB5cD3esmmPG/view?usp=share_link) You can visualize [here.](https://view.qiime2.org/)
stats table.

**Paired-end data** 
[Taxonomy table.](https://drive.google.com/file/d/1NX9-91EkOtvBuH1b1W64HzSCZzKDgzix/view?usp=share_link)
[ASV table. ](https://drive.google.com/file/d/193wXGAt7PDyQ4IZ6-YnZM3gQj_fOerSA/view?usp=share_link)

[Stats.](https://drive.google.com/file/d/1zcH4pDoT-pEPy8QMfe9creVMooVUcXaD/view?usp=share_link) you can visualize [here. ](https://view.qiime2.org/)

**My code**

```R=1
library(tidyverse)
library(qiime2R)
library(phyloseq)
library(microbiome)

#ASV table####
ASV<-read_qza("DATA/table.qza")
#Taxonomy table
taxonomyq <-read_qza("DATA/taxonomy.qza")
#Tax table transformations
taxtableq <- taxonomyq$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))
taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom)
taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum)
taxtableq$Class <- gsub("c__", "",taxtableq$Class)
taxtableq$Order <- gsub("o__", "",taxtableq$Order)
taxtableq$Family <- gsub("f__", "",taxtableq$Family)
taxtableq$Genus <- gsub("g__", "",taxtableq$Genus)
taxtableq$Species <- gsub("s__", "",taxtableq$Species)

#Clean taxa: Remove Kingdom Unassigned & Eukaryota
taxtableq <- taxtableq[grep("Eukaryota|Unassigned", taxtableq$Kingdom, invert=T),]

#Tree
tree<-read_qza("DATA/rooted-tree.qza")

#Metadata
#Here you have to create a csv file with the metadata
metadata <-read_csv("metadata.csv")

#Construct phyloseq object
ASVs <- otu_table(ASV$data, taxa_are_rows = T) 
tree <- phy_tree(tree$data) 
TAXq <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it
sample_metadata <- sample_data(metadata %>% as.data.frame()) 
sample_names(sample_metadata) <- paste(metadata$SampleID)
physeq<- merge_phyloseq(ASVs, tree, TAXq,sample_metadata)

#Remove samples with no diagnosis - Some extra samples were inlcuded by mistake-  
physeq <- subset_samples(physeq, Diagnosis!="NA")

#Agglomerate at genus level 
physeq.gen <- tax_glom(physeq, taxrank="Genus")

#Write data for sharing
data <- psmelt(physeq.gen) 
ASV_tab <- data %>% select(SampleID,Abundance,Genus) 
genus <- pivot_wider(ASV_tab, id_cols=SampleID, names_from=Genus, values_from=Abundance, values_fn = mean)
write_csv(genus,"genus.csv")

#Save phyloseq object for sharing
saveRDS(physeq.gen, "physeq.RDS")
```

:female-construction-worker: **Working Plan:**

1. Merge the new studies with the other studies in a phyloseq object. 
2. Maybe we need to use some strategy to analyze all studies togeteher an reduce batch-effect. I have used [Mmuphin](https://www.bioconductor.org/packages/release/bioc/html/MMUPHin.html) in other cases. Other tools are welcome :wink: 
3. Descriptive analysys of the data, including: 
    - Analytics of metadata
    - Alfa and beta diversities of the microbiome 
    - Whatever you think would be useful :P 

:::info 
Regarding classification, this is a topic that we have discussed (with little succes) in some meetings in the WG. I think somebody suggest just to use Healthy / Cancer categories. We have an issue here because Adenoma is something inbetween. So, first with the descriptive analysis we will see the "dimension" of the problem i.e. how many samples are in each category, and then we can decide better on this. However, I will try to ask a collegue that works with CRC genomics to know her opinion.  
But the "easy" solution will be to do two categories: Healthy(Control)/Cancer
:::


## Solving Point 1: 

```r=1
library(tidyverse)
library(qiime2R)
library(phyloseq)
library(microbiome)
library(ape)

#Single
#ASV table####
ASV<-read_qza("DATA/table-single.qza")
#Taxonomy table
taxonomyq <-read_qza("DATA/taxonomy-single.qza")
#Tax table transformations
taxtableq <- taxonomyq$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))
taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom)
taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum)
taxtableq$Class <- gsub("c__", "",taxtableq$Class)
taxtableq$Order <- gsub("o__", "",taxtableq$Order)
taxtableq$Family <- gsub("f__", "",taxtableq$Family)
taxtableq$Genus <- gsub("g__", "",taxtableq$Genus)
taxtableq$Species <- gsub("s__", "",taxtableq$Species)

#Clean taxa: Remove Kingdom Unassigned & Eukaryota
taxtableq <- taxtableq[grep("Eukaryota|Unassigned", taxtableq$Kingdom, invert=T),]

#Metadata
#Here you have to create a csv file with the metadata
metadata <-read_tsv("DATA/metadata_single.csv")

#Construct phyloseq object
ASVs <- otu_table(ASV$data, taxa_are_rows = T) 
#tree <- phy_tree(tree$data) 
TAXq <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it
sample_metadata <- sample_data(metadata %>% as.data.frame()) 
sample_names(sample_metadata) <- paste(metadata$SampleID)
sample_names(ASVs) <- paste(metadata$SampleID)
physeq.single <- merge_phyloseq(ASVs, TAXq,sample_metadata)


#Paired-end
#ASV table####
ASV<-read_qza("DATA/table-paired.qza")
#Taxonomy table
taxonomyq <-read_qza("DATA/taxonomy-paired.qza")
#Tax table transformations
taxtableq <- taxonomyq$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))
taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom)
taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum)
taxtableq$Class <- gsub("c__", "",taxtableq$Class)
taxtableq$Order <- gsub("o__", "",taxtableq$Order)
taxtableq$Family <- gsub("f__", "",taxtableq$Family)
taxtableq$Genus <- gsub("g__", "",taxtableq$Genus)
taxtableq$Species <- gsub("s__", "",taxtableq$Species)

#Clean taxa: Remove Kingdom Unassigned & Eukaryota
taxtableq <- taxtableq[grep("Eukaryota|Unassigned", taxtableq$Kingdom, invert=T),]

#Construct phyloseq object
ASVp <- otu_table(ASV$data, taxa_are_rows = T) 
TAXp <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it

#Metadata
names <- as.data.frame(sample_names(ASVp))
colnames(names) <- "SampleID"
metadataALL <- read_csv("DATA/metadataALL.csv")
metadata_paired <- metadataALL %>% right_join(names)

sample_metadatap <- sample_data(metadata_paired %>% as.data.frame()) 
sample_names(sample_metadatap) <- paste(metadata_paired$SampleID)
sample_names(ASVp) <- paste(metadata_paired$SampleID)
physeq.paired <- merge_phyloseq(ASVp, TAXp,sample_metadatap)


#Previous datasets #####
#ASV table####
ASV<-read_qza("DATA/table.qza")
#Taxonomy table
taxonomyq <-read_qza("DATA/taxonomy.qza")
#Tax table transformations
taxtableq <- taxonomyq$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))
taxtableq$Kingdom <- gsub("d__", "",taxtableq$Kingdom)
taxtableq$Phylum <- gsub("p__", "",taxtableq$Phylum)
taxtableq$Class <- gsub("c__", "",taxtableq$Class)
taxtableq$Order <- gsub("o__", "",taxtableq$Order)
taxtableq$Family <- gsub("f__", "",taxtableq$Family)
taxtableq$Genus <- gsub("g__", "",taxtableq$Genus)
taxtableq$Species <- gsub("s__", "",taxtableq$Species)

#Clean taxa: Remove Kingdom Unassigned & Eukaryota
taxtableq <- taxtableq[grep("Eukaryota|Unassigned", taxtableq$Kingdom, invert=T),]


#Metadata
metadata<-read_csv("DATA/metadata.csv")

#Construct phyloseq object
ASVpr <- otu_table(ASV$data, taxa_are_rows = T) 
TAXpr <- tax_table(as.data.frame(taxtableq) %>% select_("-Confidence") %>% column_to_rownames("Feature.ID") %>% as.matrix()) #moving the taxonomy to the way phyloseq wants it
sample_metadata <- sample_data(metadata %>% as.data.frame()) 
sample_names(sample_metadata) <- paste(metadata$SampleID)
physeq.prev <- merge_phyloseq(ASVpr, TAXpr, sample_metadata)

#Remove samples with no diagnosis - Some extra samples were inlcuded by mistake-  
physeq.prev <- subset_samples(physeq.prev, Diagnosis!="NA")


#Unir tres objetos phyloseq
physeq <- merge_phyloseq(physeq.single, physeq.paired,physeq.prev)


#Generar arbol
random_tree <-rtree(ntaxa(physeq), rooted=TRUE, tip.label=taxa_names(physeq))
physeq.tree <- merge_phyloseq(physeq,random_tree)

#Agglomerate at genus level 
physeq.gen <- tax_glom(physeq.tree, taxrank="Genus")

#Write data for sharing
data <- psmelt(physeq.gen) 
ASV_tab <- data %>% select(SampleID,Abundance,Genus) 
genus <- pivot_wider(ASV_tab, id_cols=SampleID, names_from=Genus, values_from=Abundance, values_fn = mean)
write_csv(genus,"genusALL.csv")

#Save phyloseq object for sharing
saveRDS(physeq.gen, "physeqALL.RDS")
```

:warning: **Pending work for Rajesh:**

There is a problem with the metadata, when I added the new datasets the previous metadata was deleted. You need to create a new `sample_data` object including all the clean metadata :eyes: be sure to include the correct `sample_names` and then create a new phyloseq object with the new metadata. 

You can use the previous code to see how to do it. 


[Here ](https://drive.google.com/file/d/1bBohAatEjZIYROJY7wFTM1rynMrDRMej/view?usp=share_link)is the phyloseq object with all the datasets and the BAD metadata. :cry: 

