space_genomes_apit.Rmd

---
title: "space_genomes_apit"
author: "Braden T Tierney"
date: "9/27/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(ggtree)
library(tidyverse)
library(hrbrthemes)
library(tm)
library(proustr)
library(ggrepel)
library(ggnewscale)
library(VennDiagram)
library(reshape2)
library(RColorBrewer)
library(pheatmap)
library(ComplexUpset)
library(cowplot)
theme_set(theme_cowplot())

```

```{r}
#load, plot, and annotate ggtree 
setwd('~/Dropbox (HMS)/space_genomes/ap_subproject/')

tree <- read.tree("gtdbtk_trees/gtdbtk_infer_with_bins/gtdbtk.unrooted.tree")
metadata = read.table('metadata_20211015.csv',header=T,sep=',') %>% mutate(Assembly = strsplit(ASSEMBLY,'\\.') %>% map_chr(1))
rownames(metadata)=metadata$Assembly
metadata = metadata %>% select(-Assembly,ISOLATION_SOURCE_CLEANED,HUMAN_FLIGHT_OTHER)

metadata$ISOLATION_SOURCE_CLEANED = as.factor(metadata$ISOLATION_SOURCE_CLEANED)
metadata$HUMAN_FLIGHT_OTHER = as.factor(metadata$HUMAN_FLIGHT_OTHER)

new_tip_labels = list()
for(i in seq(tree$tip.label)){
   val = tree$tip.label[[i]]
   new_tip_labels[[i]] = as.character(metadata[val,'ISOLATION_SOURCE_CLEANED'])
}

new_tip_labels=unlist(unname(new_tip_labels))

labelmap = data.frame(label = tree$tip.label, newlabel = new_tip_labels)

p = ggtree(tree,branch.length="none",layout = "circular") + geom_treescale() 

p2=gheatmap(p, metadata %>% select(HUMAN_FLIGHT_OTHER), offset=0, width=.1, colnames=FALSE, legend_title="Isolation Source")+ scale_fill_viridis_d(option="D", name="General isolation Source")  

p3 = p2 %<+% labelmap + geom_tiplab(offset = 5,aes(label=newlabel),size = 2, parse = F)

ggsave('apit_tree_316.pdf',height=10,width=10)
```

```{r}
# load in and visualize mash data

### distances generated with:
# mash sketch -s 10000 -o mash_sketch bacterial_genomes_actually_pittii/*.fna
# mash dist mash_sketch.msh mash_sketch.msh > apit_mash.tsv
setwd('~/Dropbox (Mason Lab)/space_genomes/ap_subproject/')

mashdata = read.table('apit_mash.tsv',header=F,sep='\t') %>% mutate(V1 = strsplit(V1,'/') %>% map_chr(2) %>% strsplit('\\.') %>% map_chr(1))  %>% mutate(V2 = strsplit(V2,'/') %>% map_chr(2) %>% strsplit('\\.') %>% map_chr(1)) %>% select(V1,V2,V3)
colnames(mashdata) = c('genome_1','genome_2','mash_distance')
hist(mashdata$mash_distance)

mashdata = mashdata %>% reshape2::dcast(genome_1 ~ genome_2,value.var = 'mash_distance') 

mashdata = left_join(mashdata,metadata %>% rownames_to_column('genome_1') %>% select(genome_1,HUMAN_FLIGHT_OTHER)) %>% mutate(HUMAN_FLIGHT_OTHER = as.factor(HUMAN_FLIGHT_OTHER))%>%  column_to_rownames('genome_1') 

mashdata_annotation = mashdata %>% select(HUMAN_FLIGHT_OTHER)
mashdata = mashdata %>% select(-HUMAN_FLIGHT_OTHER)

save_pheatmap_pdf <- function(x, filename, width=8, height=8) {
   stopifnot(!missing(x))
   stopifnot(!missing(filename))
   pdf(filename, width=width, height=height)
   grid::grid.newpage()
   grid::grid.draw(x$gtable)
   dev.off()
}

p = (pheatmap::pheatmap(mat = -log(mashdata+0.00001,10), show_rownames = F,show_colnames = F,annotation_row = mashdata_annotation, cellwidth = 1,cellheight = 1))

save_pheatmap_pdf(p,'mash_distance_heatmap.pdf')

space_bugs =read.csv('metadata_20211015.csv') %>% filter(grepl('MT',HUMAN_FLIGHT_OTHER)) %>% select(ASSEMBLY) %>% unlist %>% unname
space_bugs = c(space_bugs,'GCA_002143515','GCA_002136755')

mashdata_sub = mashdata[space_bugs,space_bugs]

p2 = (pheatmap::pheatmap(mat = log(mashdata_sub+0.00001,10), show_rownames = F,show_colnames = T,annotation_row = mashdata_annotation, cellwidth = 10,cellheight = 10))

save_pheatmap_pdf(p2,'mash_distance_heatmap_zoomed.pdf')

# find closest earthbound relative, for use in SNP discovery 

#mashdata_melted = mashdata %>% rownames_to_column() %>% melt
#reference_bug = mashdata_melted %>% filter(rowname %in% space_bugs, !(variable %in% space_bugs)) %>% group_by(variable) %>% mutate(average_mash_dist = mean(value)) %>% select(-value,-rowname) %>% unique %>% arrange(average_mash_dist) %>% select(variable) %>% head(1) %>% unlist %>% unname %>% as.character()
#print(reference_bug)

# pick an MT1 genome for us in SNP discovery (least similar to all the others)

#mt1_bugs = metadata %>% filter(Project=='MT1') %>% rownames %>% unlist %>% unname
#reference_mt1_bug = mashdata_melted %>% filter(rowname %in% space_bugs, variable %in% mt1_bugs) %>% group_by(variable) %>% mutate(average_mash_dist = mean(value)) %>% select(-value,-rowname) %>% unique %>% arrange(desc(average_mash_dist))%>% select(variable) %>% head(1) %>% unlist %>% unname %>% as.character()
#print(reference_mt1_bug)

```

```{r}
# build upset and cog analysis

setwd('~/Dropbox (HMS)/space_genomes/ap_subproject/gene_catalogs/')

build_upset_plots <- function(clusterfile,prokkamap,metadata,iterations){
   output_raw = list()
   output_upset = list()
   output_upset_singletons = list()
   data = read.csv(clusterfile,sep='\t',header=F)
   data = data %>% mutate(val = 1,prokka_id = strsplit(V2,'_') %>% map_chr(1))
   data_prokka = left_join(data,prokkamap,by=c('prokka_id'='V2'))
   colnames(data_prokka) = c('congene','rawgene','val','prokka_id','ASSEMBLY')
   data_prokka = data_prokka %>% select(-rawgene,-prokka_id) %>% unique
   data_prokka_metadata = left_join(data_prokka,metadata %>% select(Assembly,ISOLATION_SOURCE_CLEANED),by = c('ASSEMBLY'='Assembly'))
   tokeep = table(metadata$ISOLATION_SOURCE_CLEANED) %>% data.frame() %>% filter(Freq>3) %>% filter(Var1 != 'unknown', Var1 != 'clinical material', Var1!= 'human gut') %>% select(Var1) %>% unlist %>% unname %>% as.character
   data_prokka_metadata_filtered = data_prokka_metadata %>% filter(ISOLATION_SOURCE_CLEANED %in% tokeep) 
   for(i in seq(1,iterations)){
      assemblies_to_keep = data_prokka_metadata_filtered %>% select(-val,-congene) %>% unique %>% group_by(ISOLATION_SOURCE_CLEANED) %>% do(sample_n(.,4)) %>% select(ASSEMBLY) %>% unlist %>% unname
      
      data_prokka_metadata_filtered_sub = data_prokka_metadata_filtered %>% filter(ASSEMBLY %in% assemblies_to_keep )%>% select(-ASSEMBLY) %>% unique 
      
      data_wide = dcast(data_prokka_metadata_filtered_sub,ISOLATION_SOURCE_CLEANED ~ congene,value.var = "val" )
      data_wide[is.na(data_wide)] = 0
      data_wide = t(data_wide %>% column_to_rownames('ISOLATION_SOURCE_CLEANED'))
      data_wide = data_wide %>% data.frame %>% rownames_to_column('row') %>% select(-row)
      data_wide[data_wide>0] = 1
      output_raw[[i]] = data_wide
      output_upset[[i]] = upset(data_wide,intersect = colnames(data_wide),min_size = 10)
      output_upset_singletons[[i]] = upset(data_wide  %>% filter(rowSums(data_wide)==1),intersect = colnames(data_wide),min_size = 1)
      pdf(paste('upset',clusterfile,'_',i,'.pdf',sep=''),width=24,height=8)
      print(upset(data_wide,intersect = colnames(data_wide),min_size = 10))
      dev.off()
      pdf(paste('upset_singletons_',clusterfile,'_',i,'.pdf',sep=''),width=24,height=8)
      print(upset(data_wide  %>% filter(rowSums(data_wide)==1),intersect = colnames(data_wide),min_size = 1,keep_empty_groups= T))
      dev.off()
   }
   return(list(rawdata = output_raw,upset = output_upset))
}

# load metadata
metadata = read.table('../metadata_20211015.csv',header=T,sep=',') %>% mutate(Assembly = strsplit(ASSEMBLY,'\\.') %>% map_chr(1))

# load prokka mapping
prokkamap = read.csv('prokka_id_mapping',header=F,sep='\t')

cluster30 = build_upset_plots('clusters_30perc_cluster.tsv',prokkamap,metadata,3)
cluster50 = build_upset_plots('clusters_50perc_cluster.tsv',prokkamap,metadata,3)
cluster70 = build_upset_plots('clusters_70perc_cluster.tsv',prokkamap,metadata,3)
cluster90 = build_upset_plots('clusters_90perc_cluster.tsv',prokkamap,metadata,3)

# load cog mapping
cogmap = read.csv('cog_fx_map')

# load prokka cog data
procog = read.csv('apit_cog_data.tsv',sep='\t')
procog_sub = procog %>% select(COG,all_of(c("GCA_017166365","GCA_013449825","GCA_000369045","GCA_002143515","GCA_014168695")))
colnames(procog_sub) = c('COG','MT1','MT2','TYPE STRAIN','HUMAN INFECTION','WASTEWATER')

procog_sub_melted = melt(procog_sub) %>% filter(value!=0)

colors = brewer.pal(n=5,"Set1") 

venn.diagram(
  x = list(
    procog_sub_melted %>% filter(variable=="MT1") %>% select(COG) %>% unlist(), 
    procog_sub_melted %>% filter(variable=="MT2") %>% select(COG) %>% unlist(), 
    procog_sub_melted %>% filter(variable=="TYPE STRAIN") %>% select(COG) %>% unlist(),
    procog_sub_melted %>% filter(variable=="HUMAN INFECTION") %>% select(COG) %>% unlist(),
    procog_sub_melted %>% filter(variable=="WASTEWATER") %>% select(COG) %>% unlist()
    ),
  category.names = c('','','','',''),
  filename = 'representative_cog_venn.png',
  output = TRUE ,
          imagetype="png" ,
          height = 1500 , 
          width = 1500 , 
          resolution = 600,
          compression = "lzw",
          lwd = 1,
          col=colors,
          fill = c(alpha(colors[[1]],0.3), alpha(colors[[2]],0.3), alpha(colors[[3]],0.3), alpha(colors[[4]],0.3), alpha(colors[[5]],0.3)),
          cex = 0.5,
          fontfamily = "sans",
          cat.cex = 0.3,
          cat.default.pos = "outer",
          cat.fontfamily = "sans",
          cat.col = colors,
        )
#left_join(procog_sub,cogmap) %>% select()

#GCA_017166365 -- mt2
#GCA_013449825 -- mt1
#GCA_000369045 -- type strain
#GCA_002143515 -- earth ref
#GCA_014168695 -- wastewater


```

```{r}
# fxn-was output

setwd('~/Dropbox (HMS)/space_genomes/ap_subproject/gene_catalogs/')

earth_space = read.csv('earth_space_cleaned_fxnl_output.csv') %>% mutate(BY = p.adjust(lrt.pvalue,method = 'BY')) %>% select(af,BY,variant_h2,func,VARTYPE,beta,name) %>% unique
flight = read.csv('flight_binary_cleaned_fxnl_output.csv')  %>% mutate(BY = p.adjust(lrt.pvalue,method = 'BY')) %>% select(af,BY,variant_h2,func,VARTYPE,beta,name) %>% unique

for(v in unique(earth_space$VARTYPE)){
   if(v == 'COG'){
      earth_space_sub = earth_space %>% filter(func != '') 
      p = ggplot(data = earth_space_sub %>% filter(VARTYPE == v),aes(x=beta,y=-log(BY,10),size=variant_h2,color=af)) + geom_point(aes(alpha=.5)) + ylab('-log10(adjusted p-value)') + geom_label_repel(data = earth_space_sub%>% filter(VARTYPE == v) %>% arrange(BY) %>% filter(BY<0.05 )%>% group_by(func) %>% slice(1:3),aes(label = name),box.padding   = 0.1, point.padding = 0.1,max.overlaps=100,size=2,segment.color = 'grey50')+geom_hline(yintercept = -log(0.05,5)) + ggtitle(paste(v,'-- EARTH V SPACE')) + facet_wrap(~ func, ncol = 5)
      ggsave(paste(v,'_earth-space_volcano.pdf',sep=''),width=14,height=14)
   }
   if(v != 'COG'){
      p = ggplot(data = earth_space  %>% filter(VARTYPE == v),aes(x=beta,y=-log(BY,10),color=af,size=variant_h2)) + geom_point(aes(alpha=.5)) + ylab('-log10(adjusted p-value)') + geom_label_repel(data = earth_space%>% filter(VARTYPE == v) %>% arrange(BY) %>% filter(BY<0.05 )%>% group_by(VARTYPE) %>% slice(1:15),aes(label = name),box.padding   = 0.1, point.padding = 0.1,max.overlaps=100,size=3,segment.color = 'grey50')+geom_hline(yintercept = -log(0.05,10)) + ggtitle(paste(v,'-- EARTH V SPACE')) + xlim(-1,1)# + facet_grid(cols=vars(VARTYPE),scales= 'free')
      ggsave(paste(v,'_earth-space_volcano.pdf',sep=''),width=8,height=8)
   }
}

for(v in unique(flight$VARTYPE)){
   if(v == 'COG'){
      flight_sub = flight %>% filter(func != '') 
      p = ggplot(data = flight_sub %>% filter(VARTYPE == v),aes( x=beta,y=-log(BY,10),size=variant_h2,color=af)) + geom_point(aes(alpha=.5)) + ylab('-log10(adjusted p-value)') + geom_label_repel(data = flight_sub%>% filter(VARTYPE == v) %>% arrange(BY) %>% filter(BY<0.05 )%>% group_by(func) %>% slice(1:3),aes(label = name),box.padding   = 0.1, point.padding = 0.1,max.overlaps=100,size=2,segment.color = 'grey50')+geom_hline(yintercept = -log(0.05,5)) + ggtitle(paste(v,'-- MT1 V MT2'))  + facet_wrap(~ func, ncol = 5)
      ggsave(paste(v,'_flight_volcano.pdf',sep=''),width=14,height=14)
   }
   if(v != 'COG'){
      p = ggplot(data = flight   %>% filter(VARTYPE == v),aes(x=beta,y=-log(BY,10),size = variant_h2,color=af)) + geom_point(aes(alpha=.5)) + ylab('-log10(adjusted p-value)') + geom_label_repel(data = flight%>% filter(VARTYPE == v) %>% arrange(BY) %>% filter(BY<0.05 )%>% group_by(VARTYPE) %>% slice(1:10),aes(label = name),box.padding   = 0.1, point.padding = 0.1,max.overlaps=100,size=3,segment.color = 'grey50')+geom_hline(yintercept = -log(0.05,10)) + xlim(-1,1) + ggtitle(paste(v,'-- MT1 V MT2'))# + facet_grid(cols=vars(VARTYPE),scales= 'free')
      ggsave(paste(v,'_flight_volcano.pdf',sep=''),width=8,height=8)
   }
}

```

```{r}
# card genes by bug
setwd('~/Dropbox (Mason Lab)/space_genomes/ap_subproject/card_analysis/')

create_aro_heatmap <- function(abx_col,outputdata){
   merged_output_aro =  outputdata %>% select(Sample_ID,abx_col) %>% mutate(val = 1)
   merged_output_aro_wide = dcast(merged_output_aro,Sample_ID ~ get(abx_col), value.var = 'val') %>% column_to_rownames('Sample_ID')
   merged_output_aro_wide[merged_output_aro_wide>0]=1
   tokeep = colSums(merged_output_aro_wide) %>% data.frame %>% filter(.>1) %>%rownames %>% unlist %>% unname
   merged_output_aro_wide_filt = merged_output_aro_wide %>% select(all_of(tokeep)) %>% rownames_to_column('Sample_ID')
   merged_output_aro_wide_filt_mdat = left_join(merged_output_aro_wide_filt,metadata %>% select(Assembly,HUMAN_FLIGHT_OTHER) %>% mutate(HUMAN_FLIGHT_OTHER = as.factor(HUMAN_FLIGHT_OTHER)),by = c('Sample_ID' = 'Assembly'))
   
   pheatmap(merged_output_aro_wide_filt_mdat %>% column_to_rownames('Sample_ID') %>% select(-HUMAN_FLIGHT_OTHER),annotation_row = merged_output_aro_wide_filt_mdat %>% column_to_rownames('Sample_ID') %>% select(HUMAN_FLIGHT_OTHER),show_rownames = FALSE)
}

save_pheatmap_pdf <- function(x, filename, width=20, height=12) {
   stopifnot(!missing(x))
   stopifnot(!missing(filename))
   pdf(filename, width=width, height=height)
   grid::grid.newpage()
   grid::grid.draw(x$gtable)
   dev.off()
}

metadata = read.table('../metadata_20211015.csv',header=T,sep=',') %>% mutate(Assembly = strsplit(ASSEMBLY,'\\.') %>% map_chr(1))

card_files = list.files()
card_files = card_files[grep('.fna.card.txt',card_files)]

output = list()
for(f in card_files){
  data = read.csv(f,sep='\t',header=T) %>% select(Drug.Class,Best_Hit_ARO,AMR.Gene.Family,Resistance.Mechanism)
  sample_id = strsplit(f,'\\.') %>% map_chr(1)
  data$Sample_ID = sample_id
  output[[f]] = data
}

merged_output = bind_rows(output) 

dc = create_aro_heatmap("Drug.Class",merged_output)
aro = create_aro_heatmap("Best_Hit_ARO",merged_output)
amrf = create_aro_heatmap("AMR.Gene.Family",merged_output)
rm = create_aro_heatmap("Resistance.Mechanism",merged_output)

save_pheatmap_pdf(dc,'drug_class_presence_absence.pdf',height=20)
save_pheatmap_pdf(aro,'aro_presence_absence.pdf')
save_pheatmap_pdf(amrf,'amrf_presence_absence.pdf')
save_pheatmap_pdf(rm,'rm_presence_absence.pdf')



```

```{r}
# metapan 

setwd('~/Dropbox (Mason Lab)/space_genomes/ap_subproject/metapan_bootstrapped_alignment/')

isolate_metadata = read.csv('../metadata_20211015.csv')
metagenomic_metadata = read.csv('../metagenome_metadata_update.csv') %>% mutate(METAGENOME_SOURCE = paste(COHORT,SITE,sep='-'))
sampleidmap = read.table('sample_id_mapping',sep='\t',header=F) %>% select(V1,V2) %>% mutate(V2 = strsplit(V2, ' \\./') %>% map_chr(1) %>% gsub('_R1_clumped_bbduk_hg38removed_repaired_tadpole.fastq.gz','',.) %>% gsub('\\./','',.) %>%gsub('/','-',.))
colnames(sampleidmap) = c('NEWNAME','SAMPLE_ID')
metagenomic_metadata = left_join(metagenomic_metadata,sampleidmap) %>% mutate(SAMPLE_ID = if_else(!is.na(NEWNAME),NEWNAME,SAMPLE_ID)) %>% select(-NEWNAME)

isolate_alignment_data = read.table('genome_metagenome_alignment_data_newassemblies_nojc.tsv',sep='\t',header=F) %>% mutate(V1=strsplit(V1,'\\.') %>% map_chr(1))
colnames(isolate_alignment_data) = c('ASSEMBLY','SAMPLE_ID','raw_counts','sample_read_counts','fraction_aligned')

alignment_data_metagenomic_metadata = left_join(isolate_alignment_data,metagenomic_metadata) 
alignment_data_fullmdat = left_join(alignment_data_metagenomic_metadata,isolate_metadata %>% select(ASSEMBLY,ISOLATION_SOURCE_CLEANED) %>% unique) %>% select(METAGENOME_SOURCE,ISOLATION_SOURCE_CLEANED,PMA_TREATMENT,raw_counts,fraction_aligned,FLIGHT,ASSEMBLY,SAMPLE_ID) %>% dplyr::rename(GENOME_SOURCE=ISOLATION_SOURCE_CLEANED) %>% mutate(GENOME_SOURCE=if_else(is.na(GENOME_SOURCE),'A_PITTII_CORE_ISS_GENES',toupper(GENOME_SOURCE)))%>% mutate(GENOME_SOURCE=if_else(GENOME_SOURCE == 'CEREBROSPINAL FLUID','TYPE STRAIN',toupper(GENOME_SOURCE)))

alignment_data_fullmdat$METAGENOME_SOURCE = factor(alignment_data_fullmdat$METAGENOME_SOURCE,levels=c('MT1-ENV','MT2-ENV',"MT2-CREW-SALIVA","MT2-CREW-MOUTH","MT2-CREW-NASAL","MT2-CREW-EAR","MT2-CREW-POOL","TWINS-GUT","TWINS-SALIVA"))

# load in that coverage information

reflengths = read.table('reference_genome_lengths.tsv',sep='\t')
colnames(reflengths) = c('ASSEMBLY','GENOME_LENGTH')

coverageinfo = read.table('coverage_cleaned',sep='\t')
colnames(coverageinfo) = c('ASSEMBLY','SAMPLE_ID','COVERAGE')

coverageinfo = left_join(reflengths,coverageinfo,by='ASSEMBLY')
coverageinfo = coverageinfo %>% mutate(`PERCENT COVERED` = (COVERAGE/GENOME_LENGTH)*100)

alignment_data_fullmdat = alignment_data_fullmdat %>% mutate(PMA_TREATMENT= if_else(PMA_TREATMENT == '','NO PMA',PMA_TREATMENT))

alignment_data_fullmdat = inner_join(alignment_data_fullmdat,coverageinfo)

alignment_data_fullmdat$`PERCENT COVERED BIN` <- cut(alignment_data_fullmdat$`PERCENT COVERED`,breaks = c(0,10,20,30,40,50,100))

#isolatenames <- list(
#  'ISS-MT1'="GENOMES ISOLATED FROM MT1 MISSION",
#  'MT1-BIN'="GENOME BINS FROM MT1 METAGENOMES",
#  'ISS-MT2'="GENOMES ISOLATED FROM MT2 MISSION",
#  'TYPE STRAIN'="TYPE STRAIN GENOME"
#)

#isolate_labeller <- function(variable,value){
#  return(isolatenames[value])
#}

ggplot(data = alignment_data_fullmdat %>% filter(grepl('ENV',METAGENOME_SOURCE),GENOME_SOURCE == 'MT1-BIN' ), aes(x = METAGENOME_SOURCE,fill=METAGENOME_SOURCE,shape =ASSEMBLY, y = log(fraction_aligned+0.00001))) + geom_quasirandom(aes(color=`PERCENT COVERED BIN`),size=3) + xlab('') + ggtitle('Read alignment to A. Pittii genome bins vs the A. Pittii Type Strain')+ theme(legend.position = "bottom") + stat_compare_means(comparisons = list(c('MT1-ENV','MT2-ENV'))) + scale_color_viridis_d()
ggsave('aligned_reads.pdf',width=8,height=5)

alignment_data_fullmdat_sub = alignment_data_fullmdat %>% filter(FLIGHT != -1, !is.na(FLIGHT))

my_comparisons <- list(c(3,4),c(7, 1))

ggplot(data = alignment_data_fullmdat_sub, aes(group = FLIGHT,x=FLIGHT,y = log(fraction_aligned+0.00001))) + geom_quasirandom() + geom_line() + geom_violin(alpha=.7)+ facet_grid(rows=vars(GENOME_SOURCE)) + xlab('') + ggtitle('A. pittii abundance by MT2 flight') + stat_compare_means(comparisons=my_comparisons,method = 'wilcox.test')  + scale_x_continuous("FLIGHT", labels = as.character(alignment_data_fullmdat_sub$FLIGHT), breaks = alignment_data_fullmdat_sub$FLIGHT)+ stat_summary(fun = "mean",geom = "crossbar",color = "black") 

ggsave('flight_by_flight_succession.pdf',width=8,height=6)

forhm = alignment_data_fullmdat_sub %>% select(ASSEMBLY,SAMPLE_ID,METAGENOME_SOURCE,GENOME_SOURCE,`PERCENT COVERED`,fraction_aligned)
annocols = alignment_data_fullmdat_sub %>% select(ASSEMBLY,GENOME_SOURCE) %>% distinct %>% column_to_rownames('ASSEMBLY')
forhm_wide = dcast(data = forhm %>% select(ASSEMBLY,SAMPLE_ID,fraction_aligned,METAGENOME_SOURCE),SAMPLE_ID+METAGENOME_SOURCE~ASSEMBLY ,value.var = 'fraction_aligned') %>% column_to_rownames("SAMPLE_ID")
annorows = forhm_wide %>% select(METAGENOME_SOURCE)
r = rowAnnotation(`METAGENOME SOURCE` = annorows$METAGENOME_SOURCE)
c = columnAnnotation(`A PITTII GENOME SOURCE` = annocols$GENOME_SOURCE)

forhm_wide = forhm_wide %>% select(rownames(annocols))

Heatmap(forhm_wide,left_annotation = r,top_annotation = c,show_row_names = F, show_column_names = F, cluster_columns=F,cluster_rows=F)

```


```{r}

setwd('~/Dropbox (Mason Lab)/space_genomes/ap_subproject/metapan_bootstrapped_alignment/')

data=read.csv('merged_bracken.tsv',sep='\t',header=F)
data = data %>% filter(V1 == 'Acinetobacter pittii' | V1 =='name')
colnames(data) = gsub('\\.bracken_frac','',colnames(data))
data_t = t(as.data.frame(data))
colnames(data_t) = data_t[1,]
data_t = data_t[2:nrow(data_t),]
data_t = data_t %>% as.data.frame %>% filter(grepl('frac',name))
data_t[,2] = as.numeric(data_t[,2])
data_t$name = gsub('\\.bracken_frac','',data_t$name)
colnames(data_t) = c('SAMPLE_ID','abundance')

metagenomic_metadata = read.csv('../metagenome_metadata_update.csv') %>% mutate(METAGENOME_SOURCE = paste(COHORT,SITE,sep='-'))
sampleidmap = read.table('sample_id_mapping',sep='\t',header=F) %>% select(V1,V2) %>% mutate(V2 = strsplit(V2, ' \\./') %>% map_chr(1) %>% gsub('_R1_clumped_bbduk_hg38removed_repaired_tadpole.fastq.gz','',.) %>% gsub('\\./','',.) %>%gsub('/','-',.))
colnames(sampleidmap) = c('NEWNAME','SAMPLE_ID')
metagenomic_metadata = left_join(metagenomic_metadata,sampleidmap) %>% mutate(SAMPLE_ID = if_else(!is.na(NEWNAME),NEWNAME,SAMPLE_ID)) %>% select(-NEWNAME)

merged = left_join(data_t,metagenomic_metadata)

merged$METAGENOME_SOURCE = factor(merged$METAGENOME_SOURCE,levels=c('MT1-ENV','MT2-ENV',"MT2-CREW-SALIVA","MT2-CREW-MOUTH","MT2-CREW-NASAL","MT2-CREW-EAR","MT2-CREW-POOL","TWINS-GUT","TWINS-SALIVA"))

merged = merged %>% mutate(PMA_TREATMENT= if_else(PMA_TREATMENT == '','NO PMA',PMA_TREATMENT))

ggplot(data = merged,aes(x = METAGENOME_SOURCE,color=PMA_TREATMENT, fill = PMA_TREATMENT,y = log(abundance+0.00001))) + geom_boxplot(alpha = .5) + geom_quasirandom(dodge.width = .8)+ theme(axis.text.x = element_text(angle = 45, vjust = 1,hjust = 1,size=10)) + ylab("log10(relative abundance)") + xlab('') + ggtitle('Acinetobacter pittii relative abundance across metagenome type') + stat_compare_means(comparisons = list(c('MT1-ENV','MT2-ENV')),method = 'wilcox.test')
ggsave('ap_abundance_by_metagenome.pdf',width=8,height=5)

ggplot(data = merged%>% filter(METAGENOME_SOURCE == 'MT1-ENV' |METAGENOME_SOURCE == 'MT2-ENV' ),aes(x = METAGENOME_SOURCE,color=PMA_TREATMENT, fill = PMA_TREATMENT,y = log(abundance+0.00001))) + geom_boxplot(alpha = .5) + geom_quasirandom(dodge.width = .8)+ theme(axis.text.x = element_text(angle = 45, vjust = 1,hjust = 1,size=10)) + ylab("log10(relative abundance)") + xlab('') + ggtitle('Acinetobacter pittii relative abundance across metagenome type') + stat_compare_means()
ggsave('ap_abundance_by_metagenome_pma.pdf',width=8,height=5)

```