CrossPopGenoComp.Rmd

---
title: "Cross-Population SNP comparison"
output: 
  html_document:
    toc: true
    theme: united
    toc_depth: 3
    number_sections: true
    toc_float: true
    fig_caption: yes
---

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

<style>
p.caption {
  font-size: 1.5em;
}
</style>

```{css, echo=F}
pre code, pre, code {
  white-space: pre !important;
  overflow-x: scroll !important;
  word-break: keep-all !important;
  word-wrap: initial !important;
}
```

***

In this document we compare linkage disequilibrium (LD) patterns and minor allele frequency (MAF) values across the five super populations within the 1KG Phase 3 reference. We then evaluate polygenic scores stratified by similarity in LD and MAF across populations.

***

# Compare LD and MAF

***

## Trying out the approach

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
########
# Using bigsnpr
########
rm(list = ls())

library(bigsnpr)
library(bigreadr)
library(data.table)

NCORES <- as.integer(Sys.getenv("SLURM_JOB_CPUS_PER_NODE"))

# We need to read in the EUR genetic data, calculate LD and MAF of variants, then compare these values to values from other super populations. Do this for each chromosome seperately.

ref.bed <- bed(download_1000G("tmp-data"))
snp_readBed(ref.bed$bedfile)
ref<-snp_attach(gsub('.bed','.rds',ref.bed$bedfile))
ref_fam_2 <- bigreadr::fread2(sub_bed(ref.bed$bedfile, ".fam2"))

EUR_index <- which(ref_fam_2$`Super Population` == "EUR")

for(i in 22){
  print(i)
  ######
  # Read in EUR data and calculate LD and MAF
  ######
  
  i_index<-which(ref$map$chr ==i)[1:1000]

  CHR <- ref$map$chr[i_index]
  POS <- ref$map$physical.pos[i_index]
  POS2 <- snp_asGeneticPos(CHR, POS, dir ='/users/k1806347/brc_scratch/Data/Genetic_Map/CEU', ncores = NCORES)
  
  EUR_maf<-bed_MAF(ref.bed, ind.row = EUR_index, ind.col = i_index, ncores = NCORES)
  EUR_corr <- snp_cor(ref$genotypes, ind.row = EUR_index, ind.col = i_index, size = 3 / 1000, infos.pos=POS2, ncores = NCORES)
  
  EUR_ldscore<-Matrix::colSums(EUR_corr^2)

  ######
  # Read in data for other super populations and calculate LD and MAF  
  ######
  
  #pop<-c('AFR','AMR','EAS','SAS')
  pop<-c('EAS')
  
  for(k in 1:length(pop)){
  
    pop_k_index <- which(ref_fam_2$`Super Population` == pop[k])
  
    pop_k_maf<-bed_MAF(ref.bed, ind.row = pop_k_index, ind.col = i_index, ncores = NCORES)
    pop_k_corr <- snp_cor(ref$genotypes, ind.row = pop_k_index, ind.col = i_index, size = 3 / 1000, infos.pos=POS2, ncores = NCORES)
  
    #####
    # Calculate LD scores
    #####
    
    pop_k_ldscore<-Matrix::colSums(pop_k_corr^2)

    #####
    # Calculate the Fst
    #####
    
    fst_vec<-snp_fst(list(EUR_maf, pop_k_maf), min_maf = 0, overall = FALSE)
    
    #####
    # Calculate POPCORN C-scores (r2 covariance)
    #####
    
    colCors = function(x, y) { 
      sqr = function(x) x*x
      x   = sweep(x, 2, colMeans(x))
      y   = sweep(y, 2, colMeans(y))
      cov = colSums(x*y)
      cor = colSums(x*y) /  sqrt(colSums(sqr(x))*colSums(sqr(y)))
      return(list(cov=cov,cor=cor))
    }
    
    # Estimate correlation using sliding window approach
    window<-3
    step<-10
    slider<-(2*window)+step
    
    EUR_pop_k_cov_score<-NULL
    EUR_pop_k_cor_score<-NULL
    j<-1
    while(j){
      print(j)
    
      start<-step*(j-1)
      stop<-step*(j-1)+slider
      
      cat('Start =',start,'\n')
      cat('Stop =',stop,'\n')
      
      block<-POS2 >= start & POS2 < stop

      cat('NSNP =',sum(block),'\n')

      tmp<-colCors(as.matrix(EUR_corr[block, block]),as.matrix(pop_k_corr[block, block]))
    
      if(j == 1){
        comp_start <- 0
        comp_stop <- stop - window
      } else{
        comp_start <- start + window
        comp_stop <- stop - window
      }
  
      cat('Comp start =',comp_start,'\n')
      cat('Comp stop =',comp_stop,'\n')

      block_comp<-POS2[block] >= comp_start & POS2[block] < comp_stop

      cat('Comp NSNP =',sum(block_comp),'\n')

      EUR_pop_k_cov_score<-c(EUR_pop_k_cov_score,tmp$cov[block_comp])
      EUR_pop_k_cor_score<-c(EUR_pop_k_cor_score,tmp$cor[block_comp])
      
      if(comp_stop > max(POS2)){
        break()
      } else {
        j<-j+1
      }
    }
    
    #####
    # Tabulate data
    #####
    
    EUR_pop_k_dat<-data.frame( CHR=ref$map$chromosome[i_index],
                               SNP=ref$map$marker.ID[i_index],
                               BP=ref$map$physical.pos[i_index],
                               A1=ref$map$allele1[i_index],
                               A2=ref$map$allele2[i_index],
                               af_1=EUR_maf$af,
                               af_2=pop_k_maf$af,
                               Fst=fst_vec,
                               ldscore_1=EUR_ldscore,
                               ldscore_2=pop_k_ldscore,
                               cov_score=EUR_pop_k_cov_score,
                               cor_score=EUR_pop_k_cor_score)
    
    # Remove variants with a MAF < 0.01 in either population
    EUR_pop_k_dat<-EUR_pop_k_dat[EUR_pop_k_dat$af_1 >= 0.01 & 
                                 EUR_pop_k_dat$af_1 <= 0.99 &
                                 EUR_pop_k_dat$af_2 >= 0.01 & 
                                 EUR_pop_k_dat$af_2 <= 0.99,]
    
    EUR_pop_k_dat[,-1:-5]<-round(EUR_pop_k_dat[,-1:-5],5)

  }
}
 
# Compare to POPCORN C-scores
c_score<-fread('/scratch/groups/biomarkers-brc-mh/Reference_data/POPCORN/Ollie_Pain/CSCORES/EUR_EAS/EUR.EAS.chr22.scores.gen_effect.txt')
  
comp<-intersect(EUR_pop_k_dat$SNP, c_score$V3)

head(EUR_pop_k_dat[(EUR_pop_k_dat$SNP %in% comp),]) 
head(c_score[(c_score$V3 %in% comp),]) 

tail(EUR_pop_k_dat[(EUR_pop_k_dat$SNP %in% comp),]) 
tail(c_score[(c_score$V3 %in% comp),]) 

comp<-merge(EUR_pop_k_dat, c_score, by.x='SNP', by.y='V3')
# The output matches POPCORN nicely
# The LD scores are slightly different and the C-score are also slightly different
# This is probably due to differences in window size.

```
</details>

***

## Run for all populations

I have now created an Rscript to run this analysis so it can be split by chromosome and population and run in parallel.

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
  for chr in $(seq 1 22);do
    sbatch -p brc,shared --mem 30G -n 1 --nodes 1 /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/cross_pop/cross_pop.R \
                  --ref /scratch/groups/biomarkers-brc-mh/Reference_data/1KG_Phase3/bigsnpr/1000G_phase3_common_norel \
                  --chr ${chr} \
                  --pop1 ${pop1} \
                  --pop2 ${pop2} \
                  --n_cores 1 \
                  --output /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/chr${chr}/res
  done
done

# Create genome-wide files
pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
head -1 /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/chr1/res.chr1.${pop1}.${pop2}.txt > /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/CrossPop.${pop1}.${pop2}.GW.txt
tail -n +2 -q /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/chr*/res.chr*.${pop1}.${pop2}.txt >> /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/CrossPop.${pop1}.${pop2}.GW.txt
done

# This scripts works well enough but the memory consumption is high. I am sure bigsnpr/bigstatsr functions could make this more efficient.
# Also, the output may vary slightly if the sliding window step size is changed due to the inclusion SNPs beyond 3cM. Although these SNPs would have a correlation of 0 due to bigsnpr window size, inclusion of the zero values my slightly influence the result.
```
</details>

***

## Stratify SNPs by LD and MAF differences

Create SNP-list listing SNPs within deciles of cov_score, cor_score, and Fst. Also, output a list of SNPs with MAF > 0.01 in both samples (i.e. in the file just created), and also split SNPs into random decile. These sensitivity analyses will help understand the cause of improved prediction.

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
library(data.table)

pop1<-'EUR'

for(pop2 in c('EAS','AMR','AFR','SAS')){
  dat_all<-NULL
  for(i in 1:22){
    print(i)
    dat<-fread(paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/chr',i,'/res.chr',i,'.',pop1,'.',pop2,'.txt'))
    dat_all<-rbind(dat_all, dat)
  }
    
  dat_all$Fst_quant <- as.numeric(cut( dat_all$Fst, quantile(dat_all$Fst, prob = seq(0, 1, length = 11), type = 5), include.lowest=T))

  dat_all$cov_score_quant <- as.numeric(cut( dat_all$cov_score, quantile(dat_all$cov_score, prob = seq(0, 1, length = 11), type = 5), include.lowest=T))

  dat_all$cor_score_quant <- as.numeric(cut( dat_all$cor_score, quantile(dat_all$cor_score, prob = seq(0, 1, length = 11), type = 5), include.lowest=T))

    set.seed(1)
    dat_all$random<-rnorm(nrow(dat_all))
    dat_all$random_quant <- as.numeric(cut( dat_all$random, quantile(dat_all$random, prob = seq(0, 1, length = 11), type = 5), include.lowest=T))
  
  write.table(dat_all$SNP, paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/',pop1,'_',pop2,'.GW.MAF1.txt'), col.names=F, row.names=F, quote=F)

  for(i in 1:10){
    
    write.table(dat_all$SNP[dat_all$Fst_quant == i], paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/',pop1,'_',pop2,'.GW.Fst_quant_',i,'.txt'), col.names=F, row.names=F, quote=F)

    write.table(dat_all$SNP[dat_all$cov_score_quant == i], paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/',pop1,'_',pop2,'.GW.cov_quant_',i,'.txt'), col.names=F, row.names=F, quote=F)
    
    write.table(dat_all$SNP[dat_all$cor_score_quant == i], paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/',pop1,'_',pop2,'.GW.cor_quant_',i,'.txt'), col.names=F, row.names=F, quote=F)

    write.table(dat_all$SNP[dat_all$random_quant == i], paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_',pop1,'_',pop2,'/',pop1,'_',pop2,'.GW.random_quant_',i,'.txt'), col.names=F, row.names=F, quote=F)

  }
}

```
</details>

***

# External summary statistics

***

## Calculate polygenic scores

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
####
# Step 1. Create score files
####

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    for quant in $(seq 1 10); do
      for type in $(echo cov cor Fst random); do
      sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator.R \
        --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
        --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/EUR_samples.keep \
        --sumstats /users/k1806347/brc_scratch/Data/GWAS_sumstats/prs_comparison/cleaned/${gwas}.cleaned.gz \
        --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
        --extract /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/${pop1}_${pop2}.GW.${type}_quant_${quant}.txt \
        --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/${type}/${quant}/${pop2}_test_quant_${quant} \
        --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list
      done
    done
  done
done

# Run unstratified including EUR for comparison
# Really this code only needs to be run once as target sample does not affect the output
for pop2 in $(echo EUR EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator.R \
      --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
      --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/EUR_samples.keep \
      --sumstats /users/k1806347/brc_scratch/Data/GWAS_sumstats/prs_comparison/cleaned/${gwas}.cleaned.gz \
      --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
      --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/unstratified/${pop2}_test_unstratified \
      --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list
  done
done

# Run using unstratified but MAF > 0.01 in both pop
for pop2 in $(echo EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator.R \
      --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
      --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/EUR_samples.keep \
      --sumstats /users/k1806347/brc_scratch/Data/GWAS_sumstats/prs_comparison/cleaned/${gwas}.cleaned.gz \
      --extract /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/${pop1}_${pop2}.GW.MAF1.txt \
--plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
      --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/MAF1/${pop2}_test_MAF1 \
      --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list
  done
done

####
# Step 2. Calculate the polygenic scores in target sample
####

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    for quant in $(seq 1 10); do
      for type in $(echo cov cor Fst random); do
  
          sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
                --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
                --target_keep /users/k1806347/brc_scratch/Data/UKBB/DiverseAncestry/QC_pop/UKB.postQC.${pop2}.unrel.keep \
                --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/${type}/${quant}/${pop2}_test_quant_${quant} \
                --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/${type}/${quant}/${pop2}_test_quant_${quant}.${pop2}.scale \
                --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2}/1KGPhase3.w_hm3.${pop2}.chr \
                --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
                --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas}/${type}/${quant}/test_${pop2}_quant_${quant}.${type}
    
      done
    done
  done
done

# Run unstratified including EUR for comparison
for pop2 in $(echo EUR EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
            --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
            --target_keep /users/k1806347/brc_scratch/Data/UKBB/DiverseAncestry/QC_pop/UKB.postQC.${pop2}.unrel.keep \
            --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/unstratified/${pop2}_test_unstratified \
            --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/unstratified/${pop2}_test_unstratified.${pop2}.scale \
            --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2}/1KGPhase3.w_hm3.${pop2}.chr \
            --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
            --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas}/unstratified/test_${pop2}_unstratified
  
  done
done

# Run using MAF > 1% SNPs
for pop2 in $(echo EAS AMR AFR SAS); do
  for gwas in $(echo HEIG03 BODY04); do
    sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
            --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
            --target_keep /users/k1806347/brc_scratch/Data/UKBB/DiverseAncestry/QC_pop/UKB.postQC.${pop2}.unrel.keep \
            --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/MAF1/${pop2}_test_MAF1 \
            --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/score_files/pt_clump/${gwas}/MAF1/${pop2}_test_MAF1.${pop2}.scale \
            --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2}/1KGPhase3.w_hm3.${pop2}.chr \
            --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
            --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas}/MAF1/test_${pop2}_MAF1
  
  done
done

```
</details>

***

## Evaluate polygenic scores

<details><summary>Show code</summary>
```{r, eval=F, echo=T}
########
# Compare unstratified to prs-stratified by each metric seperately
########

gwas=$(echo HEIG03 BODY04)
pheno=$(echo Height BMI)

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 2); do
for type in $(echo cov cor Fst random); do

gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs_${type}.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/unstratified/test_${pop2}_unstratified.profiles unstratified
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/MAF1/test_${pop2}_MAF1.profiles MAF1
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/1/test_${pop2}_quant_1.${type}.profiles ${type}_1 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/2/test_${pop2}_quant_2.${type}.profiles ${type}_2 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/3/test_${pop2}_quant_3.${type}.profiles ${type}_3 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/4/test_${pop2}_quant_4.${type}.profiles ${type}_4 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/5/test_${pop2}_quant_5.${type}.profiles ${type}_5 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/6/test_${pop2}_quant_6.${type}.profiles ${type}_6 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/7/test_${pop2}_quant_7.${type}.profiles ${type}_7 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/8/test_${pop2}_quant_8.${type}.profiles ${type}_8 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/9/test_${pop2}_quant_9.${type}.profiles ${type}_9 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/${type}/10/test_${pop2}_quant_10.${type}.profiles ${type}_10
EOF

done
done
done

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 2); do
for type in $(echo cov cor Fst random); do

gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Data/UKBB/Phenotype/${pheno_i}/UKBB_${pheno_i}.score.UpdateIDs.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs_${type} \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs_${type}.txt

done
done
done

########
# Compare unstratified to all statified prs
########
# Don't include random or MAF1 prs

gwas=$(echo HEIG03 BODY04)
pheno=$(echo Height BMI)

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 2); do

gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/unstratified/test_${pop2}_unstratified.profiles unstratified
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/1/test_${pop2}_quant_1.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/2/test_${pop2}_quant_2.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/3/test_${pop2}_quant_3.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/4/test_${pop2}_quant_4.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/5/test_${pop2}_quant_5.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/6/test_${pop2}_quant_6.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/7/test_${pop2}_quant_7.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/8/test_${pop2}_quant_8.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/9/test_${pop2}_quant_9.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cov/10/test_${pop2}_quant_10.cov.profiles cov
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/1/test_${pop2}_quant_1.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/2/test_${pop2}_quant_2.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/3/test_${pop2}_quant_3.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/4/test_${pop2}_quant_4.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/5/test_${pop2}_quant_5.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/6/test_${pop2}_quant_6.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/7/test_${pop2}_quant_7.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/8/test_${pop2}_quant_8.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/9/test_${pop2}_quant_9.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/cor/10/test_${pop2}_quant_10.cor.profiles cor
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/1/test_${pop2}_quant_1.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/2/test_${pop2}_quant_2.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/3/test_${pop2}_quant_3.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/4/test_${pop2}_quant_4.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/5/test_${pop2}_quant_5.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/6/test_${pop2}_quant_6.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/7/test_${pop2}_quant_7.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/8/test_${pop2}_quant_8.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/9/test_${pop2}_quant_9.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/Fst/10/test_${pop2}_quant_10.Fst.profiles Fst
EOF

done
done

# Test association between prs and phenotype
gwas=$(echo HEIG03 BODY04)
pheno=$(echo Height BMI)

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 2); do

gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Data/UKBB/Phenotype/${pheno_i}/UKBB_${pheno_i}.score.UpdateIDs.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs.txt

done
done

########
# Run for EUR target sample for comparison
########

gwas=$(echo HEIG03 BODY04)
pheno=$(echo Height BMI)

pop1=EUR
pop2=EUR
for i in $(seq 1 2); do

gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/prs/pt_clump/${gwas_i}/unstratified/test_${pop2}_unstratified.profiles unstratified
EOF

done

# Test association between prs and phenotype
gwas=$(echo HEIG03 BODY04)
pheno=$(echo Height BMI)

pop1=EUR
pop2=EUR

for i in $(seq 1 2); do
gwas_i=$(echo ${gwas} | cut -f ${i} -d ' ')
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared --mem 40G /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Data/UKBB/Phenotype/${pheno_i}/UKBB_${pheno_i}.score.UpdateIDs.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas_i}/prs.txt
done

```
</details>

***

## Plot the results

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
library(data.table)
library(ggplot2)
library(cowplot)

type=c('cov','cor','Fst','random')
gwas=c('HEIG03','BODY04')
pheno=c('Height','BMI')
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')

# Read in European results and calculate portability improvement.
EUR_res_all<-NULL
pheno<-c('Height','BMI')
for(i in 1:length(gwas)){
EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/EUR_EUR/assoc/',gwas[i],'/prs.pred_eval.txt'))

EUR_res<-EUR_res[EUR_res$Model == 'All_group',]
EUR_res$Model<-'EUR'
EUR_res$Type<-'EUR'
EUR_res$GWAS<-gwas[i]
EUR_res$Population<-'EUR'
EUR_res$Method<-NA
EUR_res_all<-rbind(EUR_res_all, EUR_res)
}

# Read in the results
res_all<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(gwas)){
for(k in 1:length(type)){
  res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/',gwas[i],'/prs_',type[k],'.pred_eval.txt'))
  
  res_group<-res[!grepl('PredFile', res$Model),]
  res_group$Model<-gsub('_group', '', res_group$Model)
  res_group$Model<-gsub(paste0(type[k],'.'), paste0(type[k],'. Q'), res_group$Model)
  res_group$Model<-gsub('unstratified', "Unstratified", res_group$Model)
  
  res_group<-do.call(rbind, list(res_group[grepl(paste0(type[k],'. Q'), res_group$Model)], res_group[grepl('MAF1', res_group$Model)], res_group[grepl('Unstratified', res_group$Model)], res_group[grepl('All', res_group$Model)]))
  res_group$Model<-factor(res_group$Model, levels=res_group$Model)
  
  res_group$Type<-NA
  res_group$Type[grepl(paste0(type[k],'. Q'), res_group$Model)]<-paste0(type[k],'. Q')
  res_group$Type[grepl("Unstratified", res_group$Model)]<-"Unstratified"
  res_group$Type[grepl("MAF1", res_group$Model)]<-"MAF1"
  res_group$Type[grepl("All", res_group$Model)]<-"All"

  res_group_2<-res_group
  res_group_2$GWAS<-gwas[i]
  res_group_2$Population<-pop2[j]
  res_group_2$Method<-type[k]

  res_all<-rbind(res_all, res_group_2)
  
}
}
}

res_all<-rbind(res_all, EUR_res_all)
res_all$Type<-factor(res_all$Type, levels=c('cov. Q','cor. Q','Fst. Q','random. Q','MAF1','Unstratified','All','EUR'))

write.csv(res_all, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/seperate_stratification_summary.csv', row.names=F, quote=F)

# Plot the performance of stratifying by cov, cor, Fst
for(j in 1:length(pop2)){
for(i in 1:length(gwas)){
res_all_i_j<-res_all[(res_all$Population == pop2[j] | res_all$Population == 'EUR') & res_all$GWAS == gwas[i],]
res_plot<-list()
for(k in 1:length(type)){
    if(min(res_all_i_j$R - res_all_i_j$SE) < 0){
      y_lim<-c(min(res_all_i_j$R - res_all_i_j$SE),max(res_all_i_j$R + res_all_i_j$SE))
    } else {
      y_lim<-c(0,max(res_all_i_j$R + res_all_i_j$SE))
    }
  
  res_all_i_j_k<-res_all_i_j[res_all_i_j$Method == type[k] | is.na(res_all_i_j$Method),]
  res_all_i_j_k$Model<-factor(res_all_i_j_k$Model, levels=res_all_i_j_k$Model)

  res_plot[[k]]<-ggplot(res_all_i_j_k, aes(x=Model, y=R, colour=Type, shape=Population)) +
          geom_point(stat="identity", position=position_dodge( width=.25), size=5) +
          geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.25)) +
          scale_shape_manual(values=c(19, 15))+
          theme_half_open() +
          geom_hline(yintercept=0) +
          ylim(y_lim) +
          theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12), legend.position = 'none') +
          background_grid(major = 'y', minor = 'y') +
          labs(y="Correlation (SE)", title=paste0(pheno[i],': ', type[k],"-stratified (",pop2[j],")"))

}

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/',gwas[i],'/prs_res.png'), width=5000, height=2000, res=300, units='px')
  print(plot_grid(plotlist=res_plot, labels = "AUTO", nrow=1))
  dev.off()
}
}

# Plot results comparing cov, cor, Fst groups, unstratified and all results combined
combined_res<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(gwas)){
  res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/',gwas[i],'/prs.pred_eval.txt'))
  
  res_group<-res[!grepl('PredFile', res$Model),]
  res_group$Model<-gsub('_group', '', res_group$Model)
  res_group$Model<-gsub('unstratified', "Unstratified", res_group$Model)
  
  res_group$Model<-factor(res_group$Model, levels=res_group$Model)
  
  res_group$Type<-res_group$Model

  res_group_2<-res_group
  res_group_2$GWAS<-gwas[i]
  res_group_2$Population<-pop2[j]
  res_group_2$Method<-NA

  combined_res<-rbind(combined_res, res_group_2)
  
}
}

write.csv(combined_res, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/combined_stratification_summary.csv', row.names=F, quote=F)

res_all_2<-res_all[!grepl('. Q', res_all$Model),]
res_all_2<-res_all_2[!grepl('Unstratified|MAF1', res_all_2$Model),]
res_all_2<-res_all_2[!grepl('random', res_all_2$Method),]
res_all_2$Model<-as.character(res_all_2$Model)
res_all_2$Method<-as.character(res_all_2$Method)
res_all_2$Type<-as.character(res_all_2$Type)
res_all_2$Model[res_all_2$Method == 'cov']<-'Unstratified + cov'
res_all_2$Model[res_all_2$Method == 'cor']<-'Unstratified + cor'
res_all_2$Model[res_all_2$Method == 'Fst']<-'Unstratified + Fst'
res_all_2$Type[res_all_2$Method == 'cov']<-'Unstratified + cov'
res_all_2$Type[res_all_2$Method == 'cor']<-'Unstratified + cor'
res_all_2$Type[res_all_2$Method == 'Fst']<-'Unstratified + Fst'

res_all_2<-rbind(res_all_2, combined_res)

res_all_2$Type<-factor(res_all_2$Type, levels=c('cov','cor','Fst','Unstratified','Unstratified + cov','Unstratified + cor','Unstratified + Fst','All','EUR'))
res_all_2$Model<-factor(res_all_2$Model, levels=c('cov','cor','Fst','Unstratified','Unstratified + cov','Unstratified + cor','Unstratified + Fst','All','EUR'))

for(j in 1:length(pop2)){
res_plot<-list()
for(i in 1:length(gwas)){
res_all_2_i_j<-res_all_2[(res_all_2$Population == pop2[j] | res_all_2$Population == 'EUR') & res_all_2$GWAS == gwas[i],]

  if(min(res_all_2_i_j$R - res_all_2_i_j$SE) < 0){
    y_lim<-c(min(res_all_2_i_j$R - res_all_2_i_j$SE),max(res_all_2_i_j$R + res_all_2_i_j$SE))
  } else {
    y_lim<-c(0,max(res_all_2_i_j$R + res_all_2_i_j$SE))
  }
  
  res_plot[[i]]<-ggplot(res_all_2_i_j, aes(x=Model, y=R, colour=Type, shape=Population)) +
          geom_point(stat="identity", position=position_dodge( width=.25), size=5) +
          geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.25)) +
          scale_shape_manual(values=c(19, 15))+
          theme_half_open() +
          geom_hline(yintercept=0) +
          ylim(y_lim) +
          theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12), legend.position = 'none') +
          background_grid(major = 'y', minor = 'y') +
          labs(y="Correlation (SE)", title=paste0(pheno[i],': ',pop2[j]))

}

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/prs_res_summary.png'), width=2000, height=3000, res=300, units='px')
  print(plot_grid(plotlist=res_plot, labels = "AUTO", nrow=2))
  dev.off()
}

# Now make a plot comparing unstratified and all methods for all populations
res_all_3<-res_all_2[res_all_2$Model == 'All' | res_all_2$Model == 'Unstratified' | res_all_2$Model == 'EUR',]
res_all_3$Model<-as.character(res_all_3$Model)
res_all_3$Model[res_all_3$Model == 'EUR']<-'Unstratified'
res_all_3$Model[res_all_3$Model == 'All']<-'CrossPop'
res_all_3$Model<-factor(res_all_3$Model, levels=c('Unstratified','CrossPop'))
res_all_3$Type<-res_all_3$Model

write.csv(res_all_3, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/Unstratified_vs_stratified.csv', row.names=F, quote=F)

plot_list<-NULL
for(i in 1:length(gwas)){
  res_all_3_i<-res_all_3[res_all_3$GWAS == gwas[i],]
  
  if(min(res_all_3_i$R - res_all_3_i$SE) < 0){
    y_lim<-c(min(res_all_3_i$R - res_all_3_i$SE),max(res_all_3_i$R + res_all_3_i$SE))
  } else {
    y_lim<-c(0,max(res_all_3_i$R + res_all_3_i$SE))
  }

  plot_list[[i]]<-ggplot(res_all_3_i, aes(x=factor(Model), y=R, colour=Population, group=Population)) +
                            geom_point(stat="identity", position=position_dodge( width=.30)) +
                            geom_line(stat="identity", position=position_dodge( width=.30)) +
                            geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.30)) +
                            ylim(y_lim) +
                            theme_half_open() +
                            theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12)) +
                            background_grid(major = 'y', minor = 'y') +
                            labs(x='Method',y="Correlation (SE)", title=pheno[i])
  
  if(i != length(gwas)){
    plot_list[[i]]<-plot_list[[i]] + theme(legend.position = "none")
  }
}

png('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_res_summary_allpop.png', units='px', res=300, width=2000, height=900)
  plot_grid(plotlist=plot_list, nrow = 1, rel_widths = c(1.3/3, 1.7/3))
dev.off()

```
</details>

```{bash, eval=T, echo=F}
mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop

for pop1 in $(echo EUR); do
for pop2 in $(echo EAS AMR AFR SAS); do
for gwas in $(echo HEIG03 BODY04); do

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/${gwas}/prs_res.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop/${pop1}_${pop2}_${gwas}_prs_res.png

done

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/${pop1}_${pop2}/assoc/prs_res_summary.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop/${pop1}_${pop2}_prs_res_summary.png

done
done

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_res_summary_allpop.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop/
```

***

### Results for cov, cor, Fst, random seperately
<details><summary>Show figures</summary>

![EUR - AFR: BODY04](Images/CrossPop/EUR_AFR_BODY04_prs_res.png)
![EUR - EAS: BODY04](Images/CrossPop/EUR_EAS_BODY04_prs_res.png)
![EUR - AMR: BODY04](Images/CrossPop/EUR_AMR_BODY04_prs_res.png)
![EUR - SAS: BODY04](Images/CrossPop/EUR_SAS_BODY04_prs_res.png)

![EUR - AFR: HEIG03](Images/CrossPop/EUR_AFR_HEIG03_prs_res.png)
![EUR - EAS: HEIG03](Images/CrossPop/EUR_EAS_HEIG03_prs_res.png)
![EUR - AMR: HEIG03](Images/CrossPop/EUR_AMR_HEIG03_prs_res.png)
![EUR - SAS: HEIG03](Images/CrossPop/EUR_SAS_HEIG03_prs_res.png)

</details>

<details><summary>Show table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/combined_stratification_summary.csv")

res$P<-format(res$P, scientific = TRUE, digits = 3)
res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))]<-round(res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))], 3)

res<-res[,c('GWAS','Population','Model','R','SE','P','R2o','N')]

library(knitr)
kable(res, rownames = FALSE, caption='Predictive utility of stratified PRS')
```

</details>

***

* These results show that stratifying PRS by cov, cor and Fst does improve prediction over unstratified PRS.
* In most instances PRS based on SNPs with lowest cov or cor explain least variance. There is an increase in variance explained are cov and cor increase, as expected. However, the highest quantile for cov and cor often shows poor prediction. This should be investigated.
* The pattern of enrichment of Fst varies. There is a typically a decreased variance explained by high Fst PRS. BMI nicely shows this pattern except for in AFR, where high Fst increases variance explained. The pattern across Fst for height is very mixed.
* Stratify by random quantiles also improve prediction over unstratified PRS, though the improvement appears slightly smaller than when stratified by cov, cor or Fst. This suggests improvements based on cor, cov and Fst are not entirely due to these properties. This would probably be less of an issue if there were fewer PRS in the model, i.e. less quantiles, and less p-value thresholds.
* On average restricting to MAF 1% variants does not improve prediction.

***

### Results for cov, cor, and Fst combined
<details><summary>Show figures</summary>

![EUR - AFR](Images/CrossPop/EUR_AFR_prs_res_summary.png)
![EUR - EAS](Images/CrossPop/EUR_EAS_prs_res_summary.png)
![EUR - AMR](Images/CrossPop/EUR_AMR_prs_res_summary.png)
![EUR - SAS](Images/CrossPop/EUR_SAS_prs_res_summary.png)

</details>

<details><summary>Show table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/seperate_stratification_summary.csv")

res$P<-format(res$P, scientific = TRUE, digits = 3)
res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))]<-round(res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))], 3)

res<-res[,c('GWAS','Population','Method','Model','R','SE','P','R2o','N')]

library(knitr)
kable(res, rownames = FALSE, caption='Predictive utility of stratified PRS')
```

</details>

***

Results are looking promising:

* Results using cov-stratified PRS match results when using C-scores from POPCORN (Toshiki)
* Stratified PRS improve prediction in most instances
* Each type of stratification can work, but combining all approaches improves prediction the most

It is is unclear whether this improvement is due to enrichment of cross population portability metrics, or if it is merely due to splitting genetic variants into many groups. Does this matter? I suppose it would be interesting to test for cross population enrichment of cov, cor and Fst. Results across quantiles hint at enrichment, but this could be tested in a more controlled way perhaps. Perhaps estimating the correlation between genetic effects across cov, cor and Fst bins. LD clumping would be required after stratification.

***

### Summary results comparing unstratfified vs stratified
<details><summary>Show figure</summary>

![Summary](Images/CrossPop/prs_res_summary_allpop.png)

</details>

<details><summary>Show table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/Unstratified_vs_stratified.csv")

res$P<-format(res$P, scientific = TRUE, digits = 3)
res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))]<-round(res[,!grepl("Model|P|Type|GWAS|Population|Method",names(res))], 3)

res<-res[,c('GWAS','Population','Model','R','SE','P','R2o','N')]

library(knitr)
kable(res, rownames = FALSE, caption='Predictive utility of stratified PRS')
```

</details>

***

### Portabilty results

<details><summary>Show code</summary>
```{R, eval=F, echo=F}
# Read in the pred_comp results
library(data.table)
library(ggplot2)
library(cowplot)

type=c('cov','cor','Fst')
gwas=c('HEIG03','BODY04')
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')

# Read in European results and calculate portability improvement.
EUR_res_all<-NULL
pheno<-c('Height','BMI')
for(i in 1:length(gwas)){
EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/EUR_EUR/assoc/',gwas[i],'/prs.pred_eval.txt'))

EUR_res<-EUR_res[EUR_res$Model == 'All_group',]
EUR_res$GWAS<-gwas[i]
EUR_res_all<-rbind(EUR_res_all, EUR_res)
}

# Read in the cross pop results
res_comp_all<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(gwas)){
for(k in 1:length(type)){
  res_comp<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/',gwas[i],'/prs_',type[k],'.pred_comp.txt'))
  
  res_comp<-res_comp[res_comp$Model_1 == 'All' & res_comp$Model_2 == 'unstratified', ]

  res_comp$GWAS<-gwas[i]
  res_comp$Population<-pop2[j]
  res_comp$Method<-type[k]
  
  res_comp$EUR_R<-EUR_res_all$R[EUR_res_all$GWAS == gwas[i]]
  
  res_comp_all<-rbind(res_comp_all, res_comp)
  
}
res_comp<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/',pop1,'_',pop2[j],'/assoc/',gwas[i],'/prs.pred_comp.txt'))

res_comp<-res_comp[res_comp$Model_1 == 'All' & res_comp$Model_2 == 'unstratified', ]

res_comp$GWAS<-gwas[i]
res_comp$Population<-pop2[j]
res_comp$Method<-'All'

res_comp$EUR_R<-EUR_res_all$R[EUR_res_all$GWAS == gwas[i]]

res_comp_all<-rbind(res_comp_all, res_comp)

}
} 

# Calculate R2 difference
res_comp_all$Model_1_R2<-res_comp_all$Model_1_R^2
res_comp_all$Model_2_R2<-res_comp_all$Model_2_R^2
res_comp_all$R2_diff<-res_comp_all$Model_1_R2 - res_comp_all$Model_2_R2

# Calculate percentage difference
res_comp_all$R_perc_diff<-(res_comp_all$Model_1_R-res_comp_all$Model_2_R)/res_comp_all$Model_2_R
res_comp_all$R2_perc_diff<-(res_comp_all$Model_1_R2-res_comp_all$Model_2_R2)/res_comp_all$Model_2_R2

# Calculate portability
res_comp_all$Model_1_R_port<-res_comp_all$Model_1_R/res_comp_all$EUR_R
res_comp_all$Model_2_R_port<-res_comp_all$Model_2_R/res_comp_all$EUR_R
res_comp_all$R_port_perc_diff<-(res_comp_all$Model_1_R_port-res_comp_all$Model_2_R_port)/res_comp_all$Model_2_R_port
res_comp_all$R_port_diff<-(res_comp_all$Model_1_R_port-res_comp_all$Model_2_R_port)

res_comp_all$Model_1_R2_port<-res_comp_all$Model_1_R2/(res_comp_all$EUR_R^2)
res_comp_all$Model_2_R2_port<-res_comp_all$Model_2_R2/(res_comp_all$EUR_R^2)
res_comp_all$R2_port_perc_diff<-(res_comp_all$Model_1_R2_port-res_comp_all$Model_2_R2_port)/res_comp_all$Model_2_R2_port
res_comp_all$R2_port_diff<-(res_comp_all$Model_1_R2_port-res_comp_all$Model_2_R2_port)

write.csv(res_comp_all, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability.csv', row.names=F, quote=F)

# Calculate the average portability across phenotypes
average_port<-NULL
type<-c('cov','cor','Fst','All')
for(k in 1:length(type)){
  average_port<-rbind(average_port, data.frame(Method=type[k],
                                               Population='All',
                                               R_perc_diff=mean(res_comp_all$R_perc_diff[res_comp_all$Method == type[k]]),
                                               R2_perc_diff=mean(res_comp_all$R2_perc_diff[res_comp_all$Method == type[k]]),
                                               R_port_perc_diff=mean(res_comp_all$R_port_perc_diff[res_comp_all$Method == type[k]]),
                                               R2_port_perc_diff=mean(res_comp_all$R2_port_perc_diff[res_comp_all$Method == type[k]]),
                                               R_port_diff=mean(res_comp_all$R_port_diff[res_comp_all$Method == type[k]]),
                                               R2_port_diff=mean(res_comp_all$R2_port_diff[res_comp_all$Method == type[k]]),
                                               All_R_port=mean(res_comp_all$Model_1_R_port[res_comp_all$Method == type[k]]),
                                               Unstratified_R_port=mean(res_comp_all$Model_2_R_port[res_comp_all$Method == type[k]]),
                                               All_R2_port=mean(res_comp_all$Model_1_R2_port[res_comp_all$Method == type[k]]),
                                               Unstratified_R2_port=mean(res_comp_all$Model_2_R2_port[res_comp_all$Method == type[k]])))
}

# Plot the average portability across phenotypes
pop2=c('EAS','AMR','AFR','SAS')
for(k in 1:length(type)){
for(j in 1:length(pop2)){
  average_port<-rbind(average_port, data.frame(Method=type[k],
                                               Population=pop2[j],
                                               R_perc_diff=mean(res_comp_all$R_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_perc_diff=mean(res_comp_all$R2_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R_port_perc_diff=mean(res_comp_all$R_port_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_port_perc_diff=mean(res_comp_all$R2_port_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R_port_diff=mean(res_comp_all$R_port_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_port_diff=mean(res_comp_all$R2_port_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               All_R_port=mean(res_comp_all$Model_1_R_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               Unstratified_R_port=mean(res_comp_all$Model_2_R_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               All_R2_port=mean(res_comp_all$Model_1_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               Unstratified_R2_port=mean(res_comp_all$Model_2_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]])))
}
}

write.csv(average_port, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability_average.csv', row.names=F, quote=F)

# Create plot showing results
port_plot_dat<-average_port[,c('Method', 'Population','All_R_port','Unstratified_R_port','All_R2_port','Unstratified_R2_port')]

port_plot_dat<-melt(port_plot_dat)
port_plot_dat$Metric<-gsub('.*_','',gsub('_port','', port_plot_dat$variable))
port_plot_dat$Method2<-NA
port_plot_dat$Method2[grepl('All', port_plot_dat$variable)]<-paste0('Stratified: ', port_plot_dat$Method[grepl('All', port_plot_dat$variable)])
port_plot_dat$Method2[!grepl('All', port_plot_dat$variable)]<-'Unstratified'
port_plot_dat$Method<-port_plot_dat$Method2
port_plot_dat$Method2<-NULL
port_plot_dat<-port_plot_dat[!duplicated(port_plot_dat),]

port_plot_dat$Method<-factor(port_plot_dat$Method, levels=unique(c('Unstratified',port_plot_dat$Method)))

library(ggplot2)
library(cowplot)
R_plot<-ggplot(port_plot_dat[port_plot_dat$Metric == 'R',], aes(x=Population, y=value, fill=Method)) +
  geom_bar(stat="identity", position=position_dodge()) +
  labs(y='Portability', title="Average portability (R)") +
  ylim(c(0,1)) +
  theme_minimal_hgrid()

R2_plot<-ggplot(port_plot_dat[port_plot_dat$Metric == 'R2',], aes(x=Population, y=value, fill=Method)) +
  geom_bar(stat="identity", position=position_dodge()) +
  labs(y='Portability', title="Average portability (R2)") +
  ylim(c(0,1)) +
  theme_minimal_hgrid()

png('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability_summary.png', units='px', res=300, width=2000, height=1400)
plot_grid(R_plot, R2_plot, labels = c('A', 'B'), ncol=1)
dev.off()

```
</details>

```{bash, eval=T, echo=F}
mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability_summary.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop/
```

<details><summary>Show average portability figure</summary>

![Average portability](Images/CrossPop/prs_portability_summary.png)

</details>


<details><summary>Show portability table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability.csv")

res$R_diff_pval<-format(res$R_diff_pval, scientific = TRUE, digits = 3)
res[,!grepl("Model_1$|Model_2$|R_diff_pval|GWAS|Population|Method",names(res))]<-round(res[,!grepl("Model_1$|Model_2$|R_diff_pval|GWAS|Population|Method",names(res))], 3)

library(knitr)
kable(res, rownames = FALSE, caption='Portability')
```

</details>

<details><summary>Show average portability table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop/prs_portability_average.csv")

res[,!grepl("Population|Method",names(res))]<-round(res[,!grepl("Population|Method",names(res))], 3)

library(knitr)
kable(res, rownames = FALSE, caption='Average portability')
```

</details>

***

This shows portability improves when stratifying by cross population metrics, with cov, cor, Fst and All providing the large improvement respectively.

We should now expand this analysis to more phenotypes. This will require GWAS to be run in UKB.

***

# UKB summary statistics

Here we will generate GWAS summary statistics within a subset of EUR UKB participants. We will then use these for polygenic scoring in UKB. This will give us more control.

***

## Identify phenotypes

Let's use continuous phenotypes to help alleviate sample size issues in smaller populations. We need to identify phenotypes that are available in most individuals. Let use the same phenotypes as used by the Liang preprint regarding transcriptome risk scores.

<details><summary>Show code</summary>
```{r, eval=F, echo=T}
library(data.table)
liang_phen<-fread('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno.csv')
liang_phen$`UKB Field ID`<-paste0(liang_phen$`UKB Field ID`,'-0.0')

# Create vector of codes with sex and age included
codes<-c("21022-0.0","31-0.0",liang_phen$`UKB Field ID`)
  
write.table(codes, '/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_codes.txt', col.names=F, row.names=F, quote=F)

# Read in 
library(ukbkings)
bio_phen(
 project_dir = "/scratch/datasets/ukbiobank/ukb18177",
 field = "/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_codes.txt",
 out = "/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/ukb.liang_pheno"
)

dat <- readRDS("/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/ukb.liang_pheno.rds")

# Read in ancestry data
pop_dat<-bio_gen_ancestry('/scratch/datasets/ukbiobank/ukb18177')

dat<-merge(dat, pop_dat, by='eid')

# Read in PCs data
qc_dat<-bio_gen_sqc('/scratch/datasets/ukbiobank/ukb18177')
qc_dat<-qc_dat[,c('eid',paste0('pc',1:20))]
dat<-merge(dat, qc_dat, by='eid')

# Remove related individuals
rel<-bio_gen_related_remove('/scratch/datasets/ukbiobank/ukb18177', '/scratch/groups/ukbiobank/Edinburgh_Data/Software/tools/GreedyRelated/GreedyRelated')

dat<-dat[!(dat$eid %in% rel$eid),]

# Update IDs to match the genotype data I have created
fam<-fread('/users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam')[,1:2]
names(fam)[1:2]<-c('FID','IID')

dat<-merge(fam, dat, by.x='IID',by.y='eid')
dat$IID_2<-dat$IID
dat$IID<-dat$FID
dat$FID<-dat$IID_2
dat$IID_2<-NULL
names(dat)[1:2]<-c('FID','IID')

# Update column names to something meaningful
for(i in 1:length(liang_phen$`UKB Field ID`)){
  names(dat)[names(dat) == liang_phen$`UKB Field ID`[i]]<-liang_phen$Tag[i]
}
names(dat)[names(dat) == '31-0.0']<-'Sex'
names(dat)[names(dat) == '21022-0.0']<-'Age'

# Create table listing N per phenotype per population
pheno<-liang_phen$Tag
pop<-unique(dat$pop[!is.na(dat$pop)])

for(pop_i in pop){
  liang_phen$N_tmp<-NA
  names(liang_phen)[names(liang_phen) == 'N_tmp']<-paste0('N_',pop_i)
}

for(pop_i in pop){
for(pheno_i in pheno){
    liang_phen[[paste0('N_',pop_i)]][liang_phen$Tag == pheno_i]<-sum(!is.na(dat[[pheno_i]][which(dat$pop == pop_i)]))
}
}

# These phenotype look good as there are >1800 individuals for all pheno and all pop.

# Split UKB EUR individuals into training and testing samples
# This is will enable us to calculate PRS in UKB without overfitting
# Save 50k individuals for testing sample
set.seed(1)
target<-which(dat$pop == 'EUR')[sample(1:50000)]
training<-which(dat$pop == 'EUR')[!(which(dat$pop == 'EUR') %in% target)]

dat$pop[target]<-'EUR_test'
dat$pop[training]<-'EUR_train'

pop<-unique(dat$pop[!is.na(dat$pop)])

# Write out keep file for EUR_test population
fwrite(dat[dat$pop == 'EUR_test', c('FID','IID'), with=F], paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/EUR_test.keep'), sep=' ', na='NA', quote=F, col.names=F)

# For each population, normalise each phenotype and regress age, sex, and 20PCs.
covs<-names(dat)[!(names(dat) %in% c('FID','IID',pheno,'pop'))]
dat<-data.table(dat)
dat_resid<-list()
for(pop_i in pop){
  dir.create(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/',pop_i))
  dat_pop_i<-dat[which(dat$pop == pop_i),]
  dat_pop_i_resid<-dat_pop_i[,c('FID','IID'), with=F]
  for(pheno_i in pheno){
    dat_pop_i_pheno_i<-dat_pop_i[,c('FID','IID',pheno_i,covs), with=F]
    dat_pop_i_pheno_i<-dat_pop_i_pheno_i[complete.cases(dat_pop_i_pheno_i),]
    
    # Plot distribution of residuals
    png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/',pop_i,'/',pheno_i,'_raw.png'), res=300, unit='px', width=1000, height=1000)
    hist(dat_pop_i_pheno_i[[pheno_i]],xlab = pheno_i, main = NULL)
    dev.off()
  
    # Perform inverse rank based normalisation
    library(RNOmni)
    dat_pop_i_pheno_i[[pheno_i]]<-RNOmni::rankNorm(dat_pop_i_pheno_i[[pheno_i]])
    
    # Regress out covariates from pheno
    mod<-lm(as.formula(paste0(pheno_i,' ~ ',paste(covs[!(grepl('FID','IID', covs))],collapse=' + '))), data=dat_pop_i_pheno_i)
    
    # Output residuals
    dat_pop_i_pheno_i$resid<-as.numeric(scale(resid(mod)))
    
    # Plot distribution of residuals
    png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/',pop_i,'/',pheno_i,'_resid.png'), res=300, unit='px', width=1000, height=1000)
    hist(dat_pop_i_pheno_i[['resid']],xlab = pheno_i, main = NULL)
    dev.off()
    
    dat_pop_i_pheno_i_resid<-dat_pop_i_pheno_i[,c('FID','IID','resid'),with=F]
    names(dat_pop_i_pheno_i_resid)<-c('FID','IID',pheno_i)
    dat_pop_i_resid<-merge(dat_pop_i_resid, dat_pop_i_pheno_i_resid, by=c('FID','IID'), all=T)
  }
  
  dat_resid[[pop_i]]<-dat_pop_i_resid
}

# Write out the phenotype data for each population seperately.
for(pop_i in pop){
for(pheno_i in pheno){
  tmp<-dat_resid[[pop_i]][,c('FID','IID',pheno_i),with=F]
  tmp<-tmp[complete.cases(tmp),]
  fwrite(tmp, paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/',pop_i,'/',pheno_i,'_resid.pheno'), sep=' ', na='NA', quote=F)
}
}

# Write out file listing pheno
write.table(pheno, '/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', col.names=F, row.names=F, quote=F)
```
</details>

***

## Run GWAS

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
mkdir /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS

for pop in $(echo EUR_train AFR AMR EAS SAS); do
  mkdir /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop}
  for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
    mkdir /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop}/${pheno}
    for chr in $(seq 1 22); do
        sbatch -p brc,shared /users/k1806347/brc_scratch/Software/plink1.9.sh \
          --bfile /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr${chr} \
          --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop}/${pheno}_resid.pheno \
          --assoc \
          --maf 0.01 \
          --geno 0.05 \
          --allow-no-sex \
          --out /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop}/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.chr${chr}
    done
    sleep 45
  done
done

```
</details>

***

## Format sumstats

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
library(data.table)
library(qqman)

# Read in bim file
bim<-NULL
for(chr in 1:22){
  bim_chr<-fread(paste0('/users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr',chr,'.bim'))
  bim<-rbind(bim,bim_chr)
}

pheno<-fread('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

pop<-c('EUR_train','EUR_test','AFR','AMR','EAS','SAS')

for(pop_i in pop){
  for(pheno_i in pheno){
  print(pheno_i)
  gwas_i<-NULL
    for(chr in 1:22){
          gwas_i_chr<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.chr',chr,'.qassoc'))
          gwas_i_chr<-merge(gwas_i_chr, bim, by.x='SNP', by.y='V2')
          gwas_i_chr<-gwas_i_chr[,c('CHR','SNP','BP','V5','V6','BETA','SE','P','NMISS')]
          names(gwas_i_chr)<-c('CHR','SNP','ORIGBP','A1','A2','BETA','SE','P','N')
          gwas_i_chr<-gwas_i_chr[complete.cases(gwas_i_chr),]
      gwas_i<-rbind(gwas_i, gwas_i_chr)
    }
    
    fwrite(gwas_i, paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean'), sep=' ', na='NA', quote=F)
  
    # Create manhattan plot
    # Truncate p values at the minimum p-value that is not equal to 0
    gwas_i$P[gwas_i$P == 0]<-min(gwas_i$P[gwas_i$P != 0])
    bitmap(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean.Mahattan.png'), res=300, unit='px', width=4000, height=3000)
    manhattan(gwas_i[gwas_i$P < 1e-4,], bp='ORIGBP')
    dev.off()
    
    # Run GWAS sumstats cleaner
    system(paste0('/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/sumstat_cleaner/sumstat_cleaner.R --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/EUR/1KGPhase3.w_hm3.EUR.chr --output /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean.QC'))
  
  }
}

```
</details>

***

## Estimate heritability and rG using LDSC and POPCORN

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
####
# Munge sumstats
####

for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
  sbatch -p brc,shared ~/brc_scratch/Software/munge_sumstats.sh \
   --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz \
   --merge-alleles /users/k1806347/brc_scratch/Data/ldsc/w_hm3.snplist \
   --out /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC
done

####
# Estimate heritability
####

for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
   sbatch -p brc,shared ~/brc_scratch/Software/ldsc.sh \
     --h2 /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.sumstats.gz \
     --ref-ld-chr /users/k1806347/brc_scratch/Data/ldsc/eur_w_ld_chr/ \
     --w-ld-chr /users/k1806347/brc_scratch/Data/ldsc/eur_w_ld_chr/ \
     --out /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC
done

```

```{R, eval=F, echo=T}
####
# Tabulate the LDSC results
####

library(data.table)

pheno<-fread('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

res<-NULL
for(pheno_i in pheno){
      h2_i<-read.table(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean.QC.log'), header=F, sep='&')
      
      res<-rbind(res, data.frame( biomarker=pheno_i,
                                  h2_obs=gsub('Total Observed scale h2: ','',h2_i$V1[grepl('Total Observed scale h2: ',h2_i$V1)]),
                                  intercept=gsub('Intercept: ','',h2_i$V1[grepl('Intercept: ',h2_i$V1)]),
                                  lambda=gsub('Lambda GC: ','',h2_i$V1[grepl('Lambda GC: ',h2_i$V1)])))
                                  
}

# Calculate P value of heritability
res$h2_obs_est<-as.numeric(gsub(' .*','',res$h2_obs))
res$h2_obs_se<-as.numeric(gsub("\\)",'',gsub(".*\\(",'',res$h2_obs)))
res$h2_obs_p<-pnorm(-abs(res$h2_obs_est/res$h2_obs_se))
res$h2_obs_est<-NULL
res$h2_obs_se<-NULL

write.csv(res, '/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS_h2.csv', row.names=F, quote=F)

```

```{bash, eval=F, echo=T}
####
# Format for POPCORN
####

for pop in $(echo AFR AMR EAS SAS);do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
  zcat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop}/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz | sed -e '1s/BETA/beta/' -e '1s/REF.FREQ/AF/' | /users/k1806347/brc_scratch/Software/pigz > /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop}/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.popcorn.gz

done
done

for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
  zcat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz | sed -e '1s/BETA/beta/' -e '1s/REF.FREQ/AF/' | /users/k1806347/brc_scratch/Software/pigz > /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.popcorn.gz

done

####
# Run POPCORN
####

for pop2 in $(echo AFR AMR EAS SAS);do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
   sbatch -p brc,shared --mem 10G ~/brc_scratch/Software/popcorn.sh fit -v 3 \
     --cfile /scratch/groups/biomarkers-brc-mh/Reference_data/POPCORN/Ollie_Pain/CSCORES/EUR_${pop2}/EUR.${pop2}.GW.scores.gen_effect.txt \
     --sfile1 /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.popcorn.gz \
     --sfile2 /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/${pop2}/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.popcorn.gz \
     --gen_effect \
     /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/POPCORN_rG/EUR_train_vs_${pop2}_${pheno}_rG_gen_effect
done
sleep 300
done

```

```{R, eval=F, echo=T}
####
# Tabulate the POPCORN results
####

library(data.table)

pheno<-fread('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

res<-NULL
for(pop2 in c('AFR','AMR','EAS','SAS')){
  for(pheno_i in pheno){
        pop_res_i<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/POPCORN_rG/EUR_train_vs_',pop2,'_',pheno_i,'_rG_gen_effect'))
        
        res<-rbind(res, data.frame( biomarker=pheno_i,
                                    pop1='EUR_train',
                                    pop2=pop2,
                                    h2_obs_val_pop1=pop_res_i$`Val (obs)`[1],
                                    h2_obs_se_pop1=pop_res_i$SE[1],
                                    h2_obs_z_pop1=pop_res_i$Z[1],
                                    h2_obs_p_pop1=pop_res_i$`P (Z)`[1],
                                    h2_obs_val_pop2=pop_res_i$`Val (obs)`[2],
                                    h2_obs_se_pop2=pop_res_i$SE[2],
                                    h2_obs_z_pop2=pop_res_i$Z[2],
                                    h2_obs_p_pop2=pop_res_i$`P (Z)`[2],                                    
                                    pge_val=pop_res_i$`Val (obs)`[3],
                                    pge_se=pop_res_i$SE[3],
                                    pge_z_diff1=pop_res_i$Z[3],
                                    pge_p_diff1=pop_res_i$`P (Z)`[3]))
                                    
  }
}

write.csv(res, '/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS_popcorn_res.csv', row.names=F, quote=F)

```
</details>

<details><summary>Show LDSC heritability table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS_h2.csv")

res$h2_obs_p<-format(res$h2_obs_p, scientific = TRUE, digits = 3)

res<-res[,c('biomarker','h2_obs','h2_obs_p','intercept','lambda')]
names(res)<-c('Phenotype','SNP-h2','p','intercept','lambda')

library(knitr)
kable(res, rownames = FALSE, caption='Heritability')
```

</details>

<details><summary>Show POPCORN results table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS_popcorn_res.csv")

res$h2_obs_p_pop1<-format(res$h2_obs_p_pop1, scientific = TRUE, digits = 3)
res$h2_obs_p_pop2<-format(res$h2_obs_p_pop2, scientific = TRUE, digits = 3)
res$pge_p_diff1<-format(res$pge_p_diff1, scientific = TRUE, digits = 3)

library(knitr)
kable(res, rownames = FALSE, caption='POPCORN')
```

</details>

***

## Perform polygenic scoring

***

### Generate SCORE files

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
# Unstratified PRS
pop1=EUR
pop2=EUR_test
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
  sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator.R \
    --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
    --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/${pop1}_samples.keep \
    --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz \
    --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
    --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/unstratified/${pop1}_unstratified \
    --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list
done

# Stratified PRS
pop1=EUR
> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
for pop2 in $(echo EAS AMR AFR SAS); do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
for quant in $(seq 1 10); do
for type in $(echo cov cor Fst random); do

if [ ! -f /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant}.EUR.scale ]; then
echo ${pop2} ${pheno} ${quant} ${type} >> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
fi
done
done
done
done

# Create shell script to run using sbatch
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh << 'EOF'
#!/bin/sh

#SBATCH -p shared,brc
#SBATCH -J CrossPop_score

pop1=EUR
pop2=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $1}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
pheno=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $2}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
quant=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $3}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
type=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $4}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)

echo ${pop1}
echo ${pop2}
echo ${pheno}
echo ${quant} 
echo ${type} 

/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator.R \
        --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
        --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/${pop1}_samples.keep \
        --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz \
        --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
        --extract /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/${pop1}_${pop2}.GW.${type}_quant_${quant}.txt \
        --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant} \
        --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list

EOF

sbatch --array 1-$(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt | cut -d' ' -f1)%40 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh

############
# LD correlation weighted PRS
############

# Provide variant property file contaiing LD correlation for each variant across populations
# Then shrink BETA according to LD correlation, using a range of weights
# BETA*(LD_cor^weight)
# weight ranging from 0,0.5,1,1.5,2

# Stratified PRS
pop1=EUR
> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
for pop2 in $(echo EAS AMR AFR SAS); do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
for weight in 0 0.5 1 1.5 2; do
if [ ! -f /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight}.EUR.scale ]; then
echo ${pop2} ${pheno} ${weight} >> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
fi
done
done
done

# Create shell script to run using sbatch
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh << 'EOF'
#!/bin/sh

#SBATCH -p shared,brc
#SBATCH -J CrossPop_score

pop1=EUR
pop2=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $1}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
pheno=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $2}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
weight=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $3}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)

echo ${pop1}
echo ${pop2}
echo ${pheno}
echo ${weight} 

/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator_weighted.R \
        --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
        --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/${pop1}_samples.keep \
        --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz \
        --annot /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/CrossPop.${pop1}.${pop2}.GW.txt \
        --annot_col_name cor_score \
        --annot_weight ${weight} \
        --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
        --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight} \
        --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list

EOF

sbatch --array 1-$(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt | cut -d' ' -f1)%40 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh

sbatch --array 1-1 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh


############
# LD correlation weighted PRS (post-clump)
############

# Provide variant property file contaiing LD correlation for each variant across populations
# Then shrink BETA according to LD correlation, using a range of weights
# BETA*(LD_cor^weight)
# weight ranging from 0,0.5,1,1.5,2

# Stratified PRS
pop1=EUR
> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
for pop2 in $(echo EAS AMR AFR SAS); do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
for weight in 0 0.5 1 1.5 2; do
if [ ! -f /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight_postclump/${weight}/${pop2}_${weight}.EUR.scale ]; then
echo ${pop2} ${pheno} ${weight} >> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt
fi
done
done
done

# Create shell script to run using sbatch
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh << 'EOF'
#!/bin/sh

#SBATCH -p shared,brc
#SBATCH -J CrossPop_score

pop1=EUR
pop2=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $1}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
pheno=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $2}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)
weight=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $3}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt)

echo ${pop1}
echo ${pop2}
echo ${pheno}
echo ${weight} 

/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/polygenic_score_file_creator/polygenic_score_file_creator_weighted.R \
        --ref_plink_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/1KGPhase3.w_hm3.chr \
        --ref_keep /users/k1806347/brc_scratch/Data/1KG/Phase3/keep_files/${pop1}_samples.keep \
        --sumstats /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/${pheno}/UKBB.w_hm3.QCd.AllSNP.${pheno}.GW.qassoc.clean.QC.gz \
        --annot /scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_${pop1}_${pop2}/CrossPop.${pop1}.${pop2}.GW.txt \
        --annot_col_name cor_score \
        --annot_weight ${weight} \
        --postclump T \
        --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
        --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight_postclump/${weight}/${pop2}_${weight} \
        --ref_pop_scale /users/k1806347/brc_scratch/Data/1KG/Phase3/super_pop_keep.list

EOF

sbatch --array 1-$(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_todo.txt | cut -d' ' -f1)%40 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh

sbatch --array 1-1 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/score_sbatch.sh

```
</details>

***

### Characterise associations

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
library(data.table)

pheno<-fread('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

pop<-c('AFR','AMR','EAS','SAS')

mod_res_plot<-list()
mod_res_ldweighted_plot<-list()

g_cor_test_all<-NULL
g_cor_test_best_LDcor<-NULL
g_cor_test_best_Fst<-NULL
for(pheno_i in pheno){
  print(pheno_i)
  gwas<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/EUR_train/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean.QC.gz'))
  gwas<-gwas[order(gwas$CHR, gwas$BP),]

  indep<-NULL
  for(chr in 1:22){
    indep<-c(indep, fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_EUR_test/score_files/pt_clump/',pheno_i,'/unstratified/EUR_unstratified.chr',chr,'.score'))$V1)
  }
  
  gwas_indep<-gwas[(gwas$SNP %in% indep),]
  
  for(pop_i in pop){
    print(pop_i)
    # Check correlation between LD correlation, Fst and BETA
    annot<-NULL
    for(chr in 1:22){
      annot<-rbind(annot,fread(paste0('/scratch/groups/biomarkers-brc-mh/Reference_data/CrossPop/Populations_EUR_',pop_i,'/chr',chr,'/res.chr',chr,'.EUR.',pop_i,'.txt')))
    }
    
    gwas_annot_match<-merge(gwas_indep, annot[,c('SNP','A1','A2','af_1','af_2','Fst','ldscore_1','ldscore_2','cov_score','cor_score'), with=F], by=c('SNP','A1','A2'))
    gwas_annot_flip<-merge(gwas_indep, annot[,c('SNP','A1','A2','af_1','af_2','Fst','ldscore_1','ldscore_2','cov_score','cor_score'), with=F], by.x=c('SNP','A1','A2'),by.y=c('SNP','A2','A1'))
    gwas_annot_flip$af_1<-1-gwas_annot_flip$af_1
    gwas_annot_flip$af_2<-1-gwas_annot_flip$af_2
    
    gwas_annot<-rbind(gwas_annot_match, gwas_annot_flip)
    gwas_annot<-gwas_annot[order(gwas_annot$CHR, gwas_annot$BP),]
    
    # Flip effects to minor allle
    gwas_annot$A3<-NA
    gwas_annot$A3[gwas_annot$af_1 > 0.5]<-gwas_annot$A2[gwas_annot$af_1 > 0.5]
    gwas_annot$A2[gwas_annot$af_1 > 0.5]<-gwas_annot$A1[gwas_annot$af_1 > 0.5]
    gwas_annot$A1[gwas_annot$af_1 > 0.5]<-gwas_annot$A3[gwas_annot$af_1 > 0.5]
    gwas_annot$A3<-NULL
    gwas_annot$BETA[gwas_annot$af_1 > 0.5]<--gwas_annot$BETA[gwas_annot$af_1 > 0.5]
    gwas_annot$af_1[gwas_annot$af_1 > 0.5]<-1-gwas_annot$af_1[gwas_annot$af_1 > 0.5]
  
    cor(gwas_annot$af_1,gwas_annot$af_2)
    cor(abs(gwas_annot$BETA),gwas_annot$af_1)
    cor(abs(gwas_annot$BETA),gwas_annot$af_2)
    summary(lm(scale(abs(gwas_annot$BETA)) ~ scale(gwas_annot$af_1) + scale(gwas_annot$af_2)))
    cor(abs(gwas_annot$BETA),gwas_annot$Fst)
    cor(gwas_annot$ldscore_1,gwas_annot$ldscore_2)
    cor(abs(gwas_annot$BETA),gwas_annot$ldscore_1)
    cor(abs(gwas_annot$BETA),gwas_annot$ldscore_2)
    summary(lm(scale(abs(gwas_annot$BETA)) ~ scale(gwas_annot$ldscore_1) + scale(gwas_annot$ldscore_2)))
    cor(abs(gwas_annot$BETA),gwas_annot$cov_score)
    cor(abs(gwas_annot$BETA),gwas_annot$cor_score)
    
    # Low MAF SNPs have larger abs(BETA)
    # MAF correlation between ancestries is substantial
    # Controlling for EUR MAF, abs(BETA) is slightyl larger for high MAF SNPs
    # LD score highly correlated across ancestries
    # Small independent correlation between LD score in EUR and AFR ancestries and abs(BETA)
    # Therefore LD covariance is correlated with BETA
    # LD correlation and Fst show minimal correlation with abs(BETA)
    
    # Therefore stratifying SNPs by LD correlation and Fst is appropriate, as these metrics have a low correlation with abs(BETA) but may reflect the portability of SNPs across populations. Now test whether the Fst and LD_cor influence the correlation between betas across population. This can be modelled as an interaction. e.g. lm(beta_1 ~ beta_2*Fst)
    
    # Read in pop_i GWAS results
    gwas_pop_i<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/UKBB.w_hm3.QCd.AllSNP.',pheno_i,'.GW.qassoc.clean.QC.gz'))
    # Merge with gwas_annot
    names(gwas_pop_i)<-c('SNP','A1','A2','CHR','BP','BETA_2','SE_2','P_2','N_2','REF.FREQ')
    gwas_merged_match<-merge(gwas_annot, gwas_pop_i[,c('SNP','A1','A2','BETA_2','SE_2','P_2','N_2')], by=c('SNP','A1','A2'))
    gwas_merged_flip<-merge(gwas_annot, gwas_pop_i[,c('SNP','A1','A2','BETA_2','SE_2','P_2','N_2')], , by.x=c('SNP','A1','A2'),by.y=c('SNP','A2','A1'))
    
    gwas_merged_flip$BETA_2<--gwas_merged_flip$BETA
    gwas_merged<-rbind(gwas_merged_match, gwas_merged_flip)
    gwas_merged$same_sign<-F
    gwas_merged$same_sign[(gwas_merged$BETA > 0 & gwas_merged$BETA_2 >0 )|(gwas_merged$BETA < 0 & gwas_merged$BETA_2 <0)]  <-T
    
    # Perform sign test using  range of p-val thresholds 
    # Test interaction with Fst and LD correlation as well
    thresh_vec<-c(1,0.5,0.1,0.01,0.001,1e-5,5e-8)
    thresh_tab<-data.frame(p0=c(thresh_vec[-1],0),
                           p1=thresh_vec)
    g_cor_test<-NULL
    for(thresh_i in 1:nrow(thresh_tab)){
      gwas_merged_subset<-gwas_merged[gwas_merged$P > thresh_tab$p0[thresh_i] & gwas_merged$P < thresh_tab$p1[thresh_i],]
      
      basic<-summary(lm(scale(gwas_merged_subset$BETA) ~ scale(gwas_merged_subset$BETA_2)))
      Fst_int<-summary(lm(scale(gwas_merged_subset$BETA) ~ scale(gwas_merged_subset$BETA_2)*scale(gwas_merged_subset$Fst)))
      LDcor_int<-summary(lm(scale(gwas_merged_subset$BETA) ~ scale(gwas_merged_subset$BETA_2)*scale(gwas_merged_subset$cor_score)))
      
      g_cor_test<-rbind(g_cor_test, data.frame(p0=thresh_tab$p0[thresh_i],
                                               p1=thresh_tab$p1[thresh_i],
                                               thresh=paste0(thresh_tab$p0[thresh_i],' - ',thresh_tab$p1[thresh_i]),
                                               cor_est=coef(basic)[2,1],
                                               cor_se=coef(basic)[2,2],
                                               cor_p=coef(basic)[2,4],
                                               Fst_int_est=coef(Fst_int)[4,1],
                                               Fst_int_se=coef(Fst_int)[4,2],
                                               Fst_int_p=coef(Fst_int)[4,4],
                                               LDcor_int_est=coef(LDcor_int)[4,1],
                                               LDcor_int_se=coef(LDcor_int)[4,2],
                                               LDcor_int_p=coef(LDcor_int)[4,4],
                                                  sign_test_P=binom.test(sum(gwas_merged_subset$same_sign), nrow(gwas_merged_subset))$ p.value, 
                                               Nsnp=nrow(gwas_merged_subset)))
    }
    
    # ld cor strongly mediates the association between BETA and BETA2. When including the ld_cor interaction, the association between BETA and BETA2 becomes non-significant.
    # Fst does not show a strong interaction.
    
    # Format for plotting
    g_cor_test_melt<-NULL
    for(mod in c('cor','Fst_int','LDcor_int')){
      tmp<-g_cor_test[,c('thresh',paste0(mod,c('_est','_se','_p')))]
      names(tmp)<-c('pT','est','se','p')
      tmp$mod<-mod
      g_cor_test_melt<-rbind(g_cor_test_melt, tmp)
    }
    g_cor_test_melt$pT<-factor(g_cor_test_melt$pT, levels=rev(unique(g_cor_test_melt$pT)))
    
    library(ggplot2)
    library(cowplot)
    library(RColorBrewer)
    
    mod_res_plot[[paste0(pop_i,'_',pheno_i)]]<-ggplot(g_cor_test_melt, aes(x=mod, y=est, fill=pT)) +
      geom_bar(stat="identity", 
               position=position_dodge()) +
      geom_errorbar(aes(ymin=est-se, ymax=est+se), width=.2,
                     position=position_dodge(.9)) +
      scale_fill_manual(values = rev(brewer.pal(9, "Blues"))) +
      labs(x='Effect', y='Estimate', title=paste0(pop_i,': ',pheno_i)) +
  theme_half_open() +
  background_grid()
    
    png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/mod_res.png'), res=300, units = 'px', width=2000, height=1200)  
print(mod_res_plot[[paste0(pop_i,'_',pheno_i)]])
    dev.off()
    
    g_cor_test$pop1<-'EUR_train'
    g_cor_test$pop2<-pop_i
    g_cor_test$pheno<-pheno_i
    
    g_cor_test_all<-rbind(g_cor_test_all, g_cor_test)
    
    g_cor_test_best_LDcor<-rbind(g_cor_test_best_LDcor, g_cor_test[g_cor_test$LDcor_int_p == min(g_cor_test$LDcor_int_p),])
    
    g_cor_test_best_Fst<-rbind(g_cor_test_best_Fst, g_cor_test[g_cor_test$Fst_int_p == min(g_cor_test$Fst_int_p),])
    
    ## Alter EUR BETA according to LD correlatino
    ## Then retest correlation and interaction with LD correlation
    g_cor_test_2<-NULL
    for(thresh_i in 1:nrow(thresh_tab)){
        gwas_merged_subset<-gwas_merged[gwas_merged$P > thresh_tab$p0[thresh_i] & gwas_merged$P < thresh_tab$p1[thresh_i],]
      for(weight in seq(0,2,0.1)){
        gwas_merged_subset$BETA_mod<-gwas_merged_subset$BETA
        gwas_merged_subset$BETA_mod<-gwas_merged_subset$BETA*(gwas_merged_subset$cor_score^weight)
        
      basic<-summary(lm(scale(gwas_merged_subset$BETA_mod) ~ scale(gwas_merged_subset$BETA_2)))
      LDcor_int<-summary(lm(scale(gwas_merged_subset$BETA_mod) ~ scale(gwas_merged_subset$BETA_2)*scale(gwas_merged_subset$cor_score)))
      
      g_cor_test_2<-rbind(g_cor_test_2, data.frame(p0=thresh_tab$p0[thresh_i],
                                               p1=thresh_tab$p1[thresh_i],
                                               thresh=paste0(thresh_tab$p0[thresh_i],' - ',thresh_tab$p1[thresh_i]),
                                               weight=weight,
                                               cor_est=coef(basic)[2,1],
                                               cor_se=coef(basic)[2,2],
                                               cor_p=coef(basic)[2,4],
                                               LDcor_int_est=coef(LDcor_int)[4,1],
                                               LDcor_int_se=coef(LDcor_int)[4,2],
                                               LDcor_int_p=coef(LDcor_int)[4,4],                                               Nsnp=nrow(gwas_merged_subset)))        
      }
    }
    
    library(ggplot2)
    library(cowplot)
    library(RColorBrewer)
    
    g_cor_test_2$pT<-factor(g_cor_test_2$thresh, levels=rev(unique(g_cor_test_2$thresh)))

    mod_res_ldweighted_plot[[paste0(pop_i,'_',pheno_i)]]<-ggplot(g_cor_test_2, aes(x=weight, y=cor_est, fill=weight)) +
      geom_bar(stat="identity", 
               position=position_dodge()) +
      geom_errorbar(aes(ymin=cor_est-cor_se, ymax=cor_est+cor_se), width=.1,
                     position=position_dodge(.9)) +
      labs(x='Weight', y='Estimate', title=paste0(pop_i,': ',pheno_i)) +
  theme_half_open() +
  background_grid() + 
    facet_grid(. ~ pT)
    
    png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/mod_res_ldcor_weighted.png'), res=300, units = 'px', width=5000, height=1200)  
    print(mod_res_ldweighted_plot[[paste0(pop_i,'_',pheno_i)]])
    dev.off()
    # This nicely shows that the correlation between BETAs can be increased by weighting variants by LD correlation. This does not change the independent variants selected though, which it would if BETAs are modified prior to LD clumping.

    # Make a plot showing the difference in correlation between BETAs in different LD correlation deciles
    LD_cor_plot_list<-list()
    Fst_plot_list<-list()
    
    for(thresh_i in 1:nrow(thresh_tab)){
        gwas_merged_subset<-gwas_merged[gwas_merged$P > thresh_tab$p0[thresh_i] & gwas_merged$P < thresh_tab$p1[thresh_i],]
    gwas_merged_subset$cor_score_quant<-as.numeric(cut( gwas_merged_subset$cor_score, quantile(gwas_merged_subset$cor_score, prob = seq(0, 1, length = 6), type = 5), include.lowest=T))
    gwas_merged_subset$Fst_quant<-as.numeric(cut( gwas_merged_subset$Fst, quantile(gwas_merged_subset$Fst, prob = seq(0, 1, length = 6), type = 5), include.lowest=T))
    
    library(ggplot2)
    library(cowplot)
    LD_cor_plot_list[[as.character(thresh_i)]]<-ggplot(gwas_merged_subset, aes(x = BETA, y = BETA_2, colour=cor_score_quant, group=cor_score_quant)) + 
  geom_point() +
  labs(x="BETA (EUR)", y=paste0("BETA (",pop_i,")"), title=paste0(pheno_i," (pT=",thresh_i,")"), colour='LD correlation\nQuantile') +
  stat_smooth(method = "lm") +
  theme_half_open() +
  background_grid()

    Fst_plot_list[[as.character(thresh_i)]]<-ggplot(gwas_merged_subset, aes(x = BETA, y = BETA_2, colour=Fst_quant, group=Fst_quant)) + 
  geom_point() +
  labs(x="BETA (EUR)", y=paste0("BETA (",pop_i,")"), title=paste0(pheno_i," (pT=",thresh_i,")"), colour='Fst Quantile') +
  stat_smooth(method = "lm") +
  theme_half_open() +
  background_grid()

  # This fits with the findings from PRS analysis, that high LD correlation SNPs predict better across populations.
  # Now what is the best way to harness that information.      
    }
  }
}

sum(g_cor_test_best_LDcor$LDcor_int_p < 0.05/nrow(thresh_tab))
sum(g_cor_test_best_Fst$Fst_int_p < 0.05/nrow(thresh_tab))

for(pop_i in pop){
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/mod_res_all_',pop_i,'.png'), res=300, units = 'px', width=3500, height=length(which(grepl(pop_i, names(mod_res_plot))))*500)  
  print(plot_grid(plotlist=mod_res_plot[which(grepl(pop_i, names(mod_res_plot)))], ncol = 2))
  dev.off()
  
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/mod_res_all_ldcor_weighted_',pop_i,'.png'), res=300, units = 'px', width=5500, height=length(which(grepl(pop_i, names(mod_res_plot))))*500)  
  print(plot_grid(plotlist=mod_res_ldweighted_plot[which(grepl(pop_i, names(mod_res_plot)))], ncol = 2))
  dev.off()
}


  # Output plots
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/LD_cor_int.png'), res=300, units = 'px', width=2000, height=10000)  
  print(plot_grid(plotlist=LD_cor_plot_list, ncol = 1))
  dev.off()
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/GWAS/',pop_i,'/',pheno_i,'/Fst_int.png'), res=300, units = 'px', width=2000, height=10000)  
  print(plot_grid(plotlist=Fst_plot_list, ncol = 1))
  dev.off()
  }
}

  # Create miami plot comparing some EUR GWAS hits to non-EUR
  # Colour points by LD correlation
  # Show before and after adjusting BETA/P for LD correlation
  n_examp<-5
  wind<-0.05e6
  for(indep_i in gwas_indep[gwas_indep$P < 5e-8,]$SNP[1:n_examp]){
    EUR_gwas_examp<-gwas[gwas$BP > (gwas$BP[gwas$SNP == indep_i]-wind) & gwas$BP < (gwas$BP[gwas$SNP == indep_i]+wind)]
    pop_gwas_examp<-gwas_pop_i[gwas_pop_i$BP > (gwas_pop_i$BP[gwas_pop_i$SNP == indep_i]-wind) & gwas_pop_i$BP < (gwas_pop_i$BP[gwas_pop_i$SNP == indep_i]+wind)]
    
    EUR_gwas_examp$Population<-'EUR'
    pop_gwas_examp$Population<-pop_i
    names(pop_gwas_examp)<-names(EUR_gwas_examp)
    both_examp<-rbind(EUR_gwas_examp, pop_gwas_examp)
    both_examp_annot<-merge(both_examp, annot)
    both_examp_annot<-merge(both_examp, annot[,c('SNP','Fst','cor_score'), with=F], by='SNP')
    
    both_examp_annot$minuslog10P<--log10(both_examp_annot$P)
    both_examp_annot$minuslog10P[both_examp_annot$Population != 'EUR']<--both_examp_annot$minuslog10P[both_examp_annot$Population != 'EUR']
    
    print(ggplot(both_examp_annot, aes(x=BP, y=minuslog10P, colour=cor_score)) +
      geom_point() +
      scale_color_gradient2(midpoint=0, low="blue", mid="white", high="red", space ="Lab" ) +
      theme_half_open() +
      background_grid())
    
    # Modify p-value in EUR according to LD_correlation
    both_examp_annot$Z<-abs(qnorm(both_examp_annot$P/2))
    for(weight in c(0,0.5,1,2,4)){
      both_examp_annot$Z2<-both_examp_annot$Z
      both_examp_annot$Z2[both_examp_annot$Population == 'EUR']<-both_examp_annot$Z[both_examp_annot$Population == 'EUR']*(both_examp_annot$cor_score[both_examp_annot$Population == 'EUR']^weight)
      both_examp_annot$P2<-2*pnorm(-abs(both_examp_annot$Z2))
      both_examp_annot$minuslog10P2<--log10(both_examp_annot$P2)
      both_examp_annot$minuslog10P2[both_examp_annot$Population != 'EUR']<--both_examp_annot$minuslog10P2[both_examp_annot$Population != 'EUR']

      print(ggplot(both_examp_annot, aes(x=BP, y=minuslog10P2, colour=cor_score)) +
        geom_point() +
        scale_color_gradient2(midpoint=0, low="blue", mid="white", high="red", space ="Lab" ) +
        theme_half_open() +
        background_grid())
    }
    # This nicely shows how weighting variants by cor_score changes the lead variant
}
}

```
</details>

***

### Calculate polygenic scores

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
# Unstratified
pop1=EUR
for pop2 in $(echo EUR_test EAS AMR AFR SAS); do
pop2_short=${pop2}
if echo "$pop2" | grep -q "EUR_test"; then
    pop2_short=EUR
fi  

for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
    sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
            --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
            --target_keep /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno}_resid.pheno \
            --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop1}_test/score_files/pt_clump/${pheno}/unstratified/${pop1}_unstratified \
            --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop1}_test/score_files/pt_clump/${pheno}/unstratified/${pop1}_unstratified.${pop2_short}.scale \
            --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2_short}/1KGPhase3.w_hm3.${pop2_short}.chr \
            --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
            --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno}/unstratified/${pop2}_unstratified
  
  done
done

# Stratified
# This is a lot of jobs. Write into a sbatch job for better control
pop1=EUR
> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt
for pop2 in $(echo EAS AMR AFR SAS); do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
for quant in $(seq 1 10); do
for type in $(echo cov cor Fst random); do

if [ ! -f /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant}.${type}.profiles ]; then
echo ${pop2} ${pheno} ${quant} ${type} >> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt
fi
done
done
done
done

# Create shell script to run using sbatch
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_sbatch.sh << 'EOF'
#!/bin/sh

#SBATCH -p shared,brc
#SBATCH -J CrossPop_prs

pop1=EUR
pop2=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $1}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)
pheno=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $2}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)
quant=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $3}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)
type=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $4}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)

echo ${pop1}
echo ${pop2}
echo ${pheno}
echo ${quant} 
echo ${type} 

/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
                --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
                --target_keep /users/k1806347/brc_scratch/Data/UKBB/DiverseAncestry/QC_pop/UKB.postQC.${pop2}.unrel.keep \
                --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant} \
                --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant}.${pop2}.scale \
                --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2}/1KGPhase3.w_hm3.${pop2}.chr \
                --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
                --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno}/${type}/${quant}/${pop2}_quant_${quant}.${type}
    
EOF

sbatch --array 1-$(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt | cut -d' ' -f1)%100 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_sbatch.sh

########
# LD cor weighted PRS
########

pop1=EUR
> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt
for pop2 in $(echo EAS AMR AFR SAS); do
for pheno in $(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt); do
for weight in 0 0.5 1 1.5 2; do

if [ ! -f /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight}.weight.profiles ]; then
echo ${pop2} ${pheno} ${weight} >> /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt
fi
done
done
done

# Create shell script to run using sbatch
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_sbatch.sh << 'EOF'
#!/bin/sh

#SBATCH -p shared,brc
#SBATCH -J CrossPop_prs

pop1=EUR
pop2=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $1}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)
pheno=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $2}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)
weight=$(awk -v var="$SLURM_ARRAY_TASK_ID" 'NR == var {print $3}' /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt)

echo ${pop1}
echo ${pop2}
echo ${pheno}
echo ${weight} 

/users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Scaled_polygenic_scorer/Scaled_polygenic_scorer.R \
                --target_plink_chr /users/k1806347/brc_scratch/Data/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr \
                --target_keep /users/k1806347/brc_scratch/Data/UKBB/DiverseAncestry/QC_pop/UKB.postQC.${pop2}.unrel.keep \
                --ref_score /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight} \
                --ref_scale /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/score_files/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight}.${pop2}.scale \
                --ref_freq_chr /users/k1806347/brc_scratch/Data/1KG/Phase3/freq_files/${pop2}/1KGPhase3.w_hm3.${pop2}.chr \
                --plink /users/k1806347/brc_scratch/Software/plink1.9.sh \
                --output /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno}/weight/${weight}/${pop2}_${weight}.weight
    
EOF

sbatch --array 1-$(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_todo.txt | cut -d' ' -f1)%50 /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_sbatch.sh

```
</details>

***

### Evaluate polygenic scores

***

#### EUR

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
# Check prediction in EUR sample using unstratified PRS for comparison
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)
pop1=EUR
pop2=EUR_test

for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/unstratified/${pop2}_unstratified.profiles unstratified
EOF

done

# Test association between prs and phenotype
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR
pop2=EUR_test

for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do
pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared --mem 10G /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
        --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno_i}_resid.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs.txt
done

```
</details>

***

#### Non-EUR

***

#### cov, cor, Fst, random seperately

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)
pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do
for type in $(echo cov cor Fst random); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_${type}.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/unstratified/${pop2}_unstratified.profiles unstratified
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/1/${pop2}_quant_1.${type}.profiles ${type}_1 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/2/${pop2}_quant_2.${type}.profiles ${type}_2 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/3/${pop2}_quant_3.${type}.profiles ${type}_3 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/4/${pop2}_quant_4.${type}.profiles ${type}_4 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/5/${pop2}_quant_5.${type}.profiles ${type}_5 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/6/${pop2}_quant_6.${type}.profiles ${type}_6 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/7/${pop2}_quant_7.${type}.profiles ${type}_7 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/8/${pop2}_quant_8.${type}.profiles ${type}_8 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/9/${pop2}_quant_9.${type}.profiles ${type}_9 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/${type}/10/${pop2}_quant_10.${type}.profiles ${type}_10
EOF

done
done
done

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do
for type in $(echo cov cor Fst random); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno_i}_resid.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_${type} \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_${type}.txt

done
done
sleep 600
done

```
</details>

***

#### cov, cor, Fst together

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/unstratified/${pop2}_unstratified.profiles unstratified
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/1/${pop2}_quant_1.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/2/${pop2}_quant_2.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/3/${pop2}_quant_3.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/4/${pop2}_quant_4.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/5/${pop2}_quant_5.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/6/${pop2}_quant_6.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/7/${pop2}_quant_7.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/8/${pop2}_quant_8.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/9/${pop2}_quant_9.cov.profiles cov 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/10/${pop2}_quant_10.cov.profiles cov
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/1/${pop2}_quant_1.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/2/${pop2}_quant_2.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/3/${pop2}_quant_3.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/4/${pop2}_quant_4.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/5/${pop2}_quant_5.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/6/${pop2}_quant_6.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/7/${pop2}_quant_7.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/8/${pop2}_quant_8.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/9/${pop2}_quant_9.cor.profiles cor 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/10/${pop2}_quant_10.cor.profiles cor
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/1/${pop2}_quant_1.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/2/${pop2}_quant_2.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/3/${pop2}_quant_3.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/4/${pop2}_quant_4.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/5/${pop2}_quant_5.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/6/${pop2}_quant_6.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/7/${pop2}_quant_7.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/8/${pop2}_quant_8.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/9/${pop2}_quant_9.Fst.profiles Fst 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/10/${pop2}_quant_10.Fst.profiles Fst
EOF

done
done

# Test association between prs and phenotype
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno_i}_resid.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs.txt

done
done

```
</details>

***

#### cov, cor, Fst together with random

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_withRandom.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/unstratified/${pop2}_unstratified.profiles unstratified_random
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/1/${pop2}_quant_1.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/2/${pop2}_quant_2.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/3/${pop2}_quant_3.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/4/${pop2}_quant_4.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/5/${pop2}_quant_5.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/6/${pop2}_quant_6.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/7/${pop2}_quant_7.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/8/${pop2}_quant_8.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/9/${pop2}_quant_9.cov.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cov/10/${pop2}_quant_10.cov.profiles ld_freq_strat
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/1/${pop2}_quant_1.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/2/${pop2}_quant_2.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/3/${pop2}_quant_3.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/4/${pop2}_quant_4.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/5/${pop2}_quant_5.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/6/${pop2}_quant_6.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/7/${pop2}_quant_7.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/8/${pop2}_quant_8.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/9/${pop2}_quant_9.cor.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/cor/10/${pop2}_quant_10.cor.profiles ld_freq_strat
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/1/${pop2}_quant_1.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/2/${pop2}_quant_2.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/3/${pop2}_quant_3.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/4/${pop2}_quant_4.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/5/${pop2}_quant_5.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/6/${pop2}_quant_6.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/7/${pop2}_quant_7.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/8/${pop2}_quant_8.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/9/${pop2}_quant_9.Fst.profiles ld_freq_strat 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/Fst/10/${pop2}_quant_10.Fst.profiles ld_freq_strat
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/1/${pop2}_quant_1.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/2/${pop2}_quant_2.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/3/${pop2}_quant_3.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/4/${pop2}_quant_4.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/5/${pop2}_quant_5.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/6/${pop2}_quant_6.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/7/${pop2}_quant_7.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/8/${pop2}_quant_8.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/9/${pop2}_quant_9.random.profiles unstratified_random 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/random/10/${pop2}_quant_10.random.profiles unstratified_random

EOF

done
done

# Test association between prs and phenotype
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno_i}_resid.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_withRandom \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_withRandom.txt

done
done

```
</details>

#### LDcor weighted

<details><summary>Show code</summary>
```{bash, eval=F, echo=T}
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR

for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

mkdir -p /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/

# Create a file containing the name of the polygenic score file
cat > /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_LDcorWeighted.txt <<EOF
predictors group 
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/unstratified/${pop2}_unstratified.profiles unstratified
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/weight/0/${pop2}_0.weight.profiles weighted_0
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/weight/0.5/${pop2}_0.5.weight.profiles weighted_0.5
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/weight/1/${pop2}_1.weight.profiles weighted_1
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/weight/1.5/${pop2}_1.5.weight.profiles weighted_1.5
/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/prs/pt_clump/${pheno_i}/weight/2/${pop2}_2.weight.profiles weighted_2

EOF

done
done

# Test association between prs and phenotype
pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

pop1=EUR
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

sbatch -p brc,shared /users/k1806347/brc_scratch/Software/Rscript.sh /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Scripts/Model_builder/Model_builder_V2_nested.R \
      --pheno /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/Residuals/${pop2}/${pheno_i}_resid.pheno \
      --out /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_LDcorWeighted \
      --compare_predictors T \
      --assoc T \
      --predictors /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_LDcorWeighted.txt

done
done

```
</details>

***

## Plot results

### Absolute

<details><summary>Show code</summary>
```{R, eval=F, echo=T}
library(data.table)
library(ggplot2)
library(cowplot)

type=c('cov','cor','Fst','random')
pheno=read.table('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')

# Read in European results and calculate portability improvement.
EUR_res_all<-NULL
for(i in 1:length(pheno)){
EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_EUR_test/assoc/',pheno[i],'/prs.pred_eval.txt'))

EUR_res<-EUR_res[EUR_res$Model == 'All_group',]
EUR_res$Model<-'EUR'
EUR_res$Type<-'EUR'
EUR_res$GWAS<-pheno[i]
EUR_res$Population<-'EUR'
EUR_res$Method<-NA
EUR_res_all<-rbind(EUR_res_all, EUR_res)
}

# Read in the results
res_all<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(pheno)){
for(k in 1:length(type)){
  res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs_',type[k],'.pred_eval.txt'))
  
  res_group<-res[!grepl('PredFile', res$Model),]
  res_group$Model<-gsub('_group', '', res_group$Model)
  res_group$Model<-gsub(paste0(type[k],'.'), paste0(type[k],'. Q'), res_group$Model)
  res_group$Model<-gsub('unstratified', "Unstratified", res_group$Model)
  
  res_group<-rbind(res_group[grepl(paste0(type[k],'. Q'), res_group$Model)], res_group[!grepl(paste0(type[k],'. Q'), res_group$Model)])
  res_group$Model<-factor(res_group$Model, levels=res_group$Model)
  
  res_group$Type<-NA
  res_group$Type[grepl(paste0(type[k],'. Q'), res_group$Model)]<-paste0(type[k],'. Q')
  res_group$Type[grepl("Unstratified", res_group$Model)]<-"Unstratified"
  res_group$Type[grepl("All", res_group$Model)]<-"All"

  res_group_2<-res_group
  res_group_2$GWAS<-pheno[i]
  res_group_2$Population<-pop2[j]
  res_group_2$Method<-type[k]

  res_all<-rbind(res_all, res_group_2)
  
}
}
}

res_all<-rbind(res_all, EUR_res_all)
res_all$Type<-factor(res_all$Type, levels=c('cov. Q','cor. Q','Fst. Q','random. Q','Unstratified','All','EUR'))

# Plot the performance of stratifying by cov, cor, Fst
for(j in 1:length(pop2)){
for(i in 1:length(pheno)){
res_all_i_j<-res_all[(res_all$Population == pop2[j] | res_all$Population == 'EUR') & res_all$GWAS == pheno[i],]
res_plot<-list()
for(k in 1:length(type)){
    if(min(res_all_i_j$R - res_all_i_j$SE) < 0){
      y_lim<-c(min(res_all_i_j$R - res_all_i_j$SE),max(res_all_i_j$R + res_all_i_j$SE))
    } else {
      y_lim<-c(0,max(res_all_i_j$R + res_all_i_j$SE))
    }
  
  res_all_i_j_k<-res_all_i_j[res_all_i_j$Method == type[k] | is.na(res_all_i_j$Method),]
  res_all_i_j_k$Model<-factor(res_all_i_j_k$Model, levels=res_all_i_j_k$Model)

  res_plot[[k]]<-ggplot(res_all_i_j_k, aes(x=Model, y=R, colour=Type, shape=Population)) +
          geom_point(stat="identity", position=position_dodge( width=.25), size=5) +
          geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.25)) +
          scale_shape_manual(values=c(19, 15))+
          theme_half_open() +
          geom_hline(yintercept=0) +
          ylim(y_lim) +
          theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12), legend.position = 'none') +
          background_grid(major = 'y', minor = 'y') +
          labs(y="Correlation (SE)", title=paste0(pheno[i],': ', type[k],"-stratified (",pop2[j],")"))

}

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs_res.png'), width=5000, height=2000, res=300, units='px')
  print(plot_grid(plotlist=res_plot, labels = "AUTO", nrow=1))
  dev.off()
}
}

# Plot results comparing cov, cor, Fst, random groups, unstratified and all results combined
combined_res<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(pheno)){
  res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs.pred_eval.txt'))
  
  res_group<-res[!grepl('PredFile', res$Model),]
  res_group$Model<-gsub('_group', '', res_group$Model)
  res_group$Model<-gsub('unstratified', "Unstratified", res_group$Model)
  
  res_group$Model<-factor(res_group$Model, levels=res_group$Model)
  
  res_group$Type<-res_group$Model

  res_group_2<-res_group
  res_group_2$GWAS<-pheno[i]
  res_group_2$Population<-pop2[j]
  res_group_2$Method<-NA

  combined_res<-rbind(combined_res, res_group_2)
  
}
}

res_all_2<-res_all[!grepl('. Q', res_all$Model),]
res_all_2<-res_all_2[!grepl('Unstratified', res_all_2$Model),]
res_all_2$Model<-as.character(res_all_2$Model)
res_all_2$Method<-as.character(res_all_2$Method)
res_all_2$Type<-as.character(res_all_2$Type)
res_all_2$Model[res_all_2$Method == 'cov']<-'Unstratified + cov'
res_all_2$Model[res_all_2$Method == 'cor']<-'Unstratified + cor'
res_all_2$Model[res_all_2$Method == 'Fst']<-'Unstratified + Fst'
res_all_2$Model[res_all_2$Method == 'random']<-'Unstratified + random'
res_all_2$Type[res_all_2$Method == 'cov']<-'Unstratified + cov'
res_all_2$Type[res_all_2$Method == 'cor']<-'Unstratified + cor'
res_all_2$Type[res_all_2$Method == 'Fst']<-'Unstratified + Fst'
res_all_2$Type[res_all_2$Method == 'random']<-'Unstratified + random'

res_all_2<-rbind(res_all_2, combined_res)

res_all_2$Type<-factor(res_all_2$Type, levels=c('cov','cor','Fst','Unstratified','Unstratified + cov','Unstratified + cor','Unstratified + Fst','Unstratified + random','All','EUR'))
res_all_2$Model<-factor(res_all_2$Model, levels=c('cov','cor','Fst','Unstratified','Unstratified + cov','Unstratified + cor','Unstratified + Fst','Unstratified + random','All','EUR'))

for(j in 1:length(pop2)){
res_plot<-list()
for(i in 1:length(pheno)){
res_all_2_i_j<-res_all_2[(res_all_2$Population == pop2[j] | res_all_2$Population == 'EUR') & res_all_2$GWAS == pheno[i],]

  if(min(res_all_2_i_j$R - res_all_2_i_j$SE) < 0){
    y_lim<-c(min(res_all_2_i_j$R - res_all_2_i_j$SE),max(res_all_2_i_j$R + res_all_2_i_j$SE))
  } else {
    y_lim<-c(0,max(res_all_2_i_j$R + res_all_2_i_j$SE))
  }
  
  res_plot[[i]]<-ggplot(res_all_2_i_j, aes(x=Model, y=R, colour=Type, shape=Population)) +
          geom_point(stat="identity", position=position_dodge( width=.25), size=5) +
          geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.25)) +
          scale_shape_manual(values=c(19, 15))+
          theme_half_open() +
          geom_hline(yintercept=0) +
          ylim(y_lim) +
          theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12), legend.position = 'none') +
          background_grid(major = 'y', minor = 'y') +
          labs(y="Correlation (SE)", title=paste0(pheno[i],': ',pop2[j]))

}

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/prs_res_summary.png'), width=3000, height=12000, res=300, units='px')
  print(plot_grid(plotlist=res_plot, labels = "AUTO", ncol=2))
  dev.off()
}

# Now make a plot comparing unstratified and all methods for all populations
res_all_3<-res_all_2[res_all_2$Model == 'All' | res_all_2$Model == 'Unstratified' | res_all_2$Model == 'EUR',]
res_all_3$Model<-as.character(res_all_3$Model)
res_all_3$Model[res_all_3$Model == 'EUR']<-'Unstratified'
res_all_3$Model[res_all_3$Model == 'All']<-'CrossPop'
res_all_3$Model<-factor(res_all_3$Model, levels=c('Unstratified','CrossPop'))
res_all_3$Type<-res_all_3$Model

plot_list<-NULL
for(i in 1:length(pheno)){
  res_all_3_i<-res_all_3[res_all_3$GWAS == pheno[i],]
  
  if(min(res_all_3_i$R - res_all_3_i$SE) < 0){
    y_lim<-c(min(res_all_3_i$R - res_all_3_i$SE),max(res_all_3_i$R + res_all_3_i$SE))
  } else {
    y_lim<-c(0,max(res_all_3_i$R + res_all_3_i$SE))
  }

  plot_list[[i]]<-ggplot(res_all_3_i, aes(x=factor(Model), y=R, colour=Population, group=Population)) +
                            geom_point(stat="identity", position=position_dodge( width=.30)) +
                            geom_line(stat="identity", position=position_dodge( width=.30)) +
                            geom_errorbar(aes(ymin=R-SE, ymax=R+SE), width=.2, position=position_dodge(.30)) +
                            ylim(y_lim) +
                            theme_half_open() +
                            theme(axis.text.x = element_text(angle = 55, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.4, size=12)) +
                            background_grid(major = 'y', minor = 'y') +
                            labs(x='Method',y="Correlation (SE)", title=pheno[i])
  
}

png('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_res_summary_allpop.png', units='px', res=300, width=2000, height=8000)
  plot_grid(plotlist=plot_list, ncol = 2)
dev.off()

```
</details>

```{bash, eval=T, echo=F}
mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop_intern

pheno=$(cat /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt)

for pop1 in $(echo EUR); do
for pop2 in $(echo EAS AMR AFR SAS); do
for i in $(seq 1 $(wc -l /users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt | cut -d' ' -f1)); do

pheno_i=$(echo ${pheno} | cut -f ${i} -d ' ')

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/${pheno_i}/prs_res.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop_intern/${pop1}_${pop2}_${pheno_i}_prs_res.png

done

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/${pop1}_${pop2}/assoc/prs_res_summary.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop_intern/${pop1}_${pop2}_prs_res_summary.png

done
done

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_res_summary_allpop.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop_intern/
```

<details><summary>Show detailed figures</summary>

![EUR - AFR](Images/CrossPop_intern/EUR_AFR_prs_res_summary.png)
![EUR - EAS](Images/CrossPop_intern/EUR_EAS_prs_res_summary.png)
![EUR - AMR](Images/CrossPop_intern/EUR_AMR_prs_res_summary.png)
![EUR - SAS](Images/CrossPop_intern/EUR_SAS_prs_res_summary.png)

</details>

<details><summary>Show figure</summary>

![Summary](Images/CrossPop_intern/prs_res_summary_allpop.png)

</details>

***

### Portability

<details><summary>Show code</summary>
```{R, eval=F, echo=F}
# Read in the pred_comp results
library(data.table)
library(ggplot2)
library(cowplot)

type=c('cov','cor','Fst')
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')

# Read in European results and calculate portability improvement.
EUR_res_all<-NULL
pheno=read.table('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1
for(i in 1:length(pheno)){
EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_EUR_test/assoc/',pheno[i],'/prs.pred_eval.txt'))

EUR_res<-EUR_res[EUR_res$Model == 'All_group',]
EUR_res$GWAS<-pheno[i]
EUR_res_all<-rbind(EUR_res_all, EUR_res)
}

# Read in the cross pop results
res_comp_all<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(pheno)){
for(k in 1:length(type)){
  res_comp<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs_',type[k],'.pred_comp.txt'))
  
  res_comp<-res_comp[res_comp$Model_1 == 'All' & res_comp$Model_2 == 'unstratified', ]

  res_comp$GWAS<-pheno[i]
  res_comp$Population<-pop2[j]
  res_comp$Method<-type[k]
  
  res_comp$EUR_R<-EUR_res_all$R[EUR_res_all$GWAS == pheno[i]]
  
  res_comp_all<-rbind(res_comp_all, res_comp)
  
}
res_comp<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs.pred_comp.txt'))

res_comp<-res_comp[res_comp$Model_1 == 'All' & res_comp$Model_2 == 'unstratified', ]

res_comp$GWAS<-pheno[i]
res_comp$Population<-pop2[j]
res_comp$Method<-'All'

res_comp$EUR_R<-EUR_res_all$R[EUR_res_all$GWAS == pheno[i]]

res_comp_all<-rbind(res_comp_all, res_comp)

}
} 

# Calculate R2 difference
res_comp_all$Model_1_R2<-res_comp_all$Model_1_R^2
res_comp_all$Model_2_R2<-res_comp_all$Model_2_R^2
res_comp_all$R2_diff<-res_comp_all$Model_1_R2 - res_comp_all$Model_2_R2

# Calculate percentage difference
res_comp_all$R_perc_diff<-(res_comp_all$Model_1_R-res_comp_all$Model_2_R)/res_comp_all$Model_2_R
res_comp_all$R2_perc_diff<-(res_comp_all$Model_1_R2-res_comp_all$Model_2_R2)/res_comp_all$Model_2_R2

# Calculate portability
res_comp_all$Model_1_R_port<-res_comp_all$Model_1_R/res_comp_all$EUR_R
res_comp_all$Model_2_R_port<-res_comp_all$Model_2_R/res_comp_all$EUR_R
res_comp_all$R_port_perc_diff<-(res_comp_all$Model_1_R_port-res_comp_all$Model_2_R_port)/res_comp_all$Model_2_R_port
res_comp_all$R_port_diff<-(res_comp_all$Model_1_R_port-res_comp_all$Model_2_R_port)

res_comp_all$Model_1_R2_port<-res_comp_all$Model_1_R2/(res_comp_all$EUR_R^2)
res_comp_all$Model_2_R2_port<-res_comp_all$Model_2_R2/(res_comp_all$EUR_R^2)
res_comp_all$R2_port_perc_diff<-(res_comp_all$Model_1_R2_port-res_comp_all$Model_2_R2_port)/res_comp_all$Model_2_R2_port
res_comp_all$R2_port_diff<-(res_comp_all$Model_1_R2_port-res_comp_all$Model_2_R2_port)

write.csv(res_comp_all, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability.csv', row.names=F, quote=F)

res_comp_all_nocov<-res_comp_all[res_comp_all$Method != 'cov',]
sum(res_comp_all$R_diff_pval < 0.05 & res_comp_all$Method == 'All')
sum(res_comp_all$Method == 'All')

sum(res_comp_all$R_diff_pval < 0.05 & res_comp_all$Method == 'cor')
sum(res_comp_all$Method == 'cor')

sum(res_comp_all$R_diff_pval < 0.05 & res_comp_all$Method == 'Fst')
sum(res_comp_all$Method == 'Fst')

library(ggplot2)
library(cowplot)

gg_qqplot(res_comp_all$R_diff_pval) +
  theme_bw(base_size = 24) +
  theme(
    axis.ticks = element_line(size = 0.5),
    panel.grid = element_blank()
    # panel.grid = element_line(size = 0.5, color = "grey80")
  )

ggplot(res_comp_all, aes(sample=R_diff_pval))+stat_qq()

# Calculate the average portability across phenotypes
average_port<-NULL
type<-c('cov','cor','Fst','All')
for(k in 1:length(type)){
  average_port<-rbind(average_port, data.frame(Method=type[k],
                                               Population='All',
                                               R_perc_diff=mean(res_comp_all$R_perc_diff[res_comp_all$Method == type[k]]),
                                               R2_perc_diff=mean(res_comp_all$R2_perc_diff[res_comp_all$Method == type[k]]),
                                               R_port_perc_diff=mean(res_comp_all$R_port_perc_diff[res_comp_all$Method == type[k]]),
                                               R2_port_perc_diff=mean(res_comp_all$R2_port_perc_diff[res_comp_all$Method == type[k]]),
                                               R_port_diff=mean(res_comp_all$R_port_diff[res_comp_all$Method == type[k]]),
                                               R2_port_diff=mean(res_comp_all$R2_port_diff[res_comp_all$Method == type[k]]),
                                               All_R_port=mean(res_comp_all$Model_1_R_port[res_comp_all$Method == type[k]]),
                                               Unstratified_R_port=mean(res_comp_all$Model_2_R_port[res_comp_all$Method == type[k]]),
                                               All_R2_port=mean(res_comp_all$Model_1_R2_port[res_comp_all$Method == type[k]]),
                                               Unstratified_R2_port=mean(res_comp_all$Model_2_R2_port[res_comp_all$Method == type[k]])))
}

# Plot the average portability across phenotypes
pop2=c('EAS','AMR','AFR','SAS')
for(k in 1:length(type)){
for(j in 1:length(pop2)){
  average_port<-rbind(average_port, data.frame(Method=type[k],
                                               Population=pop2[j],
                                               R_perc_diff=mean(res_comp_all$R_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_perc_diff=mean(res_comp_all$R2_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R_port_perc_diff=mean(res_comp_all$R_port_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_port_perc_diff=mean(res_comp_all$R2_port_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R_port_diff=mean(res_comp_all$R_port_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               R2_port_diff=mean(res_comp_all$R2_port_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               All_R_port=mean(res_comp_all$Model_1_R_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               Unstratified_R_port=mean(res_comp_all$Model_2_R_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               All_R2_port=mean(res_comp_all$Model_1_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
                                               Unstratified_R2_port=mean(res_comp_all$Model_2_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]])))
}
}

write.csv(average_port, '/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability_average.csv', row.names=F, quote=F)

# Create plot showing results
port_plot_dat<-average_port[,c('Method', 'Population','All_R_port','Unstratified_R_port','All_R2_port','Unstratified_R2_port')]

port_plot_dat<-melt(port_plot_dat)
port_plot_dat$Metric<-gsub('.*_','',gsub('_port','', port_plot_dat$variable))
port_plot_dat$Method2<-NA
port_plot_dat$Method2[grepl('All', port_plot_dat$variable)]<-paste0(port_plot_dat$Method[grepl('All', port_plot_dat$variable)],'-stratified')
port_plot_dat$Method2[!grepl('All', port_plot_dat$variable)]<-'Unstratified'
port_plot_dat$Method<-port_plot_dat$Method2
port_plot_dat$Method2<-NULL
port_plot_dat<-port_plot_dat[!duplicated(port_plot_dat),]
port_plot_dat$Method<-gsub('cor','LDcor', port_plot_dat$Method)
port_plot_dat$Method<-gsub('All-stratified','All', port_plot_dat$Method)
  
port_plot_dat$Method<-factor(port_plot_dat$Method, levels=unique(c('Unstratified',port_plot_dat$Method)))
port_plot_dat$Population<-factor(port_plot_dat$Population, levels=c('AFR','AMR','EAS','SAS','All'))

library(ggplot2)
library(cowplot)
R_plot<-ggplot(port_plot_dat[port_plot_dat$Metric == 'R',], aes(x=Population, y=value, fill=Method)) +
  geom_bar(stat="identity", position=position_dodge()) +
  labs(y='Portability', title="Average portability (R)") +
  ylim(c(0,1)) +
  theme_minimal_hgrid()

R2_plot<-ggplot(port_plot_dat[port_plot_dat$Metric == 'R2',], aes(x=Population, y=value, fill=Method)) +
  geom_bar(stat="identity", position=position_dodge()) +
  labs(y='Portability', title="Average portability (R2)") +
  ylim(c(0,1)) +
  theme_minimal_hgrid()

png('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability_summary.png', units='px', res=300, width=2000, height=1400)
plot_grid(R_plot, R2_plot, labels = c('A', 'B'), ncol=1)
dev.off()

R2_plot_poster<-ggplot(port_plot_dat[port_plot_dat$Metric == 'R2' & (port_plot_dat$Method == 'LDcor-stratified' | port_plot_dat$Method == 'Fst-stratified' | port_plot_dat$Method == 'Unstratified'),], aes(x=Population, y=value, fill=Method)) +
  geom_bar(stat="identity", position=position_dodge()) +
  labs(y='Portability') +
  ylim(c(0,1)) +
  theme_minimal_hgrid()

png('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability_summary_poster.png', units='px', res=300, width=1450, height=800)
R2_plot_poster
dev.off()

```
</details>

```{bash, eval=T, echo=F}
mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop

cp /users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability_summary.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/CrossPop_intern/
```

<details><summary>Show average portability figure</summary>

![Average portability](Images/CrossPop_intern/prs_portability_summary.png)

</details>

<details><summary>Show portability table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability.csv")

res$R_diff_pval<-format(res$R_diff_pval, scientific = TRUE, digits = 3)
res[,!grepl("Model_1$|Model_2$|R_diff_pval|GWAS|Population|Method",names(res))]<-round(res[,!grepl("Model_1$|Model_2$|R_diff_pval|GWAS|Population|Method",names(res))], 3)

library(knitr)
kable(res, rownames = FALSE, caption='Portability')
```

</details>

<details><summary>Show average portability table</summary>

```{r, echo=F, eval=T, results='asis'}
res<-read.csv("/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/prs_portability_average.csv")

res[,!grepl("Population|Method",names(res))]<-round(res[,!grepl("Population|Method",names(res))], 3)

library(knitr)
kable(res, rownames = FALSE, caption='Average portability')
```

</details>

### LDweighted

```{R, eval=F, echo=T}
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')
pheno=read.table('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

for(pop_i in pop2){
  per_pT_plot<-list()
  best_pT_plot<-list()
  best_pT_abs_diff_plot<-list()
  best_pT_rel_diff_plot<-list()
  for(pheno_i in pheno){
    # Read in European results and calculate portability improvement.
    EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_EUR_test/assoc/',pheno[i],'/prs.assoc.txt'))
    
    # Read in pop2 results
    pop2_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_',pop_i,'/assoc/',pheno_i,'/prs_LDcorWeighted.assoc.txt'))

    pop2_res$Group<-gsub('.Pred.*','',gsub('Group_','',pop2_res$Predictor))
    pop2_res$pT<-as.numeric(gsub('e.','e-',gsub('.*SCORE_','',pop2_res$Predictor)))
    pop2_res$pT<-factor(pop2_res$pT, levels=unique(pop2_res$pT))
    
    library(ggplot2)
    library(cowplot)
    library(RColorBrewer)

    per_pT_plot[[pheno_i]]<-ggplot(pop2_res, aes(x=pT, y=BETA, fill=pT)) +
      geom_bar(stat="identity", position=position_dodge()) +
      geom_errorbar(aes(ymin=BETA-SE, ymax=BETA+SE), width=.2, position=position_dodge(.9)) +
      labs(x='pT', y="BETA (SE)", title=paste0(pop_i,': ',pheno_i)) +
      theme_half_open() +
      background_grid() +
      facet_wrap(~Group, nrow=1)
    
    best_pT<-NULL
    for(group_i in unique(pop2_res$Group)){
      best_pT<-rbind(best_pT, pop2_res[pop2_res$Group == group_i & pop2_res$BETA == max(pop2_res$BETA[pop2_res$Group == group_i]),])
    }
    
    best_pT_plot[[pheno_i]]<-ggplot(best_pT, aes(x=Group, y=BETA, fill=Group)) +
      geom_bar(stat="identity", position=position_dodge()) +
      geom_errorbar(aes(ymin=BETA-SE, ymax=BETA+SE), width=.2, position=position_dodge(.9)) +
      labs(x='pT', y="BETA (SE)", title=paste0(pop_i,': ',pheno_i)) +
      theme_half_open() +
      background_grid() +
      theme(axis.text.x = element_text(angle = 45,hjust=1))
    
    best_pT$Unstrat_abs_diff_R2<-best_pT$Obs_R2 - best_pT$Obs_R2[best_pT$Group == 'unstratified']
    best_pT$Unstrat_rel_diff_R2<-best_pT$Unstrat_abs_diff_R2/best_pT$Obs_R2[best_pT$Group == 'unstratified']
    
    best_pT_abs_diff_plot[[pheno_i]]<-ggplot(best_pT, aes(x=Group, y=Unstrat_abs_diff_R2, fill=Group)) +
      geom_bar(stat="identity", position=position_dodge()) +
      labs(x='pT', y="Absolute improvement R2", title=paste0(pop_i,': ',pheno_i)) +
      theme_half_open() +
      background_grid() +
      theme(axis.text.x = element_text(angle = 45,hjust=1))

    best_pT_rel_diff_plot[[pheno_i]]<-ggplot(best_pT, aes(x=Group, y=Unstrat_rel_diff_R2, fill=Group)) +
      geom_bar(stat="identity", position=position_dodge()) +
      labs(x='pT', y="Relative improvement R2", title=paste0(pop_i,': ',pheno_i)) +
      theme_half_open() +
      background_grid() +
      theme(axis.text.x = element_text(angle = 45,hjust=1))

  }

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_',pop_i,'/assoc/EUR_',pop_i,'_per_pT.png'), res=300, units='px', width=5000, height=12000)
    print(plot_grid(plotlist=per_pT_plot, ncol = 1))
  dev.off()

  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_',pop_i,'/assoc/EUR_',pop_i,'_best_pT.png'), res=300, units='px', width=4000, height=6000)
    print(plot_grid(plotlist=best_pT_plot, ncol = 2))
  dev.off()
  
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_',pop_i,'/assoc/EUR_',pop_i,'_best_pT_abs.png'), res=300, units='px', width=4000, height=6000)
    print(plot_grid(plotlist=best_pT_abs_diff_plot, ncol = 2))
  dev.off()
  
  png(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_',pop_i,'/assoc/EUR_',pop_i,'_best_pT_rel.png'), res=300, units='px', width=4000, height=6000)
    print(plot_grid(plotlist=best_pT_rel_diff_plot, ncol = 2))
  dev.off()
  
}

# The LD weighted approach isn't working well. Looking back at the evidence of interaction between LDcor and rG between ancestries, I am wondering whether this approach should even work.

```

***

# Discussion

The results using internal GWAS summary statistics are far less encouraging. I am not sure why, though the GIANT BMI and Height GWAS havve behaved strangely in the past. This, in combination with the observation that random stratification works almost as well, makes me want to rethink the approach. 

The use of the pT+clump approach, may increase the likelihood of random stratification leading to improved prediction by including many predictors in the elastic net. Perhaps use of a single score method, such as PRScs, SBayesR, or DBSLMM would be a cleaner approach.

Also, use of nested stratification by cov, cor and Fst may be a better approach, as this will retain 'good' SNPs, whilst allowing 'less good' SNPs to be assigned a large effect size. I.e, if a locus has a strong effect that has a high LD cor across populations, it will be used in the PRS, even when allowing SNPs with lower LD cor to be considered. Alternatively, a continuous weight could be applied, so BETAs are shrunk when low cross population LD cor. A cross-validation procedure could be used to determine optimal shrinkage.

Something that would be useful would be a comparison of GW-sig loci across populations. Does LD correlation across populations correlate with GWAS effect size? Examining LD correlation structures would also be useful. Are there LD correlation blocks. Perhaps low LD cor blocks should use the top SNP, assumming it is most likely to be causal. This reminds me that weights by likelihood of causality (PIP) is also a good idea for cross pop prediction. Perhaps LD cor and PIP can be used together.



***
***
***

# OLD
<details><summary>Show code</summary>
```{R, eval=F, echo=F}
# Read in the pred_comp results
library(data.table)
library(ggplot2)
library(cowplot)

type=c('cov','cor','Fst','random')
pop1='EUR'
pop2=c('EAS','AMR','AFR','SAS')
pheno=read.table('/users/k1806347/brc_scratch/Misc/Toshiki/Phenotypes/liang_pheno_names.txt', header=F)$V1

# Read in European results and calculate portability improvement.
EUR_res_all<-NULL
for(i in 1:length(pheno)){
EUR_res<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/EUR_EUR_test/assoc/',pheno[i],'/prs.pred_eval.txt'))

EUR_res<-EUR_res[EUR_res$Model == 'All_group',]
EUR_res$Model<-'EUR'
EUR_res$Type<-'EUR'
EUR_res$GWAS<-pheno[i]
EUR_res$Population<-'EUR'
EUR_res$Method<-NA
EUR_res_all<-rbind(EUR_res_all, EUR_res)
}

# Read in the cross pop results
res_comp_all<-NULL

for(j in 1:length(pop2)){
for(i in 1:length(pheno)){
for(k in 1:length(type)){
  res_comp<-fread(paste0('/users/k1806347/brc_scratch/Misc/Toshiki/CrossPop_intern/',pop1,'_',pop2[j],'/assoc/',pheno[i],'/prs_',type[k],'.pred_comp.txt'))
  
  res_comp<-res_comp[res_comp$Model_1 == 'All' & res_comp$Model_2 == 'unstratified', ]

  res_comp$GWAS<-pheno[i]
  res_comp$Population<-pop2[j]
  res_comp$Method<-type[k]
  
  res_comp$EUR_R<-EUR_res_all$R[EUR_res_all$GWAS == pheno[i]]
  
  res_comp_all<-rbind(res_comp_all, res_comp)
  
}
}
} 

# Calculate R2 difference
res_comp_all$Model_1_R2<-res_comp_all$Model_1_R^2
res_comp_all$Model_2_R2<-res_comp_all$Model_2_R^2
res_comp_all$R2_diff<-res_comp_all$Model_1_R2 - res_comp_all$Model_2_R2

# Calculate percentage difference
res_comp_all$R2_perc_diff<-(res_comp_all$Model_1_R2-res_comp_all$Model_2_R2)/res_comp_all$Model_2_R2

# Calculate portability
res_comp_all$Model_1_R2_port<-res_comp_all$Model_1_R2/(res_comp_all$EUR_R^2)
res_comp_all$Model_2_R2_port<-res_comp_all$Model_2_R2/(res_comp_all$EUR_R^2)
res_comp_all$R2_port_diff<-(res_comp_all$Model_1_R2_port-res_comp_all$Model_2_R2_port)

# Make table showing the mean portability, improvement etc
meta_res<-NULL
for(j in 1:length(pop2)){
for(k in 1:length(type)){
  meta_res<-rbind(meta_res, data.frame(Population=pop2[j],
                                       Type=type[k],
                                       Unstratified_R2=mean(res_comp_all$Model_2_R2[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           CrossPop_R2=mean(res_comp_all$Model_1_R2[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           R2_diff=mean(res_comp_all$R2_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           R2_perc_diff=mean(res_comp_all$R2_perc_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           Unstratified_port=mean(res_comp_all$Model_2_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           CrossPop_port=mean(res_comp_all$Model_1_R2_port[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]]),
           port_diff=mean(res_comp_all$R2_port_diff[res_comp_all$Method == type[k] & res_comp_all$Population == pop2[j]])))
}
}

res_comp_all_port_mod1<-res_comp_all[,c('Population','Method','Model_1_R2_port','R2_diff','R2_perc_diff')]
names(res_comp_all_port_mod1)<-c('Population','Method','Portability','R2_diff','R2_perc_diff')
res_comp_all_port_mod2<-res_comp_all[,c('Population','Method','Model_2_R2_port')]
names(res_comp_all_port_mod2)<-c('Population','Method','Portability')
res_comp_all_port_mod2$Method<-'Unstratified'
res_comp_all_port_mod2$R2_diff<-0
res_comp_all_port_mod2$R2_perc_diff<-0
res_comp_all_port<-rbind(res_comp_all_port_mod1, res_comp_all_port_mod2)

library(ggplot2)
library(cowplot)

ggplot(res_comp_all_port, aes(x=Population, y=Portability, fill=Method)) +
  geom_hline(yintercept=1) +
  geom_violin(position=position_dodge(0.7)) +
  geom_boxplot(width=0.1, position=position_dodge(0.7)) +
  theme_half_open() +
  background_grid(major = 'y', minor = 'y')

ggplot(res_comp_all_port, aes(x=Population, y=R2_diff, fill=Method)) +
  geom_hline(yintercept=0) +
  geom_violin(position=position_dodge(0.7)) +
  geom_boxplot(width=0.1, position=position_dodge(0.7)) +
  theme_half_open() +
  background_grid(major = 'y', minor = 'y') +
  labs(y="R-squared difference")

ggplot(res_comp_all_port, aes(x=Population, y=R2_perc_diff*100, fill=Method)) +
  geom_hline(yintercept=0) +
  geom_violin(width=0.5, position=position_dodge(0.7)) +
  geom_boxplot(width=0.5, position=position_dodge(0.7)) +
  theme_half_open() +
  background_grid(major = 'y', minor = 'y') +
  labs(y="R-squared difference (%)")

# The improvement is largest in AFR 18-36%. Improvements are generall small. This should be run including random stratification for comparison.

```
</details>

***