new code for best reference species selection

[taylorreiter]

• ignores euks
• based on abundance

library(readr)
library(dplyr)

gtdb_lineages <- read_csv(snakemake@input[['gtdb_lineages']])
#gtdb_lineages <- read_csv("inputs/gtdb-rs207/gtdb-rs207.taxonomy.csv.gz")

#gather <- read_csv("outputs/rnaseq_sourmash_gather/SRR11413720_gtdb_k31.csv") %>%
gather <- read_csv(snakemake@input[['gather']]) %>%
  mutate(genome_accession = gsub(" .*", "", name)) %>%   # get accession of match
  left_join(gtdb_lineages, by = c("genome_accession" = "ident")) %>% # join to GTDB taxonomy lineages
  filter(!is.na(superkingdom)) %>% # filter to only GTDB results (removes eukaryote matches; all euks should be contamination)
  slice_max(f_unique_weighted) %>% # select reference species by GTDB match with highest abundance-weighted unique overlap with query
  select(name, query_name, genome_accession, superkingdom, phylum, class, order, family, genus, species) %>%
  mutate(species_no_space = gsub(' ', '_', species)) %>%
  mutate(sra_to_ref_species = paste0(species_no_space, "-",  query_name))

write_csv(gather, snakemake@output[['sra_to_ref_species']])