<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Create a flotilla study from the Illumina Bodymap 2.0 data
By
Step1:  What does this look like? Let's look at the top of the file with head.
Step2:  We'll use the pandas data analysis library to read the data. But we'll need a few caveats. 
The first three lines of the data don't match the tabular format, so we'll need to ignore them with skiprows=3.
The row names indicating the gene name are the first two columns, so we'll need to specify that with index_col=[0, 1], which says that the index columns/row names are the 0th and 1st columns, since we're counting from zero in the computer world.
Let's import pandas and all the other packages we'll need.
Step3:  Interesting side note
Step4:  Now, let's also replace all our NAs with 0s. This makes sense because if a gene is not detected, then its expression value will be 0. From the header file of the original data, it said
Step5:  Cool, pretty close to 0.5. Now let's replace all NAs with 0, with fillna(0).
Step6:  Finally, we will add 1 and log-transform the data, so it's closer to normally distributed. Gene expression data is known to be log-normal.
Step7:  Create metadata about the expression features
The other thing that will make this data simpler to work with is making the columns of the data just use the crazy unique ENSEMBL id like "ENSG00000280433". First, let's check to see if the common names are actually unique. We'll expand out the column names to a list of tuples using expression.columns.tolist()
Step8:  So there's fewer unique gene names, meaning we should use the ENSEMBL ids for the unique IDs. We'll do this by resetting the columns of expression, and creating metadata about the expression features, stored as expression_feature_data.
First, let's reassign the columns as the ensembl_ids we created from before.
Step9:  Now let's create the expression_feature_data DataFrame, and add a column of 'gene_name' for the renamed feature.
Step10:  Preparing the metadata
For flotilla, every project is required to have metadata. We'll create one from scratch using the sample names from the expression data.
Step11:  The first category we'll use is pretty straightforward, it'll just be the name of the tissue. We can just use the index as the 'phenotype'.
Step12:  Next, let's add some categories on the data, grouping different tissue types together that have the same structure or function. My awesome MD/PhD friend Cynthia Hsu (you have to search for her name on the webpage) came up with these categories.
Step13:  Assign colors to the phenotypes
Now we need to choose colors for the data. Since we have 16 samples, none of the usual ColorBrewer colors apply because the maximum color set there is Paired, which is 12 colors. So we'll use a human-friendly version of the hue-lightness-saturation (hls) scale, called "husl". Read more at the seaborn color tutorial.
Step14:  Let's create an iterator so we can easily loop over this list of colors without having to reference indices.
Step15:  Finally, let's make the phenotype to color mapping as a dictionary.
Step16:  Make a flotilla study!
<ASSISTANT_TASK:>
Python Code:
! curl http://www.ebi.ac.uk/gxa/experiments/E-MTAB-513.tsv > E-MTAB-513.tsv
! curl http://www.ebi.ac.uk/arrayexpress/files/E-MTAB-513/E-MTAB-513.sdrf.txt>  E-MTAB-513.sdrf.txt
Explanation: Create a flotilla study from the Illumina Bodymap 2.0 data
By: Olga Botvinnik
In this example, we will create a flotilla study from the Illumina Bodymap 2.0 data. This is a nice dataset because it is simple yet complex. It has several nice features:
Already processed to publically available gene expression values
16 human tissue types, so its relatable to people across biological disciplines, and to non-scientists (everyone knows what a brain is)
Can group the tissue types as meta-types
These data are available from the European Bioinformatics Institute (EBI) website here (you can browse their other datasets on their Data library website). There's a link on the "Expression Atlas" view of the data to download the processed data, which is the link I use below. So let's get this data!
Download the data and supplement
End of explanation
! head E-MTAB-513.tsv
! head E-MTAB-513.sdrf.txt
Explanation: What does this look like? Let's look at the top of the file with head.
End of explanation
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import flotilla
expression = pd.read_table('E-MTAB-513.tsv', skiprows=3, index_col=[0, 1])
expression.head()
Explanation: We'll use the pandas data analysis library to read the data. But we'll need a few caveats. 
The first three lines of the data don't match the tabular format, so we'll need to ignore them with skiprows=3.
The row names indicating the gene name are the first two columns, so we'll need to specify that with index_col=[0, 1], which says that the index columns/row names are the 0th and 1st columns, since we're counting from zero in the computer world.
Let's import pandas and all the other packages we'll need.
End of explanation
expression = expression.T
expression.head()
Explanation: Interesting side note: you'll notice that "animal ovary" is the only thing indicated as "animal." This is because in the Experimental Factor Ontology (EFO), they are defined as "animal ovary." This is because the EFO defines plant ovary separately. Apparently, in the venn diagram of the organs in animals versus the organs in plants, ovary is one of the few things that overlap! Cool! Thanks to Nick Semenkovich for figuring this out.
For flotilla, we follow the machine learning standard of using matrices in the format $(\text{samples}) \times (\text{features})$, like this
Since this data is $(\text{features}) \times (\text{samples})$, we simply need to transpose the matrix with .T
End of explanation
expression.min().min()
Explanation: Now, let's also replace all our NAs with 0s. This makes sense because if a gene is not detected, then its expression value will be 0. From the header file of the original data, it said:
Genes matching: 'protein_coding' exactly, specifically expressed in any Organism part above the expression level cutoff: 0.5 in experiment E-MTAB-513
So we're forcing everything with expression less than 0.5 down to 0. Out of curiousity, let's see what was the minimum value left over after they did this 0.5 filter:
End of explanation
expression = expression.fillna(0)
Explanation: Cool, pretty close to 0.5. Now let's replace all NAs with 0, with fillna(0).
End of explanation
expression = np.log2(expression + 1)
Explanation: Finally, we will add 1 and log-transform the data, so it's closer to normally distributed. Gene expression data is known to be log-normal.
End of explanation
ensembl_ids = pd.Index([a for a, b in expression.columns.tolist()])
gene_names = pd.Index([b for a, b in expression.columns.tolist()])
len(ensembl_ids.unique())
len(gene_names.unique())
Explanation: Create metadata about the expression features
The other thing that will make this data simpler to work with is making the columns of the data just use the crazy unique ENSEMBL id like "ENSG00000280433". First, let's check to see if the common names are actually unique. We'll expand out the column names to a list of tuples using expression.columns.tolist()
End of explanation
expression.columns = ensembl_ids
Explanation: So there's fewer unique gene names, meaning we should use the ENSEMBL ids for the unique IDs. We'll do this by resetting the columns of expression, and creating metadata about the expression features, stored as expression_feature_data.
First, let's reassign the columns as the ensembl_ids we created from before.
End of explanation
expression_feature_data = pd.DataFrame(index=ensembl_ids)
expression_feature_data['gene_name'] = gene_names
expression_feature_data.head()
Explanation: Now let's create the expression_feature_data DataFrame, and add a column of 'gene_name' for the renamed feature.
End of explanation
metadata = pd.DataFrame(index=expression.index)
metadata.head()
Explanation: Preparing the metadata
For flotilla, every project is required to have metadata. We'll create one from scratch using the sample names from the expression data.
End of explanation
metadata['phenotype'] = metadata.index
metadata.head()
Explanation: The first category we'll use is pretty straightforward, it'll just be the name of the tissue. We can just use the index as the 'phenotype'.
End of explanation
# All of these tissue types are part of the reproductive system
metadata['reproductive'] = metadata.phenotype.isin(['animal ovary', 'testis'])
# ALl of these tissue types generate hormones
metadata['hormonal'] = metadata.phenotype.isin(['animal ovary', 'testis', 'adrenal gland', 'thyroid'])
# These tissues are part of in the immune system
metadata['immune'] = metadata.phenotype.isin(['leukocyte', 'thyroid', 'lymph node'])
# These tissues are fatty
metadata['fatty'] = metadata.phenotype.isin(['adipose tissue', 'brain', 'breast'])
# These tissues contain either smooth (involuntary) or skeletal (voluntary) muscle
metadata['muscle'] = metadata.phenotype.isin(['colon', 'heart', 'prostate', 'skeletal muscle'])
# These tissues' main function is to filter blood in some way
metadata['filtration'] = metadata.phenotype.isin(['colon', 'kidney', 'liver'])
# These tissues have high blood flow to them, compared to other tissues
metadata['high_blood_flow'] = metadata.phenotype.isin(['brain', 'colon', 'kidney', 'liver', 'lung'])
Explanation: Next, let's add some categories on the data, grouping different tissue types together that have the same structure or function. My awesome MD/PhD friend Cynthia Hsu (you have to search for her name on the webpage) came up with these categories.
End of explanation
colors = sns.color_palette('husl', len(expression.index))
sns.palplot(colors)
Explanation: Assign colors to the phenotypes
Now we need to choose colors for the data. Since we have 16 samples, none of the usual ColorBrewer colors apply because the maximum color set there is Paired, which is 12 colors. So we'll use a human-friendly version of the hue-lightness-saturation (hls) scale, called "husl". Read more at the seaborn color tutorial.
End of explanation
colors_iter = iter(colors)
Explanation: Let's create an iterator so we can easily loop over this list of colors without having to reference indices.
End of explanation
phenotype_to_color = {phenotype: colors_iter.next() for phenotype in metadata.phenotype}
phenotype_to_color
Explanation: Finally, let's make the phenotype to color mapping as a dictionary.
End of explanation
study = flotilla.Study(metadata, expression_data=expression, 
                       metadata_phenotype_to_color=phenotype_to_color, 
                       expression_feature_data=expression_feature_data, 
                       expression_feature_rename_col='gene_name',
                       species='hg19')
study.expression.feature_data.head()
study.interactive_pca()
study.plot_gene('RBFOX1')
study.plot_gene('ADIPOQ')
study.plot_gene('OCA2')
study.plot_gene('MAPT')
era_metadata = pd.read_table("E-MTAB-513.sdrf.txt")
from pprint import pprint
pprint(sorted(era_metadata.columns.tolist()))
era_metadata['Comment[biosource provider]']
era_metadata['Scan Name']
era_metadata['Material Type']
era_metadata['Material Type.1']
era_metadata.Description
study.interactive_clustermap()
study.save('bodymap2')
Explanation: Make a flotilla study!
End of explanation
<END_TASK>