---
title: "Updating population labels and geographical coordinates in metadata"
author: "Mary T. Yohannes"
date: "2023-12-05"
output: html_document
---


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

# import needed libraries 
library(tidyverse)
library(RColorBrewer)
library(plotly)
library(cowplot)
library(grid)
library(gridExtra)
library(maptools)
library(ggpubr)
```


## Update population labels in metadata 
### before getting into the PCA plots, we first have to update the population labels for the HGDP samples to reflect changes made in the 2020 paper which are more appropriate than those traditionally used 
### Got sample metadata with new popualtion labels from Anders Bergstrom 
### Updated three fields: hgdp_tgp_meta.Population, hgdp_tgp_meta.Latitude and hgdp_tgp_meta.Longitude 
```{r}
# read-in metadata from gnomAD 
hgdp_tgp_unfil <- read.delim('gnomad_meta_v1.tsv', header=T, sep='\t') 

# some sample IDs have 'v3.1::' prefix as done in gnomAD so stripping that in order to match correctly with the pc score files 
hgdp_tgp_unfil$s <- gsub('v3.1::', '', hgdp_tgp_unfil$s)

# exclude to only HGDP samples 
hgdp_meta <- hgdp_tgp_unfil[hgdp_tgp_unfil$hgdp_tgp_meta.Project == "HGDP",]

# read in file with the updated labels 
new_labels <- read.delim('hgdp_wgs.20190516.metadata.txt', header=T, sep='\t') 

# split "library" column on "." since there are some sample IDs that are the second half of those row values 
new_labels2 <- new_labels %>% separate(library, into = c("First", "Second"), "[.]", fill="left")

# samples found in the new file but not in old metadata 
not_old <- new_labels2[!(new_labels2$sample %in% hgdp_meta$s),]

# now check if there are any overlaps between the second half of the split column in the new file and the sample IDs from the old metadata 
not_old2 <- not_old[!(not_old$Second %in% hgdp_meta$s), ] # NOT FOUND IN OUR DATASET - 8 SAMPLES 

# samples found in the old metadata but not in new file 
not_new <- hgdp_meta[!(hgdp_meta$s %in% new_labels2$sample),] 

# now check if there are any overlaps between the second half of the split column in the new file and the sample IDs from the old metadata 
not_new2 <- not_new[!(not_new$s %in% new_labels2$Second),] # FOUND IN OURS BUT NOT THEIRS - 27 SAMPLES

# finding ones that are overlapped between the old metadata and new table
in_both1 <- new_labels2[(new_labels2$sample %in% hgdp_meta$s),] # compare with "sample" column
in_both1 <- in_both1[, -c(2:3)] # delete columns 2 & 3 to make joining easier later on 
in_both2 <- new_labels2[(new_labels2$Second %in% hgdp_meta$s),] # compare with "Second" column
in_both2 <- in_both2[, -c(1:2)] # delete columns 1 & 2
colnames(in_both2)[1] <- "sample" # rename "Second" (column name) to "sample" so we can join it with the "in_both1" table 
in_both <- rbind(in_both1, in_both2) # join two tables by rows 

# update old metadata file with new population labels from new table 
final_metadata <- hgdp_tgp_unfil %>% left_join(in_both, by = c("s"="sample")) 
final_metadata <- final_metadata[,-1]
final_metadata <- final_metadata %>% mutate(population = coalesce(population,hgdp_tgp_meta.Population),
                          latitude = coalesce(latitude,hgdp_tgp_meta.Latitude),
                          longitude = coalesce(longitude,hgdp_tgp_meta.Longitude))

# update the population labels for the samples that are not in the new table but are found in the old metadata 
# checked beforehand to see which population labels were still present that needed to be changed to the new ones - used table(final_metadata$population)
final_metadata <- final_metadata %>% mutate(population=ifelse(population=="BiakaPygmy", "Biaka", population),
                                            population=ifelse(population=="MbutiPygmy", "Mbuti", population),
                                            population=ifelse(population=="Mongola", "Mongolian", population),
                                            population=ifelse(population=="Melanesian", "Bougainville", population))

# in the new table, the geographical coordinates reported for Mongolian samples are (48.5, 119). In the old metadata, they were (45,111) so make sure the geographical coordinate are updated for all Mongolian samples
final_metadata <- final_metadata %>% mutate(latitude = case_when(population=="Mongolian" ~ 48.5,
                                                                 TRUE ~ latitude),
                                            longitude = case_when(population=="Mongolian" ~ 119,
                                                                  TRUE ~ longitude))

# do the same for the one Han sample that is found in the old metadata but not in the new table - new latitude coordinate is 32.3
final_metadata <- final_metadata %>% mutate(latitude = case_when(population=="Han" ~ 32.3,
                                                                 TRUE ~ latitude))

# do the same for the one Daur sample that is found in the old metadata but not in the new table - new latitude coordinate is 48.5
final_metadata <- final_metadata %>% mutate(latitude = case_when(population=="Daur" ~ 48.5,
                                                                 TRUE ~ latitude))

# do the same for the one Hazara sample that is found in the old metadata but not in the new table - new latitude coordinate is 33.5
final_metadata <- final_metadata %>% mutate(latitude = case_when(population=="Hazara" ~ 33.5,
                                                                 TRUE ~ latitude))

# do the same for the one Kalash sample that is found in the old metadata but not in the new table - new latitude coordinate is 36
final_metadata <- final_metadata %>% mutate(latitude = case_when(population=="Kalash" ~ 36,
                                                                 TRUE ~ latitude))

# sanity check 
# which rows/samples had population label difference between the old metadata and new table 
pop_diff <- final_metadata[!(final_metadata$hgdp_tgp_meta.Population == final_metadata$population), ]
# check that the number of samples that had changes to their population labels add up and are correct/consistent: 
table(pop_diff$hgdp_tgp_meta.Population)
table(pop_diff$population)

# select the columns needed for plotting 
plot_metadata <- final_metadata %>% select(s, project_meta.title, hgdp_tgp_meta.Genetic.region, population, hgdp_tgp_meta.Continent.colors, hgdp_tgp_meta.Pop.colors, latitude, longitude)

# rename column names to refer to them easily
colnames(plot_metadata) <- c("s", "project", "genetic.region", "population", "continent.colors", "pop.colors", "latitude", "longitude")    

# write out files
write.table(final_metadata, 'gnomad_meta_updated.tsv', row.names=FALSE, sep="\t", quote = FALSE) # updated metadata - everything merged 
write.table(plot_metadata, 'gnomad_meta_for_plots.tsv', row.names=FALSE, sep="\t", quote = FALSE) # for plotting - subsetted metadata 
```

