<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Basic data analysis with Pandas
Today we will learn how to
Step1:  Loading Data
Pandas can load data from multiple sources
Step2:  Multiple types of data interspersed together.
Let's separate them.
Accessing the DataFrame
Two main methods to access elements
Step3:  For example, 'RPS11' or 'ELMO2' are the indices, or 'names' of two rows while 'Cy72_CD45_H02_S758_comb' is the name of the first column. 
However, they are also the fifth and sixth rows or the first column, respectively. 
Let's select the 'RPS11' row
Step4:  Now let's select the fifth row.
Keep in mind that Python is zero-indexed (i.e. the index of the first element is 0) thus, we need to pick the '4th' element.
Step5:  How can we compare these? 
Let's retrieve the values as an array and see if they are all equal to each other.
Step6:  We get a lot of values boolean values...
Step7:  4645 values in fact---one for each column (i.e. cell).
The first and last three look identical, but there are thousands of others. If just one of them is different, however, the equality check should fail. How can we check them all?
Step8:  All the entries are identical.
What if the numbers were just so slightly off, but not enough that we'd mind.
Example
Step9:  This is our way of checking if all values are close to each other.
'Close' is defined as within some fixed relative or absolute tolerance.
Pandas Series
What if we wanted to sort the cells by their expression of a particular gene?
First we need to get the expression of that gene in every cell.
Let's do this for an example gene, CD4.
Step10:  Let's find the cells with the highest expression
Step11:  Let's find the index of the 10 cells with the highest CD4 expression
Step12:  How about the 10 highest and 10 lowest?
First, a distinction.
Step13:  These are disjoint, so let's unite them.
Step14:  Let's find the average expression over all the cells
Step15:  Alternatively
Step16:  or
Step17:  Let's see if they are all close
Step18:  Something interesting happened there... 
Broadcasting
Broadcasting works with Pandas Series or DataFrame objects as well.
Let's see some examples
Step19:  DataFrame level summary operations
Thus far, we have mostly summarized information for the Series object.
This is easier because Series is like a simple list.
What if we want to operate on the whole expression dataset?
Example goal
Step20:  Let's strip the information in the first three rows
Step21:  We now have celltype stored in two different ways
Step22:  Retrieving values
To rank normalize, first we need all the non-zero expression values.
Step23:  Note
Step24:  Now let's normalize each expression value
Step25:  Can we do this in one-line?
Step26:  Is the output equivalent to the other method?
Step27:  Yes, these are equivalent.
What would the distribution be?
Step28:  It's distributed identically, except now it is normalized to the 0-1 interval.
What if we want to analyze groups of data?
Example
Step29:  Let's get a cell count using the value_counts() method of a Series
Step30:  Now let's group...
Grouping
Now, let's use groupby
This useful functionality allows us to group elements. 
First we'll transpose the data, and then apply the grouping.
Step31:  What do we do with this?
Handling groups
The most straightforward way is to get a group and work with it.
Step32:  That's the right size.
Step33:  A straightforward extension of this to calculating per-group average is
Step34:  A more elegant way to achieve this
Step35:  What if we wanted to see a few different marker genes?
Here's a rough (and inaccurate) guideline
Step36:  What if we wanted to do something more complicated?
For example, all these marker genes have zero expression at minimum in their respective cell populations. 
How many cells actually have zero expression of their markers?
Step37:  Output
Finally, let's see how we can write Pandas data frames to files.
For example, let's write the CD3D expression of T cells into an Excel file called T-cells.xlsx in the appropriate format.
Any ideas?
<ASSISTANT_TASK:>
Python Code:
import sys,os
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
import pandas as pd
import numpy as np
import scipy
sns.set_style('white')
if not os.path.isdir('Figures'):
    os.mkdir("Figures")
Explanation: Basic data analysis with Pandas
Today we will learn how to:
- load a single cell RNA-seq expression dataset
- access the data
- group the data
- compute summary statistics
We will use a real scRNA-seq dataset from Tirosh et al, Science, 2016. 
Please first download the dataset from this link.
End of explanation
file_loc = "/Users/murat/Data/tirosh/GSE72056_melanoma_single_cell_revised_v2.txt.gz"
sc_data = pd.read_csv(file_loc, "\t", index_col=0)
sc_data
Explanation: Loading Data
Pandas can load data from multiple sources:
- Delimited files ("TSV", "CSV", etc)
- Excel files
- JSON
- HTML
- and everything else listed here...
Delimited file format is highly common, and the example dataset we will use has this format. The others are similar.
End of explanation
sc_data
Explanation: Multiple types of data interspersed together.
Let's separate them.
Accessing the DataFrame
Two main methods to access elements:
- indices
- numbers
End of explanation
rps11 = sc_data.loc['RPS11']
rps11
Explanation: For example, 'RPS11' or 'ELMO2' are the indices, or 'names' of two rows while 'Cy72_CD45_H02_S758_comb' is the name of the first column. 
However, they are also the fifth and sixth rows or the first column, respectively. 
Let's select the 'RPS11' row:
End of explanation
fifth = sc_data.iloc[4,:]
fifth
Explanation: Now let's select the fifth row.
Keep in mind that Python is zero-indexed (i.e. the index of the first element is 0) thus, we need to pick the '4th' element.
End of explanation
rps11.values
fifth.values
fifth.values == rps11.values
Explanation: How can we compare these? 
Let's retrieve the values as an array and see if they are all equal to each other.
End of explanation
len(fifth.values)
Explanation: We get a lot of values boolean values...
End of explanation
np.all(fifth.values == rps11.values)
Explanation: 4645 values in fact---one for each column (i.e. cell).
The first and last three look identical, but there are thousands of others. If just one of them is different, however, the equality check should fail. How can we check them all?
End of explanation
fifth_full = fifth.values
fifth_full
fifth_full.dtype
fifth_small = np.float16(fifth.values)
fifth_small
fifth_small.dtype
np.all(fifth_small == fifth_full)
np.any(fifth_small == fifth_full)
np.allclose(fifth_full, fifth_small, rtol=1e-3)
Explanation: All the entries are identical.
What if the numbers were just so slightly off, but not enough that we'd mind.
Example: same data but with a different (less memory consuming) data structure
End of explanation
cd4 = sc_data.loc['CD4']
cd4
type(sc_data)
type(cd4)
Explanation: This is our way of checking if all values are close to each other.
'Close' is defined as within some fixed relative or absolute tolerance.
Pandas Series
What if we wanted to sort the cells by their expression of a particular gene?
First we need to get the expression of that gene in every cell.
Let's do this for an example gene, CD4.
End of explanation
cd4.sort_values(ascending=False)
Explanation: Let's find the cells with the highest expression
End of explanation
cd4.sort_values(ascending=False).index[:10]
Explanation: Let's find the index of the 10 cells with the highest CD4 expression
End of explanation
type(cd4.sort_values(ascending=False).index[:10])
type(cd4.sort_values(ascending=False).index[:10].values)
[cd4.sort_values(ascending=False).index[:10].values, 
                cd4.sort_values(ascending=False).index[-10:].values]
Explanation: How about the 10 highest and 10 lowest?
First, a distinction.
End of explanation
np.concatenate([cd4.sort_values(ascending=False).index[:10].values, 
                cd4.sort_values(ascending=False).index[-10:].values])
Explanation: These are disjoint, so let's unite them.
End of explanation
cd4.mean()
Explanation: Let's find the average expression over all the cells
End of explanation
np.mean(cd4)
Explanation: Alternatively
End of explanation
sum(cd4.values)/len(cd4)
Explanation: or
End of explanation
np.allclose(sum(cd4.values)/len(cd4),[cd4.mean(),np.mean(cd4.values)])
Explanation: Let's see if they are all close:
End of explanation
cd4+10
list(zip(cd4, cd4+10, cd4*10, np.power(cd4,2)))
Explanation: Something interesting happened there... 
Broadcasting
Broadcasting works with Pandas Series or DataFrame objects as well.
Let's see some examples:
End of explanation
sc_data
Explanation: DataFrame level summary operations
Thus far, we have mostly summarized information for the Series object.
This is easier because Series is like a simple list.
What if we want to operate on the whole expression dataset?
Example goal: normalize expression
Let's replace expression values with values in the range 0 to 1, such that 0 and 1 respectively correspond to the least and highest expression values in the whole dataset.
We can't use the data frame as it stands because there is non-expression data. We need to strip them away.
End of explanation
patient = sc_data.iloc[0,:]
malignant = sc_data.iloc[1,:]
celltype = sc_data.iloc[2,:]
exp = sc_data.iloc[3:,:]
Explanation: Let's strip the information in the first three rows
End of explanation
celltype[malignant==2] = 7
Explanation: We now have celltype stored in two different ways: 
1. malignant
2. cell type
Let's unite them under celltype
End of explanation
vals = exp.values.flatten()
plt.figure(figsize=(9,6))
_=sns.distplot(vals[vals>0])
plt.savefig("Figures"+os.sep+"tirosh_exp_hist.png")
Explanation: Retrieving values
To rank normalize, first we need all the non-zero expression values.
End of explanation
min_val = np.min(vals)
max_val = np.max(vals)
print('Min: {0:.2f}, Max: {1:.2f}'.format(min_val,max_val))
Explanation: Note: We'll cover plots later. For now, take this code for granted. 
Now we need to normalize each expression value to the 0-1 interval.
We need the minimum and maximum expression values.
End of explanation
norm_exp = (exp-min_val) / max_val
norm_exp
print('Normalized min: {0:.1f}, max: {1:.1f}'.format(np.min(norm_exp.values), np.amax(norm_exp.values)))
Explanation: Now let's normalize each expression value:
End of explanation
norm_exp2 = (exp - exp.values.flatten().min()) / exp.values.flatten().max()
Explanation: Can we do this in one-line?
End of explanation
np.all(norm_exp == norm_exp2)
Explanation: Is the output equivalent to the other method?
End of explanation
norm_vals = norm_exp.values.flatten()
plt.figure(figsize=(9,6))
_=sns.distplot(norm_vals[norm_vals>0])
plt.savefig("Figures"+os.sep+"tirosh_normexp_hist.png")
Explanation: Yes, these are equivalent.
What would the distribution be?
End of explanation
celltype = sc_data.iloc[2,:]
celltype[malignant==2] = 7
celltype_dict = {1:"T",2:"B",3:"Macro",4:"Endo",5:"CAF",6:"NK",7:"Mal",0:"unspec"}
celltype.rename('celltype',inplace=True)
celltype = celltype.apply(lambda x:celltype_dict[x])
Explanation: It's distributed identically, except now it is normalized to the 0-1 interval.
What if we want to analyze groups of data?
Example: calculate average expression of T-cell marker CD3D in each cell type
Our expression dataset does not have the cell type information any more. 
First, let's replace cell type indices with human readable names.
End of explanation
celltype.value_counts()
Explanation: Let's get a cell count using the value_counts() method of a Series
End of explanation
exp = pd.DataFrame(sc_data.iloc[3:,:], dtype=np.float)
print(exp.shape)
expT = exp.transpose()
print(expT.shape)
exp_g = expT.groupby(celltype.values)
print(type(exp_g))
exp_g.groups
Explanation: Now let's group...
Grouping
Now, let's use groupby
This useful functionality allows us to group elements. 
First we'll transpose the data, and then apply the grouping.
End of explanation
t_cells = exp_g.get_group('T')
print(t_cells.shape)
Explanation: What do we do with this?
Handling groups
The most straightforward way is to get a group and work with it.
End of explanation
print('Average CD3D expression in T cells: {0:.2f}'.format(t_cells['CD3D'].mean()))
Explanation: That's the right size.
End of explanation
for g_ind in exp_g.groups.keys():
    group = exp_g.get_group(g_ind)
    print('Average CD3D expression in {2:s} cells is {0:.2f} (n={1:d} cells)'.format(
        group['CD3D'].mean(), group.shape[0], g_ind))
Explanation: A straightforward extension of this to calculating per-group average is:
End of explanation
exp_g.aggregate({'CD3D': ['count', 'mean', 'min', 'max']})
Explanation: A more elegant way to achieve this:
End of explanation
exp_g.aggregate({
    'CD3D': ['count', 'mean', 'min', 'max'],
    'CD14': ['mean', 'min', 'max'],
    'MS4A1': ['mean', 'min', 'max'],
    'NCAM1': ['mean', 'min', 'max'],})
Explanation: What if we wanted to see a few different marker genes?
Here's a rough (and inaccurate) guideline:
* CD14: macrophages
* CD20: B cells (aka MS4A1)
* CD56: NK
End of explanation
f = lambda x: sum(x==0)
exp_g.aggregate({
    'CD3D': ['count',f],
    'CD14': f,
    'MS4A1': f,
    'NCAM1': f,})
Explanation: What if we wanted to do something more complicated?
For example, all these marker genes have zero expression at minimum in their respective cell populations. 
How many cells actually have zero expression of their markers?
End of explanation
t_cells['CD3D'].to_excel('T-cells.xlsx')
!open T-cells.xlsx
Explanation: Output
Finally, let's see how we can write Pandas data frames to files.
For example, let's write the CD3D expression of T cells into an Excel file called T-cells.xlsx in the appropriate format.
Any ideas?
End of explanation
<END_TASK>