Lab 6: GWAS

This is an R Markdown Notebook. When you execute a code within the notebook, the results appear beneath the code.

Place your cursor inside and press the green arrow on the right side of each chunk to execute it.

A. GWAS Catalog

4a. Click on the first SNPs on chromosomes 1, 2, and 3. Summarize the information in the following table.

Chr.	1st SNP	p-value	EFO mapping	GWAS trait	Study type
1
2
3

4b. Repeat steps 3 and 4 for Cancer.

Chr.	1st SNP	p-value	EFO mapping	GWAS trait	Study type
1
2
3

You will run PLINK on a unix shell (bash) and produce plots using R codes.

Note: if you get error during knitting (converting R markdown document into the formatted document such as Word or HTML), delete the code that produced the error and run Knit again.

Set CRAN repository: Go to Tools –> Global Options –> Packages –> Package Management –> Primary CRAN repository –> USA (IA) –> OK. You can also use this code:

options(repos=c(CRAN="https://mirror.las.iastate.edu/CRAN/"))

Part B. Software

10. By default, the working directory for R code chunks is the directory that contains the Rmd document. This could change when running different code chunks. For consistency, set the directory using the function knitr. Use the function getwd (get working directory) to check the directory.

knitr::opts_knit$set(root.dir = "R:/Biology/BIOL-4129-5129/hickljaa/Lab6_GWAS")
getwd()

## [1] "R:/Biology/BIOL-4129-5129/hickljaa/Lab6_GWAS"

Load the packages knitr and qqman. knitr is used to generate the markdown report and qqman to create Q-Q and Manhattan plots

# Can also load them from the packages list. 
library(knitr)
library(qqman)

## 

## For example usage please run: vignette('qqman')

## 

## Citation appreciated but not required:

## Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731, https://doi.org/10.21105/joss.00731.

## 

Part C. Data

14. Check the path of the working directory for PLINK

cd "R:/Biology/BIOL-4129-5129/hickljaa/Lab6_GWAS"
pwd = print working directory

## /r/Biology/BIOL-4129-5129/hickljaa/Lab6_GWAS

15. Check the working directory contains the files HapMap_3_r3_1.bed, HapMap_3_r3_1.bim, and HapMap_3_r3_1.fam.

find . -type f \( -name "*.bed" -o -name ".bim" -o -name "*.fam" \)

## ./HapMap_3_r3_1.bed
## ./HapMap_3_r3_1.fam
## ./HapMap_3_r3_10.bed
## ./HapMap_3_r3_10.fam
## ./HapMap_3_r3_11.bed
## ./HapMap_3_r3_11.fam
## ./HapMap_3_r3_12.bed
## ./HapMap_3_r3_12.fam
## ./HapMap_3_r3_13.bed
## ./HapMap_3_r3_13.fam
## ./HapMap_3_r3_4.bed
## ./HapMap_3_r3_4.fam
## ./HapMap_3_r3_5.bed
## ./HapMap_3_r3_5.fam
## ./HapMap_3_r3_6.bed
## ./HapMap_3_r3_6.fam
## ./HapMap_3_r3_7.bed
## ./HapMap_3_r3_7.fam
## ./HapMap_3_r3_8.bed
## ./HapMap_3_r3_8.fam
## ./HapMap_3_r3_9.bed
## ./HapMap_3_r3_9.fam
## ./HapMap_hwe_filter_step1.bed
## ./HapMap_hwe_filter_step1.fam

16. View the first few lines of the files using the commands head and -n followed by the number of lines (rows).

head -n 10 HapMap_3_r3_1.bim    # head --lines 10 HapMap_3_r3_1.fam works as well
head -n 10 HapMap_3_r3_1.fam

## 1    rs2185539   0   556738  T   C
## 1    rs11510103  0   557616  G   A
## 1    rs11240767  0   718814  T   C
## 1    rs3131972   0   742584  A   G
## 1    rs3131969   0   744045  A   G
## 1    rs1048488   0   750775  C   T
## 1    rs12562034  0   758311  A   G
## 1    rs12124819  0   766409  G   A
## 1    rs4040617   0   769185  G   A
## 1    rs2905036   0   782343  C   T
## 1328 NA06989 0 0 2 2
## 1377 NA11891 0 0 1 2
## 1349 NA11843 0 0 1 1
## 1330 NA12341 0 0 2 2
## 1444 NA12739 NA12748 NA12749 1 -9
## 1344 NA10850 0 NA12058 2 -9
## 1328 NA06984 0 0 1 2
## 1463 NA12877 NA12889 NA12890 1 -9
## 1418 NA12275 0 0 2 1
## 13291 NA06986 0 0 1 1

Part D. GWAS quality control

18. Check which SNPs are missing in a large proportion of samples (–geno) and which individuals have high rates of missing genotypes

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_1 -missing

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to plink.log.
## Options in effect:
##   --bfile HapMap_3_r3_1
##   --missing
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1457897 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997378.
## --missing: Sample missing data report written to plink.imiss, and variant-based
## missing data report written to plink.lmiss.
## Warning: 225 het. haploid genotypes present (see plink.hh ); many commands
## treat these as missing.

19. Load the plink.imiss (individuals) and plink.lmiss (SNPs) files into RStudio and check the number of individuals missing SNPs and the number of genotypes missing, respectively.

plink.imiss <- read.table(file="plink.imiss", header=TRUE)
plink.lmiss <- read.table(file="plink.lmiss", header=TRUE)

19a. Have a quick look at the plink.imiss data to make sure it contains the correct columns.

str(plink.imiss)

## 'data.frame':    165 obs. of  6 variables:
##  $ FID       : int  1328 1377 1349 1330 1444 1344 1328 1463 1418 13291 ...
##  $ IID       : chr  "NA06989" "NA11891" "NA11843" "NA12341" ...
##  $ MISS_PHENO: chr  "N" "N" "N" "N" ...
##  $ N_MISS    : int  4203 20787 1564 6218 29584 2631 9638 3788 5349 1758 ...
##  $ N_GENO    : int  1457897 1457897 1457897 1457897 1457897 1457897 1457897 1457897 1457897 1457897 ...
##  $ F_MISS    : num  0.00288 0.01426 0.00107 0.00426 0.02029 ...

19b. Check the plink.lmiss data as well.

str(plink.lmiss)

## 'data.frame':    1457897 obs. of  5 variables:
##  $ CHR   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ SNP   : chr  "rs2185539" "rs11510103" "rs11240767" "rs3131972" ...
##  $ N_MISS: int  0 4 0 0 0 1 0 1 0 0 ...
##  $ N_GENO: int  165 165 165 165 165 165 165 165 165 165 ...
##  $ F_MISS: num  0 0.0242 0 0 0 ...

20. Count the number of missing phenotypes (MISS_PHENO), missing SNPs (N_MISS), missing genotypes (N_GENO), and the minimum and maximum values for the proportion of missing SNPs (F_MISS) in the plink.imiss file.

20a. Use the table() function in R to count the number of missing phenotypes (MISS_PHENO).

table(plink.imiss['MISS_PHENO']=="Y") # MISS_PHENO = Y/N

## 
## FALSE  TRUE 
##   112    53

20b. Add the number of SNPs missing in each individual (N_MISS) using the colSums() function in R.

colSums(plink.imiss['N_MISS'])

## N_MISS 
## 630620

20c. Add the number of non-obligatory missing genotypes (N_GENO) in each individual using the colSums() function.

colSums(plink.imiss['N_GENO'])

##    N_GENO 
## 240553005

20d. Find the minimum and maximum values for the frequency of missing SNPs (F_MISS) using the min() and max() functions in R, respectively.

min(plink.imiss['F_MISS'])

## [1] 0.0004198

max(plink.imiss['F_MISS'])

## [1] 0.02029

- Summarize the information for plink.imiss data in this table.

Missingness	Header	No. missing
Per individual	FID
(.imiss)	IID
	MISS_PHENO	53
	N_MISS	630620
	N_GENO	240553005
	F_MISS	0.0004198-0.02029

21. Count the number of missing phenotypes (MISS_PHENO), missing SNPs (N_MISS), missing genotypes (N_GENO), and the minimum and maximum values for the proportion of missing SNPs (F_MISS) in the plink.imiss data. Record this in the table.

21a. Use the colSums() function in R to count the number of missing SNPs (N_MISS).

colSums(plink.lmiss['N_MISS'])

## N_MISS 
## 630620

21b. Add the number of non-obligatory missing genotypes (N_GENO) using the colSums() function.

colSums(plink.lmiss['N_GENO'])

##    N_GENO 
## 240553005

21c. Find the minimum and maximum values for the frequency of missing SNPs (F_MISS) using the min() and max() functions in R, respectively.

min(plink.lmiss['F_MISS'])

## [1] 0

max(plink.lmiss['F_MISS'])

## [1] 0.04848

- Summarize the information for plink.lmiss data in this table.

|——————|———-|—————————————-|—————–| | Missingness | Header | Description | No. missing | |==================|==========|========================================|=================| | Per SNP (.lmiss) | CHR | Chromosome number | | | | SNP | SNP identifier | | | | N_MISS | Number of individuals missing this SNP | 630620 | | | | (add the rows in this column) | | | | N_GENO | Number of missing genotypes | 240553005 | | | | (add the rows in this column) | | | | F_MISS | Proportion of sample missing for | | | this SNP (give a range) | 0-0.04848 | ——————————————————————————————

####22. Plot a histogram in R to visualize the distribution of individuals missing SNPs in the plink.imiss data.

hist(plink.imiss[,6],main="Histogram of individual missingness") # selects column 6, F_MISS

23. Plot a histogram in R to visualize the distribution of genotypes missing in the plink.lmiss data.

hist(plink.lmiss[,5],main="Histogram of SNP missingness")  # selects column 5, F_MISS

####24. Remove variants (SNPs) with missing call rates using –Geno with a threshold >2%

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_1 --geno 0.02 --make-bed --out HapMap_3_r3_4 

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_4.log.
## Options in effect:
##   --bfile HapMap_3_r3_1
##   --geno 0.02
##   --make-bed
##   --out HapMap_3_r3_4
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1457897 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997378.
## 27454 variants removed due to missing genotype data (--geno).
## 1430443 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_4.bed + HapMap_3_r3_4.bim + HapMap_3_r3_4.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.
## Warning: 225 het. haploid genotypes present (see HapMap_3_r3_4.hh ); many
## commands treat these as missing.

####25. Remove individuals with missingness using –Geno with a stringent threshold >0.02

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_4 --mind 0.02 --make-bed --out HapMap_3_r3_5

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_5.log.
## Options in effect:
##   --bfile HapMap_3_r3_4
##   --make-bed
##   --mind 0.02
##   --out HapMap_3_r3_5
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1430443 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## 0 people removed due to missing genotype data (--mind).
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997899.
## 1430443 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_5.bed + HapMap_3_r3_5.bim + HapMap_3_r3_5.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.
## Warning: 179 het. haploid genotypes present (see HapMap_3_r3_5.hh ); many
## commands treat these as missing.

– Read the information printed on the screen and summarize the number of SNPs and individuals removed due to missingness in this table.

Data filtration based on missingness Number of initial variants loaded	| ==========+ 1457897 |
Number of people in the original file (both cases and controls)	165 |
Number of variants that passed the filter	27454 |
Variants removed due to missing genotype	1430443
Number of phenotypes (cases and controls) that passed the filter	112 |
Number of phenotypes (cases and controls) missing	53 |

26. Check discrepancies (discordance) between assigned sex and the genetic sex of individuals estimated from the F coefficient.

26a. Check sex discrepancy using the funtion –check-sex. By default, males have F estimates > 0.8 and females < 0.2.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_5 --check-sex

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to plink.log.
## Options in effect:
##   --bfile HapMap_3_r3_5
##   --check-sex
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1430443 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997899.
## 1430443 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --check-sex: 23424 Xchr and 0 Ychr variant(s) scanned, 1 problem detected.
## Report written to plink.sexcheck .
## Warning: 179 het. haploid genotypes present (see plink.hh ); many commands
## treat these as missing.

26c. Plot histograms to visualize the sex-check results.

gender <- read.table("plink.sexcheck", header=T,as.is=T)
hist(gender[,6],main="Gender", xlab="F")

male=subset(gender, gender$PEDSEX==1)
hist(male[,6],main="Men",xlab="F")

female=subset(gender, gender$PEDSEX==2)
hist(female[,6],main="Women",xlab="F")

How many individuals have sex discrepancy? What was the frequency??

Only one individual, about 0.9 frequency #### 26d. Remove individuals with sex discrepancy

grep "PROBLEM" plink.sexcheck| awk '{print$1,$2}'> sex_discrepancy.txt
"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_6.log.
## Options in effect:
##   --bfile HapMap_3_r3_5
##   --make-bed
##   --out HapMap_3_r3_6
##   --remove sex_discrepancy.txt
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1430443 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## --remove: 164 people remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 52 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate in remaining samples is 0.99798.
## 1430443 variants and 164 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (52 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_6.bed + HapMap_3_r3_6.bim + HapMap_3_r3_6.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.
## Warning: 179 het. haploid genotypes present (see HapMap_3_r3_6.hh ); many
## commands treat these as missing.

26e. Impute (replace) sex assignments based on the SNP data using the function –impute-sex.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_5 --impute-sex --make-bed --out HapMap_3_r3_6

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_6.log.
## Options in effect:
##   --bfile HapMap_3_r3_5
##   --impute-sex
##   --make-bed
##   --out HapMap_3_r3_6
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1430443 variants loaded from .bim file.
## 165 people (80 males, 85 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997899.
## 1430443 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --impute-sex: 23424 Xchr and 0 Ychr variant(s) scanned, all sexes imputed.
## Report written to HapMap_3_r3_6.sexcheck .
## --make-bed to HapMap_3_r3_6.bed + HapMap_3_r3_6.bim + HapMap_3_r3_6.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.
## Warning: 179 het. haploid genotypes present (see HapMap_3_r3_6.hh ); many
## commands treat these as missing.

27. Remove SNPs with low minor allele frequency (MAF) because they are rare and prone to genotyping error.

27a. Select autosomal SNPs (i.e., chromosomes 1 to 22).

awk '{ if ($1 >= 1 && $1 <= 22) print $2 }' HapMap_3_r3_4.bim > snp_1_22.txt
"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_7.log.
## Options in effect:
##   --bfile HapMap_3_r3_6
##   --extract snp_1_22.txt
##   --make-bed
##   --out HapMap_3_r3_7
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1430443 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## --extract: 1398544 variants remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99806.
## 1398544 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_7.bed + HapMap_3_r3_7.bim + HapMap_3_r3_7.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

27b. Calculate minor allele frequencies (MAFs) for each SNP using the function –freq.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_7 --freq --out MAF_check

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to MAF_check.log.
## Options in effect:
##   --bfile HapMap_3_r3_7
##   --freq
##   --out MAF_check
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1398544 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99806.
## --freq: Allele frequencies (founders only) written to MAF_check.frq .

27c. Look at the head of the MAF_check file.

head MAF_check.frq

##  CHR         SNP   A1   A2          MAF  NCHROBS
##    1   rs2185539    T    C            0      224
##    1  rs11240767    T    C            0      224
##    1   rs3131972    A    G       0.1652      224
##    1   rs3131969    A    G       0.1339      224
##    1   rs1048488    C    T       0.1667      222
##    1  rs12562034    A    G       0.1027      224
##    1  rs12124819    G    A       0.2902      224
##    1   rs4040617    G    A       0.1295      224
##    1   rs2905036    C    T            0      224

27d. d. Plot the distribution of minor allele frequencies (MAFs).

maf_freq <- read.table("MAF_check.frq", header =TRUE, as.is=T)
hist(maf_freq[,5],main = "MAF distribution", xlab = "MAF")

27e. Remove SNPs with MAF <0.05.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_8.log.
## Options in effect:
##   --bfile HapMap_3_r3_7
##   --maf 0.05
##   --make-bed
##   --out HapMap_3_r3_8
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1398544 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99806.
## 325318 variants removed due to minor allele threshold(s)
## (--maf/--max-maf/--mac/--max-mac).
## 1073226 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_8.bed + HapMap_3_r3_8.bim + HapMap_3_r3_8.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

- Read the information printed on the screen and record the number of autosomal SNPs, individuals, cases, and controls that passed filtration based on the MAF threshold.

Number of variants loaded	1398544
Number of people loaded (both cases and controls)	165 |
Number of variants removed due to MAF threshold	325318 |
Number of variants that passed the filter	10732226|
Number of cases that passed the filter	56 |
Number of controls that passed the filter	56 |

28. Remove alleles that are not in Hardy–Weinberg equilibrium (HWE).

28a. Generate a list of genotype counts and HW tests for each SNP using the option –hardy.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_8 --hardy

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to plink.log.
## Options in effect:
##   --bfile HapMap_3_r3_8
##   --hardy
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998047.
## --hardy: Writing Hardy-Weinberg report (founders only) to plink.hwe ... 0%0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

28b. Plot the distribution (histogram) of HWE p-values.

hwe<-read.table (file="plink.hwe", header=TRUE)
hist(hwe[,9],main="Histogram HWE")

28c. Select SNPs with HWE p-value <1e−5. We will use this to zoom in on SNPs that deviate strongly.

awk '{ if ($9 <0.00001) print $0 }' plink.hwe > plinkzoomhwe.hwe

####28d. Look closely at the SNPs that deviate strongly.

hwe_zoom<-read.table (file="plinkzoomhwe.hwe", header=TRUE)
hist(hwe_zoom[,9],main="Histogram HWE: SNPS that deviate strongly")

29. Filter controls and cases based on the HWE test.

29a. First, filter controls using a stringent HWE threshold (p-value cut off of 1e-6).

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_8 --hwe 1e-6 --make-bed --out HapMap_hwe_filter_step1

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_hwe_filter_step1.log.
## Options in effect:
##   --bfile HapMap_3_r3_8
##   --hwe 1e-6
##   --make-bed
##   --out HapMap_hwe_filter_step1
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998047.
## --hwe: 0 variants removed due to Hardy-Weinberg exact test.
## 1073226 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_hwe_filter_step1.bed + HapMap_hwe_filter_step1.bim +
## HapMap_hwe_filter_step1.fam ... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

29b. Then, filter the cases using a stringent threshold p-value of 1e-6.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_hwe_filter_step1 --hwe 1e-10 --hwe-all --make-bed --out HapMap_3_r3_9

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_9.log.
## Options in effect:
##   --bfile HapMap_hwe_filter_step1
##   --hwe 1e-10
##   --hwe-all
##   --make-bed
##   --out HapMap_3_r3_9
## 
## Note: --hwe-all flag deprecated.  Use "--hwe include-nonctrl".
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998047.
## --hwe: 0 variants removed due to Hardy-Weinberg exact test.
## 1073226 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_9.bed + HapMap_3_r3_9.bim + HapMap_3_r3_9.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

30. Exclude individuals that deviate from the mean heterozygosity rate.

30a. Check heterozygosity on SNPs that are not highly correlated. A file called “inversion.txt” is needed for thsi step. Check your folder for for this file.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_9 --exclude inversion.txt --range --indep-pairwise 50 5 0.2 --out indepSNP

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to indepSNP.log.
## Options in effect:
##   --bfile HapMap_3_r3_9
##   --exclude inversion.txt
##   --indep-pairwise 50 5 0.2
##   --out indepSNP
##   --range
## 
## Note: --range flag deprecated.  Use e.g. "--extract range <filename>".
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## --exclude range: 9893 variants excluded.
## --exclude range: 1063333 variants remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998042.
## 1063333 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## 1%2%3%4%5%6%6%7%8%9%10%11%12%13%13%14%15%15%16%17%18%19%20%21%22%22%23%24%25%26%27%28%29%30%31%31%32%33%34%35%36%37%38%38%39%40%40%41%42%43%44%45%46%47%47%48%49%50%51%52%53%54%55%56%56%57%58%59%60%61%62%63%63%64%65%65%66%67%68%69%70%71%72%72%73%74%75%76%77%78%79%80%81%81%82%83%84%85%86%87%88%88%89%90%90%91%92%93%94%95%96%97%97%98%99%Pruned 79323 variants from chromosome 1, leaving 8397.
## 1%1%2%3%4%5%6%6%7%8%9%10%11%11%12%13%14%15%16%16%17%18%19%20%21%21%22%23%24%25%26%26%27%28%29%30%31%31%32%33%34%35%36%36%37%38%39%40%41%41%42%43%44%45%46%46%47%48%49%50%51%52%53%53%54%55%56%57%58%58%59%60%61%62%63%63%64%65%66%67%68%68%69%70%71%72%73%73%74%75%76%77%78%78%79%80%81%82%83%83%84%85%86%87%88%88%89%90%91%92%93%93%94%95%96%97%98%98%99%Pruned 81833 variants from chromosome 2, leaving 7965.
## 1%2%2%3%4%5%6%6%7%8%9%10%11%12%13%13%14%15%16%17%17%18%19%20%21%22%23%24%24%25%26%27%28%29%30%31%31%32%33%34%35%35%36%37%38%39%40%41%42%42%43%44%45%46%46%47%48%49%50%51%52%53%53%54%55%56%57%58%59%60%60%61%62%63%64%64%65%66%67%68%69%70%71%71%72%73%74%75%75%76%77%78%78%79%80%81%82%82%83%84%85%86%87%88%89%89%90%91%92%93%93%94%95%96%97%98%99%Pruned 67905 variants from chromosome 3, leaving 6957.
## 1%2%3%3%4%5%6%7%7%8%9%10%11%11%12%13%14%15%15%16%17%18%19%19%20%21%22%23%23%24%25%26%27%28%29%30%31%32%33%34%34%35%36%37%38%38%39%40%41%42%42%43%44%45%46%46%47%48%49%50%50%51%52%53%54%54%55%56%57%58%59%60%61%62%63%64%65%65%66%67%68%69%69%70%71%72%73%73%74%75%76%77%77%78%79%80%81%81%82%83%84%85%85%86%87%88%89%90%91%92%93%94%95%96%96%97%98%99%Pruned 59914 variants from chromosome 4, leaving 6215.
## 0%1%1%2%2%3%3%4%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%29%30%30%31%31%32%32%33%33%34%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%58%59%59%60%60%61%61%62%62%63%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%88%89%89%90%90%91%91%92%92%93%93%94%95%96%97%98%99%Pruned 62681 variants from chromosome 5, leaving 6336.
## 1%1%2%3%3%4%5%5%6%7%7%8%9%9%10%11%11%12%13%13%14%15%15%16%17%17%18%26%27%28%29%30%31%32%33%33%34%35%35%36%37%37%38%39%39%40%41%41%42%43%43%44%45%45%46%47%47%48%49%49%50%51%51%52%53%53%54%55%55%56%57%57%58%59%59%60%61%61%62%63%63%64%65%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%Pruned 59440 variants from chromosome 6, leaving 6092.
## 1%2%3%4%4%5%6%7%8%9%9%10%11%12%13%14%14%15%16%17%18%19%19%20%21%22%23%24%25%26%27%28%28%29%30%31%32%33%33%34%35%36%37%38%38%39%40%41%42%43%43%44%45%46%47%48%48%49%50%51%52%52%53%54%55%56%57%57%58%59%60%61%62%62%63%64%65%66%67%67%68%69%70%71%72%72%73%74%75%76%76%77%78%79%80%81%81%82%83%84%85%86%86%87%88%89%90%91%91%92%93%94%95%96%96%97%98%99%Pruned 53506 variants from chromosome 7, leaving 5598.
## 1%2%3%3%4%5%6%6%7%8%9%13%14%14%15%16%17%18%19%20%21%21%22%23%24%25%26%27%28%28%29%30%31%32%33%34%35%35%36%37%38%39%40%41%42%42%43%44%45%46%47%48%49%49%50%51%52%53%54%55%56%56%57%58%59%60%61%62%63%63%64%65%66%67%68%69%70%70%71%72%73%74%75%76%77%77%78%79%80%81%82%83%84%84%85%86%87%88%89%90%91%91%92%93%94%95%96%97%98%98%99%Pruned 51227 variants from chromosome 8, leaving 5025.
## 1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%52%53%54%55%55%56%57%58%58%59%60%61%61%62%63%64%64%65%66%67%67%68%69%70%70%71%72%73%73%74%75%76%76%77%78%79%79%80%81%82%82%83%84%85%85%86%87%88%88%89%90%91%91%92%93%94%94%95%96%97%97%98%99%Pruned 45349 variants from chromosome 9, leaving 4983.
## 1%2%3%4%5%6%6%7%8%8%9%10%11%12%13%14%15%16%17%17%18%19%20%21%22%23%24%25%26%26%27%28%28%29%30%31%32%33%34%35%36%37%37%38%39%40%41%42%43%44%45%46%46%47%48%48%49%50%51%52%53%54%55%56%57%57%58%59%60%61%62%63%64%65%66%66%67%68%68%69%70%71%72%73%74%75%76%77%77%78%79%80%81%82%83%84%85%86%86%87%88%88%89%90%91%92%93%94%95%96%97%97%98%99%Pruned 51843 variants from chromosome 10, leaving 5382.
## 1%1%2%3%4%5%6%7%8%8%9%10%11%12%13%14%15%15%16%17%18%19%20%21%22%22%23%24%25%26%27%28%29%29%30%31%31%32%33%34%35%36%36%37%38%38%39%40%41%42%43%44%45%45%46%47%48%49%50%51%52%52%53%54%55%56%57%58%59%59%60%61%62%63%64%65%66%66%67%68%68%69%70%71%72%73%73%74%75%75%76%77%78%79%80%80%81%82%82%83%84%85%86%87%88%89%89%90%91%92%93%94%95%96%96%97%98%99%Pruned 50263 variants from chromosome 11, leaving 5021.
## 0%1%2%3%4%5%6%7%7%8%8%9%10%11%12%13%14%15%15%16%17%18%19%20%21%22%23%23%24%25%26%27%28%29%30%30%31%31%32%33%34%35%36%37%38%38%39%40%41%42%43%44%45%46%46%47%48%49%50%51%52%53%53%54%54%55%56%57%58%59%60%61%61%62%62%63%64%65%66%67%68%69%69%70%71%72%73%74%75%76%76%77%77%78%79%80%81%82%83%84%84%85%85%86%87%88%89%90%91%92%92%93%94%95%96%97%98%99%Pruned 47315 variants from chromosome 12, leaving 5250.
## 1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%40%41%41%42%42%43%43%44%44%45%45%46%46%47%47%48%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%90%91%91%92%92%93%93%94%94%95%95%96%96%97%97%98%98%99%Pruned 36460 variants from chromosome 13, leaving 4030.
## 1%2%3%3%4%5%6%6%7%8%9%10%10%11%12%13%13%14%15%16%17%17%18%19%20%20%21%22%23%24%25%26%27%27%28%29%30%31%32%33%34%34%35%36%37%38%39%40%41%41%42%43%44%45%46%47%48%48%49%50%51%52%53%54%55%55%56%57%58%59%60%61%62%62%63%64%65%65%66%67%68%69%69%70%71%72%72%73%74%75%76%76%77%78%79%79%80%81%82%83%83%84%85%86%86%87%88%89%90%91%92%93%93%94%95%96%97%98%99%Pruned 31645 variants from chromosome 14, leaving 3499.
## 0%1%1%2%2%3%3%4%4%5%5%6%6%7%7%8%8%9%9%10%10%11%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%62%63%63%64%64%65%65%66%66%67%67%68%68%69%69%70%70%71%71%72%72%73%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%Pruned 28604 variants from chromosome 15, leaving 3404.
## 1%2%3%4%4%5%6%7%7%8%9%10%10%11%12%13%14%15%16%17%18%18%19%20%21%21%22%23%24%25%26%27%28%29%29%30%31%32%32%33%34%35%36%37%38%39%40%40%41%42%43%43%44%45%46%47%48%49%50%51%51%52%53%54%54%55%56%57%57%58%59%60%61%62%62%63%64%65%65%66%67%68%68%69%70%71%72%73%74%75%76%76%77%78%79%79%80%81%82%83%84%85%86%87%87%88%89%90%90%91%92%93%94%95%96%97%98%98%99%Pruned 30134 variants from chromosome 16, leaving 3685.
## 1%2%3%4%4%5%6%7%8%8%9%10%11%12%12%13%14%15%16%16%17%18%19%20%21%21%22%23%24%25%25%26%27%28%29%29%30%31%32%33%33%34%35%36%37%38%39%40%41%42%42%43%44%45%46%46%47%48%49%50%50%51%55%55%56%57%58%59%60%60%61%62%63%64%64%65%66%67%68%68%69%70%71%72%72%73%74%75%76%77%78%79%80%81%81%82%83%84%85%85%86%87%88%89%89%90%91%92%93%94%95%96%97%98%98%99%Pruned 24471 variants from chromosome 17, leaving 3366.
## 1%2%3%4%5%6%6%7%8%9%9%10%11%12%12%13%14%15%16%17%18%19%20%21%22%22%23%24%25%25%26%27%28%29%30%31%32%33%34%35%35%36%37%38%38%39%40%41%41%42%43%44%45%46%47%48%48%49%50%51%51%52%53%54%54%55%56%57%58%59%60%61%62%63%64%64%65%66%67%67%68%69%70%70%71%72%73%74%75%76%77%77%78%79%80%80%81%82%83%83%84%85%86%87%88%89%90%91%92%93%93%94%95%96%96%97%98%99%Pruned 28242 variants from chromosome 18, leaving 3413.
## 1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%15%16%17%18%18%19%20%21%21%22%23%24%24%25%26%27%27%28%29%30%30%31%32%33%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%52%53%54%55%55%56%57%58%58%59%60%61%61%62%63%64%64%65%66%67%67%68%69%70%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%86%87%88%89%89%90%91%92%92%93%94%95%95%96%97%98%98%99%Pruned 16854 variants from chromosome 19, leaving 2808.
## 1%2%3%4%5%6%7%7%8%8%9%10%11%12%13%14%15%16%17%17%18%18%19%20%21%22%23%24%25%26%26%27%27%28%29%30%31%32%33%34%35%36%36%37%37%38%39%40%41%42%43%44%45%46%46%47%47%48%49%50%51%52%53%54%55%55%56%56%57%58%59%60%61%62%63%64%65%65%66%66%67%68%69%70%71%72%73%74%75%75%76%77%78%79%80%81%82%83%84%84%85%85%86%87%88%89%90%91%92%93%94%94%95%95%96%97%98%99%Pruned 24897 variants from chromosome 20, leaving 3051.
## 1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%22%23%23%24%24%25%25%26%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%50%51%51%52%52%53%53%54%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%77%78%78%79%79%80%80%81%81%82%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%Pruned 13769 variants from chromosome 21, leaving 1713.
## 1%2%3%4%5%6%7%8%9%10%11%12%12%13%13%14%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%28%29%29%30%30%31%32%33%34%35%36%37%38%39%40%41%42%43%43%44%44%45%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%59%60%60%61%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%75%76%76%77%77%78%79%80%81%82%83%84%85%86%87%88%89%90%90%91%91%92%92%93%94%95%96%97%98%99%Pruned 13514 variants from chromosome 22, leaving 1954.
## Pruning complete.  959189 of 1063333 variants removed.
## Writing...Marker lists written to indepSNP.prune.in and indepSNP.prune.out .

30b. Prune the data set. Pruning is reducing the size of the data by removing irrelevant, redundant, or low-value data. In this exercise, we prune the data so that the SNPs that remain are not associated with each other, i.e., they are in linkage equilibrium.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_9 --extract indepSNP.prune.in --het --out R_check

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to R_check.log.
## Options in effect:
##   --bfile HapMap_3_r3_9
##   --extract indepSNP.prune.in
##   --het
##   --out R_check
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## --extract: 104144 variants remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 112 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998036.
## 104144 variants and 165 people pass filters and QC.
## Among remaining phenotypes, 56 are cases and 56 are controls.  (53 phenotypes
## are missing.)
## --het: 104144 variants scanned, report written to R_check.het .

30c. Plot the distribution of the rate of heterozygosity in R.

het <- read.table("R_check.het", head=TRUE)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
hist(het$HET_RATE, xlab="Heterozygosity Rate", ylab="Frequency", main= "Heterozygosity Rate")

30d. Produce a list of individuals that deviate more than 3 standard deviations from the mean heterozygosity rate in R. The output of this script is fail-het-qc.txt.

het <- read.table("R_check.het", head=TRUE)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
het_fail = subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE)));
het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE);
write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)

30e. Format the fail-het-qc.txt file at the bash to make it compatible with PLINK by removing all quotation marks and selecting only the first two columns.

sed 's/"// g' fail-het-qc.txt | awk '{print$1, $2}'> het_fail_ind.txt

30f. Now, remove heterozygosity rate outliers, i.e., those that deviate >3 sd from the mean heterozygosity rate.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_10.log.
## Options in effect:
##   --bfile HapMap_3_r3_9
##   --make-bed
##   --out HapMap_3_r3_10
##   --remove het_fail_ind.txt
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 165 people (81 males, 84 females) loaded from .fam.
## 112 phenotype values loaded from .fam.
## --remove: 163 people remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate in remaining samples is 0.998071.
## 1073226 variants and 163 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.  (53 phenotypes
## are missing.)
## --make-bed to HapMap_3_r3_10.bed + HapMap_3_r3_10.bim + HapMap_3_r3_10.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

####- Read the information printed above and record the number of autosomal SNPs, individuals, cases and controls that passed filters and QC.

Table. Filtration based on standard deviation.

|——————————————-|—–| | Filtration based on standard deviation | | |——————————————-|—–| | Number of variants that passed the filter | 1073226 | | Number of cases that passed the filter | 55 | | Number of controls that passed the filter | 55 | |——————————————-|—–|

31. Remove genetically related individuals (parent-offspring relationship). PLINK estimates genetic relatedness based on identity by descent (IBD). IBD is a situation where two genomes inherited identical alleles from a common ancestor in the past as opposed to identity by state (IBS) where they have identical alleles that arise independently due to similar mutations. PLINK produces IBD based on 𝜋̂ (pi-hat), which is the proportion of IBD calculated as P(IBD=2) + 0.5*P(IBD=1), and the probabilities of sharing 0, 1, or 2 alleles (Z0, Z1, Z2). The function –genome calculates IBD of all sample pairs based on a subset of autosomal SNPs that are approximately in linkage equilibrium (independent). Another function –min creates a list of individuals with relatedness above the chosen 𝜋̂.

31b. Exclude all individuals above pi-hat 0.2, which indicates they are related, using the function –genome.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_10 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to pihat_min0.2.log.
## Options in effect:
##   --bfile HapMap_3_r3_10
##   --extract indepSNP.prune.in
##   --genome
##   --min 0.2
##   --out pihat_min0.2
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 163 people (79 males, 84 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## --extract: 104144 variants remaining.
## Using up to 8 threads (change this with --threads).
## Before main variant filters, 110 founders and 53 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998051.
## 104144 variants and 163 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.  (53 phenotypes
## are missing.)
## 1152 markers complete.2304 markers complete.3456 markers complete.4608 markers complete.5760 markers complete.6912 markers complete.8064 markers complete.9216 markers complete.10368 markers complete.11520 markers complete.12672 markers complete.13824 markers complete.14976 markers complete.16128 markers complete.17280 markers complete.18432 markers complete.19584 markers complete.20736 markers complete.21888 markers complete.23040 markers complete.24192 markers complete.25344 markers complete.26496 markers complete.27648 markers complete.28800 markers complete.29952 markers complete.31104 markers complete.32256 markers complete.33408 markers complete.34560 markers complete.35712 markers complete.36864 markers complete.38016 markers complete.39168 markers complete.40320 markers complete.41472 markers complete.42624 markers complete.43776 markers complete.44928 markers complete.46080 markers complete.47232 markers complete.48384 markers complete.49536 markers complete.50688 markers complete.51840 markers complete.52992 markers complete.54144 markers complete.55296 markers complete.56448 markers complete.57600 markers complete.58752 markers complete.59904 markers complete.61056 markers complete.62208 markers complete.63360 markers complete.64512 markers complete.65664 markers complete.66816 markers complete.67968 markers complete.69120 markers complete.70272 markers complete.71424 markers complete.72576 markers complete.73728 markers complete.74880 markers complete.76032 markers complete.77184 markers complete.78336 markers complete.79488 markers complete.80640 markers complete.81792 markers complete.82944 markers complete.84096 markers complete.85248 markers complete.86400 markers complete.87552 markers complete.88704 markers complete.89856 markers complete.91008 markers complete.92160 markers complete.93312 markers complete.94464 markers complete.95616 markers complete.96768 markers complete.97920 markers complete.99072 markers complete.100224 markers complete.101376 markers complete.102528 markers complete.103680 markers complete.104144 markers complete.IBD calculations complete.  
## Writing... 99%Finished writing pihat_min0.2.genome .

31c. Determine parent-offspring relations in the dataset based on the probability of sharing 1 allele (Z1) is > 0.9.

awk '{ if ($8 >0.9) print $0 }' pihat_min0.2.genome > zoom_pihat.genome

31d. Plot a histogram to visualize the parent-offspring relationship.

# PO = parent-offspring pairs have Z0 and z2 close to 0, Z1 close to 1.
# UN = unrelated individuals typically have Z0 close to 1 and Z1 and Z2 close to 0.
# Full siblings have Z0 ~ 0.25, Z1 ~ 0.5, and Z2 ~ 0.25
# Subset data for parent-offspring pairs using the subset function with the relationship type (RT):

# Using scale x = 0 to 0.1, y = 0 to 1
relatedness = read.table("pihat_min0.2.genome", header=T)
par(pch=16, cex=1)
with(relatedness,plot(Z0,Z1, xlim=c(0,1), ylim=c(0,1), type="n"))
with(subset(relatedness,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness,RT=="UN") , points(Z0,Z1,col=3))

# Zooming in on x = 0 to 0.02, y = 0.98 to 1
relatedness_zoom = read.table("zoom_pihat.genome", header=T)
par(pch=16, cex=1)
with(relatedness_zoom,plot(Z0,Z1, xlim=c(0,0.02), ylim=c(0.98,1), type="n"))
with(subset(relatedness_zoom,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness_zoom,RT=="UN") , points(Z0,Z1,col=3))

# Histogram
relatedness = read.table("pihat_min0.2.genome", header=T)
hist(relatedness[,10],main="Histogram relatedness", xlab= "Pihat")

- Explain the relationship in the data.

The graph shows the proportion of identity by descent and their recent ancestery. The high Pihat indicates that those points are related. There is some data with a low pihat, but it is much less frequent. #### 31e. Delete the parent-child relationship by keeping founders (individuals without parents in the dataset), i.e., individuals for whom the paternal and maternal individual codes and both 0.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_11.log.
## Options in effect:
##   --bfile HapMap_3_r3_10
##   --filter-founders
##   --make-bed
##   --out HapMap_3_r3_11
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 163 people (79 males, 84 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## 53 people removed due to founder status (--filter-founders).
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate in remaining samples is 0.998025.
## 1073226 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## --make-bed to HapMap_3_r3_11.bed + HapMap_3_r3_11.bim + HapMap_3_r3_11.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

31f. Let’s look again for individuals with pi-hat >0.2.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_11 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2_in_founders

# Use the cat (for concatenate) command to display the contents of pihat_min0.2_in_founders.genome
cat pihat_min0.2_in_founders.genome

# Use the following commands to create a text file with related pairs (PI_HAT>0.2)

awk '{ if ($10 >0.2) print $0 }' pihat_min0.2_in_founders.genome >0.2_low_call_rate_pihat.txt

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to pihat_min0.2_in_founders.log.
## Options in effect:
##   --bfile HapMap_3_r3_11
##   --extract indepSNP.prune.in
##   --genome
##   --min 0.2
##   --out pihat_min0.2_in_founders
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## --extract: 104144 variants remaining.
## Using up to 8 threads (change this with --threads).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997978.
## 104144 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## 1152 markers complete.2304 markers complete.3456 markers complete.4608 markers complete.5760 markers complete.6912 markers complete.8064 markers complete.9216 markers complete.10368 markers complete.11520 markers complete.12672 markers complete.13824 markers complete.14976 markers complete.16128 markers complete.17280 markers complete.18432 markers complete.19584 markers complete.20736 markers complete.21888 markers complete.23040 markers complete.24192 markers complete.25344 markers complete.26496 markers complete.27648 markers complete.28800 markers complete.29952 markers complete.31104 markers complete.32256 markers complete.33408 markers complete.34560 markers complete.35712 markers complete.36864 markers complete.38016 markers complete.39168 markers complete.40320 markers complete.41472 markers complete.42624 markers complete.43776 markers complete.44928 markers complete.46080 markers complete.47232 markers complete.48384 markers complete.49536 markers complete.50688 markers complete.51840 markers complete.52992 markers complete.54144 markers complete.55296 markers complete.56448 markers complete.57600 markers complete.58752 markers complete.59904 markers complete.61056 markers complete.62208 markers complete.63360 markers complete.64512 markers complete.65664 markers complete.66816 markers complete.67968 markers complete.69120 markers complete.70272 markers complete.71424 markers complete.72576 markers complete.73728 markers complete.74880 markers complete.76032 markers complete.77184 markers complete.78336 markers complete.79488 markers complete.80640 markers complete.81792 markers complete.82944 markers complete.84096 markers complete.85248 markers complete.86400 markers complete.87552 markers complete.88704 markers complete.89856 markers complete.91008 markers complete.92160 markers complete.93312 markers complete.94464 markers complete.95616 markers complete.96768 markers complete.97920 markers complete.99072 markers complete.100224 markers complete.101376 markers complete.102528 markers complete.103680 markers complete.104144 markers complete.IBD calculations complete.  
## Writing... 99%Finished writing pihat_min0.2_in_founders.genome .
##    FID1     IID1   FID2     IID2 RT    EZ      Z0      Z1      Z2  PI_HAT PHE       DST     PPC   RATIO
##   13291  NA07045   1454  NA12813 UN    NA  0.2572  0.5007  0.2421  0.4924   0  0.839777  1.0000  9.7022

31g. Open the file pihat_min0.2_in_founders.genome in Notepad or Excel. Look at the first four columns (FID1 = family ID for first sample, IID1 = individual ID for first sample, FID2 = family ID for second sample, IID2 = individual ID for second sample). These are likely to be full sib or DZ twin pairs.

####- How many pairs of individuals with a pi-hat >0.2 remained in the data and what is the pi-hat value? These are likely to be full sibs or DZ twin pairs based on the Z values. 110 people passed the filters. In the table there are 94 rows. The Pi-hat value is around 0.5. The Z values were about z0=0, z1=1 and z2=0. #### 31h. For each pair of related individuals with a pi-hat >0.2, remove the individual with the lowest call rate.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_11 --missing

# Check the total number of individuals missing SNPs using the wc (word count) command
wc -l plink.imiss # prints the number of lines in a file
head -n 5 plink.imiss
wc -l 0.2_low_call_rate_pihat.txt

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to plink.log.
## Options in effect:
##   --bfile HapMap_3_r3_11
##   --missing
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.998025.
## --missing: Sample missing data report written to plink.imiss, and variant-based
## missing data report written to plink.lmiss.
## 111 plink.imiss
##     FID       IID MISS_PHENO   N_MISS   N_GENO   F_MISS
##    1328   NA06989          N     2181  1073226 0.002032
##    1377   NA11891          N    13586  1073226  0.01266
##    1349   NA11843          N      813  1073226 0.0007575
##    1330   NA12341          N     3697  1073226 0.003445
## 2 0.2_low_call_rate_pihat.txt

31i. Delete individuals with the lowest call rate in related pairs with a pi-hat >0.2.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_11 --remove 0.2_low_call_rate_pihat.txt --make-bed --out HapMap_3_r3_12

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_12.log.
## Options in effect:
##   --bfile HapMap_3_r3_11
##   --make-bed
##   --out HapMap_3_r3_12
##   --remove 0.2_low_call_rate_pihat.txt
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1073226 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## --remove: 109 people remaining.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 109 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate in remaining samples is 0.998028.
## 1073226 variants and 109 people pass filters and QC.
## Among remaining phenotypes, 54 are cases and 55 are controls.
## --make-bed to HapMap_3_r3_12.bed + HapMap_3_r3_12.bim + HapMap_3_r3_12.fam ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

32. Population stratification

32b. Create an MDS plot with k=10 dimensions

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_13 --extract indepSNP.prune.in --genome --out HapMap_3_r3_13
"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_13 --read-genome HapMap_3_r3_13.genome --cluster --mds-plot 10 --out HapMap_3_r3_13_mds

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_13.log.
## Options in effect:
##   --bfile HapMap_3_r3_13
##   --extract indepSNP.prune.in
##   --genome
##   --out HapMap_3_r3_13
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1092489 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## --extract: 104144 variants remaining.
## Using up to 8 threads (change this with --threads).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.997978.
## 104144 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## 1152 markers complete.2304 markers complete.3456 markers complete.4608 markers complete.5760 markers complete.6912 markers complete.8064 markers complete.9216 markers complete.10368 markers complete.11520 markers complete.12672 markers complete.13824 markers complete.14976 markers complete.16128 markers complete.17280 markers complete.18432 markers complete.19584 markers complete.20736 markers complete.21888 markers complete.23040 markers complete.24192 markers complete.25344 markers complete.26496 markers complete.27648 markers complete.28800 markers complete.29952 markers complete.31104 markers complete.32256 markers complete.33408 markers complete.34560 markers complete.35712 markers complete.36864 markers complete.38016 markers complete.39168 markers complete.40320 markers complete.41472 markers complete.42624 markers complete.43776 markers complete.44928 markers complete.46080 markers complete.47232 markers complete.48384 markers complete.49536 markers complete.50688 markers complete.51840 markers complete.52992 markers complete.54144 markers complete.55296 markers complete.56448 markers complete.57600 markers complete.58752 markers complete.59904 markers complete.61056 markers complete.62208 markers complete.63360 markers complete.64512 markers complete.65664 markers complete.66816 markers complete.67968 markers complete.69120 markers complete.70272 markers complete.71424 markers complete.72576 markers complete.73728 markers complete.74880 markers complete.76032 markers complete.77184 markers complete.78336 markers complete.79488 markers complete.80640 markers complete.81792 markers complete.82944 markers complete.84096 markers complete.85248 markers complete.86400 markers complete.87552 markers complete.88704 markers complete.89856 markers complete.91008 markers complete.92160 markers complete.93312 markers complete.94464 markers complete.95616 markers complete.96768 markers complete.97920 markers complete.99072 markers complete.100224 markers complete.101376 markers complete.102528 markers complete.103680 markers complete.104144 markers complete.IBD calculations complete.  
## Writing... 20%Writing... 41%Writing... 61%Writing... 83%Writing... 99%Finished writing HapMap_3_r3_13.genome .
## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to HapMap_3_r3_13_mds.log.
## Options in effect:
##   --bfile HapMap_3_r3_13
##   --cluster
##   --mds-plot 10
##   --out HapMap_3_r3_13_mds
##   --read-genome HapMap_3_r3_13.genome
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1092489 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99752.
## 1092489 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## Clustering... [sorting IBS values]Clustering... [100 merges performed]Clustering... done.                        
## Writing cluster solution... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%100%Cluster solution written to HapMap_3_r3_13_mds.cluster1 ,
## HapMap_3_r3_13_mds.cluster2 , and HapMap_3_r3_13_mds.cluster3 .
## Performing multidimensional scaling analysis (SVD algorithm, 10
## dimensions)... done.
## MDS solution written to HapMap_3_r3_13_mds.mds .

# Load the MDS results into R and make a pair-wise MDS plot using dimension 1 (C1) and dimension 2 (C2).

mds_data <- read.table("HapMap_3_r3_13_mds.mds", header=TRUE)

# Create the plot
plot(mds_data$C1, mds_data$C2, xlab="Dimension 1", ylab="Dimension 2", main="MDS Plot")

32c. Change the format of the .mds file into a plink covariate file.

awk '{print$1, $2, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13}' HapMap_3_r3_13_mds.mds > covar_mds.txt

Part E. Association analysis and visualization

34. Test association using the –assoc option.

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_13 --assoc --out assoc_results

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to assoc_results.log.
## Options in effect:
##   --assoc
##   --bfile HapMap_3_r3_13
##   --out assoc_results
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1092489 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99752.
## 1092489 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## Writing C/C --assoc report to assoc_results.assoc ... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.

36. Multiple testing correction

"C:/Program Files/PLINK/plink.exe" --bfile HapMap_3_r3_13 --assoc --adjust --out adjusted_assoc_results

## PLINK v1.90b6.26 64-bit (2 Apr 2022)           www.cog-genomics.org/plink/1.9/
## (C) 2005-2022 Shaun Purcell, Christopher Chang   GNU General Public License v3
## Logging to adjusted_assoc_results.log.
## Options in effect:
##   --adjust
##   --assoc
##   --bfile HapMap_3_r3_13
##   --out adjusted_assoc_results
## 
## 16106 MB RAM detected; reserving 8053 MB for main workspace.
## 1092489 variants loaded from .bim file.
## 110 people (55 males, 55 females) loaded from .fam.
## 110 phenotype values loaded from .fam.
## Using 1 thread (no multithreaded calculations invoked).
## Before main variant filters, 110 founders and 0 nonfounders present.
## Calculating allele frequencies... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
## Total genotyping rate is 0.99752.
## 1092489 variants and 110 people pass filters and QC.
## Among remaining phenotypes, 55 are cases and 55 are controls.
## Writing C/C --assoc report to adjusted_assoc_results.assoc ...
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%done.
## --adjust: Genomic inflation est. lambda (based on median chisq) = 1.0157.
## 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99%--adjust values (1092489 variants) written to
## adjusted_assoc_results.assoc.adjusted .

37. Visualize SNP-phenotype associations with Manhattan and Q-Q plots.

37a. Make a Manhattan plot with highlighted SNPs.

# Load the GWAS results into R

results_as <- read.table("assoc_results.assoc", head=TRUE)

# Define the significance threshold

significance_threshold <- 5e-5

# Create the Manhattan plot
manhattan(results_as,
          chr="CHR",
          bp="BP",
          p="P",
          snp="SNP", 
          col = c("tomato3", "blue"),
          suggestiveline = -log10(1e-5), 
          genomewideline = -log10(significance_threshold), 
          highlight = results_as$SNP[results_as$P < significance_threshold], 
          main = "Manhattan plot with significant SNPs")

###Annotate the name of SNPs with a high p-value using the annotatePval argument.

manhattan(results_as, annotatePval = 0.01)

- Filter significant SNPs (p-value <5×10-8) using the bash function awk and save.

awk '$9 < 5e-5' adjusted_assoc_results.assoc > significant_snps.txt

####- According to the Manhattan plot, how many SNPs are associated with the phenotype and on which chromosomes are they located? 22 SNP’s were left after filtering the significant SNPs and they are on chromosomes. ####- Does a significant association between genetic variants and a phenotype definitively predict that the SNP is responsible for expression of the phenotype? Yes because the results are below the genome wide significance value and thus it follows the suggested significance. ####- The majority of genetic changes in GWAS lie outside of the coding region. Explain how changes in non-coding regions could affect expression of the phenotype. Changes outside the coding regions can mess with regulatory functions and thus alter the gene, in turn impacting the genotype. The alteration of gene expression can mess with the control. #### 37b. Make a Q-Q plot.

qq(results_as$P, main = "Q-Q plot of GWAS p-values: log")

- Interpret the results of the Q-Q plot.

There isn’t population stratification seen. The plot hugs the linear line, thus fitting the expected p values. There is no early separation seen, rather separation only occurs further along the x-axis, where the observed P values are slightly greated than the expected. #### 38. Compile the report (knit). Click on Knit above the Editor window. The knit function takes an input file, extracts the R codes, evaluates the codes, writes the compiled document into a file in a format you specified in the YAML.

- You can achieve the same thing by rendering the R Markdown document (.Rmd) to the specified output format using the function package.render(“file_name”).

- Note: if for some reason you get error, delete the code that produces the error and run Knit.

40. When you are done, click on the wheel next to Knit –> Clear all output.

41. Close the .Rmd document in the Editor window. Clear objects in the Console and Environment (use the brush icon). Then, type quit() at the Console to close R.

42. Upload the report to the Labs folder in Blackboard for grading.