Load libraries, data and metadata

################## 
# FASTA FUNCTION #
##################
write.fasta <- function(df, file.out, open = "w"){
  outfile <- file(description = file.out, open = "w")
  write.oneseq <- function(sequence, name){
    writeLines(paste(">", name, sep = ""), outfile)
    writeLines(sequence, outfile)
    return(name)
  }
  map2_chr(df$Sequence, df$Seq_ID, ~write.oneseq(.x, .y))
  close(outfile)
}

Reshape seqtab

seqtab.nochim.reshaped <- map(seqtabs, ~.x %>%
                            rownames_to_column(var = "Sample_ID") %>%
                            pivot_longer(., -Sample_ID) %>%
                            dplyr::rename(Sequence="name", Reads="value")
)
seqtab.taxa <- seqtab.nochim.reshaped %>%
                map2(., taxtabs, ~left_join(.x, .y, by="Sequence"))

Creat Short ID For Sequences across runs

Create Query files

# Identify sequences not assigned by RDP at species level
unnasigned_0 <- seqtab.taxa.all %>%
  filter(., is.na(Species)) #%>%
  #filter(., Reads != 0) 

# Create Query files
write.fasta(Seq_ID_key, paste0(result_dir,"BVAB_Query.fasta"))
write.fasta(unnasigned_0, paste0(result_dir,"optivag_Query.fasta"))

BLAST

BLAST seraches

# Make uniqe identifyers in the BVAB file
fa = readDNAStringSet(paste0(db_dir,"BVAB_rRNA_database.fa"))
ID <-  make.unique(paste0(str_match(names(fa), "^BVAB._"),"0"), sep = "")
names(fa) = paste(ID, names(fa))
writeXStringSet(fa, paste0(db_dir,"BVAB_db.fasta"))

for(name in c("BVAB", "optivag")) {
  if (isFALSE(is_empty(list.files(path = db_dir, pattern = paste0(name,".+\\.nhr"))))) {
    message("using database provided in path:")
    message(paste0(db_dir, name, "_reference.db"))
  } else {
    system2(command = "makeblastdb", 
                          args = c("-in", paste0(db_dir, name, "_db.fasta"), 
                                  "-out", paste0(db_dir, name, "_reference.db"), 
                                  "-parse_seqids", 
                                  "-dbtype", "nucl"),
                         wait = TRUE,
                         stdout = TRUE)
  }
  blast_out <- system2(command = "blastn", 
                       args = c("-db", paste0(db_dir, name, "_reference.db"), 
                                "-query", paste0(result_dir, name, "_Query.fasta"), 
                                #"-out", paste0(result_dir, name, "_hits.tab"),
                                "-outfmt", 6, 
                                "-perc_identity", 99.5,
                                "-qcov_hsp_perc", 90),
                       wait = TRUE,
                       stdout = TRUE) %>%
  as_tibble() %>% 
  separate(col = value, sep = "\t", convert = TRUE,
           into = c("qseqid","sseqid","pident","length","mismatch","gapopen",
                    "qstart","qend","sstart","send","evalue","bitscore") )
  
  write_tsv(blast_out, paste0(result_dir, name, "_hits.tsv"))
}

Load BLAST results

blast_OptiVag <- read_table2(paste0(result_dir,"optivag_hits.tsv"))
blast_BVAB <- read_table2(paste0(result_dir, "BVAB_hits.tsv"))

Retrive Genus and Species hits from BLAST search

# Resolve ambiguously asigned genus and Lacto spp.
OptiVag_hits <- blast_OptiVag %>%
                left_join(select(OptiVag, Full_taxonomy, Species, Seq_ID), 
                          by=c(sseqid='Seq_ID')) %>%
                select(Seq_ID=qseqid, Full_taxonomy, Species, Ref_iD=sseqid, 
                       Identity=pident, Mismatch=mismatch) %>%
                group_by(Seq_ID) %>% 
# removes hits with mismatch=1, for sequences with multiple hits where on of the hits has mismatch=0 
                filter(!(n()>1 & any(grepl('0', Mismatch)) & Mismatch == 1)) %>%
                ungroup() %>%
                separate(Species, into=c('Genus', 'Species'), sep = '_') %>% 
                separate(Full_taxonomy, into=c("Full_taxonomy", NA), sep = '(;)(?:.(?!;))+$') %>%
                select(Seq_ID, Full_taxonomy, Genus, Species) %>%
                unique() %>%
# the last group is dropped after summarize (.groups argument)
                group_by(Seq_ID, Full_taxonomy, Genus) %>% 
                summarise(Species = paste(Species, collapse="/"), 
                          .groups='drop_last') %>%
                mutate(Genus = paste(Genus, collapse="/")) %>%
                left_join(select(unnasigned_0, Seq_ID, Reads_per_taxa), 
                          by='Seq_ID') %>%
                ungroup() 

BVAB_hits <- blast_BVAB %>%
              select(qseqid, sseqid) %>%
              left_join(select(seqtab.taxa.all, Seq_ID, Reads_per_taxa),
                          by=c('qseqid'='Seq_ID')) %>%
              separate(sseqid, c('Genus', 'Species')) %>%
              dplyr::rename(Seq_ID=qseqid) %>%
              select(Seq_ID, Genus, Reads_per_taxa) %>%
              unique() 

all_hits <- bind_rows(BVAB_hits, OptiVag_hits) %>%
              unique()

BLAST result info

print(paste0("Number of sequences annotated by OptiVag:  ", length(unique(OptiVag_hits$Seq_ID))))

print(paste0("Number of sequences annotated by BVAB:     ", length(unique(BVAB_hits$Seq_ID))))
 
overlap <- intersect(BVAB_hits$Seq_ID, OptiVag_hits$Seq_ID)
ifelse(length(unique(all_hits$Seq_ID)) == length(all_hits$Seq_ID), "There is no duplicate species annotation", paste0(message("Resolving overlapping annotation between BVAB and OptiVag:"), paste(overlap,  collapse=', ')))

doNextChunk <- isFALSE(length(unique(all_hits$Seq_ID)) == length(all_hits$Seq_ID))

## [1] "Number of sequences annotated by OptiVag:  416"
## [1] "Number of sequences annotated by BVAB:     24"
## [1] "Seq_00601, Seq_00693, Seq_00694, Seq_01193, Seq_01198, Seq_01202, Seq_01203, Seq_01212, Seq_01214, Seq_01216, Seq_02430, Seq_02442, Seq_02728, Seq_02729, Seq_02732, Seq_02737, Seq_02740, Seq_02741, Seq_02745, Seq_02746, Seq_02773, Seq_02774, Seq_02776, Seq_02791"

Resolve Taxonomy

Resolve overlapping annotation between BVAB and OptiVag

# Identify which annotation that overlap for BVAB sequences 
duplicate <- all_hits %>%
  group_by(Seq_ID) %>% 
  filter(n()>1) %>%
  filter((grepl('BVAB', Genus))) %>%
  ungroup()

# Filter out duplicated annotation for BVAB
filtered_hits <- all_hits %>%
  filter(!(Seq_ID %in% duplicate$Seq_ID & !(grepl('BVAB', Genus)))) 

# is there duplicates other than BVAB?
d <- filtered_hits %>% group_by(Seq_ID) %>% summarize(n(), .groups='drop_last') %>% filter(`n()`>1)  %>% .$Seq_ID

# prints message about duplicate status
ifelse(length(unique(filtered_hits$Seq_ID)) == length(filtered_hits$Seq_ID), "Overlaping annotation resolved", paste0(message("Resolving duplicate taxonomy asignment:"), paste(d,  collapse=', ')))

# If duplicates are identified, run next chunk
doNextChunk <- isTRUE(length(d)>=1)

## [1] "Seq_00354, Seq_00355, Seq_00359, Seq_00504, Seq_01449, Seq_01575, Seq_01781, Seq_01783, Seq_01786, Seq_01787, Seq_01788, Seq_01791, Seq_01792, Seq_02265, Seq_02434"

Resolve duplicate taxonomy

# Collapse species annotation for sequences with several hits
filtered_2 <- aggregate(. ~Seq_ID, data = filtered_hits, na.action=NULL, FUN=function(x) { paste(unique(x), collapse = '/') }) 

# Sequence ids that are duplicate:
dd <- filtered_2 %>% group_by(Seq_ID) %>% summarize(n(), .groups='drop_last') %>% filter(`n()`>1)  %>% .$Seq_ID

# Message to terminal
ifelse(length(unique(filtered_hits$Seq_ID)) == length(filtered_2$Seq_ID), "Duplicate annotation resolved", paste0(message("WARNING! still duplicate annotation:"), paste(dd,  collapse=', ')))

filtered_hits <- filtered_2

## [1] "Duplicate annotation resolved"

Merge New Annotation With Tax Table

#### Update taxa ####
update_taxa <- seqtab.taxa.all %>%
  left_join(., select(filtered_hits, Seq_ID, Genus_u="Genus",
                      Species_u="Species", Full_taxonomy), by='Seq_ID') %>%
  mutate(identical = ifelse(str_detect(.$Genus_u, "BVAB")|
                                     .$Genus_u == .$Genus|
                                           is.na(.$Genus), "yes", "no"))

new_taxtable <- update_taxa %>%
  mutate(Genus_f = ifelse(.$identical=="no"|is.na(.$identical), 
                          .$Genus, .$Genus_u)) %>%
  mutate(Species_f = ifelse(.$identical=="no"|is.na(.$identical), 
                            .$Species, .$Species_u)) %>%
  select(-Genus, -Species, -Genus_u, -Species_u, -identical, -Sequencing_run) %>%
  dplyr::rename(Genus="Genus_f", Species="Species_f") %>%
  select(Seq_ID, everything(), -Full_taxonomy, -Reads_per_taxa)

write.csv(new_taxtable, file=paste0(result_dir,"SeqID_to_taxa.csv"), quote=FALSE, row.names=FALSE)

Report of unasigned sequences

#### Taxa that still have no Genus or Lacto spp. annotation ####
unnasigned_final <- new_taxtable %>% 
  left_join(select(seqtab.taxa.all, Sequence, Reads_per_taxa, Sequencing_run), by="Sequence") %>%
  filter(is.na(.$Genus)) %>%
  arrange(desc(Reads_per_taxa)) %>%
  arrange(desc(Family)) %>%
  filter(Reads_per_taxa >100)

write.csv(unnasigned_final, 
          file=paste0(result_dir,"Unnasigned_tax_report.csv"), quote=FALSE, row.names=FALSE)

unnasigned_final

## # A tibble: 39 x 11
##    Seq_ID Sequence Kingdom Phylum Class Order Family Genus Species
##    <chr>  <chr>    <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>  
##  1 Seq_0… GACAGAG… Bacter… Cyano… Chlo… Chlo… Strep… <NA>  <NA>   
##  2 Seq_0… GCGAGCG… Bacter… Firmi… Clos… Clos… Rumin… <NA>  <NA>   
##  3 Seq_0… TACGGAA… Bacter… Bacte… Bact… Bact… Prevo… <NA>  <NA>   
##  4 Seq_0… GCGAGCG… Bacter… Bacte… Bact… Bact… Porph… <NA>  <NA>   
##  5 Seq_0… GCGAGCG… Bacter… Bacte… Bact… Bact… Porph… <NA>  <NA>   
##  6 Seq_0… GCTAGCG… Bacter… Prote… Alph… Rhiz… Phyll… <NA>  <NA>   
##  7 Seq_0… GCAAGCG… Bacter… Firmi… Clos… Clos… Lachn… <NA>  <NA>   
##  8 Seq_0… GCAAGCG… Bacter… Firmi… Clos… Clos… Lachn… <NA>  <NA>   
##  9 Seq_0… GCAAGCG… Bacter… Firmi… Clos… Clos… Lachn… <NA>  <NA>   
## 10 Seq_0… TACGTAG… Bacter… Firmi… Clos… Clos… Lachn… <NA>  <NA>   
## # … with 29 more rows, and 2 more variables: Reads_per_taxa <int>,
## #   Sequencing_run <chr>

Save results

merge seqtab with new taxonomy

# Add sample info and updated taxonomy 
ASV_wCtrl <- seqtab.nochim.reshaped %>% 
            map2(., sample_info, 
                 ~left_join(.x, select(.y, SampleID, ID, Sample_type, Visit, Place_sequenced,
                                            Description, Include, Include_final),
                            by=c("Sample_ID"="SampleID")) ) %>%
                  {. ->> ASV } %>%
            map(., ~ .x %>%
                  select(Sequence, Sample_ID, Reads) %>%
                  pivot_wider(., names_from=Sample_ID, values_from=Reads) %>%
                  left_join(., new_taxtable, by="Sequence")
                  )

#### ASV with all controls and replicate samples removed ####
# ASV_names <- c("ASV_CVL_V2_B1", "ASV_CVL_V3_B1", "ASV_tissue_V3_B1", 
#                "ASV_CVL_V2_B2", "ASV_tissue_V3_B2" )
ASV_names <- c("ASV_CVL_V3_B1", "ASV_tissue_V3_B1", "ASV_tissue_V3_B2" )
ASV_list <- map(ASV, ~ .x %>%
                  filter(grepl("yes", .$Include)) %>%
                  group_split(Sample_type, Visit, .keep = T)
            ) %>%
            flatten() %>%
            set_names(ASV_names) %>%
            map(., ~ .x %>%
                select(Sequence, ID, Reads) %>%
                pivot_wider(., ,names_from=ID, values_from=Reads) %>%
                left_join(., new_taxtable, by="Sequence") 
                )

Save csv

imap(ASV_list, 
     ~write.csv(.x, file=paste0(result_dir,.y,".csv"), quote=FALSE, row.names=F) )

saveRDS(ASV, paste0(result_dir,"ASV_wCtrl_long_format.RDS"))

## $ASV_CVL_V3_B1
## NULL
## 
## $ASV_tissue_V3_B1
## NULL
## 
## $ASV_tissue_V3_B2
## NULL

Improve Taxonomy Assignment

23 Februari, 2022

Load libraries, data and metadata

Reshape seqtab

Creat Short ID For Sequences across runs

Create Query files

BLAST

BLAST seraches

Load BLAST results

Retrive Genus and Species hits from BLAST search

BLAST result info

Resolve Taxonomy

Resolve overlapping annotation between BVAB and OptiVag

Resolve duplicate taxonomy

Merge New Annotation With Tax Table

Report of unasigned sequences

Save results

merge seqtab with new taxonomy

Save csv