<jupyter_start><jupyter_text>Excercise GWASThe goal of this excercise is to reproduce a published GWAS study.We will be using the ["<NAME> (2013). Genome-wide association study using cellular traits identifies a new regulator of root development in Arabidopsis. Nature Genetics"](https://www.nature.com/articles/ng.2824)You can access the study and interactive Manhattan plots on [GWA-Portal](https://gwas.gmi.oeaw.ac.at//study/739/phenotypes). We also added the static Manhattan plots below generated by the [PyGWAS library](https://github.com/timeu/PyGWAS) (another GWAS library)|   |   ||---|---||![phenotype.png](images/meijon/phenotype.png)  | ![phenotype.png](images/meijon/paper.png)  |First we initialize the inline plotting and disable some warnings<jupyter_code>%matplotlib inline
from warnings import simplefilter 
simplefilter(action='ignore', category=FutureWarning)
import pylab as pl<jupyter_output><empty_output><jupyter_text>Load the librariesYou will need to load the necessary libraries: - h5py - pandas - numpy - scipy<jupyter_code># load the libraries above<jupyter_output><empty_output><jupyter_text>Parse and convert phenotype dataFirst we need to parse the 2 phenotypes that are stored in *data/meijon/MC_length.csv* and *data/meijon/MZ_length.csv* into one panda dataframe and store the panda dataframe in *data/meijon/gwas.h5* using the panda hdf5 format You should merge the two datafarmes before saving them. Make sure to properly set the index and check the panda documentation how to merge 2 dataframes<jupyter_code># load MC_length.csv

# load MZ_length.csv

# merge both dataframes

# store merged panda dataframe into hdf5 file<jupyter_output><empty_output><jupyter_text>Load and convert the genotype data Download the data from [here](https://drive.google.com/drive/folders/0BwYlrRGOfVARS01aMnFVVWotbXc), extract it and load the genotype data that is stored in an hdf5 file and create a Panda dataframe out if it. Then store the panda dataframe in *data/meijon/gwas.h5* using the panda hdf5 format The genotype is stored slightly differently inside this HDF5 compared to our previous notebooks. The chromosomes for example are not stored in an array but the chromosome regions (indices) are stored as an array in the attributes of the "positions" dataset. You need to generate the chromosome array using those indices.Also don't subset the SNPs. Use all of them for the analysis below Store the chromosomes along with the positions as a Panda Multiindex. This will make it possible to filter by both indices (chromosome and position)<jupyter_code># load the genotype from the hdf5 file and create a dataframe that contains the SNP matrix. 
# The columns should be the sample_IDs and the index should be a MultiIndex of chromosomes and positions

# YOUR CODE

# Quiz:
# how much memory does the dataframe consume ? 
# are there ways to make it more efficient ?<jupyter_output><class 'pandas.core.frame.DataFrame'>
MultiIndex: 206088 entries, (1.0, 657) to (5.0, 26974904)
Columns: 1386 entries, 9381 to 9127
dtypes: float64(1386)
memory usage: 2.1 GB<jupyter_text>When using the panda function to save to an hdf5 file, you can specify a 'key' under which the panda dataframe is stored inside the hdf5 file, allowing multiple different dataframes to be stored in the same hdf5 file<jupyter_code># store the genotype dataframe in the gwas.h5 dataset side by side with the phenotype

# use bash command to check the size of the gwas.h5 file on the filesystem


# Quiz:
# does the file system on disk match the size in memory ?
# use HDF5 command line tools to check the hierarchical structure of the gwas.h5 file


# Quiz:
# what do you notice in the structure ?<jupyter_output>/                        Group
/genotype                Group
/genotype/axis0          Dataset {1386}
/genotype/axis1_label0   Dataset {206088}
/genotype/axis1_label1   Dataset {206088}
/genotype/axis1_level0   Dataset {5}
/genotype/axis1_level1   Dataset {205494}
/genotype/block0_items   Dataset {1386}
/genotype/block0_values  Dataset {206088, 1386}
/phenotype               Group
/phenotype/axis0         Dataset {2}
/phenotype/axis1         Dataset {201}
/phenotype/block0_items  Dataset {2}
/phenotype/block0_values Dataset {201, 2}<jupyter_text>Cleanup phenotypes and align genotype dataframe with phenotype dataframeRemove N/A from the phenotype dataframe and make sure that the genotype and phenotype dataframes only contain accessions that occur in both<jupyter_code>#Remove NaNs from the phenotype dataset

#Filter the SNP matrix to only include the accessions for which we have genotypes


# Quiz
# have all accessions that have been phenotyped also genotyped ?<jupyter_output><empty_output><jupyter_text>Check model assumption. Next we will check the distrubtion of the two phenotypes and perform the Shapiro-Wilk test for normality.  To calculate the Shapiro-Wilk statistics, you can use the scipy.stats package.<jupyter_code># Use limix to plot the phenotypic distrubition of the phenotype and add the Shapiro-Wilk pvalue to the title of the plot
# Don't forget to import the plot function from the limix library.

# YOUR CODE

# Quiz:
# which phenotype is more normal ?<jupyter_output>/Users/uemit.seren/miniconda3/envs/gwas-lecture/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6571: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "<jupyter_text>Transformation of the phenotypePerfom box-cox and log normal transformation, plot the histogram and calculcate the shapiro wilk score For the log-normal transformation you can use following formular: $$Y_{i}=\log( Y_{i} - \min(Y) + 0.1 * var(Y)) $$<jupyter_code># Plot the histogram and calculate the Shapiro-Wilk pvalue for box-cox and log transformed phenotypes

# YOUR CODE

# Quiz:
# Does the boxcox transformation help ?<jupyter_output><empty_output><jupyter_text>Standard linear regressionRun the standard linear regression for all 3 phenotypes (raw, box-cox and log transformed) and genetate the Manhattan plots<jupyter_code># Run the linear regression and create a dataframe of the pvalues


# Plot the manhattan plots for both phenotypes<jupyter_output><empty_output><jupyter_text>Compare your plots to these: MC length![MC_length](images/meijon/MC_lm.png) MZ length![MC_length](images/meijon/MZ_lm.png)  QQ-Plotsgenerate qq-plots and check the inflation<jupyter_code>## plot the QQ plots for both phenotypes

# YOUR CODE
    
# Quiz:
# can we see any inflation ?<jupyter_output><empty_output><jupyter_text>Compare your plots to these: MC length![MC_length](images/meijon/MC_qq_lm.png) MZ length![MC_length](images/meijon/MZ_qq_lm.png)  Linear Mixel ModelRun the same analysis with the Linear Mixed Model Make sure to filter out monophormic SNPs from the genotype dataset before you calculate the kinship matrix and run the linear mixed model<jupyter_code># Filter out monomorphic SNP
# one way to do it is to sum up the SNPs (they are 0,1) and remove them if they are equal 0 or the number of accessions
# see the Numpy Lecture for how to use masks to filter arrays

# YOUR CODE

# Quiz:
# how many SNP have been filtered out ?
# Estimate sample relatedness (kinship) using the monomorphic SNPs from above

#normalize it

# Plot the kinship

# YOUR CODE

# Quiz:
# Do you think it is possible to calculate the kinship for all accession once and use this even for a subset or 
# do you always have to re-calculate the kinship for the accessions that you are running the GWAS for ?
# Run linear mixed model

# Retrieve pvalues and create a dataframe for both phenotypes
# retrieve the P-values, the Betas, and the standard errors.

# Plot the manhattan plots for both phenotypes

# YOUR CODE

# Quiz:
# What makes these two 2 phenotypes/GWAS results special ?
# What furhter analysis would you do ?<jupyter_output><empty_output><jupyter_text>Compare your plots to these: MC length![MC_length](images/meijon/MC_lmm.png) MZ length![MC_length](images/meijon/MZ_lmm.png)<jupyter_code>## plot the QQ plots for both phenotypes

# YOUR CODE
    
# Quiz:
# does the QQ plot improve compared to the LM ?<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06607554961739341, 'fraction_numerical': 0.017976436292967326, 'mean_word_length': 4.819081272084806, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 48, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '634974', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 2642, 'n_tokens_neox': 2414, 'n_words': 1025, 'path': 'limix1/Exercise-1-GWAS.ipynb', 'repo_name': 'mahort/gwas-lecture', 'size': '8240'}