Target_sample_phenotype_prep.Rmd

---
title: "Training Sample Phenotype Preparation"
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, error=TRUE)
```

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

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

***

To evaluate genotype-based prediction approaches, I am using a range of phenotypes measured in the UK Biobank and TEDS samples.

# UKBB

Phenotype:

* Depression
* Intelligence
* BMI
* Height
* CAD (Coronary artery disease)
* T2D (Type II diabetes)
* Inflammatory Bowel Disorder (IBD)
* Rheumatoid arthritis (RheumArth)
* Multiple Slerosis (MultiScler)
* Prostate Cancer
* Breast Cancer

Below I provide a description of how each phenotype was derived and some descriptive statistics.

```{bash, echo=F, eval=F}
######
# Update file listing IDs after QC to match genetic data
######

module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)
keep<-fread(qc_indiv_list)
fam<-fread(paste0(UKBB_output,'/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam'))

keep_fam<-merge(keep, fam, by=c('V2'))
keep_new<-keep_fam[,c('V1.y','V2'), with=F]

fwrite(keep_new, paste0(Results_dir,'/ukb18177_glanville_post_qc_id_list.UpdateIDs.fam'), sep=' ', na='NA', col.names=F)

q()
n
```

<br/>

***

## Depression
The depression phenotype was shared with me by Kylie Glanville.

<details><summary>Preparation of Depression phenotype</summary>
```{bash, echo=T, eval=F}
######
# Modify depression phenotype file so IDs match genotypic data
######

module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)
pheno<-fread(paste0(UKBB_output, '/Phenotype/Depression/ever_depressed_pheno_final.txt'))
fam<-fread(paste0(UKBB_output,'/UKBB/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam')
names(fam)[1:2]<-c('FID','IID')

pheno_fam<-merge(pheno, fam, by=c('IID'))
pheno_new<-pheno_fam[,c('FID.y','IID','Depressed.Ever.Clean'), with=F]
names(pheno_new)[1]<-'FID'

fwrite(pheno_new, paste0(UKBB_output,'/Phenotype/Depression/ever_depressed_pheno_final.UpdateIDs.txt'), sep=' ', na='NA')

q()
n
```
</details>

<br/>

***

## Intelligence
The Intelligence phenotype is derived from the Fluid.intelligence.score.0.0 (f.20191.0.0) variable.

<details><summary>Preparation of Intelligence phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Extract intelligence variable, update IDs and calculate descriptives
#####
module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')

# Read in the field finder to find he code for the phenotype of interest.
field_finder<-read.table(paste0(ukbb_pheno_file,'/ukb18177_glanville_field_finder.txt'), stringsAsFactors=F)
field_finder[grepl('Fluid.intelligence', field_finder$V3),]
# The Fluid intelligence variable has a code of f.20191.0.0 and is called Fluid.intelligence.score.0.0

# Identify which column in the phenotype file contains the relevant code
library(data.table)
pheno_header<-names(data.frame(fread(cmd=paste0('head -n 1 ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt'))))
rel_cols<-c(1,which(grepl('f.20191.0.0',pheno_header)))

# Read in the IID and fluid intelligence columns
pheno<-data.frame(fread(cmd=paste0('cut -f ',paste(rel_cols,collapse=','),' ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt')))
names(pheno)<-c('IID','Fluid.intelligence.score')

# Remove individuals with missing phenotype data
pheno<-pheno[complete.cases(pheno),]

# Calculate descriptive statistics for the phenotype
library(e1071)
descriptives<-data.frame(Phenotype='Fluid.intelligence.score',
                        Code='f.20127.0.0',
                        Mean=mean(pheno$Fluid.intelligence.score),
                        SD=sd(pheno$Fluid.intelligence.score),
                        Min=min(pheno$Fluid.intelligence.score),
                        Max=max(pheno$Fluid.intelligence.score),
                        Skewness=skewness(pheno$Fluid.intelligence.score),
                        N=dim(pheno)[1])

system(paste0('mkdir ',UKBB_output,'/Phenotype/Intelligence'))
write.table(descriptives,paste0(UKBB_output,'/Phenotype/Intelligence/Descriptives.txt'), row.names=F, quote=F)

# Plot the distribution
png(paste0(UKBB_output,'/Phenotype/Intelligence/Histogram.png'))
hist(pheno$Fluid.intelligence.score)
dev.off()

# Write out the phenotype in PLINK format.
pheno_plink<-data.frame(FID=pheno$IID,
                        IID=pheno$IID,
                        Fluid.intelligence.score=pheno$Fluid.intelligence.score)
                        
fwrite(pheno_plink, paste0(UKBB_output,'/Phenotype/Intelligence/UKBB_Fluid.intelligence.score.pheno'), sep='\t')

q()
n
```
</details>

```{bash, echo=F, eval=T}
# Set variables
. /users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config

mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Intelligence

cp ${UKBB_output}/Phenotype/Intelligence/Histogram.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Intelligence/
```

<details><summary>Descriptives of Intelligence phenotype</summary>

<center>

![Figure 1: Distribution of Intelligence in UK Biobank](/users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Intelligence/Histogram.png)

\center

<br/>

***

```{r, echo=F, eval=T, results='asis'}
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
desc<-read.table(paste0(UKBB_output,'/Phenotype/Intelligence/Descriptives.txt'), header=T)
desc$Mean<-round(desc$Mean,3)
desc$SD<-round(desc$SD,3)
desc$Min<-round(desc$Min,3)
desc$Max<-round(desc$Max,3)
desc$Skewness<-round(desc$Skewness,3)
desc$N<-round(desc$N,3)

library(knitr)
kable(desc, row.names = FALSE, caption='Descriptives of Intelligence in UK Biobank')
```

</details>

<br/>

***

## Height
The Height phenotype is derived from the Standing.height.0.0 (f.50.0.0) variable.

<details><summary>Preparation of Height phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Extract Height variable, update IDs and calculate descriptives
#####
module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')

# Read in the field finder to find he code for the phenotype of interest.
field_finder<-read.table(paste0(ukbb_pheno_file,'/ukb18177_glanville_field_finder.txt'), stringsAsFactors=F)
field_finder[grepl('height', field_finder$V3),]
# The best height variable has a code of f.50.0.0 and is called Standing.height.0.0

# Identify which column in the phenotype file contains the relevant code
library(data.table)
pheno_header<-names(data.frame(fread(cmd=paste0('head -n 1 ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt'))))
rel_cols<-c(1,which(grepl('f.50.0.0',pheno_header)))

# Read in the IID and fluid intelligence columns
pheno<-data.frame(fread(cmd=paste0('cut -f ',paste(rel_cols,collapse=','),' ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt')))
names(pheno)<-c('IID','Height')

# Remove individuals with missing phenotype data
pheno<-pheno[complete.cases(pheno),]

# Calculate descriptive statistics for the phenotype
library(e1071)
descriptives<-data.frame(Phenotype='Standing.height.0.0',
                        Code='f.50.0.0',
                        Mean=mean(pheno$Height),
                        SD=sd(pheno$Height),
                        Min=min(pheno$Height),
                        Max=max(pheno$Height),
                        Skewness=skewness(pheno$Height),
                        N=dim(pheno)[1])

system(paste0('mkdir ',UKBB_output,'/Phenotype/Height'))
write.table(descriptives,paste0(UKBB_output,'/Phenotype/Height/Descriptives.txt'), row.names=F, quote=F)

# Plot the distribution
png(paste0(UKBB_output,'/Phenotype/Height/Histogram.png'))
hist(pheno$Height)
dev.off()

# Remove outliers
pheno_no_outlier<-pheno[(pheno$Height < mean(pheno$Height)+(3*sd(pheno$Height))) & (pheno$Height > mean(pheno$Height)-(3*sd(pheno$Height))),]

# Plot the distribution
png(paste0(UKBB_output,'/Phenotype/Height/Histogram_noOutlier.png'))
hist(pheno_no_outlier$Height)
dev.off()

# Write out the phenotype in PLINK format.
pheno_plink<-data.frame(FID=pheno_no_outlier$IID,
                        IID=pheno_no_outlier$IID,
                        Height=pheno_no_outlier$Height)
                        
fwrite(pheno_plink, paste0(UKBB_output,'/Phenotype/Height/UKBB_Height.pheno'), sep='\t')

q()
n
```
</details>

```{bash, echo=F, eval=T}
# Set variables
. /users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config

mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Height

cp ${UKBB_output}/Phenotype/Height/Histogram.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Height/
cp ${UKBB_output}/Phenotype/Height/Histogram_noOutlier.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Height/
```

<details><summary>Descriptives of Height phenotype</summary>

<center>

![Distribution of Height in UK Biobank](/users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Height/Histogram.png)

![Distribution of Height in UK Biobank after removal of outliers](/users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/Height/Histogram_noOutlier.png)

\center

<br/>

***

```{r, echo=F, eval=T, results='asis'}
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
desc<-read.table(paste0(UKBB_output,'/Phenotype/Height/Descriptives.txt'), header=T)
desc$Mean<-round(desc$Mean,3)
desc$SD<-round(desc$SD,3)
desc$Min<-round(desc$Min,3)
desc$Max<-round(desc$Max,3)
desc$Skewness<-round(desc$Skewness,3)
desc$N<-round(desc$N,3)

library(knitr)
kable(desc, row.names = FALSE, caption='Descriptives of Height in UK Biobank')
```

</details>

<br/>

***

## BMI
The BMI phenotype is derived from the Body.mass.index.BMI.0.0 (f.21001.0.0) variable.

<details><summary>Preparation of BMI phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Extract BMI variable, update IDs and calculate descriptives
#####

module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')

# Read in the field finder to find he code for the phenotype of interest.
field_finder<-read.table(paste0(ukbb_pheno_file,'/ukb18177_glanville_field_finder.txt'), stringsAsFactors=F)
field_finder[grepl('bmi', field_finder$V3),]
# The best height variable has a code of f.21001.0.0 and is called Body.mass.index.BMI.0.0

# Identify which column in the phenotype file contains the relevant code
library(data.table)
pheno_header<-names(data.frame(fread(cmd=paste0('head -n 1 ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt'))))
rel_cols<-c(1,which(grepl('f.21001.0.0',pheno_header)))

# Read in the IID and fluid intelligence columns
pheno<-data.frame(fread(cmd=paste0('cut -f ',paste(rel_cols,collapse=','),' ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt')))
names(pheno)<-c('IID','BMI')

# Remove individuals with missing phenotype data
pheno<-pheno[complete.cases(pheno),]

# Calculate descriptive statistics for the phenotype
library(e1071)
descriptives<-data.frame(Phenotype='Body.mass.index.BMI.0.0',
                        Code='f.50.0.0',
                        Mean=mean(pheno$BMI),
                        SD=sd(pheno$BMI),
                        Min=min(pheno$BMI),
                        Max=max(pheno$BMI),
                        Skewness=skewness(pheno$BMI),
                        N=dim(pheno)[1])

system(paste0('mkdir ',UKBB_output,'/Phenotype/BMI'))
write.table(descriptives,paste0(UKBB_output,'/Phenotype/BMI/Descriptives.txt'), row.names=F, quote=F)

# Plot the distribution
png(paste0(UKBB_output,'/Phenotype/BMI/Histogram.png'))
hist(pheno$BMI)
dev.off()

# Remove outliers
pheno_no_outlier<-pheno[(pheno$BMI < mean(pheno$BMI)+(3*sd(pheno$BMI))) & (pheno$BMI > mean(pheno$BMI)-(3*sd(pheno$BMI))),]

# Plot the distribution
png(paste0(UKBB_output,'/Phenotype/BMI/Histogram_noOutlier.png'))
hist(pheno_no_outlier$BMI)
dev.off()

# Write out the phenotype in PLINK format.
pheno_plink<-data.frame(FID=pheno_no_outlier$IID,
                        IID=pheno_no_outlier$IID,
                        BMI=pheno_no_outlier$BMI)
                        
fwrite(pheno_plink, paste0(UKBB_output,'/Phenotype/BMI/UKBB_BMI.pheno'), sep='\t')

q()
n
```

</details>

```{bash, echo=F, eval=T}
# Set variables
. /users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config

mkdir -p /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/BMI

cp ${UKBB_output}/Phenotype/BMI/Histogram.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/BMI/
cp ${UKBB_output}/Phenotype/BMI/Histogram_noOutlier.png /users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/BMI/
```

<details><summary>Descriptives of BMI phenotype</summary>

<center>

![Figure 4: Distribution of BMI in UK Biobank](/users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/BMI/Histogram.png)

![Distribution of BMI in UK Biobank after removal of outliers](/users/k1806347/brc_scratch/Software/MyGit/GenoPred/Images/Target_sample_phenotype_prep/BMI/Histogram_noOutlier.png)

\center

<br/>

***

```{r, echo=F, eval=T, results='asis'}
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
desc<-read.table(paste0(UKBB_output,'/Phenotype/BMI/Descriptives.txt'), header=T)
desc$Mean<-round(desc$Mean,3)
desc$SD<-round(desc$SD,3)
desc$Min<-round(desc$Min,3)
desc$Max<-round(desc$Max,3)
desc$Skewness<-round(desc$Skewness,3)
desc$N<-round(desc$N,3)

library(knitr)
kable(desc, row.names = FALSE, caption='Descriptives of BMI in UK Biobank')
```
</details>

<br/>

## CAD
T2D phenotypes were shared by Saskia Hagenaars and Anna F&uuml;rtjes.

<details><summary>Preparation of CAD phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Modify CAD phenotype file so IDs match genotypic data
#####

module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)
pheno<-fread(paste0(UKBB_output, '/Phenotype/CAD/cad_only_111119.txt'))
names(pheno)[1]<-'IID'
fam<-fread(paste0(UKBB_output,'/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam'))
names(fam)[1:2]<-c('FID','IID')

pheno_fam<-merge(pheno, fam, by=c('IID'))
pheno_new<-pheno_fam[,c('FID','IID','CAD_all'), with=F]

fwrite(pheno_new, paste0(UKBB_output,'/Phenotype/CAD/cad_only_111119.UpdateIDs.txt'), sep=' ', na='NA')

q()
n
```
</details>

<br/>

## Type-2 Diabetes (T2D)
T2D phenotypes were shared by Saskia Hagenaars and Anna F&uuml;rtjes.

<details><summary>Preparation of T2D phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Modify T2D phenotype file so IDs match genotypic data
#####

module add apps/R/3.6.0
R
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)
pheno<-fread(paste0(UKBB_output, '/Phenotype/T2D/t2d_only_111119.txt'))
names(pheno)[1]<-'IID'
fam<-fread(paste0(UKBB_output,'/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam'))
names(fam)[1:2]<-c('FID','IID')

pheno_fam<-merge(pheno, fam, by=c('IID'))
pheno_new<-pheno_fam[,c('FID','IID','t2d_all'), with=F]

fwrite(pheno_new, paste0(UKBB_output,'/Phenotype/T2D/t2d_only_111119.UpdateIDs.txt'), sep=' ', na='NA')

q()
n
```
</details>

<br/>

## Autoimmune disorders
Phenotypes were shared by Kylie Glanville.

<details><summary>Preparation of Autoimmune disorder phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Modify Autoimmune Disorder phenotype file so IDs match genotypic data
#####

module add apps/R/3.6.0
R
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)
pheno<-fread(paste0(UKBB_output, '/Phenotype/AutoImmune/new_autoimmune_all'))
pheno$FID<-NULL
fam<-fread(paste0(UKBB_output,'/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam'))
names(fam)[1:2]<-c('FID','IID')

pheno_fam<-merge(pheno, fam, by=c('IID'))

# Extract probable phnotypes only
pheno_new<-pheno_fam[,c('FID','IID',names(pheno_fam)[grepl('_prob', names(pheno_fam))]), with=F]

# Identify individuals that don't have any autoimmune disorders
pheno_new_control<-pheno_new[which(rowSums(pheno_new[,-1:-2]) == 0),]

# Create a pheno file for Crohns MultiScler RheuArth
pheno_name<-c('IBD','MultiScler','RheuArth')
pheno_col<-c('ibd_prob',"ms_prob","arthritis_prob")
for(i in 1:3){
  pheno_new_i_case<-pheno_new[pheno_new[[pheno_col[i]]] == 1,]
  pheno_clean_i<-rbind(pheno_new_i_case, pheno_new_control)
  pheno_clean_i<-pheno_clean_i[,c('FID','IID',pheno_col[i]), with=F]
  system(paste0('mkdir ',UKBB_output,'/Phenotype/',pheno_name[i]))
  fwrite(pheno_clean_i, paste0(UKBB_output,'/Phenotype/',pheno_name[i],'/UKBB.',pheno_name[i],'.txt'), sep=' ', na='NA', quote=F)

}

q()
n
```
</details>

<br/>

***

## Prostate and Breast cancer

<details><summary>Preparation of cancer phenotype</summary>
```{bash, echo=T, eval=F}
#####
# Extract Height variable, update IDs and calculate descriptives
#####
module add apps/R/3.6.0
R

source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')

# Read in the field finder to find the code for the phenotype of interest.
field_finder<-read.table(paste0(ukbb_pheno_file,'/ukb18177_glanville_field_finder.txt'), stringsAsFactors=F)
field_finder[grepl('Cancer|cancer', field_finder$V3),]
# The self reported cancer variable code is f.20001
# The sex variable code is f.31.0.0

# Identify which column in the phenotype file contains the relevant code
library(data.table)
pheno_header<-names(data.frame(fread(cmd=paste0('head -n 1 ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt'))))
rel_cols<-c(1,which(grepl('f.31.0.0|f.20001',pheno_header)))

# Read in the IID, sex, and cancer columns
pheno<-data.frame(fread(cmd=paste0('cut -f ',paste(rel_cols,collapse=','),' ',ukbb_pheno_file,'/ukb18177_glanville_phenotypes.txt')))
names(pheno)[1:2]<-c('IID','sex')
names(pheno)<-gsub('X..','',names(pheno))
names(pheno)<-gsub('\\.\\.','',names(pheno))

# Update IDs to match genetic data
fam<-fread(paste0(UKBB_output,'/Genotype/Harmonised/UKBB.w_hm3.QCd.AllSNP.chr22.fam'))
names(fam)[1:2]<-c('FID','IID')

pheno<-merge(fam[,c('FID','IID')],pheno, by=c('IID'))

# Identify individuals with prostate cancer ID '1044' in any column
prostate<-NULL
for(i in names(pheno)[-1:-3]){
  prostate<-c(prostate, which(pheno[[i]] == '1044'))
  prostate<-unique(prostate)
}

pheno$prostate_cancer<-NA
pheno$prostate_cancer[pheno$sex == 1]<-0
pheno$prostate_cancer[prostate]<-1

# Identify individuals with breast cancer ID '1002' in any column
breast<-NULL
for(i in names(pheno)[-1:-3]){
  breast<-c(breast, which(pheno[[i]] == '1002'))
  breast<-unique(breast)
}

pheno$breast_cancer<-NA
pheno$breast_cancer[pheno$sex == 0]<-0
pheno$breast_cancer[breast]<-1

# Write out the phenotype in PLINK format.
pheno_plink_prostate<-data.frame(FID=pheno$FID,
                        IID=pheno$IID,
                        prostate_cancer=pheno$prostate_cancer)
                        
pheno_plink_prostate<-pheno_plink_prostate[complete.cases(pheno_plink_prostate),]

pheno_plink_breast<-data.frame(FID=pheno$FID,
                        IID=pheno$IID,
                        breast_cancer=pheno$breast_cancer)

pheno_plink_breast<-pheno_plink_breast[complete.cases(pheno_plink_breast),]

dir.create(paste0(UKBB_output,'/Phenotype/Prostate_Cancer/'))
dir.create(paste0(UKBB_output,'/Phenotype/Breast_Cancer/'))

fwrite(pheno_plink_prostate, paste0(UKBB_output,'/Phenotype/Prostate_Cancer/UKBB_Prostate_Cancer.pheno'), sep='\t')
fwrite(pheno_plink_breast, paste0(UKBB_output,'/Phenotype/Breast_Cancer/UKBB_Breast_Cancer.pheno'), sep='\t')

q()
n
```
</details>

<br/>

***

## Create subset of UK Biobank individuals
The UK Biobank sample is very large making the computation of the various scores and derivation of the models computationally intensive. We do not need such large samples to compare PRS methods, and it would also be useful to have a non-overlapping subset of UKBiobank which can be used as a reference sample.

<details><summary>Subsetting phenotypic data</summary>
```{bash, echo=T, eval=F}
#####
# Idenitfy subset of UK biobank with phenotypic data and European ancestry
#####

module add apps/R/3.6.0
R
source('/users/k1806347/brc_scratch/Software/MyGit/GenoPred/config_used/Phenotype_prep.config')
library(data.table)

pheno<-c('Depression','Intelligence','BMI','Height','T2D','CAD','IBD','MultiScler','RheuArth','Breast_Cancer','Prostate_Cancer')
pheno_file<-c('ever_depressed_pheno_final.UpdateIDs.txt','UKBB_Fluid.intelligence.score.UpdateIDs.pheno','UKBB_BMI.score.UpdateIDs.pheno','UKBB_Height.score.UpdateIDs.pheno','t2d_only_111119.UpdateIDs.txt','cad_only_111119.UpdateIDs.txt','UKBB.IBD.txt','UKBB.MultiScler.txt','UKBB.RheuArth.txt','UKBB_Breast_Cancer.pheno','UKBB_Prostate_Cancer.pheno')

pheno_list<-list()
for(i in 1:length(pheno)){
  pheno_list[[pheno[i]]]<-fread(paste0(UKBB_output,'/Phenotype/',pheno[i],'/',pheno_file[i]))
}

pheno_all<-Reduce(function(...) merge(..., all=T, by=c('FID','IID')), pheno_list)

# We want to retain all individuals that are cases, and those that have complete data.
pheno_all_2<-pheno_all[complete.cases(pheno_all[,!grepl('cancer',names(pheno_all)),with=F]) | pheno_all$Depressed.Ever.Clean == 1 | pheno_all$t2d_all == 1 | pheno_all$CAD_all == 1 | pheno_all$ibd_prob == 1 | pheno_all$ms_prob == 1 | pheno_all$arthritis_prob == 1 | pheno_all$breast_cancer == 1 | pheno_all$prostate_cancer == 1,]

# Extract individuals with EUR ancestry
EUR<-fread(paste0(UKBB_output,'/Projected_PCs/AllAncestry/UKBB.w_hm3.AllAncestry.EUR.keep'))
pheno_all_2_EUR<-pheno_all_2[(pheno_all_2$IID %in% EUR$V2),]

# Extract individuals that survive QC
keep<-fread('/users/k1806347/brc_scratch/Analyses/PRS_comparison/UKBB_outcomes_for_prediction/ukb18177_glanville_post_qc_id_list.UpdateIDs.fam')
pheno_all_2_EUR<-pheno_all_2_EUR[(pheno_all_2_EUR$IID %in% keep$V2),]

# Calculate the number of individuals surviving genotype QC and ancestry check
sum(EUR$V2 %in% keep$V2) # 379437

# Remove individuals in the UKB reference
UKB_ref<-fread('/users/k1806347/brc_scratch/Data/UKBB/UKBB_ref/genotype/UKBB.noPheno.EUR.10K.chr22.fam')
pheno_all_2_EUR<-pheno_all_2_EUR[!(pheno_all_2_EUR$IID %in% UKB_ref$V2),]
# This seems circular as the phenotype file was used to select UKB individuals in the reference. This is due to the addition of breast cancer and prostate cancer phenotype. Using the original set of reference individuals losses <150 cases for each phenotype and saves recreating UKB reference files. 

# There are still >100k individuals. Write out files for each phenotype, retaining individuals with complete data, with a maximum of 50000 individuals for quantitative outcomes.
system(paste0('mkdir ',UKBB_output,'/Phenotype/PRS_comp_subset'))
pheno_N<-NULL
pheno_any<-NULL
for(i in 1:length(pheno)){
  pheno_all_2_EUR_i<-pheno_all_2_EUR[,c('FID','IID',names(pheno_all_2_EUR)[2+i]),with=F]
  pheno_all_2_EUR_i<-pheno_all_2_EUR_i[complete.cases(pheno_all_2_EUR_i),]
  if(sum(pheno_all_2_EUR_i[[3]] == 0 | pheno_all_2_EUR_i[[3]] == 1) == length(pheno_all_2_EUR_i[[3]])){
    pheno_all_2_EUR_i_case<-pheno_all_2_EUR_i[pheno_all_2_EUR_i[[3]] == 1,]
    pheno_all_2_EUR_i_case<-pheno_all_2_EUR_i_case[1:25000,]
    pheno_all_2_EUR_i_case<-pheno_all_2_EUR_i_case[complete.cases(pheno_all_2_EUR_i_case),]
    pheno_all_2_EUR_i_con<-pheno_all_2_EUR_i[pheno_all_2_EUR_i[[3]] == 0,]
    pheno_all_2_EUR_i_con<-pheno_all_2_EUR_i_con[1:(50000 - dim(pheno_all_2_EUR_i_case)[1]),]
    pheno_all_2_EUR_i<-rbind(pheno_all_2_EUR_i_case, pheno_all_2_EUR_i_con)
    
    pheno_N<-rbind(pheno_N,data.frame(Pheno=pheno[i], 
                                      N=length(pheno_all_2_EUR_i[[3]]),
                                      N_case=sum(pheno_all_2_EUR_i[[3]] == 1,na.rm=T),
                                      N_con=sum(pheno_all_2_EUR_i[[3]] == 0,na.rm=T)))
  } else {
    pheno_all_2_EUR_i<-pheno_all_2_EUR_i[1:50000,]
    pheno_N<-rbind(pheno_N,data.frame(Pheno=pheno[i], 
                                      N=length(pheno_all_2_EUR_i[[3]]),
                                      N_case=NA,
                                      N_con=NA))
  }
  
  pheno_any<-rbind(pheno_any,pheno_all_2_EUR_i[,c('FID','IID'),with=F])
  pheno_any<-pheno_any[!duplicated(pheno_any),]
  
  write.table(pheno_all_2_EUR_i, paste0(UKBB_output,'/Phenotype/PRS_comp_subset/UKBB.',pheno[i],'.txt'), col.names=T, row.names=F, quote=F)
}

write.csv(pheno_N, paste0(UKBB_output,'/Phenotype/PRS_comp_subset/Pheno_N.csv'), row.names=F, quote=F)

# Calculate the number of individuals across all samples
dim(pheno_any) # 104278

```
</details>

<br/>