<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 CS446/546 - Class Session 20 - Partial correlation network
In this class session we will continue with analyzing the tumor gene expression dataset from the NIH human bladder cancer cohort (M=414 tumors), building on what we learned in Class Session 19 (Correlation Network).  In order to keep the analysis simple&ast;, in this notebook will restrict our analysis to a set of N=164 genes that are very highly expressed in bladder cancer. Using the 164 x 414 matrix of transcript abundance measurements, we will construct a network based on gene-gene partial correlation coefficients. We will also compare the distribution of partial correlation coefficients to the distribution of Pearson correlation coefficients. Do you think they will be different?  In what way would you expect them to be different?
&ast; Here, "simple" means that the covariance matrix will be nonsingular, so that we can obtain the partial correlation matrix by inversion.
We'll import all of the python modules that we will need for this exercise
Step1:  Read the tab-deliminted text file of gene expression measurements (rows correspond to genes, columns correspond to bladder tumor samples). (use Pandas, pandas.read.csv, and as_matrix). As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using shape.
Step2:  Convert your data frame to a numpy matrix, using the pandas.DataFrame.as_matrix method.
Step3:  As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using shape.
Step4:  Compute the median expression level for each row of your matrix
Step5:  Filter the matrix to include only rows for which the gene's median expression is > 12; matrix should now be 164 x 414; this will enable us to easily compute the partial correlation matrix using the inverse of the covariance matrix. Print the size of the filtered matrix, as a sanity check.
Step6:  Print the shape of your filtered matrix, as a sanity check.  It should be 164x414.
Step7:  Compute the 164 x 164 matrix of gene-gene Pearson correlation coefficients, using numpy.corrcoef (this function treats each row as a random variable, so you don't have to do any transposing of the matrix, unlike the situation in R).
Step8:  Compute the covariance matrix using sklearn.covariance.empirical_covariance (from the sklearn.covariance package, . Make sure you take the transpose of the matrix_filt matrix before passing it to the empirical_covariance function!
Step9:  Use numpy.linalg.inv to get the inverse matrix.
Step10:  Use a double for loop to "scale" the negative of the precision matrix, which will give you the partial correlation. Print the dimensions of the matrix you get back, as a sanity check.
Step11:  Get the correlation coefficients and the partial correlation coefficients of the lower triangle of the matrix (not including the diagonal), as two vectors cor_values and pcor_values; your resulting vectors should each have length 13,366.  You will want to use numpy.tri and numpy.where (see class session 19 exercise)
Step12:  plot the histograms of the correlation coefficients (upper triangle only) and the partial correlation coefficients, on the same plot using alpha blending (refer to class session 17 exercise)
Step13:  Fisher transform the partial correlation values, using numpy.log
Step14:  Compute a p-value for each gene pair (upper triangle only), using the fact that sqrt(M-N-5) times the fisher Z sore should be approximately univariate normal (with zero mean) under the null hypothesis that a given gene pair's measurements (conditioned on the measurements for all the other 162 genes) are independent. You will want to use scipy.stats.norm.cdf, numpy.abs, and math.sqrt function (see class session 19 exercise).
Step15:  How many gene pairs have a P value less than 0.01? (use which and length)
Step16:  What are the sizes of the components in the undirected graph whose edges have P &lt; 0.05 in the statistical test that you did?  You will need to use zip, tolist, list, and igraph.Graph.TupleList (see class session 19 exercise)
Step17:  Plot the graph degree distribution on log-log scale
<ASSISTANT_TASK:>
Python Code:
import pandas  ## data file loading
import numpy
import sklearn.covariance  ## for covariance matrix calculation
import matplotlib.pyplot
import matplotlib
import pylab
import scipy.stats  ## for calculating the CDF of normal distribution
import igraph  ## for network visualization and finding components
import math
Explanation: CS446/546 - Class Session 20 - Partial correlation network
In this class session we will continue with analyzing the tumor gene expression dataset from the NIH human bladder cancer cohort (M=414 tumors), building on what we learned in Class Session 19 (Correlation Network).  In order to keep the analysis simple&ast;, in this notebook will restrict our analysis to a set of N=164 genes that are very highly expressed in bladder cancer. Using the 164 x 414 matrix of transcript abundance measurements, we will construct a network based on gene-gene partial correlation coefficients. We will also compare the distribution of partial correlation coefficients to the distribution of Pearson correlation coefficients. Do you think they will be different?  In what way would you expect them to be different?
&ast; Here, "simple" means that the covariance matrix will be nonsingular, so that we can obtain the partial correlation matrix by inversion.
We'll import all of the python modules that we will need for this exercise
End of explanation
gene_matrix_for_network_df = pandas.read_csv("shared/bladder_cancer_genes_tcga.txt", sep="\t")
Explanation: Read the tab-deliminted text file of gene expression measurements (rows correspond to genes, columns correspond to bladder tumor samples). (use Pandas, pandas.read.csv, and as_matrix). As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using shape.
End of explanation
gene_matrix_for_network = gene_matrix_for_network_df.as_matrix()
Explanation: Convert your data frame to a numpy matrix, using the pandas.DataFrame.as_matrix method.
End of explanation
print(gene_matrix_for_network.shape)
Explanation: As always, sanity check that the file that you loaded has the expected dimensions (4,473 x 414) using shape.
End of explanation
genes_median_expression = numpy.median(gene_matrix_for_network, axis=1)
Explanation: Compute the median expression level for each row of your matrix
End of explanation
gene_matrix_np = numpy.array(gene_matrix_for_network)
genes_keep = numpy.where(genes_median_expression > 12)
matrix_filt = gene_matrix_np[genes_keep, ][0]
matrix_filt.shape
N = matrix_filt.shape[0]
Explanation: Filter the matrix to include only rows for which the gene's median expression is > 12; matrix should now be 164 x 414; this will enable us to easily compute the partial correlation matrix using the inverse of the covariance matrix. Print the size of the filtered matrix, as a sanity check.
End of explanation
matrix_filt.shape
Explanation: Print the shape of your filtered matrix, as a sanity check.  It should be 164x414.
End of explanation
matrix_cor = numpy.corrcoef(matrix_filt)
Explanation: Compute the 164 x 164 matrix of gene-gene Pearson correlation coefficients, using numpy.corrcoef (this function treats each row as a random variable, so you don't have to do any transposing of the matrix, unlike the situation in R).
End of explanation
matrix_cov = sklearn.covariance.empirical_covariance(numpy.matrix.transpose(matrix_filt))
Explanation: Compute the covariance matrix using sklearn.covariance.empirical_covariance (from the sklearn.covariance package, . Make sure you take the transpose of the matrix_filt matrix before passing it to the empirical_covariance function!
End of explanation
matrix_cov_inv = numpy.linalg.inv(matrix_cov)
Explanation: Use numpy.linalg.inv to get the inverse matrix.
End of explanation
matrix_pcor = -matrix_cov_inv
for i in range(N):
    for j in range(N):
        matrix_pcor[i,j] /= numpy.sqrt(matrix_cov_inv[i,i]*matrix_cov_inv[j,j])
print(matrix_pcor.shape)
Explanation: Use a double for loop to "scale" the negative of the precision matrix, which will give you the partial correlation. Print the dimensions of the matrix you get back, as a sanity check.
End of explanation
cor_values = matrix_cor[numpy.where(numpy.tri(*matrix_cor.shape, k=-1))]
pcor_values = matrix_pcor[numpy.where(numpy.tri(*matrix_pcor.shape, k=-1))]
print(len(cor_values))
print(len(pcor_values))
Explanation: Get the correlation coefficients and the partial correlation coefficients of the lower triangle of the matrix (not including the diagonal), as two vectors cor_values and pcor_values; your resulting vectors should each have length 13,366.  You will want to use numpy.tri and numpy.where (see class session 19 exercise)
End of explanation
matplotlib.pyplot.hist(cor_values, normed=1, alpha=0.5, label="cor")
matplotlib.pyplot.hist(pcor_values, normed=1, alpha=0.5, label="pcor")
matplotlib.pyplot.legend(loc="upper left")
matplotlib.pyplot.xlabel("R")
matplotlib.pyplot.ylabel("frequency")
matplotlib.pyplot.show()
Explanation: plot the histograms of the correlation coefficients (upper triangle only) and the partial correlation coefficients, on the same plot using alpha blending (refer to class session 17 exercise)
End of explanation
z_scores = 0.5*numpy.log((1+pcor_values)/
                         (1-pcor_values))
Explanation: Fisher transform the partial correlation values, using numpy.log:
End of explanation
M = gene_matrix_for_network_df.shape[1]
P_values = 2*scipy.stats.norm.cdf(-numpy.abs(z_scores)*(math.sqrt((M-N-5))))
Explanation: Compute a p-value for each gene pair (upper triangle only), using the fact that sqrt(M-N-5) times the fisher Z sore should be approximately univariate normal (with zero mean) under the null hypothesis that a given gene pair's measurements (conditioned on the measurements for all the other 162 genes) are independent. You will want to use scipy.stats.norm.cdf, numpy.abs, and math.sqrt function (see class session 19 exercise).
End of explanation
len(numpy.where(P_values < 0.01)[0])
Explanation: How many gene pairs have a P value less than 0.01? (use which and length)
End of explanation
inds_tri = numpy.where(numpy.tri(*matrix_pcor.shape, k=-1))
inds_sig = numpy.where(P_values < 0.01)
graph_edge_list = list(zip(inds_tri[1][inds_sig].tolist(),
                                                inds_tri[0][inds_sig].tolist()))
final_network = igraph.Graph.TupleList(graph_edge_list, directed=False)
final_network.summary()
Explanation: What are the sizes of the components in the undirected graph whose edges have P &lt; 0.05 in the statistical test that you did?  You will need to use zip, tolist, list, and igraph.Graph.TupleList (see class session 19 exercise)
End of explanation
degree_dist = final_network.degree_distribution()
xs, ys = zip(*[(left, count) for left, _, count in degree_dist.bins()])
matplotlib.pyplot.loglog(xs, ys, marker="o")
pylab.xlabel("k")
pylab.ylabel("N(k)")
pylab.show()
Explanation: Plot the graph degree distribution on log-log scale
End of explanation
<END_TASK>