In [1]:
# para plotar diretamente no browser usando ipython notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import math
import operator
import scipy.spatial.distance as distance
import sklearn.preprocessing as preprocessing
import scipy.cluster.hierarchy as sch
import matplotlib.gridspec as gridspec
import sys
from sklearn.neighbors import NearestNeighbors
####Confirmar que a versão da biblioteca pandas é superior ou igual a 0.15
#print pd.__version__
In [2]:
####Carregar o arquivo com os dados
Data = pd.read_csv('http://www.molbiolcell.org/content/suppl/2003/09/24/9.12.3273.DC1/CDCDATA.txt', sep='\t', skip_blank_lines=True)
#Data = pd.read_csv('CDCDATA.txt', sep='\t', skip_blank_lines=True)
In [3]:
####Remoção dos NaN####
####Nas linhas#####
Data.dropna(how="all", inplace=True)
In [4]:
####Filtrar as linhas com menos que 2 NaN's
data_clean = Data[Data.isnull().sum(axis=1) <= 2]
In [5]:
####Renomeando primeira coluna para Nome Gene####
Data = data_clean.rename(columns={'Unnamed: 0': 'Nome Gene'})
In [6]:
#####Separando Nome dos genes e os valores numéricos####
D_GNames=Data['Nome Gene']
del Data['Nome Gene']
D_Matrix=Data
#print 'Número de linhas (Matriz Numérica): {0}'.format(D_Matrix.shape[0])
#print 'Número de colunas (Matriz Numérica): {0}'.format(D_Matrix.shape[1])
#print 'Tamanho da Matriz numérica: {0}'.format(D_Matrix.shape)
#print 'Número de genes no array Nome gene: {0}'.format(D_GNames.shape[0])
In [7]:
# não mais usado
def distancia(gene1_index, gene2_index):
distance = 0
# length = 0
for col in range(len(D_Matrix.columns)):
if not (math.isnan(D_Matrix.iloc[gene1_index,col]) or math.isnan(D_Matrix.iloc[gene2_index,col])):
distance += math.pow(math.log(D_Matrix.iloc[gene1_index,col]+5.0) - math.log(D_Matrix.iloc[gene2_index,col]+5.0), 2.0)
# length += 1
return math.sqrt(distance)
In [8]:
# não mais usado
def KNN (gene1_index, k):
distancias = []
for row in range(len(D_Matrix.index)):
if (row != gene1_index):
dist = distancia(gene1_index,row)
distancias.append((row,dist))
distancias.sort(key=operator.itemgetter(1))
vizinhos = []
for i in range(k):
vizinhos.append((distancias[i][0],-distancias[i][1]))
return vizinhos
#print [ y for x, y in vizinhos]
In [9]:
# este trecho está levando tempo demais para executar. Embora não seja um primor de código otimizado,
# o código equivalente em C++ rodou em menos de 2 minutos - aqui levou 8 horas e não tinha passado de 80 linhas executadas
# Python é lento demais da conta!!!!!!!!!!!!!
#for row in range(len(D_Matrix.index)):
# # para acompanhar a execução
# print "Linha {linha:03d}".format(linha=row)
# sys.stdout.flush()
# for col in range(len(D_Matrix.columns)):
# # para acompanhar a execução
# #print " Coluna {coluna:03d}".format(coluna=col)
# if math.isnan(D_Matrix.iloc[row,col]):
# vizinhos = KNN(row,15)
# vizinhos_valor = map(D_Matrix.iloc[operator.itemgetter(0),col], vizinhos)
# vizinhos_peso = map(D_Matrix.iloc[operator.itemgetter(1),col], vizinhos)
# novo_valor = np.average(vizinhos_valor,weights=vizinhos_peso)
# #novo_valor = sum([ D_Matrix.iloc[x,col] * y for x, y in vizinhos if not math.isnan(D_Matrix.iloc[x,col])]) / \
# # sum([ y for x, y in vizinhos if not math.isnan(D_Matrix.iloc[x,col])])
# D_Matrix.iloc[row,col] = novo_valor
In [10]:
# tentando ver se o que fizeram é menos lento...
X = D_Matrix.copy() #se nao fizer o copy, ao mudar X o DD_Matrix muda tbm. n sei o porquê
# evitando erros por tentar fazer log de número negativo
for i in range(len(D_Matrix.index)):
sys.stdout.flush()
for j in range(len(D_Matrix.columns)):
sys.stdout.flush()
if X.iloc[i,j] != np.nan:
X.iloc[i,j] += 5.0
# agora sim, tudo certo
X = np.log(X)
In [11]:
# Método sugerido por outro aluno
# rodou mais rápido, manter isto. (nunca vi como python é cheio de gambiarras para conseguir rodar mais rápido.)
## ============== ##
## imputar os KNN ##
## ============== ##
Y = X.dropna()
for col1 in range(len(D_Matrix.columns)):
#print "Coluna1 {coluna:03d}".format(coluna=col1)
#sys.stdout.flush()
for col2 in range(col1+1,len(D_Matrix.columns)):
#print "Coluna2 {coluna:03d}".format(coluna=col2)
#sys.stdout.flush()
pair = (col1, col2)
## listar com todas as colunas, exceto col1 e col2
filtraExp = list(set(range(len(D_Matrix.columns))) - set(pair))
## genes alvo = todos os genes que NÂO tenham NA fora da col1 ou col2
genesAlvo = list(X[X[filtraExp].isnull().sum(axis=1) == 0].index)
## dentre os genes alvo, selecionar os que possuem NA
e1NaN = X.ix[genesAlvo,][X.ix[genesAlvo,col1].isnull()].index
e2NaN = X.ix[genesAlvo,][X.ix[genesAlvo,col2].isnull()].index
## se houver algum gene,
if len(e1NaN)+len(e2NaN) > 0:
neigh = NearestNeighbors(15, p=2).fit(Y[filtraExp])
distances, indices = neigh.kneighbors(X.ix[e1NaN | e2NaN,filtraExp])
for i,gene in enumerate(X.ix[e1NaN | e2NaN, filtraExp].index):
if gene in e1NaN:
#print (i,gene)
neighbors = [Y.iloc[x,col1] for x in indices[i,]]
similarity = np.multiply(distances[i,],-1)
weightedMean = np.average(neighbors,weights=similarity) ## calcula média ponderada da expressao usando as disntancias como peso
X.ix[gene, col1] = weightedMean
if gene in e2NaN:
#print (i,gene)
neighbors = [Y.iloc[x,col2] for x in indices[i,]]
similarity = np.multiply(distances[i,],-1)
weightedMean = np.average(neighbors,weights=similarity)
X.ix[gene, col2] = weightedMean
In [12]:
# normalizando
D_Matrix_scaled = []
for row in range(len(X.index)):
D_Matrix_scaled.append(preprocessing.scale(np.asarray(X.iloc[row]), axis=0, with_mean=True, with_std=True, copy=True))
D_Matrix_scaled = pd.DataFrame(D_Matrix_scaled)
#print D_Matrix_scaled
In [13]:
pairwise_dists = distance.squareform(distance.pdist(D_Matrix_scaled,metric='correlation'))
pairwise_dists_columns = distance.squareform(distance.pdist(D_Matrix_scaled.T,metric='correlation'))
In [14]:
# cluster
clusters = sch.linkage(pairwise_dists, method='average')
In [15]:
# helper for cleaning up axes by removing ticks, tick labels, frame, etc.
def clean_axis(ax):
"""Remove ticks, tick labels, and frame from axis"""
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
for sp in ax.spines.values():
sp.set_visible(False)
In [16]:
testDF = pd.DataFrame(np.array(D_Matrix_scaled, dtype='float'))
fig = plt.figure(figsize=(10, 100))
heatmapGS = gridspec.GridSpec(1,2,wspace=0.0,hspace=0.0,width_ratios=[0.25,1])
### row dendrogram ###
denAX = fig.add_subplot(heatmapGS[0,0])
denD = sch.dendrogram(clusters,orientation='right')
clean_axis(denAX)
### heatmap ###
heatmapAX = fig.add_subplot(heatmapGS[0,1])
axi = heatmapAX.imshow(testDF.ix[denD['leaves']],interpolation='nearest',aspect='auto',origin='lower',cmap=plt.cm.RdBu)
clean_axis(heatmapAX)
In [ ]: