<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Imputation using UNCURL
Here, we will test imputation of missing data.
Step1:  Loading data
There are two datasets
Step2:  Artificially inducing missing values
One problem is that sometimes, all values for a gene are set to zero, which completely breaks uncurl. So we'll try not including these genes?
Step3:  Perform UNCURL for imputation
Step5:  Imputation with MW
Step6:  Imputation test with MAGIC
MAGIC is a scRNA-seq preprocessing method.
TODO
Step7:  Try with varying missing value fractions
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import random
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.io
from scipy import sparse
import uncurl
Explanation: Imputation using UNCURL
Here, we will test imputation of missing data.
End of explanation
data_z = scipy.io.loadmat('../data/GSE60361_dat.mat')
data_10x = scipy.io.loadmat('../data/10x_pooled_400.mat')
data_z_mat = data_z['Dat']
genes_z = uncurl.max_variance_genes(data_z_mat, 1, 0.2)
data_z_mat = data_z_mat[genes_z, :]
data_10x_mat = data_10x['data']
genes_10x = uncurl.max_variance_genes(data_10x_mat, 1, 0.2)
data_10x_mat = data_10x_mat[genes_10x, :]
data_z_mat.shape
data_10x_mat.shape
Explanation: Loading data
There are two datasets: 'z' is the Zeisel (2015) dataset, with 7 cell types. '10x' is a subset of pooled data from 8 different PBMC types, 50 cells from each type.
End of explanation
p_impute = 0.2
data_z_mat = sparse.coo_matrix(data_z_mat)
data_10x_mat = sparse.coo_matrix(data_10x_mat)
indices_z = random.sample(range(len(data_z_mat.data)), int(p_impute*len(data_z_mat.data)))
indices_10x = random.sample(range(len(data_10x_mat.data)), int(p_impute*len(data_10x_mat.data)))
# save deleted data for comparison
true_values_z = data_z_mat.data[indices_z]
sampled_rows_z = data_z_mat.row[indices_z]
sampled_cols_z = data_z_mat.col[indices_z]
true_values_10x = data_10x_mat.data[indices_10x]
sampled_rows_10x = data_10x_mat.row[indices_10x]
sampled_cols_10x = data_10x_mat.col[indices_10x]
# do the data deletion
data_z_mat.data[indices_z] = 0
data_z_mat.eliminate_zeros()
data_z_mat = sparse.csc_matrix(data_z_mat)
data_10x_mat.data[indices_10x] = 0
data_10x_mat.eliminate_zeros()
data_10x_mat = sparse.csc_matrix(data_10x_mat)
(data_z_mat.max(1).toarray() == 0).sum()
(data_10x_mat.max(1).toarray() == 0).sum()
Explanation: Artificially inducing missing values
One problem is that sometimes, all values for a gene are set to zero, which completely breaks uncurl. So we'll try not including these genes?
End of explanation
k_z = 7
%time M_z, W_z, ll = uncurl.poisson_estimate_state(data_z_mat, k_z, threads=4, disp=False)
k_10x = 8
%time M_10x, W_10x, ll = uncurl.poisson_estimate_state(data_10x_mat, k_10x, threads=4, disp=False)
Explanation: Perform UNCURL for imputation
End of explanation
def calculate_imputation_error(sampled_rows, sampled_cols, indices, true_values, new_matrix):
    
    Returns imputed values, rmse, correlation, spearman.
    
    imputed_values = []
    for i in range(len(indices)):
        imputed_value = new_matrix[sampled_rows[i], sampled_cols[i]]
        imputed_values.append(imputed_value)
    imputed_values = np.array(imputed_values)
    rmse = np.sqrt(np.mean(((imputed_values - true_values)/true_values)**2))
    print('rmse:', rmse)
    correlation = np.corrcoef(imputed_values, true_values)[0,1]
    print('correlation:', correlation)
    spearman = scipy.stats.spearmanr(imputed_values, true_values).correlation
    print('spearman:', spearman)
    return imputed_values, rmse, correlation, spearman
print('Results for 10x imputation with UNCURL:')
MW = M_10x.dot(W_10x)
imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_10x,
                                               sampled_cols_10x,
                                               indices_10x,
                                               true_values_10x, MW)
print()
print('Results for Zeisel imputation with UNCURL:')
MW = M_z.dot(W_z)
imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_z,
                                               sampled_cols_z,
                                               indices_z,
                                               true_values_z, MW)
Explanation: Imputation with MW
End of explanation
import magic
import pandas as pd
table_z = pd.DataFrame(data_z_mat.toarray().T)
table_z.columns = table_z.columns.astype(str)
scdata = magic.mg.SCData(table_z, data_type='sc-seq')
scdata = scdata.normalize_scseq_data()
scdata.run_magic(n_pca_components=k_z, random_pca=True,
        t=6, k=30, ka=10, epsilon=1, rescale_percent=99)
magic_output_z = scdata.magic.data.as_matrix().T
table_10x = pd.DataFrame(data_10x_mat.toarray().T)
table_10x.columns = table_10x.columns.astype(str)
scdata = magic.mg.SCData(table_10x, data_type='sc-seq')
scdata = scdata.normalize_scseq_data()
scdata.run_magic(n_pca_components=k_10x, random_pca=True,
        t=6, k=30, ka=10, epsilon=1, rescale_percent=99)
magic_output_10x = scdata.magic.data.as_matrix().T
print('Results for 10x imputation with Magic:')
MW = magic_output_10x
imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_10x,
                                               sampled_cols_10x,
                                               indices_10x,
                                               true_values_10x, MW)
print('Results for Zeisel imputation with Magic:')
imputed_values = []
MW = magic_output_z
imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_z,
                                               sampled_cols_z,
                                               indices_z,
                                               true_values_z, MW)
Explanation: Imputation test with MAGIC
MAGIC is a scRNA-seq preprocessing method.
TODO: perhaps include other methods such as scImpute?
End of explanation
p_impute_vals = [0.2, 0.1, 0.05, 0.01, 0.005, 0.001]
results = pd.DataFrame(columns=['dataset', 'method', 'metric', 'p_impute'])
data_points = []
for p_impute in p_impute_vals:
    print()
    print('p_impute:', str(p_impute))
    data_z_mat = data_z['Dat']
    genes_z = uncurl.max_variance_genes(data_z_mat, 1, 0.2)
    data_z_mat = data_z_mat[genes_z, :]
    data_10x_mat = data_10x['data']
    genes_10x = uncurl.max_variance_genes(data_10x_mat, 1, 0.2)
    data_10x_mat = data_10x_mat[genes_10x, :]
    
    data_z_mat = sparse.coo_matrix(data_z_mat)
    data_10x_mat = sparse.coo_matrix(data_10x_mat)
    indices_z = random.sample(range(len(data_z_mat.data)), int(p_impute*len(data_z_mat.data)))
    indices_10x = random.sample(range(len(data_10x_mat.data)), int(p_impute*len(data_10x_mat.data)))
    # save deleted data for comparison
    true_values_z = data_z_mat.data[indices_z]
    sampled_rows_z = data_z_mat.row[indices_z]
    sampled_cols_z = data_z_mat.col[indices_z]
    true_values_10x = data_10x_mat.data[indices_10x]
    sampled_rows_10x = data_10x_mat.row[indices_10x]
    sampled_cols_10x = data_10x_mat.col[indices_10x]
    # do the data deletion
    data_z_mat.data[indices_z] = 0
    data_z_mat.eliminate_zeros()
    data_z_mat = sparse.csc_matrix(data_z_mat)
    data_10x_mat.data[indices_10x] = 0
    data_10x_mat.eliminate_zeros()
    data_10x_mat = sparse.csc_matrix(data_10x_mat)
    
    # run uncurl
    M_10x, W_10x, ll = uncurl.poisson_estimate_state(data_10x_mat, k_10x, threads=4, disp=False)
    M_z, W_z, ll = uncurl.poisson_estimate_state(data_z_mat, k_z, threads=4, disp=False)
    
    print('Results for 10x imputation with UNCURL:')
    MW = M_10x.dot(W_10x)
    imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_10x,
                                               sampled_cols_10x,
                                               indices_10x,
                                               true_values_10x, MW)
    data_points.append(['10x', 'uncurl', 'rmse', rmse])
    data_points.append(['10x', 'uncurl', 'corr', corr])
    data_points.append(['10x', 'uncurl', 'sp', sp])
    print()
    print('Results for Zeisel imputation with UNCURL:')
    MW = M_z.dot(W_z)
    imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_z,
                                               sampled_cols_z,
                                               indices_z,
                                               true_values_z, MW)
    data_points.append(['z', 'uncurl', 'rmse', rmse])
    data_points.append(['z', 'uncurl', 'corr', corr])
    data_points.append(['z', 'uncurl', 'sp', sp])
    print()
    
    # run MAGIC
    table_z = pd.DataFrame(data_z_mat.toarray().T)
    table_z.columns = table_z.columns.astype(str)
    scdata = magic.mg.SCData(table_z, data_type='sc-seq')
    scdata = scdata.normalize_scseq_data()
    scdata.run_magic(n_pca_components=k_z, random_pca=True,
        t=6, k=30, ka=10, epsilon=1, rescale_percent=99)
    magic_output_z = scdata.magic.data.as_matrix().T
    
    table_10x = pd.DataFrame(data_10x_mat.toarray().T)
    table_10x.columns = table_10x.columns.astype(str)
    scdata = magic.mg.SCData(table_10x, data_type='sc-seq')
    scdata = scdata.normalize_scseq_data()
    scdata.run_magic(n_pca_components=k_10x, random_pca=True,
        t=6, k=30, ka=10, epsilon=1, rescale_percent=99)
    magic_output_10x = scdata.magic.data.as_matrix().T
    
    print()
    print('Results for 10x imputation with Magic:')
    MW = magic_output_10x
    imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_10x,
                                               sampled_cols_10x,
                                               indices_10x,
                                               true_values_10x, MW)
    data_points.append(['10x', 'magic', 'rmse', rmse])
    data_points.append(['10x', 'magic', 'corr', corr])
    data_points.append(['10x', 'magic', 'sp', sp])
    
    print()
    print('Results for Zeisel imputation with Magic:')
    imputed_values = []
    MW = magic_output_z
    imputed_vals, rmse, corr, sp = calculate_imputation_error(sampled_rows_z,
                                               sampled_cols_z,
                                               indices_z,
                                               true_values_z, MW)
    data_points.append(['z', 'magic', 'rmse', rmse])
    data_points.append(['z', 'magic', 'corr', corr])
    data_points.append(['z', 'magic', 'sp', sp])
    print()
results = pd.DataFrame(data_points, columns=['dataset', 'method', 'metric', 'p_impute'])
results.head()
# okay so this was an error... we have to correct the results.
results['value'] = results['p_impute']
results['p_impute'] = [x for y in [[p]*12 for p in p_impute_vals] for x in y]
results.head()
uncurl_rmse_10x = results[(results.method=='uncurl') &\
                          (results.dataset=='10x') &\
                          (results.metric=='rmse')]
uncurl_rmse_10x = uncurl_rmse_10x[['p_impute', 'value']]
uncurl_corr_10x = results[(results.method=='uncurl') &\
                          (results.dataset=='10x') &\
                          (results.metric=='corr')]
uncurl_corr_10x = uncurl_corr_10x[['p_impute', 'value']]
uncurl_sp_10x = results[(results.method=='uncurl') &\
                          (results.dataset=='10x') &\
                          (results.metric=='sp')]
uncurl_sp_10x = uncurl_sp_10x[['p_impute', 'value']]
uncurl_rmse_z = results[(results.method=='uncurl') &\
                          (results.dataset=='z') &\
                          (results.metric=='rmse')]
uncurl_rmse_z = uncurl_rmse_z[['p_impute', 'value']]
uncurl_corr_z = results[(results.method=='uncurl') &\
                          (results.dataset=='z') &\
                          (results.metric=='corr')]
uncurl_corr_z = uncurl_corr_z[['p_impute', 'value']]
uncurl_sp_z = results[(results.method=='uncurl') &\
                          (results.dataset=='z') &\
                          (results.metric=='sp')]
uncurl_sp_z = uncurl_sp_z[['p_impute', 'value']]
magic_rmse_10x = results[(results.method=='magic') &\
                          (results.dataset=='10x') &\
                          (results.metric=='rmse')]
magic_rmse_10x = magic_rmse_10x[['p_impute', 'value']]
magic_corr_10x = results[(results.method=='magic') &\
                          (results.dataset=='10x') &\
                          (results.metric=='corr')]
magic_corr_10x = magic_corr_10x[['p_impute', 'value']]
magic_sp_10x = results[(results.method=='magic') &\
                          (results.dataset=='10x') &\
                          (results.metric=='sp')]
magic_sp_10x = magic_sp_10x[['p_impute', 'value']]
magic_rmse_z = results[(results.method=='magic') &\
                          (results.dataset=='z') &\
                          (results.metric=='rmse')]
magic_rmse_z = magic_rmse_z[['p_impute', 'value']]
magic_corr_z = results[(results.method=='magic') &\
                          (results.dataset=='z') &\
                          (results.metric=='corr')]
magic_corr_z = magic_corr_z[['p_impute', 'value']]
magic_sp_z = results[(results.method=='magic') &\
                          (results.dataset=='z') &\
                          (results.metric=='sp')]
magic_sp_z = magic_sp_z[['p_impute', 'value']]
plt.figure(dpi=100)
plt.semilogx(uncurl_rmse_10x.p_impute, uncurl_rmse_10x.value, '--o', label='uncurl')
plt.semilogx(magic_rmse_10x.p_impute, magic_rmse_10x.value, '--o', label='magic')
plt.legend()
plt.title('RMSE of imputed values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('RMSE')
plt.figure(dpi=100)
plt.semilogx(uncurl_corr_10x.p_impute, uncurl_corr_10x.value, '--o', label='uncurl')
plt.semilogx(magic_corr_10x.p_impute, magic_corr_10x.value, '--o', label='magic')
plt.legend()
plt.title('Correlation of imputed values vs actual values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('Correlation')
plt.figure(dpi=100)
plt.semilogx(uncurl_sp_10x.p_impute, uncurl_sp_10x.value, '--o', label='uncurl')
plt.semilogx(magic_sp_10x.p_impute, magic_sp_10x.value, '--o', label='magic')
plt.legend()
plt.title('Spearman Correlation of imputed values vs actual values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('Spearman Correlation')
plt.figure(dpi=100)
plt.semilogx(uncurl_rmse_z.p_impute, uncurl_rmse_z.value, '--o', label='uncurl')
plt.semilogx(magic_rmse_z.p_impute, magic_rmse_z.value, '--o', label='magic')
plt.legend()
plt.title('RMSE of imputed values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('RMSE')
plt.figure(dpi=100)
plt.semilogx(uncurl_corr_z.p_impute, uncurl_corr_z.value, '--o', label='uncurl')
plt.semilogx(magic_corr_z.p_impute, magic_corr_z.value, '--o', label='magic')
plt.legend()
plt.title('Correlation of imputed values vs actual values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('Correlation')
plt.figure(dpi=100)
plt.semilogx(uncurl_sp_z.p_impute, uncurl_sp_z.value, '--o', label='uncurl')
plt.semilogx(magic_sp_z.p_impute, magic_sp_z.value, '--o', label='magic')
plt.legend()
plt.title('Spearman Correlation of imputed values vs actual values for UNCURL and Magic')
plt.xlabel('Fraction removed')
plt.ylabel('Spearman Correlation')
Explanation: Try with varying missing value fractions
End of explanation
<END_TASK>