filter the files for analysis

library(tidyverse)
library(xlsx)
library(rotemfuns)
library(vegan)
library(funrar)
library(subSeq)

dbmap = read.csv("/pita/pub/data/16S_DBs/maps/DB1-12_premap_v2.csv")
dbmap$reads_number = dbmap$reads_number %>% as.character() %>% as.numeric()

NAs introduced by coercion

done = 
dbmap %>% 
  filter(Cohort == "eyes") %>% 
  # select(-Family) %>% 
  filter(Organism == "Human") %>% 
  arrange(desc(reads_number)) %>% 
  distinct(sample_ID, .keep_all = T)

prev_families = read.csv("DB3_control_families.csv")    

done = 
done %>% 
  full_join(dbmap %>% filter(X.SampleID %in% prev_families$X.SampleID))

Joining, by = c("X.SampleID", "BarcodeSequence", "LinkerPrimerSequence", "Barcode", "Database", "Cohort", "sample_ID", "Country", "Organism", "Dx", "FFPE_fresh", "Source", "Enzyme", "remark", "remarks2", "DUPLICATE_code", "Date", "Swab", "pn_ID", "Time", "Time2", "amount", "Age", "Gender", "Family", "position", "box", "Description", "reads_number")

done = filter(done, reads_number > 1000)
done$sample_ID = gsub("-", "_", done$sample_ID %>% as.character())
done = distinct(done, sample_ID, .keep_all = T)
names(done)[1] = "sampleid"
# write.table(done, "human/16s_map.tsv", row.names = F, sep = "\t", quote = F)

# Old
# map = read.xlsx("human/Table for rotem- Families 22.6.20.xlsx", 1)
map = read.xlsx("Table for rotem- Families RH_6Aug20.xlsx",1)
map$Sample.ID = str_replace_all(map$Sample.ID, "-", replacement = "_")
# map$Sample.ID = paste0("E_", map$Sample.ID)
done$sample_ID = gsub("E_", "", done$sample_ID)
map$pnid = paste0(map$Family,"-", map$Subject_.)
# Verify that all samples that already run fount Hadas table
setdiff(done$sample_ID , map$Sample.ID)

 [1] "HY_fresh" "HR_fresh" "HP"       "HD"       "DA_fresh" "CU_fresh" "CR"      
 [8] "Cya"      "Cyo"      "FS"       "FO"       "CK"       "CO"       "CS"      
[15] "BM"       "BS"       "NE"       "NN"       "NS"      

intersect(done$sample_ID , map$Sample.ID) %>% length() == length(done$sample_ID)

[1] FALSE

names(map)[1] = "sample_ID"

#  I made some manual cahnges for this file, including: unify age, family pnID
# done %>%
#   select(-Dx, -FFPE_fresh, -Source, -Enzyme, -remark, remarks2) %>% 
#   left_join(map %>% select(-Family), by = 'sample_ID') %>% 
#   write.table("16s_map.tsv", row.names = F, sep = "\t", quote = F)

Samples cross QC (Reads above 2000)

data = read.csv("16s_map.tsv",sep = "\t")
sum(data$reads_number > 2000)

[1] 183

Total number of families that passed QC

data %>% 
  filter(reads_number > 2000) %>% 
  select(Family) %>% 
  unique() %>% 
  nrow()

[1] 26

PCoA not included the families from old DB3

pca = read_core_pcoa("res/core-metrics-results/unweighted_unifrac_pcoa_results.qza")
names(pca) = paste0("PCoA_", 1:length(names(pca)))

data = 
data %>% 
  left_join(pca %>% rownames_to_column("sampleid")) 

Joining, by = "sampleid"Column `sampleid` joining factor and character vector, coercing into character vector

pl = 
ggplot(data, aes(x = PCoA_1, y = PCoA_2, color = Age, name = sample_ID)) + 
  geom_point() + 
  theme_rh
plotly::ggplotly(pl)

pl = 
ggplot(data, aes(x = PCoA_1, y = PCoA_2, color = Family, name = sample_ID)) + 
  geom_point() + 
  theme_rh
plotly::ggplotly(pl)

pl = 
ggplot(data, aes(x = PCoA_1, y = PCoA_2, color = Gender, name = sample_ID)) + 
  geom_point() + 
  theme_rh
plotly::ggplotly(pl)

biom_path = "/pita/pub/data/16S_DBs/processed/16S_DB1-12_merged/feature-table.biom"
biom = read_biom(biom_path)
samples_to_filter = done$sampleid

colnums = match(as.character(samples_to_filter), names(biom)) %>% na.omit()
biom = 
  biom %>% 
    select(1, colnums) %>% 
    column_to_rownames(names(biom)[1])

biom =
  biom %>% 
  distinct(.keep_all = T)

biom %>% dim()

[1] 2660  185