Clusterização de um conjunto de dados de expressão gênica e viabilidade celular - Desafio 3º Imersão Dados da Alura¶
Meu portfolio: https://edugvs.github.io/
1. Importando bibliotecas necessárias¶
# Para manipulação dos dados
import pandas as pd
import numpy as np
# Para visualização dos dados
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import pylab
%matplotlib inline
# Para modelagem preditiva
from sklearn.cluster import KMeans
from sklearn.metrics import homogeneity_completeness_v_measure
from sklearn.model_selection import train_test_split
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import scale
# Redução de dimensionalidade
from sklearn.decomposition import PCA
# Para o cálculo de algumas estatísticas
from scipy.stats import kurtosis, skew
from scipy.spatial.distance import cdist, pdist
from pyclustertend import hopkins
# Para ignorar os avisos
import warnings
warnings.filterwarnings("ignore")
2. Dicionário de dados¶
| Variável | Descrição |
|---|---|
| id | ID da amostra |
| tratamento | Indica o grupo controle e o grupo tratado |
| tempo | Tempo em horas (h) de exposição a droga |
| dose | Indica quantas doses foram utilizadas no experimento (D1/D2) |
| droga | ID da droga utilizada no tratamento |
| g-0 : g-771 | Expressão gênica da amostra, quanto maior, mais expresso o gene |
| c-0 : c-99 | Viabilidade celular da amostra, quanto menor, menor é a viabilidade celular |
3. Carregando a base de dados¶
# Carregando os dados
url_dados = 'https://github.com/alura-cursos/imersaodados3/blob/main/dados/dados_experimentos.zip?raw=true'
dados = pd.read_csv(url_dados, compression = 'zip')
display(dados)
# Criando uma cópia do dataset
df = dados
É sempre uma boa prática criarmos uma cópia do dataset para termos a base de dados original em caso de necessidade.
# Visualizando o tipo de objeto dos dados
type(df)
# Checando se há dados duplicados
df.duplicated().sum()
# Checando se há valores ausentes
sns.heatmap(df.isnull(), cbar=False)
Não encontramos nenhum valor NA na base de dados.
# Verificando as dimensões do objeto
df.shape
O dataset contém 23814 registros (linhas) e 877 variáveis (colunas).
# Contando valores únicos
info = df[['id', 'tratamento', 'tempo', 'dose', 'droga', 'g-0', 'c-0']].nunique().sort_values()
# Determinando o tipo de dados para cada variável
info = pd.DataFrame(info.values, index = info.index, columns = ['NUniques'])
# Atribuindo as informações sobre o tipo de dados das variáveis em um DataFrame
info['dtypes'] = df.dtypes
# Exibe o dataframe
display(info)
Aparentemente os dados foram carregados corretamente, com o tipo de dado correspondendo ao que cada variável representa.
4. Análise exploratória dos dados¶
Criando uma função para nos auxiliar no cálculo de algumas estatísticas fundamentais.
# Definindo uma função para gerar um dataframe com estatísticas de variáveis numéricas.
def varStats(col, data, target = ''):
if target == '':
stats = pd.DataFrame({
'Min' : data[col].min(),
'Q1' : data[col].quantile(.25),
'Median': data[col].median(),
'Mean' : data[col].mean(),
'Q3' : data[col].quantile(.75),
'Max' : data[col].max(),
'SD' : data[col].std(),
'SK' : skew(data[col]),
'KU' : kurtosis(data[col])
}, index = [col])
else:
stats = pd.concat([
df[[col, target]].groupby(target).min(),
df[[col, target]].groupby(target).quantile(.25),
df[[col, target]].groupby(target).median(),
df[[col, target]].groupby(target).mean(),
df[[col, target]].groupby(target).quantile(.75),
df[[col, target]].groupby(target).max(),
df[[col, target]].groupby(target).std(),
df[[col, target]].groupby(target).skew(),
df[[col, target]].groupby(target).apply(lambda group: kurtosis(group)[0])
], axis = 1)
stats.columns = ['Min', 'Q1', 'Median', 'Mean', 'Q3', 'Max', 'SD', 'SK', 'KU']
return stats
O coeficiente de assimetria (Skewness) indica como os dados são distribuídos e para interpretar seu resultado podemos olhar a seguinte tabela:
| Skewness | Descrição |
|---|---|
| SK ≈ 0 | Os dados são simétricos. Tanto a cauda direita quanto a esquerda da função de densidade de probabilidade são iguais |
| SK < 0 | A assimetria é negativa. A cauda do lado esquerdo da função de densidade de probabilidade é maior do que a cauda à direita. |
| SK > 0 | A assimetria é positiva. A cauda do lado direito da função de densidade de probabilidade é maior do que a cauda do lado esquerdo. |
from IPython.display import Image
Image("img/skew.png")
O coeficiente de curtose (curtose) é uma medida que caracteriza o achatamento da curva da função de distribuição e para interpretar o seu resultado podemos observar a seguinte tabela:
| Kurtosis | Descrição |
|---|---|
| KU ≈ 0 | A distribuição é normal e é chamada de curtose mesocúrtica |
| KU < 0 | A curva é mais plana do que o normal. Para um coeficiente de curtose negativo existe uma curtose platicúrtica. |
| KU > 0 | A curva é mais proeminente do que o normal. Para um coeficiente de curtose positivo, existe uma curtose leptocúrtica. |
from IPython.display import Image
Image("img/kurtosis.png")
# Resumo das variáveis categóricas
df.describe(include=object)
Com essa tabela vemos que temos 2 valores posssíveis para dose (D1/D2) e para tratamento (Tratado/Controle). é interessante observar que a droga 'cacb2b860' é o valor mais frequente, com um total de 1866 registros. Devemos investigar isso mais detalhadamente.
# Resumo das variáveis numéricas
df.describe()
É praticamente impossivel analisar cada gene e viabilidade celular olhando somente essa tabela. Entretanto conseguimos tirar uma informação valiosa ao observar os valores mínimos e máximos das variáveis que representam a expressão gênica (g) e viabilidade celular (c). Os valorem oscilam entre -10.0 e 10.0, o que nos mostra que esses dados já estão normalizados.
No caso da análise da expressão gênica, o que os pesquisadores geralmente fazem é converter os resultados da amplificação gênica pelo RT-PCR (ou qualquer outro método de biologia molecular para amplificação de ácidos nucleicos) em um valor de Log na base 2, o que pode explicar o porquê dos dados estarem nessa escala.
Começaremos nossa análise exploratória olhando para a variável 'droga', descobrindo quantas drogas distintas foram utilizadas nesses experimentos. Nesse conjunto de dados foram utilizadas mais de 3000 drogas, será que há alguma diferença na frequência que elas aparecem?
# Contando os valores da coluna 'droga' e retornando as 10 mais utilizadas
df['droga'].value_counts().head(10)
# Retornando os dados em porcentagem
dados['droga'].value_counts(normalize = True).map('{:.2%}'.format).head(10)
Podemos observar que apenas 9 drogas, dentre as mais de 3000, aparecem com muito mais frequência do que a média.
# Criando uma variável ordenada com as drogas mais utilizadas
cod_drogas = df['droga'].value_counts().index[0:10]
# Criando um grafico de barras
sns.set()
plt.figure(figsize=(15, 5))
ax = sns.countplot(x = 'droga', data=df.query('droga in @cod_drogas'), order=cod_drogas, hue='tratamento')
ax.set_title('Top 10 drogas utilizadas')
plt.show()
A droga "cacb2b860" na verdade é um código para demonstrar que a amostra não recebeu tratamento nenhum, por isso está com mais registros do que as outras.
# Verificando a proporção entre os dados na var 'tratamento'
df['tratamento'].value_counts(normalize = True).map('{:.2%}'.format)
Aparentemente os dados estão desbalanceados entre os grupos tratado e controle. Uma hipótese para explicar essa diferença é que não precisamos utilizar um controle para cada amostra, geralmente os pesquisadores comparam vários grupos tratado com um controle, desde que a única diferença entre os grupos seja a variável de interesse, que no caso é o tratamento com a droga específica.
# Verificando a proporção entre os dados na var 'dose'
df['dose'].value_counts(normalize = True).map('{:.2%}'.format)
# Verificando a proporção entre os dados na var 'tempo'
df['tempo'].value_counts(normalize = True).map('{:.2%}'.format)
Para as outras variáveis os dados estão bem distribuidos.
pd.crosstab([df['dose'], df['tempo']], df['tratamento'])
Com o crosstab podemos ver melhor como os dados estão organizados. O experimento foi desenhado para utilizar duas doses em três janelas de tempo diferentes (24h, 48h e 72h), sempre comparando o grupo tratado com o grupo controle. Ao observar as alterações nos níveis de expressão gênica e na viabilidade celular, podemos inferir se alguma droga está modulando determinados genes e como isso está afetando a viabilidade das células.
4.1 Análise das distribuições¶
# Criando um histograma para verificar a distribuição das médias das expressões gênicas de g-0 até g-771
df.loc[:,'g-0':'g-771'].describe().T['mean'].hist(bins=30)
# Criando uma variável para receber os valores das médias das expressões gênicas de g-0 até g-771
g_mean_values = df.loc[:,'g-0':'g-771'].describe().T['mean']
# Transformando em uma série do pandas
pd.Series(g_mean_values)
# Ordenando a série em ordem decrescente
g_mean_values.sort_values(ascending=False)
Podemos perceber que o gene 'g-707' é o que apresenta a maior média de expressão, seguidos dos genes 'g-100', 'g-744' e 'g-392'. Enquanto os genes 'g-370', 'g-508' e 'g-37' apresentaram a menor média de expressão.
# Transformando a série em um df do pandas
g_mean_values_df = g_mean_values.to_frame()
# Definindo a coluna para calcular as estatísticas
col = 'mean'
# Aplicando a função varStats no conjunto de dados
varStats(col, data = g_mean_values_df)
A distribuição dos dados se assemelha muito a uma curva normal. Podemos observar uma assimetria levemente negativa, indicando que cauda do lado esquerdo da distribuição é maior do que a cauda à direita. A kurtosis levemente positiva indica que a curva é mais proeminente do que o normal, os dados estão levemente mais concentrados em torno da média da distribuição com um desvio padrão baixo.
df.loc[:,'c-0':'c-99'].describe().T['mean'].hist(bins=30)
# Criando uma variável para receber os valores das médias das expressões gênicas de g-0 até g-771
c_mean_values = df.loc[:,'c-0':'c-99'].describe().T['mean']
# Transformando em uma série do pandas
pd.Series(c_mean_values)
# Ordenando a série em ordem decrescente
c_mean_values.sort_values(ascending=False)
Observamos que as variáveis que representam a viabilidade celular que apresentaram a maior média foram as 'c-74, 'c-37', 'c-7', 'c-76' e 'c-69'. Enquanto as menores médias foram as variáveis 'c-65', 'c-18', 'c-63', 'c-38' e 'c-26' respectivamente. É interessante observar que todas as médias foram negativas, o que pode indicar baixa viabilidade celular para todas as variáveis do dataset.
# Transformando a série em um df do pandas
c_mean_values_df = c_mean_values.to_frame()
# Definindo a coluna para calcular as estatísticas
col = 'mean'
# Aplicando a função varStats no conjunto de dados
varStats(col, data = c_mean_values_df)
A distribuição dos dados se assemelha muito a uma curva normal. Podemos observar uma assimetria bem positiva, indicando que cauda do lado direito da distribuição é maior do que a cauda à direita, os dados estão concentrados no lado esquerdo da distribuição. A kurtosis fortemente positiva indica que a curva é bem mais proeminente do que o normal, os dados estão muito mais concentrados em torno da média da distribuição (-0,432231) com um desvio padrão baixo (0,093).
Para esse etapa, iremos selecionar apenas 10 genes e 10 dados de viabilidade celular por motivos de tempo hábil para executar o projeto. Os critérios serão baseados nas top 5 maiores e menores médias de expressão para as variáveis de expressão gênica e o mesmo critério também será aplicado para as variáveis de viabilidade celular.
# Criando duas lista com as variáveis selecionadas
# Para as variáveis de viabilidade celular
c_list = ['c-74',
'c-37',
'c-7',
'c-76',
'c-69',
'c-26',
'c-38',
'c-63',
'c-18',
'c-65']
# Para as variáveis de expressão gênica
g_list = ['g-707',
'g-100',
'g-744',
'g-392',
'g-38',
'g-672',
'g-50',
'g-37',
'g-508',
'g-370']
# Criando plot para as distribuições dos dados para cada variável
plt.figure(figsize=(16,16))
sns.set_style('white')
genes = np.random.choice(len(g_list),16)
for i,col in enumerate(genes):
plt.subplot(4,4,i+1)
plt.hist(df.loc[:,g_list[col]],bins=100,color='blue')
plt.title(g_list[col])
Podemos verificar com esses histogramas como estão distribuidos os dados de expressão gênica para cada gene dentro dos top 10 selecionados.
# Criando plot para as distribuições dos dados para cada variável
plt.figure(figsize=(16,16))
sns.set_style('white')
viabilidade = np.random.choice(len(c_list),16)
for i,col in enumerate(viabilidade):
plt.subplot(4,4,i+1)
plt.hist(df.loc[:,c_list[col]],bins=100,color='green')
plt.title(c_list[col])
Podemos verificar com esses histogramas como estão distribuidos os dados de viabilidade celular para cada variável dentro das top 10 selecionados.
4.2 Análise de correlação¶
Analisaremos a correlação entre as features selecionadas anteriormente.
# Criando um dataframe somente com as variáveis de interesse
var_df = df.iloc[:,[42,43,55,105, 375, 397, 513, 677, 712, 749, 784, 795, 803, 815, 814, 840, 846, 842, 851, 853]]
# Criando um heatmap para a análise correlação
plt.figure(figsize=(20, 20))
heatmap = sns.heatmap(var_df.corr(), vmin=-1, vmax=1, annot=True)
heatmap.set_title('Heatmap de Correlação', fontdict={'fontsize':22}, pad=12);
A partir desse mapa de correlação podemos inferir que alguns genes se correlacionam tanto positivamente como negativamente com a viabilidade celular. Em outras palavras, existem situações em que a expressão gênica aumenta e a viabilidade celular também (correlação positiva) e casos em que a expressão gênica diminui/aumenta e a viabilidade celular aumenta/diminui (correlação negativa). Vamos pegar, apenas como exemplo, o gene 'g-38' que está correlacionado negativamente a praticamente todas as variáveis de viabilidade celular e o gene 'g-672' que está correlacionado positivamente. Para fins de demonstração, vamos selecionar a var 'c-26', será que a expressão aumenta quando a viabilidade diminui ou o contrário?
O ideal seria analisar todos os genes do dataset, com todas as combinações possíveis, vamos escolher esses dois apenas para ilustrar o conceito.
4.2.1 Agrupando os dados de correlação por clusters¶
# Criando um dataframe somente com as variáveis de interesse
var_df = df.iloc[:,[42,43,55,105, 375, 397, 513, 677, 712, 749, 784, 795, 803, 815, 814, 840, 846, 842, 851, 853]]
# Criando um clustermap para a análise correlação em grupos
clustermap = sns.clustermap(var_df.corr(), vmin=-1, vmax=1, annot=False, figsize=(15, 15), method='centroid', metric='euclidean')
Podemos ver dois grandes clusters dentro do mapa de correlação, os genes 'g-392', 'g-100', 'g-38' e o 'g-744' se correlacionam positivamente entre si e negativamente com as outras variáveis. Os genes 'g-370' e 'g-707' não se correlacionam com nenhuma outra variável, um possível explicação para isso é que esses genes possam ser controle. Esse tipo de visualização do mapa de correlação pode ser muito útil para identificar grupos que se correlacionam entre si.
# Definindo a variável a ser calculada
col = 'c-26'
# Aplicando a função
varStats(col, target = 'tempo', data = df)
Somente olhando essas estatísticas poderíamos inferir que quanto maior o tempo e exposição a droga, menor a viabilidade celular. Contudo, precisaríamos rodar testes estatísticos específicos (como um ANOVA) para verificar se essa diferença é significante.
# Definindo a variável a ser calculada
col = 'g-38'
# Aplicando a função
varStats(col, target = 'tempo', data = df)
Somente olhando essas estatísticas poderíamos inferir que quanto maior o tempo e exposição a droga, maior a expressão gênica. Contudo, precisaríamos rodar testes estatísticos específicos (como um ANOVA) para verificar se essa diferença é significante.
# Configurando o grid e o scatterplot
plt.figure(figsize=(10, 5))
a_plot = sns.lmplot(data=df, x='g-38', y='c-26', line_kws={'color': 'red'}, col='dose', row='tempo', hue="tratamento", palette="mako")
# Ajustando a escala
a_plot.set(xlim=(-10, 10))
a_plot.set(ylim=(-10, 10))
Podemos afirmar que quando a expressão gênica de g-38 aumenta a viabilidade tende a diminuir. Pensando nos efeitos biológios, a droga está afetando a viabilidade celular modulando a expressão gênica desse gene específico.
# # Definindo a variável a ser calculada
col = 'g-672'
# Aplicando a função
varStats(col, target = 'tempo', data = df)
Somente olhando essas estatísticas poderíamos inferir que quanto maior o tempo e exposição a droga, menor a expressão gênica. Contudo, precisaríamos rodar testes estatísticos específicos (como um ANOVA) para verificar se essa diferença é significante.
# Configurando o grid e o scatterplot
plt.figure(figsize=(10, 5))
b_plot = sns.lmplot(data=df, x='g-672', y='c-26', line_kws={'color': 'red'}, col='dose', row='tempo', hue="tratamento", palette="mako")
# Ajustando a escala
b_plot.set(xlim=(-10, 10))
b_plot.set(ylim=(-10, 10))
Podemos afirmar que quando a expressão gênica de g-672 aumenta a viabilidade tende a aumentar. Pensando nos efeitos biológios, a droga está afetando a viabilidade celular modulando a expressão gênica desse gene específico.
4.3 Analisando a distribuição entre os grupos com boxplots¶
# boxplot controle
plt.figure(figsize=(12,6))
sns.set_style("whitegrid")
sns.boxplot(data = df.query('droga == "cacb2b860"'), y= 'g-38', x='dose', hue='tempo', linewidth=0.7, palette='mako')
plt.title('Boxplot controle (cacb2b860)')
Perceba como não há diferenças significativas no grupo controle. Precisamos sempre utilizar esse grupo ao fazer comparações entre os grupos.
# Criando boxplots para cada droga dentro das top 10 mais utilizadas par avaliar a expressão do gene g-38 entre as doses
# e os diferentes intervalos de tempo.
for droga in cod_drogas[1:]: # Para retirar o código do controle 'cacb2b860'
fig, axs = plt.subplots(1,2,figsize=(15,8))
sns.set_style("whitegrid")
sns.boxplot(data = df.query('droga == "cacb2b860"'), y= 'g-38', x='dose', hue='tempo', linewidth=0.7, palette='mako', ax=axs[0])
axs[0].set_title('Boxplot controle')
sns.boxplot(data = df.query('droga == @droga'), y= 'g-38', x='dose', hue='tempo', linewidth=0.7, palette='mako', ax=axs[1])
axs[1].set_title('Boxplot droga ' + droga )
plt.show()
Ao comparar os grupos podemos inferir que cada uma desses drogas está modulando a expressão de g-38. Ainda assim, seriam necessários alguns testes estatísticos para comprovar essa hipótese.
# Criando boxplots para cada droga dentro das top 10 mais utilizadas par avaliar a viabilidade celular na var 'c-26' entre as
# doses e os diferentes intervalos de tempo.
for droga in cod_drogas[1:]: # Para retirar o código do controle 'cacb2b860'
fig, axs = plt.subplots(1,2,figsize=(15,8))
sns.set_style("whitegrid")
sns.boxplot(data = df.query('droga == "cacb2b860"'), y= 'c-26', x='dose', hue='tempo', linewidth=0.7, palette='mako', ax=axs[0])
axs[0].set_title('Boxplot controle')
sns.boxplot(data = df.query('droga == @droga'), y= 'c-26', x='dose', hue='tempo', linewidth=0.7, palette='mako', ax=axs[1])
axs[1].set_title('Boxplot droga ' + droga )
plt.show()
Com essa série de boxplots, podemos perceber como diferentes drogas alteram o nível de expressão gênica, sempre comparadas com o grupo controle. Os intervalos de tempo também são importantes para modular a expressão, no caso do gene 'g-38' e da viabilidade celular 'c-26' que foram analisados, além disso, podemos comparar também a variável dose e verificar se existe alguma diferença.
Está análise é só um exemplo do que pode ser feito, o ideal seria analisar todos os genes e realizar testes estatísticos para verificar se as diferenças são significativas e depois descobrir que droga produz que tipo de efeito nos genes: superexpressão ou subexpressão gênica e como altera a viabilidade celular.
5. Carregando o segundo dataset, com os dados dos mecanismos de ação das drogas.¶
Conceitos importantes:
| Mecanismo de ação | Descrição |
|---|---|
| Agonista | Droga capaz de se ligar a um receptor celular e ativá-lo para provocar uma resposta biológica. |
| Antagonista | Droga capaz de se ligar a um receptor celular e bloqueá-lo para provocar uma resposta biológica. |
resultados_df = pd.read_csv('https://github.com/alura-cursos/imersaodados3/blob/main/dados/dados_resultados.csv?raw=true')
resultados_df.head()
# Contando valores únicos
info = resultados_df.drop('id', axis=1).sum().sort_values(ascending=False)
# Determinando o tipo de dados para cada variável
info = pd.DataFrame(info.values, index = info.index, columns = ['Count'])
# Atribuindo as informações sobre o tipo de dados das variáveis em um DataFrame
info['dtypes'] = resultados_df.dtypes
# Exibe o dataframe
display(info)
Mecanismos de ação que possuem maior quantidade de resultados positivos.
Neste momento, vamos mesclar os dois dataframes para criarmos algumas visualizações interessantes, mas antes disso vamos criar uma nova coluna chamada 'n_moa', que representa o número de mecanismos de ação ativos.
# Criando a var n_moa somando todos os valores em cada linha do df
resultados_df['n_moa'] = resultados_df.drop('id', axis=1).sum(axis=1)
# Verificando se a coluna apresenta resultados diferentes de 0. Caso algum mecanismo de ação tenha sido ativado o resultado
# sera > 0
resultados_df['n_moa'] != 0
# Add a coluna 'ativo_moa' no dataframe
resultados_df['ativo_moa'] = (resultados_df['n_moa'] != 0)
resultados_df.head()
# Unindo os dois dataframes pelo ID
df_merged = pd.merge(df, resultados_df[['id','n_moa', 'ativo_moa']], on='id')
df_merged.head()
Agora temos em um mesmo dataframe todas as variáveis de expressão gênica, viabilidade celular, dose, tempo, droga, além dos dados dos mecanismos de ação envolvidos, se algum foi ativado e total ativo.
# Verificando se no grupo controle existe algum mecanismo de ação ativado
df_merged.query('tratamento == "com_controle"' )['ativo_moa'].value_counts()
Como o total de registros no grupo controle é de 1866, nenhum mecanismo de ação foi ativado no grupo controle, o que é ideal.
6. Como diferentes drogas modulam a expressão de diferentes genes e a viabilidade celular?¶
# Criando uma variável para receber o nome das top 5 drogas mais utilizadas, excluindo o controle
cod_drogas = df_merged['droga'].value_counts().index[1:6]
# Plotando alguns boxplots das top 5 drogas utilizadas para a lista de genes criada anteriormente e comparado com o grupo controle
for gene in g_list:
fig, axs = plt.subplots(1,2,figsize=(15,8))
sns.set_style("whitegrid")
sns.boxplot(data = df_merged.query('droga == "cacb2b860"'), y= gene, x='droga', hue='ativo_moa', linewidth=0.7, palette='magma', ax=axs[0])
axs[0].set_title('Boxplot controle')
sns.boxplot(data = df_merged.query('droga in @cod_drogas'), y= gene, x='droga', hue='ativo_moa', linewidth=0.7, palette='mako', ax=axs[1])
axs[1].set_title('Boxplot droga ' + droga )
plt.show()
Todas as drogas selecionadas ativaram algum mecanismo de ação (lado direito). Com essa série de gráficos conseguimos ver como diferentes drogas modulam a expressão de diferentes genes.
# Plotando alguns boxplots das top 5 drogas utilizadas para a lista de viabilidade celular criada anteriormente e comparado
# com o grupo controle
for cells in c_list:
fig, axs = plt.subplots(1,2,figsize=(15,8))
sns.set_style("whitegrid")
sns.boxplot(data = df_merged.query('droga == "cacb2b860"'), y= cells, x='droga', hue='ativo_moa', linewidth=0.7, palette='rocket_r', ax=axs[0])
axs[0].set_title('Boxplot controle')
sns.boxplot(data = df_merged.query('droga in @cod_drogas'), y= cells, x='droga', hue='ativo_moa', linewidth=0.7, palette='viridis', ax=axs[1])
axs[1].set_title('Boxplot droga ' + droga )
plt.show()
Todas as drogas selecionadas ativaram algum mecanismo de ação (lado direito). Com essa série de gráficos conseguimos ver como diferentes drogas alteram a viabilidade celular.
7. Modelagem preditiva - Criando um modelo de aprendizagem Não-Supervisionada¶
Neste projeto vamos utilizar o algoritmo KMeans para agrupar nosso conjunto de dados e compreender como eles estão organizados. Em um trabalho publicado na revista Nature por Patrik D'haeseleer (referência do artigo no Readme do repositório no GitHub), o pesquisador diz que método de agrupamento utilizando o KMeans é bem popular. Esse algoritmo de clusterização utiliza um método de particionamento, que subdivide todo o conjunto de dados, os genes no nosso caso, em um número predeterminado (k) de agrupamentos.
O algoritmo é inicializado com k centróides do cluster escolhidos aleatoriamente, e cada gene é atribuído ao cluster com o centróide do armário. Em seguida, os centróides são redefinidos para a média dos genes em cada cluster. Este processo continua até que nenhum outro gene mude de cluster. Posições de centróide iniciais diferentes podem produzir resultados de cluster diferentes, e é importante executar o algoritmo várias vezes com sementes aleatórias diferentes.
7.1 Avaliando a Tendência de Cluster¶
Um teste estatístico que permite adivinhar se os dados seguem uma distribuição uniforme. Se o teste for positivo (uma pontuação de hopkins que tende a 0), significa que os dados não estão uniformemente distribuídos. Portanto, o agrupamento pode ser útil para classificar as observações. No entanto, se a pontuação for muito alta (acima de 0,5 por exemplo); os dados são distribuídos uniformemente e o clustering pode não ser realmente útil para o problema em questão.
| Estatística Hopkins para o conjunto de dados | Descrição |
|---|---|
| Valores > 0.5 | significam que o dataset não é "clusterizável" |
| Valores < 0.5 | significam que o dataset é "clusterizável" |
| Valor ≅ 0 | Quanto mais próximo de zero melhor. |
# Vamos transformar a variável dose em numérica, substituindo D1 e D2 por 0 e 1, pois ela pode ser relevante para a clusterização
df_merged_clustering = df_merged.replace(["D1", "D2"], [0,1]); df_merged_clustering
# Selecionando apenas variáveis numéricas
df_numeric = df_merged_clustering.select_dtypes(exclude=['object'])
# Aplicando o teste hopkins
X = scale(df_numeric)
hopkins(X, 1000) # número de amostras
Valor = 0.148 indica que o dataset é fortemente "clusterizável".
7.2 Redução de dimensionalidade¶
# Aplica redução de dimensionalidade
pca = PCA(n_components = 2).fit_transform(df_numeric)
Ao invés de trabalhar com mais de 800 variáveis vamos utilizar apenas 2 componentes principais. A matemática por trás do Principal Component Analysis utiliza auto valores e auto vetores para compressão dos dados, mudando o conjunto de dados de um espaço de características para outro. Com essa mudança nós conseguimos representar a mesma informação, os números podem mudar mas a informação continua a mesma.
7.3 Determinando valor ideial de k (número ideial de clusters)¶
# Determinando um range de K
k_range = range(1,15)
# Aplicando o modelo K-Means para cada valor de K (esta célula pode levar bastante tempo para ser executada)
k_means_var = [KMeans(n_clusters = k).fit(pca) for k in k_range]
# Ajustando o centróide do cluster para cada modelo criado anteriormente
centroids = [X.cluster_centers_ for X in k_means_var]
Agora vamos executar todos os passos matemáticos para calcular a distância euclidiana entre os pontos e depois plotar a curva de elbow.
Links úteis:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html
# Calculando a distância euclidiana de cada ponto de dado para o centróide
k_euclid = [cdist(pca, cent, 'euclidean') for cent in centroids]
dist = [np.min(ke, axis = 1) for ke in k_euclid]
# Soma dos quadrados das distâncias dentro do cluster
soma_quadrados_intra_cluster = [sum(d**2) for d in dist]
# Soma total dos quadrados
soma_total = sum(pdist(pca)**2)/pca.shape[0]
# Soma dos quadrados entre clusters
soma_quadrados_inter_cluster = soma_total - soma_quadrados_intra_cluster
# Curva de Elbow
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(k_range, soma_quadrados_inter_cluster/soma_total * 100, 'b*-')
ax.set_ylim((0,100))
plt.grid(True)
plt.xlabel('Número de Clusters')
plt.ylabel('Percentual de Variância Explicada')
plt.title('Variância Explicada x Valor de K')
Olhando a curva de elbow podemos observar um platô a partir do valor de k = 8, a partir de 8 cluster conseguimos explicar a variância nos dados. O valor ideal de k que utilizeramos no Kmeans e algo acima desse valor. Portanto vamos clusterizar nosso dataset em 8 e 12 grupos e depois comparar os resultados. Nem sempre um maior número de cluster pode melhorar a acurácia do modelo pois isso pode levar a sobreposição dos dados.
7.4 Criando o modelo utilizando o algoritmo KMeans¶
# Criando um modelo com K = 8
modelo_v1 = KMeans(n_clusters = 8)
modelo_v1.fit(pca)
# Obtém os valores mínimos e máximos e organiza o shape
x_min, x_max = pca[:, 0].min() - 5, pca[:, 0].max() - 1
y_min, y_max = pca[:, 1].min() + 1, pca[:, 1].max() + 5
# obs: Esses subtrações e adições são necessárias para encontrar os valores ideais
xx, yy = np.meshgrid(np.arange(x_min, x_max, .02), np.arange(y_min, y_max, .02))
Z = modelo_v1.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
# Plot das áreas dos clusters
plt.figure(1)
plt.clf()
plt.imshow(Z,
interpolation = 'nearest',
extent = (xx.min(), xx.max(), yy.min(), yy.max()),
cmap = plt.cm.Paired,
aspect = 'auto',
origin = 'lower')
# Plot dos centróides
plt.plot(pca[:, 0], pca[:, 1], 'k.', markersize = 4)
centroids = modelo_v1.cluster_centers_
inert = modelo_v1.inertia_
plt.scatter(centroids[:, 0], centroids[:, 1], marker = 'x', s = 169, linewidths = 3, color = 'r', zorder = 8)
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.xticks(())
plt.yticks(())
plt.show()
Vamo utilizar uma métrica de clusterização que calcula o coeficiente de silhueta médio de todas as amostras. O melhor valor é 1 e o pior valor é -1. Valores próximos a 0 indicam clusters sobrepostos. Valores negativos geralmente indicam que uma amostra foi atribuída ao cluster errado, pois um cluster diferente é mais semelhante.
# Silhouette Score
labels = modelo_v1.labels_
silhouette_score(pca, labels, metric = 'euclidean')
# Criando um modelo com K = 12
modelo_v2 = KMeans(n_clusters = 12)
modelo_v2.fit(pca)
# Obtém os valores mínimos e máximos e organiza o shape
x_min, x_max = pca[:, 0].min() - 5, pca[:, 0].max() - 1
y_min, y_max = pca[:, 1].min() + 1, pca[:, 1].max() + 5
xx, yy = np.meshgrid(np.arange(x_min, x_max, .02), np.arange(y_min, y_max, .02))
Z = modelo_v2.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
# Plot das áreas dos clusters
plt.figure(1)
plt.clf()
plt.imshow(Z,
interpolation = 'nearest',
extent = (xx.min(), xx.max(), yy.min(), yy.max()),
cmap = plt.cm.Paired,
aspect = 'auto',
origin = 'lower')
# Plot dos centróides
plt.plot(pca[:, 0], pca[:, 1], 'k.', markersize = 4)
centroids = modelo_v2.cluster_centers_
inert = modelo_v2.inertia_
plt.scatter(centroids[:, 0], centroids[:, 1], marker = 'x', s = 169, linewidths = 3, color = 'r', zorder = 8)
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.xticks(())
plt.yticks(())
plt.show()
# Silhouette Score
labels = modelo_v2.labels_
silhouette_score(pca, labels, metric = 'euclidean')
A segunda versão do modelo apresentou um score inferior ao primeiro 0.72 vs 0.77. O valor ideal de cluster neste caso é 8.
Concluindo este projeto, vamos criar um clustermap com a média do número de mecanismos de ação ativos por cluster.
# Lista com nomes das colunas
nomes_colunas = ['tratamento', 'tempo', 'dose', 'droga', 'n_moa', 'ativo_moa']
# Cria o cluster map
cluster_map = pd.DataFrame(df_merged, columns = nomes_colunas)
cluster_map['cluster'] = modelo_v1.labels_
display(cluster_map)
# Calcula a média de consumo de energia por cluster
n_moa_mean = cluster_map.groupby('cluster')['n_moa'].mean().sort_values(ascending=False); display(n_moa_mean)
O cluster 2 apresentou a maior média (1,78) de mecanismos de ação ativos.
8. Conclusões¶
Não há valores duplicados.
Não há valores ausentes.
A droga ‘cacb2b860’ aparece com mais frequência do que as outras e representa que a amostra não recebeu tratamento nenhum.
Apenas 9 drogas aparecem com muito mais frequência do que as demais dentro da base de dados.
Os dados nas variáveis ‘dose’ e ‘tempo’ estão bem distribuídos entre os grupos.
O desenho experimental foi feito para utilizar duas doses em três janelas distintas de tempo (24h, 48 e 72h) para verificar a modulação da expressão gênica em mais de 700 genes e na viabilidade celular por diferentes drogas com diferentes mecanismos de ação.
O gene 'g-707’ é o que apresenta a maior média de expressão, seguidos dos genes 'g-100', 'g-744' e 'g-392'. Enquanto os genes 'g-370', 'g-508' e 'g-37' apresentaram a menor média de expressão.
A distribuição dos dados para a média das expressões gênicas se assemelha muito a uma curva normal, tendo assimetria levemente negativa e kurtosis levemente positiva.
As variáveis que representam a viabilidade celular que apresentaram a maior média foram as 'c-74', 'c-37', 'c-7', 'c-76' e 'c-69'. Enquanto as menores médias foram as variáveis 'c-65', 'c-18', 'c-63', 'c-38' e 'c-26' respectivamente. É interessante observar que todas as médias foram negativas, o que pode indicar baixa viabilidade celular para todas as variáveis do conjunto de dados.
A distribuição dos dados para a média de viabilidade celular não se assemelha a uma curva normal, tendo assimetria positiva e kurtosis fortemente positiva. Os dados estão muito mais concentrados em torno da média da distribuição (-0,432231) com um desvio padrão baixo (0,093).
Podemos inferir que alguns genes se correlacionam tanto positivamente como negativamente com a viabilidade celular. Em outras palavras, existem situações em que a expressão gênica aumenta e a viabilidade celular também (correlação positiva) e casos em que a expressão gênica diminui/aumenta e a viabilidade celular aumenta/diminui respectivamente (correlação negativa).
Analisando as features selecionadas para fins didáticos, podemos afirmar que os genes 'g-38' e 'g-672' estão correlacionados negativamente, o que significa que quando a expressão gênica de um aumenta a do outro diminui. Quando a expressão do 'g-38' aumenta a viabilidade (c-26) tende a diminuir. Quando a expressão de 'g-672' aumenta a viabilidade celular também tende a aumentar. Esses resultados podem indicar que esses genes possam estar intrinsicamente conectados por alguma via que foi modulada pela droga em questão. O tempo de exposição é relevante para os efeitos observados superexpressão gênica de g-38 e subexpressão do g-672 foram mais intensas no tempo de 72h.
Podemos inferir que dentre as drogas analisadas há uma diferença na expressão gênica de 'g-38' quando comparados os grupos tratado e controle. Testes estatísticos específicos poderiam confirmar essa hipótese, que por sinal é bem forte.
Existem mais de 200 mecanismos de ação no segundo conjunto de dados.
Os tipos de mecanismo de ação mais frequentes foram ‘nfkb_inhibitor’,‘proteasome-inhibitor’, ‘cyclooxygenase_inhibitor’,‘dopamine_receptor_antagonist’, e ‘serotonin_receptor_antagonist’.
Nenhum mecanismo de ação foi ativado no grupo controle, o que é esperado e ideal.
O conjunto de dados é altamente “clusterizável” obtendo um resultado de 0,14 no teste de Hopkins.
O número ideal de cluster para se trabalhar com essa base de dados foi de 8 grupos, com um score de silhueta de 0,77.
O grupo que apresentou a maior média de mecanismos de ação ativos foi o grupo 2.
Curiosidade: Este notebook levou 05:09,87ms para ser executado.
Apesar de conseguido fazer algumas análises interessantes, o curto prazo para a execução deste projeto teve forte influência dentro do que foi possível ser feito. Concentrei meus esforços no que foi mais relevante dentro do conjunto de dados, como as drogas mais utilizadas, os genes que tiveram a média de expressão mais alterada e as variáveis de viabilidade celular que mais divergiram em relação ao grupo controle. Todas as análises foram feitas utilizando apenas essas 20 variáveis selecionadas de acordo com esse critério. O ideal seria trabalhar com todos os genes, o que levaria muito mais tempo, entretanto, creio que o resultado ficou excepcional dentro do prazo de 3 dias proposto para executar este projeto.
Se você chegou até aqui, meus sinceros agradecimentos, sei que o notebook ficou imenso e é preciso muita resiliência para ver com calma os resultados e chegar em alguma conclusão. Esse desafio foi uma oportunidade incrível de testar meus conhecimentos e de trabalhar dentro de uma área incrível que é a biologia. Por fim, gostaria de parabenizar todos os membros da equipe da Alura pelo excelente treinamento e suporte durante toda a 3 Imersão dados.
Versionamento¶
| Versão | Editor | Data | Observação |
|---|---|---|---|
| 1.0 | Eduardo Gonçalves | 09/05/2021 | Primeira versão do projeto. |