Bioinfo delicinha II
Salve!
Há alguns meses, eu compartilhei aqui uma introdução sobre Bioinformática, contando como podemos usar R para carregar uma tabela com todos os genes encontrados em seres humanos, e um link para o Tabula Muris Senis (TMS), um atlas com os níveis de cada gene em camundongos com difrentes idades entre 1 e 27 meses.
Hoje, nós vamos mergulhar no atlas TMS e investigar genes associados ao envelhecimento dos nossos queridos roedores. Especificamente, focaremos em encontrar genes expressos no cérebro que estão associados ao envelhecimento.
Para começar, precisamos baixar os dados disponíveis no seguinte link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132040
Ele leva ao Gene Expression Omnibus, uma base de dados em que cientistas publicam os resultados de seus experimentos. Os dados em que estamos interessados estão na parte inferior da tela, em uma tabela chamada "Supplementary file". O arquivo GSE132040_190214_A00111_0269_AHH3J3DSXX_190214_A00111_0270_BHHMFWDSXX.csv.gz contém os dados de expressão dos genes em si - e de todos os tecidos dos camundongos. Já o arquivo GSE132040_MACA_Bulk_metadata.csv.gz contém uma tabela com os metadados do experimento, e nos permite encontrar várias informações sobre cada amostra, como o tecido de origem e a idade do indivíduo.
Agora que temos nossos dados, a diversão começa! No R, podemos importar as duas tabelas. Felizmente, elas pode ser importadas mesmo estando comprimidas, usando a função read.csv.
dataFolder = 'coloqueAquiOCaminhoDoDownloadDoGeo'
geneExpression = read.csv(paste0(dataFolder, 'GSE132040_190214_A00111_0269_AHH3J3DSXX_190214_A00111_0270_BHHMFWDSXX.csv.gz'))
geneMetadataAll = read.csv(paste0(dataFolder, 'GSE132040_MACA_Bulk_metadata.csv.gz'))
Como a tabela de metadados tem muitas informações (primeira linha no bloco abaixo - os comentários iniciando em '#' logo abaixo mostram a saída do comando), vamos selecionar somente nossas colunas de interesse: tecido da amostra e idade no momento da coleta.
colnames(geneMetadataAll)
# [1] "Sample.name" "title" "source.name"
# [4] "organism" "characteristics..age" "characteristics..developmental.stage"
# [7] "characteristics..sex" "molecule" "description"
# [10] "processed.data.file" "raw.file" "BioSample"
# [13] "Instrument.Model"
geneMetadata = geneMetadataAll[,c('Sample.name', 'source.name', 'characteristics..age')]
head(geneMetadata)
# Sample.name source.name characteristics..age
# 1 A1_384Bulk_Plate1_S1 BAT_24 6
# 2 A1_384Bulk_Plate3_S1 SCAT_43 3
# 3 A10_384Bulk_Plate1_S10 Brain_47 18
# 4 A10_384Bulk_Plate2_S10 GAT_39 1
# 5 A10_384Bulk_Plate3_S10 Lung_6 12
# 6 A11_384Bulk_Plate1_S11 Limb_Muscle_11 3
Agora, vamos extrair o nome das amostras coletadas do sistema nervoso central da tabela de metadados e selecionar somente estas amostras na tabela com os níveis de cada gene. É interessante notar que, nesta segunda tabela, cada amostra corresponde a uma coluna; enquanto, na primeira, cada amostra corresponde a uma linha.
samplesSourceBrain = grep('Brain', geneMetadata$source.name)
geneMetadataBrain = geneMetadata[samplesSourceBrain,]
sampleNameBrain = paste0(geneMetadataBrain$Sample.name, '.gencode.vM19')
geneExpressionBrain = geneExpression[,sampleNameBrain]
nrow(geneExpressionBrain)
# [1] 54357
ncol(geneExpressionBrain)
# [1] 56
No final, temos aproximadamente 54 mil genes, medidos em 56 amostras diferentes. Vamos inspecionar quantas amostras nós temos em relação a cada idade. O commando "table" nos mostra quantas vezes cada elemento se repete. A primeira linha da saída mostra a idade em meses de cada indivíduo, ordenada em ordem lexicográfica. Cada grupo tem seis amostras, exceto pelos camundongos com 24 e 27 meses, que contém quatro indivíduos por grupo.
table(geneMetadataBrain$characteristics..age)
# 1 12 15 18 21 24 27 3 6 9
# 6 6 6 6 6 4 4 6 6 6
Antes de partirmos para a ação, é sempre saudável ver como nossas amostras são comparáveis e normalizá-las, se preciso. Normalização é um assunto que pode ir looooonge, e a técnica que eu ilustro aqui é bem rudimentar. Como experimentos de genômica podem ter eficiências diferentes, é possível observar que algumas amostras possuem uma discrepâncias nos níveis totais de expressão - por exemplo, uma amostra contém mais de 55 milhões de fragmentos de RNA, enquanto a amostra com menor eficiência apresenta apenas 3 milhões. Neste caso, um gene com 500 fragmentos na primeira amostra tem, aproximadamente, a mesma proporção em relação ao total se ele apresentar 30 fragmentos na segunda. Para tornar as duas amostras comparáveis (apesar da graaande discrepância), vamos multiplicar os níveis de expressão em cada amostra de forma que elas tenham 26 milhões de fragmentos, o valor médio observado (novamente, esta prática é bem simples, e métodos mais eficientes são comumente utilizados).
summary(colSums(geneExpressionBrain))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 3288190 21984764 26572072 26338253 31085600 55313948
sampleNormalizationFactors = mean(colSums(geneExpressionBrain))/colSums(geneExpressionBrain)
normalizedGeneExpressionBrain = sweep(geneExpressionBrain, 2, sampleNormalizationFactors, '*')
summary(colSums(normalizedGeneExpressionBrain))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 26338253 26338253 26338253 26338253 26338253 26338253
Agora é hora de ver quais genes estão associados ao envelhecimento! Existem milhares de formas de definir esta relação. Neste post, vamos ver genes que apresentam uma correlação (positiva ou negativa) com a idade. Como a correlação não faz sentido para variáveis que não mudam, nós vamos analisar somente genes com desvio padrão (standard deviation) diferente de zero. Ainda assim, nós temos uns 50 mil genes: o próximo bloco pode demorar alguns minutos para executar.
sdGeneBrain = apply(normalizedGeneExpressionBrain, 1, sd)
sum(sdGeneBrain == 0)
# [1] 2408
correlationTestBrain = apply(normalizedGeneExpressionBrain[sdGeneBrain != 0,], 1, function(geneLevels){
cor.test(geneLevels, as.numeric(geneMetadataBrain$characteristics..age))
})
Extraindo os resultados (campo estimate do resultado da correlação) e fazendo um histograma da distribuição dos valores de correlação de atividade genética com a idade.
correlationEstimateBrain = sapply(correlationTestBrain, function(ctb) ctb$estimate)
summary(correlationEstimateBrain)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# -0.6773914 -0.1056527 -0.0014906 -0.0008069 0.1017325 0.8564320
hist(correlationEstimateBrain)
Podemos notar que a maioria dos genes tem uma correlação baixa com a idade, próxima de zero. Entretanto, também identificamos genes que tanto diminuem com a idade, quanto genes que se tornam mais ativos. Vamos identificar quantos genes efetivamente apresentam uma correlação entre tempo de vida ao observar um parâmetro essencial em bioinformática: o valor p. Este valor vem da estatística, e serve como um indicador da probabilidade de observar a variância obtida experimentalmente caso nossas amostras venham da mesma fonte. Valores pequenos de p indicam que nossas amostras de fato têm efeito sobre a variável de interesse, e podemos dizer que ela responde às condições do experimento. Aqui, vamos adotar a convenção de que valores significativos são indicados por p menor que 5 %.
correlationPValueBrain = sapply(correlationTestBrain, function(ctb) ctb$p.value)
hist(correlationPValueBrain)
sum(correlationPValueBrain < 0.05)
# [1] 4176
Aproximadamente 4000 mil genes do cérebro têm uma correlação singificante com a idade. Para destacar tais genes, vamos ordenar as correlações baseadas em seu valor absoluto (em outras palavras, correlações negativas, mas com grande amplitude, vão aparecer no topo também). Nós também temos de voltar à tabela original com níveis de expressão, onde os nomes de cada gene estão armazenados na primeira coluna. Ao final do próximo bloco, nós construímos nossa própria tabela ligando os nomes dos genes às suas respectivas correlações.
sortedCorrelationEstimateBrain = sort.int(abs(correlationEstimateBrain), index.return = TRUE, decreasing = TRUE)
geneNamesBrain = geneExpression[, 'gene']
geneNamesSdBrain = geneNamesBrain[sdGeneBrain != 0]
sortedGenesWithNames = data.frame(geneName = geneNamesSdBrain[sortedCorrelationEstimateBrain$ix],
ageCorrelation = correlationEstimateBrain[sortedCorrelationEstimateBrain$ix])
Na lista abaixo, podemos inspecionar os 10 genes com maior correlação. Os genes C4b e C3 estão associados a respostas do sistema imune - que de fato é mais ativo em indivíduos com mais idade. Já o gene Neat1 não codifica uma proteína, mas está ligado à regulação da atividade de outros genes e está associado ao câncer.
sortedGenesWithNames[seq(10),]
# geneName ageCorrelation
# 5531.cor C4b 0.8564320
# 5509.cor C3 0.7875731
# 42574.cor Neat1 0.7754526
# 45044.cor Pcdhb8 0.7501735
# 39548.cor Lgals3bp 0.7375911
# 39944.cor Lyz2 0.7339193
# 39547.cor Lgals3 0.7325344
# 53518.cor Zc3hav1 0.7291223
# 49582.cor Spag6 0.7288461
# 38719.cor Itgb2 0.7247033
Por hoje, é só tudo isso! Pra vc que chegou até aqui, muito obrigado pela sua atenção e espero que esta pequena caminhada possa ter despertado seu interesse pelos processos que ocorrem dentro de nós enquanto vc lê este texto - e até quando vc pensa sobre eles! Pra quem quiser começar a flexionar seus músculos de bioinformata, recomendo tentar o mesmo procedimento descrito aqui, para um outro órgão. Quem sabe, até organizar esse código em funções e executar para todos os órgão coletados - e comparar as listas obtidas em cada caso para uma descrição sistêmica do envelhecimento :)