+ Page Not Found
+ +Try searching the whole site for the content you want:
+ + +diff --git a/preview/pr-27/404.html b/preview/pr-27/404.html new file mode 100644 index 0000000..2fdfc17 --- /dev/null +++ b/preview/pr-27/404.html @@ -0,0 +1,515 @@ + + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Try searching the whole site for the content you want:
+ + +Our lab is part of the Federal University of Rio Grande do Norte’s Bioinformatics Multidisciplinary Environment. Our team is headquartered on the 2nd Floor of the Brain Institute (ICe).
+ + + + + + +The Dalmolin Systems Biology Group advances biomedical knowledge through Bioinformatics and Systems Biology. Our mission is to unravel the intricate complexities of biological systems through innovative approaches and advanced methodologies. Committed to excellence, our team is dedicated to addressing a diverse array of biological challenges using the powerful lens of Systems Biology. Our research group is led by professor Rodrigo Dalmolin, Ph.D., and it is associated with the Graduate Program in Bioinformatics from the Federal University of Rio Grande do Norte (UFRN).
+ + +Our Research
+ + +Our research focuses on the evolution of biological systems, analysis of biological networks, transcriptional analysis, metagenomic analyses, and the development of novel bioinformatics tools.
+ + + + +Tools by researchers, for researchers
+ + +Streamline your investigation with our tools and pipelines. All open source and well documented. Don’t know how to use them? Get in touch with us! See the contact information.
+ + + + +Building excellence in bioinformatics
+ + +We are a team of interdisciplinary researchers working in the various fields of bioinformatics. We strive for excellence in research and diversity and inclusivity are our values.
+ + + + +Some Markdown content
+Metagenomics data, or data containing the full genetic content of an environmental sample, constitutes one of our current research topics at the Dalmolin Group. Metagenomics data provides a great potential for studying the diversity and functional activity of complex microbial environments.
+ +Our studies of metagenomics data began with the development of our in-house pipeline for processing these data, MEDUSA (Morais et al, 2022). With MEDUSA, we aimed to build a sensitive, flexible and reproducible pipeline that produces taxonomic classification and functional annotation of metagenomics reads. Below, you can see a diagram detailing the different steps of the original MEDUSA pipeline. We later reimplemented MEDUSA in another workflow orchestration language, Nextflow, this time being called EURYALE.
+ + + +Some of our work exploring metagenomics data is highlighted by “Metagenomic Analyses Reveal the Influence of Depth Layers on Marine Biodiversity on Tropical and Subtropical Regions” (Santiago et al, 2023). In this work, using MEDUSA, we observed differences in diversity between three ocean layers in the Tropical and Subtropical regions and also noted different functional annotations between these layers.
+ +And now, with MEDUSA/EURYALE and all of the software in the metagenomics ecosystem, we aim to keep exploring metagenomics data, with a focus on functional analysis.
+ + +The development of sequencing technologies in the past two decades allowed the study of the phylogenetic relationship of several model organisms. An interesting question relies on the mechanisms behind how new genes arise and how they interact and evolve together to build the intricate network of genes observed in the organism’s biochemical pathways. Gene duplication is highlighted as a crucial event in genome evolution, offering a prime source for genetic material and allowing evolutionary forces to generate novelty. In the Dalmolin et al, 2011 paper, the importance of duplication events in the evolutionary history of orthologous groups is assessed by the ratio of components to the number of organisms with items from the group. The study extends to analyze Eukaryotic Clusters of Orthologous Groups (KOG) in STRING, evaluating the evolutionary plasticity and conservation of groups based on component distribution across eukaryotic genomes, proposing an equation for plasticity, and identifying correlations between evolutionary distance and plasticity. The algorithm used in this paper to infer the evolutionary root of orthologous groups, based on a phylogenetic tree from eukaryotes, is implemented as an R/Bioconductor package called geneplast.
+ + + +In the Viscardi et al. 2020 paper, we aimed to reconstruct the evolutionary history of the human neurotransmission gene network by analyzing genes in major neurotransmitter systems. The findings suggest that the emergence of receptor families, particularly ionotropic receptors, at the Human–Cnidaria Last Common Ancestor played a crucial role in the evolution of synapses, even before their establishment in Ctenophores.
+ + + +In De Souza et al, 2021 paper, essential genes from five model eukaryotes were analyzed to determine if they differ in age from non-essential genes. The study compared network properties and investigated the biological functions of essential genes in each species, revealing that essential genes are rooted in different evolutionary periods. Despite varying origins, essential genes tend to be older and more connected than other genes across all evaluated model species.
+ + +Some Markdown content
+Some Markdown content
+Graduated in Information Technology (Bachelor) with minor certification in Bioinformatics (2019) and Master in Bioinformatics from PPg-Bioinfo (2022), both from the Federal University of Rio Grande do Norte (UFRN). She is a Computer Technician at the Integrated Technical Course in Computer Science at the Federal Institute of Education, Science and Technology of Rio Grande do Norte (IFRN)-2015. Recently she contributed to the Fiocruz Genomic Network working on the development of pipelines for the analysis of viromas. She is currently a PhD student in Bioinformatics at PPg-Bioinfo at UFRN, where she is part of the Bioinformatics Multidisciplinary Environment (BioME-IMD, UFRN) working at the Dalmolin Systems Biology Group (DSBG) laboratory in the areas of Systems Biology and Metagenomics.
+ + + + + + ++ + Search for Bianca Santiago's papers on the Research page + +
+ + + + +Pursuing a master’s degree in the Postgraduate Program in Bioinformatics at IMD-UFRN. Bachelor in Biological Sciences from the Federal University of Rio Grande do Norte (2021). During the undergraduate period, I conducted research in the field of Physiology, with an emphasis on Animal Behavior in fish (Danio rerio). I hold a technical degree in Computer Science from the Federal Institute of Education, Science and Technology of Rio Grande do Norte (2015).
+ + + + + + ++ + Search for Bruno Fernandes's papers on the Research page + +
+ + + + +Undergraduate at Infantaria from Academia Militar das Agulhas Negras (1985), Master’s at Programa de Pós-Graduação em Informática - PPGI from Universidade Federal da Paraíba (2015) and Doctorate in Bioinformatics from Universidade Federal do Rio Grande do Norte (2019). Has experience in Computer Science, acting on the following subjects: impulsive detection methods, partial discharge, synonymous mutation, beijerinckia and molecular evolution.
+ + + + + + ++ + Search for Clóvis Reis's papers on the Research page + +
+ + + + +I am a Master’s student in Bioinformatics and BSc in Biomedicine at UFRN, where I carry out research in systems biology and deal with data integration, biological networks analysis, the evolution of biological systems and the development of computational tools. I also work as a bioinformatician at the genomics lab Mendelics, developing analysis pipelines for genomic data as well as developing LIMS for multiple sequencing platforms. I am an IT technician experienced in R, data visualization, full-stack development and databases.
+ + + + + + ++ + Search for Danilo Imparato's papers on the Research page + +
+ + + + +I am a computer science technician from the Federal Institute of Brasília (IFB). Currently, I am pursuing a bachelor’s degree in Biological Sciences at the Federal University of Rio Grande do Norte (UFRN). I am also a participant in the Institutional Program for Scientific Initiation Scholarships (PIBIC) funded by CNPq, working in the bioinformatics laboratory at the Institute Metrópole Digital (IMD).
+ + + + + + ++ + Search for Dante von Zuben's papers on the Research page + +
+ + + + +I have a degree in Computer Science from the State University of Rio Grande do Norte (2011), a specialization in Computer Forensics from the Universidade Potiguar (2016), and a master’s and a doctorate in Bioinformatics from the Federal University of Rio Grande do Norte (2018 and 2022). My experience lies in the fields of Computer Science and Bioinformatics, with a focus on Computer Networks and Systems Biology.
+ + + + + + ++ + Search for Diego Morais's papers on the Research page + +
+ + + + +B.Sc. in Science and Technology (2021) and in Biomedical Engineering (2022), M.Sc, in Bioinformatics (2023), and is currently a Ph.D. student in Bioinformatics, all at the Federal University of Rio Grande do Norte (UFRN). His training in Bioinformatics focused on Systems Biology, working with transcriptomics and genomics analyses. He has experience in the analysis of biological networks through Systems Biology methodologies associated with machine learning techniques. He is currently developing a research project at the Multiuser Center of Bioinformatics (BioME - IMD/UFRN) in the areas of oncology, non-coding RNAs and machine learning.
+ + + + + + ++ + Search for Epitácio Farias's papers on the Research page + +
+ + + + +I hold a degree in Information Systems from UNI-RN (2007) and a degree in Biomedicine from UFRN (2021). I worked at Petrobras S.A. in the IT and TCOM (Telecommunications) fields. In the development areas, I utilized the Pascal/Delphi programming language. I managed network infrastructure based on TCP/IP protocols with support for Linux servers and applications. I also managed TCOM infrastructure with support for routing protocols such as OSPF, RIP, and BGP, involving both local and wireless network protocols.
+ + + + + + ++ + Search for Fábio Ferreira Silva's papers on the Research page + +
+ + + + +Graduating in Biomedicine from the Federal University of Rio Grande do Norte since 2018.1, expected to be completed in 2024.2. My focus of interest is in the area of Bioinformatics, with an emphasis on RNA-seq data analysis and Metagenomics. I am currently part of the BioME team, where I work as a Scientific Initiation Fellow. My participation in BioME involves the application of Bioinformatics tools to analyze different types of data and biological processes. At the same time, I work in the area of Microbiology as a monitor in the Medical Bacteriology discipline. In this context, I collaborate in the development of the bacteria identification manual aimed at academic teaching, in addition to providing assistance in practical classes. This initiative has provided improvements in bacterial identification techniques and the opportunity to contribute to the teaching of microbiology. Previous experiences include monitoring the subjects of Biochemistry and organic chemistry, offering support to students and helping them learn fundamental concepts. Furthermore, I participated in another Scientific Initiation project also focused on the area of Bioinformatics, in which I discovered my interest in the area. My goal is to pursue a research career in Bioinformatics, I am open to research and collaboration opportunities in order to expand its applications in the biomedical field and contribute to scientific advances.
+ + + + + + ++ + Search for Gleison Medeiros's papers on the Research page + +
+ + + + +I have a bachelor’s degree in Biotechnology - Bioinformatics from the Federal University of Rio Grande do Sul (2021) and a Master’s degree in Bioinformatics at the Federal University of Rio Grande do Norte. Currently, I am a PhD student in Bioinformatics at the Federal University of Rio Grande do Norte.
+ + + + + + ++ + Search for Gustavo Michaelsen's papers on the Research page + +
+ + + + +B.Sc. in Biomedical Sciences (2016) and Statistics (2023), M.Sc. in Bioinformatics (2017), and Ph.D. in Bioinformatics (2022), all from the Federal University of Rio Grande do Norte (UFRN). His training in bioinformatics focused on Systems Biology, working with transcriptomics and genomics analyses. She has experience in the analysis of biological networks through systems biology methodologies and biochemical systems evolution. In addition, she is interested in data science and in the application of statistical methods in bioinformatics. She is the cofounder of the Lexanomics startup, a company focused on bioinformatics solutions. She is currently developing her research projects at the Bioinformatics Multidisciplinary Environment (BioME - IMD/UFRN) in the areas of toxicogenomics, neuropsychiatric, and molecular evolution.
+ + + + + + ++ + Search for Iara Dantas de Souza's papers on the Research page + +
+ + + + +Master in Bioinformatics at the Federal University of Rio Grande do Norte (UFRN) / Biome (2019), he develops research within the scope of Systems Biology. Specialized in Information Technology Applied to the Legal Area from UFRN/Tribunal de Justiça do Rio Grande do Norte (TJRN), graduated in Information Technology from UFRN (2018), Systems Analysis and Development from Fatec Carapicuíba (2012) and Business Administration Companies by Universidade Presbiteriana Mackenzie (2013).
+ + + + + + ++ + Search for Igor Brandão's papers on the Research page + +
+ + + + ++ + Search for João Cavalcante's papers on the Research page + +
+ + + + +Pharmacy student at the Federal University of Rio Grande do Norte (UFRN).
+ + + + + + ++ + Search for João Vitor Almeida's papers on the Research page + +
+ + + + +I have a degree in Computational Engineering, with an emphasis in Systems and Computing, and a Master’s in Computational and Electric Engineering. As a researcher, I have started designing intelligent systems using a reinforcement learning technique for routing packets in simulated networks. I have also been working as a developer since 2004 in the Federal University of Rio Grande do Norte - UFRN, obtaining vast experience in Java programming language, web/mobile systems, relational and non-relational databases, project management, and agile processes. I worked on many projects which yielded software patents to the university. I’m currently a Ph.D. student at Bioinformatics Multidisciplinary Environment’s Graduate Program from the Federal University of Rio Grande do Norte - UFRN. I have experience in the areas of Artificial Intelligence and Machine Learning, Data Science, Computational Biology, and Systems Biology.
+ + + + + + ++ + Search for Leonardo Campos's papers on the Research page + +
+ + + + +Marcel holds a Computer Engineering degree from the Federal University of Rio Grande do Norte, a graduate degree in Big Data from the Instituto Metrópole Digital at UFRN, another graduate degree in Health Informatics from UFRN, a Master of Science in Bioinformatics from the Graduate Program in Bioinformatics at IMD/UFRN and a doctorate degree from the Sorbonne Université in Paris. During his PhD, he worked as a researcher at the Institut Curie, part of the Université de recherche Paris Sciences et Lettres (PSL Research University), in the UMR168 unit in the laboratory “Evolution of Biomolecular Networks, RNA Dynamics”, under the supervision and doctoral guidance of Prof. Dr. Hervé Isambert.Co-founder of the Laboratory of Technological Innovation in Health (LAIS) at Hospital Universitário Onofre Lopes (HUOL-UFRN), he participated in research activities there for more than 7 years in the areas of health informatics, medical devices, hospital networks, real time, patient telemonitoring, teleradiology, health human resources systems and artificial intelligence. He also participated in research activities resulting from international cooperation, such as with Harvard University and MIT, as well as management, internationalization and fundraising activities for research.He currently works as a Developer Advocate at Seqera, a multinational leader in the development of technology for orchestrating data pipelines, and as a graduate professor at Universidade Potiguar. He is interested in the following topics: causal inference, systems biology, biological networks, bioinformatics, data engineering, data science, machine learning and artificial intelligence.
+ + + + + + ++ + Search for Marcel Ribeiro-Dantas's papers on the Research page + +
+ + + + +PhD in Biochemistry and Molecular Biology from the Federal University of Rio Grande do Norte (UFRN), Specialist in Molecular Biology-Applied to Human Health-UFRN, Master in Biochemistry from UFRN (IMT), Graduate in Biological Sciences, Universidade Potiguar, Laureate International Universities - UnP, worked with bioactive principles of plant compounds, under the guidance of Matheus Augusto Bittencourt Pasquali-UFCG and co-supervision of Rodrigo J. S. Dalmolin-UFRN/BIOME. Currently a businesswoman in the field of Clinical Analysis and health (NatLab).
+ + + + + + ++ + Search for Natália Souza's papers on the Research page + +
+ + + + +He has a degree in Mathematics from the Federal University of Rio Grande do Norte (2007), a Bachelor’s degree in Information Technology from the Federal University of Rio Grande do Norte (2020) and a master’s degree in Mathematics - ProfMat/IMPA from the Federal University of Rio Grande do Norte (2014). He has experience in the areas of Mathematics, Bioinformatics and Information Technology, with main interest in the following topics: R Language, Graphs, Systems Biology, Transcriptomics and Forensic Genetics. He is a master’s student in the Postgraduate Program in Bioinformatics at UFRN.
+ + + + + + ++ + Search for Odilon Santos's papers on the Research page + +
+ + + + +I graduated in Biological Sciences (Bachelor’s degree) from the Federal University of Maranhão (2018). I obtained a master’s degree from the Graduate Program in Oncology and Medical Sciences at the Federal University of Pará (2021). I have experience in Oncogenetics, Immunogenetics, and Molecular Biology, where I conduct research in these areas. I am a member of the following research groups: Human Population Genomics Research Network and the Autoimmune Diseases group.
+ + + + + + ++ + Search for Rafaella Ferraz's papers on the Research page + +
+ + + + +I hold a Bachelor’s degree in Biomedicine with a specialization in Clinical Analysis/Pathology from the Federal University of Rio Grande do Norte - UFRN (2014), where I conducted research with a focus on bioinformatics applied to nucleic acids. I used this approach to analyze economically significant viruses. I completed a Master’s degree in Biochemistry from the Graduate Program in Biochemistry at UFRN (2017), with a project also centered on bioinformatics, this time employing systems biology approaches to analyze the evolution of gene systems in Arabidopsis thaliana. I earned a Ph.D. in Bioinformatics from the Graduate Program in Bioinformatics at UFRN (2021), where I studied transcription factors as potential master regulators in the development of complex inflammatory diseases.
+ + + + + + ++ + Search for Raffael Oliveira's papers on the Research page + +
+ + + + +Prof. Dalmolin is a Ph.D. in Biochemistry (2012) by the Federal University of Rio Grande do Sul (UFRGS), with a master’s degree also in Biochemistry (2008) and a bachelor’s degree in Biological Sciences with an emphasis in Molecular, Cellular, and Functional Biology (2005) by the same institution. Currently, he is an Adjunct Professor and the head of the Biochemistry Department of the Federal University of Rio Grande do Norte (UFRN). He is a permanent faculty member of the graduate program in Biochemistry and the graduate program in Bioinformatics, developing his research activities in the Bioinformatics Multidisciplinary Environment at UFRN (BioME-IMD, UFRN). His research interests rely on the evaluation of biological networks, especially in the study of protein-protein interactions and regulatory networks, the evolution of biochemical systems, and the development of bioinformatics tools.
+ + + + + + ++ + Search for Rodrigo J. S. Dalmolin's papers on the Research page + +
+ + + + +Computer Engineer graduated from the Federal University of Paraíba, with a sandwich period at the Milwaukee School of Engineering (USA), funded by the Coordination for the Improvement of Higher Education Personnel (CAPES), through the Science Without Borders program. Master’s degree in Bioinformatics from the Federal University of Rio Grande do Norte.
+ + + + + + ++ + Search for Tayrone Monteiro's papers on the Research page + +
+ + + + +I am a PhD in Physics from the Federal University of Rio Grande do Norte. During my bachelor’s degree, I conducted research in the areas of Classical Field Theory and Population Dynamics. I established collaborations with researchers at the Faculty of Sciences of the University of Amsterdam, where I developed numerical simulations of population biodynamics. In my master’s degree, I delved into the field of Gravitation; I studied the 3+1 Formalism of General Relativity and Numerical Relativity, and I conducted simulations of black hole collisions. In my doctoral research, I worked on gravitational waves in collaboration with LIGO (Laser Interferometer Gravitational Wave Observatory) and carried out a project involving machine learning and the quality of numerical gravitational waves. I enjoy theoretical work in the areas of Cosmology and Gravitation, numerical work using C/C++, Python, and Julia, and Machine Learning for any interdisciplinary project in the field of Exact Sciences.
+ + + + + + ++ + Search for Tibério Azevedo Pereira's papers on the Research page + +
+ + + + +BSc in Biomedical Science at the Federal University of Rio Grande do Norte (UFRN), with a scientific initiation scholarship from the Bioinformatics Multi-User Center (BioME/IMD). Possesses skills in molecular biology and bioinformatics. Engages in work in the field of bioinformatics, with a special interest in neuropsychiatric diseases. Specifically studies sex-specific differences in major depressive disorder from a systems biology perspective. Proficient in RNA-seq analysis, Microarray, differential gene expression, splicing variants, and biological networks. Has a strong affinity for the field of data science, primarily using R and Bash, but also has experience with Python and JavaScript. Also interested in machine learning, development of computational tools, and mobile applications.
+ + + + + + ++ + Search for Vítor Saldanha's papers on the Research page + +
+ + + + +Below you can find the papers published by our group. Found something interesting, want to collaborate with us or can’t access a paper? Get in contact with us. In case you want to know more about our lines of research, check out the research page.
+ + + + + + + +Here you can see the lines of research in bioinformatics and systems biology that we work with. Clicking on a card below will take you to a brief explanation of the research we do within this field and a few examples.
+The Dalmolin Group routinely organizes courses in Bioinformatics as part of BioME. Below, you can find the teaching materials used in these courses, which encompass both the undergrad and the post-graduate levels. The courses are in Brazilian Portuguese.
+O interpretador da linguagem R pode ser instalado em Linux, Mac e +Windows1, +e encontra-se disponível gratuitamente no Comprehensive R Archive Network +(CRAN).
+O RStudio é um ambiente de desenvolvimento integrado que inclui +console, editor ciente de sintaxe e diversas outras ferramentas, que +visam o aumento da produtividade do desenvolvedor. Possui edições +gratuitas e comerciais, que podem ser obtidas em RStudio.com.
+# Adição
+2 + 5
## [1] 7
+# Subtração
+5 - 2
## [1] 3
+# Multiplicação
+2 * 5
## [1] 10
+# Divisão
+8 / 2
## [1] 4
+# Exponenciação
+2 ^ 5
## [1] 32
+2 ** 5
## [1] 32
+# Resto da divisão
+5 %% 2
## [1] 1
+# Igual
+3 == 5
## [1] FALSE
+# Diferente
+3 != 5
## [1] TRUE
+# Maior que
+3 > 5
## [1] FALSE
+# Menor que
+3 < 5
## [1] TRUE
+# Maior ou igual
+3 >= 5
## [1] FALSE
+# Menor ou igual
+3 <= 5
## [1] TRUE
+Operações podem ser concatenadas:
+((2 + 5 - 3) * 10) ^ 4 / 7 ^ 4
## [1] 1066.222
+Atribuição de valores:
+x <- 1
+# Sobrescreve o conteúdo anterior da variável x
+x <- 5
+y <- "gol do Grêmio"
Exibindo conteúdo de variáveis:
+x
## [1] 5
+y
## [1] "gol do Grêmio"
+Armazenando o resultado de operações:
+x <- 2 + 5
+y = 5 - 2
+2 * 5 -> w
+z <- 8 / 2
+
+resultado <- (((x - y) * w) ^ z)/(x ^ z)
+resultado
## [1] 1066.222
+Chamando funções:
+sum(1, 3, 5)
## [1] 9
+a <- rep("Aluno", times = 3)
+a
## [1] "Aluno" "Aluno" "Aluno"
+Estas funções buscam e exibem a documentação de funções:
+help(sum)
+?sd
+??plot
Estas funções manipulam o diretório de trabalho:
+# Verifica o caminho para o diretório de trabalho
+getwd()
+
+# Define o diretório de trabalho
+setwd()
+
+# Lista os arquivos presentes no diretório de trabalho
+list.files()
+
+# Carrega um arquivo binário do diretório de trabalho para o ambiente
+load()
+
+# Salva o conteúdo de uma variável no diretório de trabalho
+save()
# O comando abaixo cria um objeto o qual iremos salvar.
+esportes <- data.frame(nomes = c("Carlos", "Roberto", "Olivio", "Jomar"),
+ esportes = c("futebol", "remo", "sumo", "maratona"))
+
+# Aqui nós salvamos o objeto "esportes" no HD
+save(esportes, file = "./dados_curso/outputs/intro_r/esportes.RData")
+
+# O comando "rm" remove o objeto do ambiente global
+rm(esportes)
+
+# Aqui carregamos o objeto "esportes" para o ambiente global
+load("./dados_curso/outputs/intro_r/esportes.RData")
+
+# Aqui deletamos o arquivo do HD
+file.remove("./dados_curso/outputs/intro_r/esportes.RData")
## [1] TRUE
+Função de concatenação c()
:
number <- c(1, 2, 3, 4, 5)
+letter <- c("x", "y", "z", "w", "j")
+logico <- c(TRUE, FALSE, FALSE, TRUE, FALSE)
+seq <- 1:10
+complexo <- 4i
A função class()
pode ser usada para acessar a classe de
+um determinado objeto:
class(number)
## [1] "numeric"
+Vetores comportam apenas uma classe de elementos. Quando um vetor é
+criado com valores pertecentes a classes distintas, é feita uma
+conversão implícita. Um valor logical
é convertido para
+numeric
, e um valor numeric
é convertido para
+character
:
class(c(1, 2, 3))
## [1] "numeric"
+class(c("1", "2", "3"))
## [1] "character"
+class(c(TRUE, FALSE, FALSE))
## [1] "logical"
+class(c("TRUE", "FALSE", "FALSE"))
## [1] "character"
+class(c(1, "a", TRUE))
## [1] "character"
+class(c(1, "a"))
## [1] "character"
+class(c(1, T))
## [1] "numeric"
+class(c("a", T))
## [1] "character"
+Com esta hierarquia, é possível somar valores lógicos, sendo
+TRUE
equivalente a 12, e FALSE
equivalente a 0:
logical <- c(TRUE, FALSE, FALSE, TRUE, FALSE)
+sum(logico)
## [1] 2
+Uma conversão explícita pode ser feita com as funções
+as.<nome da classe>
:
x <- 0:10
+x
## [1] 0 1 2 3 4 5 6 7 8 9 10
+class(x)
## [1] "integer"
+a <- as.numeric(x)
+a
## [1] 0 1 2 3 4 5 6 7 8 9 10
+class(a)
## [1] "numeric"
+b <- as.character(x)
+b
## [1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
+class(b)
## [1] "character"
+c <- as.logical(x)
+c
## [1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
+class(c)
## [1] "logical"
+Valores não disponíveis são representados por NA
+(Not Available), e valores impossíveis, como o resultado de uma
+divisão por 0, são representados por NaN
(Not a
+Number).
x <- c(1, 2, 3, NA)
+y <- c("a", "b", "c", NA)
+is.na(x)
## [1] FALSE FALSE FALSE TRUE
+w <- rep(NA, 10)
+w
## [1] NA NA NA NA NA NA NA NA NA NA
+class(w)
## [1] "logical"
+z <- rep(NA_integer_, 10)
+z
## [1] NA NA NA NA NA NA NA NA NA NA
+class(z)
## [1] "integer"
+a <- c(1, 3, NA, 7, 9)
+sum(a)
## [1] NA
+sum(a, na.rm = TRUE)
## [1] 20
+Todos os objetos possuem atributos:
+x <- 1:5
+x
## [1] 1 2 3 4 5
+length(x)
## [1] 5
+dim(x)
## NULL
+attributes(x)
## NULL
+names(x) <- c("a", "b", "c", "d", "e")
+x
## a b c d e
+## 1 2 3 4 5
+attributes(x)
## $names
+## [1] "a" "b" "c" "d" "e"
+Um vetor da classe factor
é um vetor categórico que
+possui o atributo levels
:
x <- factor(c("s", "n", "n", "s", "s"))
+z <- factor(c("alto", "baixo", "medio"))
+x
## [1] s n n s s
+## Levels: n s
+z
## [1] alto baixo medio
+## Levels: alto baixo medio
+Usamos []
para acessar elementos de vetores:
letter <- c("x", "y", "z", "w", "j")
+
+# Acessa o segundo elemento do vetor
+letter[2]
## [1] "y"
+# Podemos usar sequências de valores
+letter[2:4]
## [1] "y" "z" "w"
+# Usamos a função c() para valores não contíguos
+letter[c(1, 4)]
## [1] "x" "w"
+# Usamos números negativos para excluir um ou mais valores
+letter[-2]
## [1] "x" "z" "w" "j"
+letter[c(-2, -5)]
## [1] "x" "z" "w"
+# Podemos criar índices numéricos
+idx <- c(1, 4)
+letter[idx]
## [1] "x" "w"
+x <- 1:10
+
+# Podemos usar operadores relacionais como filtros
+x[x > 7]
## [1] 8 9 10
+# Também funciona com caracteres, levando em consideração a ordem lexicográfica
+letter[letter > "k"]
## [1] "x" "y" "z" "w"
+letter[letter < "k"]
## [1] "j"
+letter == "z"
## [1] FALSE FALSE TRUE FALSE FALSE
+Funções para identificar valores extremos:
+# Definindo uma semente para a geração de valores aleatórios
+set.seed(1)
+s <- sample(-1000:1000, 200)
+
+# Procura a posição do maior valor
+which.max(s)
## [1] 126
+# Exibe o maior valor
+max(s)
## [1] 997
+# Exibe o menor valor
+min(s)
## [1] -982
+# Exibe o intervalo dos valores do vetor
+range(s)
## [1] -982 997
+# Cria um vetor lógico
+s > 0
## [1] TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE
+## [13] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE TRUE FALSE
+## [25] FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE
+## [37] FALSE TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE
+## [49] FALSE TRUE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
+## [61] TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE
+## [73] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE
+## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE
+## [97] FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE
+## [109] FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE TRUE TRUE
+## [121] FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE
+## [133] TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE
+## [145] TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE
+## [157] FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE
+## [169] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE
+## [181] TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE
+## [193] TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE
+# Cria um vetor com as posições que satisfazem o comando
+which(s > 0)
## [1] 1 2 6 10 11 13 14 15 16 18 19 20 21 23 27 28 31 32
+## [19] 38 40 42 43 45 50 52 56 61 62 63 66 67 69 70 72 74 75
+## [37] 77 78 79 80 81 85 86 87 88 89 90 91 93 94 95 99 100 102
+## [55] 105 110 111 113 117 118 119 120 122 123 124 126 130 131 133 134 136 138
+## [73] 142 143 145 146 147 148 149 151 153 154 156 161 163 166 168 169 170 177
+## [91] 178 181 182 185 187 190 191 192 193 194 198
+Funções de ordenamento:
+x <- c(3, 8, 2, 1, 5, 9, 7, 7, 3)
+x
## [1] 3 8 2 1 5 9 7 7 3
+# Ordena um vetor
+sort(x)
## [1] 1 2 3 3 5 7 7 8 9
+sort(x, decreasing = T)
## [1] 9 8 7 7 5 3 3 2 1
+# Informa a ordem na qual cada elemento deve ser acessado para exibir o conteúdo do vetor em ordem crescente
+order(x)
## [1] 4 3 1 9 5 7 8 2 6
+# Exibe o conteúdo do vetor de forma aleatória, e uma única vez, cada posição
+sample(x)
## [1] 2 5 3 3 8 1 7 7 9
+# Elimina as replicatas
+unique(x)
## [1] 3 8 2 1 5 9 7
+# Exibe um vetor lógico referente à posição das replicatas
+duplicated(x)
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
+Matrizes são vetores bidimensionais que possuem o atributo dimensão. +Por serem vetores, comportam apenas uma classe de elementos:
+x <- 1:20
+x
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
+attributes(x)
## NULL
+m <- matrix(x, 4, 5)
+m
## [,1] [,2] [,3] [,4] [,5]
+## [1,] 1 5 9 13 17
+## [2,] 2 6 10 14 18
+## [3,] 3 7 11 15 19
+## [4,] 4 8 12 16 20
+attributes(m)
## $dim
+## [1] 4 5
+dim(x) <- c(4,5)
+x
## [,1] [,2] [,3] [,4] [,5]
+## [1,] 1 5 9 13 17
+## [2,] 2 6 10 14 18
+## [3,] 3 7 11 15 19
+## [4,] 4 8 12 16 20
+identical(x, m)
## [1] TRUE
+a <- 1:5
+b <- -1:-5
+c <- c(3, 6, 4, 9, 1)
+
+# A função cbind() concatena colunas
+m <- cbind(a, b, c)
+m
## a b c
+## [1,] 1 -1 3
+## [2,] 2 -2 6
+## [3,] 3 -3 4
+## [4,] 4 -4 9
+## [5,] 5 -5 1
+# A função rbind() concatena linhas
+m1 <- rbind(a, b, c)
+m1
## [,1] [,2] [,3] [,4] [,5]
+## a 1 2 3 4 5
+## b -1 -2 -3 -4 -5
+## c 3 6 4 9 1
+# Elementos são acessados pelos índices das duas dimenções [linha, coluna]
+m[1,3]
## c
+## 3
+# Toda a linha
+m[1, ]
## a b c
+## 1 -1 3
+m[2:3, ]
## a b c
+## [1,] 2 -2 6
+## [2,] 3 -3 4
+# Atribuição
+m[1,] <- NA
+m
## a b c
+## [1,] NA NA NA
+## [2,] 2 -2 6
+## [3,] 3 -3 4
+## [4,] 4 -4 9
+## [5,] 5 -5 1
+Um array
é um vetor que possui mais de duas
+dimensões:
# Criando um vetor multidimensional com 4 matrizes de 5 linhas e 10 colunas
+ar <- array(1:200, c(5, 10, 4))
+ar
## , , 1
+##
+## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 1 6 11 16 21 26 31 36 41 46
+## [2,] 2 7 12 17 22 27 32 37 42 47
+## [3,] 3 8 13 18 23 28 33 38 43 48
+## [4,] 4 9 14 19 24 29 34 39 44 49
+## [5,] 5 10 15 20 25 30 35 40 45 50
+##
+## , , 2
+##
+## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 51 56 61 66 71 76 81 86 91 96
+## [2,] 52 57 62 67 72 77 82 87 92 97
+## [3,] 53 58 63 68 73 78 83 88 93 98
+## [4,] 54 59 64 69 74 79 84 89 94 99
+## [5,] 55 60 65 70 75 80 85 90 95 100
+##
+## , , 3
+##
+## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 101 106 111 116 121 126 131 136 141 146
+## [2,] 102 107 112 117 122 127 132 137 142 147
+## [3,] 103 108 113 118 123 128 133 138 143 148
+## [4,] 104 109 114 119 124 129 134 139 144 149
+## [5,] 105 110 115 120 125 130 135 140 145 150
+##
+## , , 4
+##
+## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 151 156 161 166 171 176 181 186 191 196
+## [2,] 152 157 162 167 172 177 182 187 192 197
+## [3,] 153 158 163 168 173 178 183 188 193 198
+## [4,] 154 159 164 169 174 179 184 189 194 199
+## [5,] 155 160 165 170 175 180 185 190 195 200
+# Acessando a primeira matriz [linha, coluna, matriz]
+ar[,,1]
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 1 6 11 16 21 26 31 36 41 46
+## [2,] 2 7 12 17 22 27 32 37 42 47
+## [3,] 3 8 13 18 23 28 33 38 43 48
+## [4,] 4 9 14 19 24 29 34 39 44 49
+## [5,] 5 10 15 20 25 30 35 40 45 50
+Listas são tipos especiais de vetores que comportam elementos de +diferentes classes:
+a <- c(1, 3, NA, 7, 9)
+b <- matrix(1:200, 20,10)
+c <- "Gol do Gremio"
+z <- factor(c("alto", "baixo", "medio"))
+ls <- list(a, b, c, z)
+
+# Cada elemento da lista aparece com [[]]
+ls
## [[1]]
+## [1] 1 3 NA 7 9
+##
+## [[2]]
+## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 1 21 41 61 81 101 121 141 161 181
+## [2,] 2 22 42 62 82 102 122 142 162 182
+## [3,] 3 23 43 63 83 103 123 143 163 183
+## [4,] 4 24 44 64 84 104 124 144 164 184
+## [5,] 5 25 45 65 85 105 125 145 165 185
+## [6,] 6 26 46 66 86 106 126 146 166 186
+## [7,] 7 27 47 67 87 107 127 147 167 187
+## [8,] 8 28 48 68 88 108 128 148 168 188
+## [9,] 9 29 49 69 89 109 129 149 169 189
+## [10,] 10 30 50 70 90 110 130 150 170 190
+## [11,] 11 31 51 71 91 111 131 151 171 191
+## [12,] 12 32 52 72 92 112 132 152 172 192
+## [13,] 13 33 53 73 93 113 133 153 173 193
+## [14,] 14 34 54 74 94 114 134 154 174 194
+## [15,] 15 35 55 75 95 115 135 155 175 195
+## [16,] 16 36 56 76 96 116 136 156 176 196
+## [17,] 17 37 57 77 97 117 137 157 177 197
+## [18,] 18 38 58 78 98 118 138 158 178 198
+## [19,] 19 39 59 79 99 119 139 159 179 199
+## [20,] 20 40 60 80 100 120 140 160 180 200
+##
+## [[3]]
+## [1] "Gol do Gremio"
+##
+## [[4]]
+## [1] alto baixo medio
+## Levels: alto baixo medio
+# A função vector() pode criar listas vazias
+ls1 <- vector("list", 5)
+ls1
## [[1]]
+## NULL
+##
+## [[2]]
+## NULL
+##
+## [[3]]
+## NULL
+##
+## [[4]]
+## NULL
+##
+## [[5]]
+## NULL
+Listas podem ser acessadas com os operadores []
,
+[[]]
e $
(para listas nomeadas):
# [] extrai uma lista
+ls[1]
## [[1]]
+## [1] 1 3 NA 7 9
+# [[]] extrai o objeto interno
+ls[[1]]
## [1] 1 3 NA 7 9
+class(ls[1])
## [1] "list"
+class(ls[[1]])
## [1] "numeric"
+# Posição na lista e posição no elemento
+ls[[c(1,2)]]
## [1] 3
+ls[[2]][2,]
## [1] 2 22 42 62 82 102 122 142 162 182
+names(ls) <- c("Arilson", "Roger", "Paulo Nunes", "Jardel")
+ls$Roger
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
+## [1,] 1 21 41 61 81 101 121 141 161 181
+## [2,] 2 22 42 62 82 102 122 142 162 182
+## [3,] 3 23 43 63 83 103 123 143 163 183
+## [4,] 4 24 44 64 84 104 124 144 164 184
+## [5,] 5 25 45 65 85 105 125 145 165 185
+## [6,] 6 26 46 66 86 106 126 146 166 186
+## [7,] 7 27 47 67 87 107 127 147 167 187
+## [8,] 8 28 48 68 88 108 128 148 168 188
+## [9,] 9 29 49 69 89 109 129 149 169 189
+## [10,] 10 30 50 70 90 110 130 150 170 190
+## [11,] 11 31 51 71 91 111 131 151 171 191
+## [12,] 12 32 52 72 92 112 132 152 172 192
+## [13,] 13 33 53 73 93 113 133 153 173 193
+## [14,] 14 34 54 74 94 114 134 154 174 194
+## [15,] 15 35 55 75 95 115 135 155 175 195
+## [16,] 16 36 56 76 96 116 136 156 176 196
+## [17,] 17 37 57 77 97 117 137 157 177 197
+## [18,] 18 38 58 78 98 118 138 158 178 198
+## [19,] 19 39 59 79 99 119 139 159 179 199
+## [20,] 20 40 60 80 100 120 140 160 180 200
+Um data.frame
é um tipo especial de lista, onde todos os
+elementos devem possuir o mesmo length. Por ser uma lista, cada posição
+comporta elementos de diferentes classes. Do ponto de vista prático, o
+data.frame
funciona como uma planilha bidimensional formado
+por vetores de mesmo tamanho, sendo cada vetor uma coluna:
number <- c(1, 2, 3, 4, 5)
+letter <- c("x", "y", "z", "w", "j")
+logical <- c(TRUE, FALSE, FALSE, TRUE, FALSE)
+seq <- 1:10
+dt <- data.frame(number, letter, logical)
+class(dt)
## [1] "data.frame"
+# Usamos $ para acessar as colunas de um data.frame
+dt$letter
## [1] "x" "y" "z" "w" "j"
+# Vetores de caracteres são interpretados como fatores
+class(dt$letter)
## [1] "character"
+# Argumento stringsAsFactors = F altera este comportamento padrão
+dt <- data.frame(number, letter, logical, stringsAsFactors = F)
+dt$letter
## [1] "x" "y" "z" "w" "j"
+class(dt$letter)
## [1] "character"
+# Data.frames possuem colnames e rownames como atributos
+attributes(dt)
## $names
+## [1] "number" "letter" "logical"
+##
+## $class
+## [1] "data.frame"
+##
+## $row.names
+## [1] 1 2 3 4 5
+colnames(dt)
## [1] "number" "letter" "logical"
+row.names(dt)
## [1] "1" "2" "3" "4" "5"
+# Acessamos data.frames da mesma forma que matrizes
+dt[5,2]
## [1] "j"
+Para acessar data.frames podemos usar os operadores []
,
+[[]]
e $
:
dt <- data.frame(number=c(1, 2, 3, 4, 5),
+ letter = c("x", "y", "z", "w", "j"),
+ logical = c(TRUE, FALSE, FALSE, TRUE, FALSE))
+
+# [[ ]] acessa cada coluna por posição
+dt[[1]]
## [1] 1 2 3 4 5
+# [ ] acessa as coordenadas [linha, coluna]
+dt[,1]
## [1] 1 2 3 4 5
+# $ acessa a coluna pelo nome
+dt$number
## [1] 1 2 3 4 5
+# Carrega o data.frame mtcars
+cars<-mtcars
+
+# Mostra as 6 primeiras linhas
+head(cars)
## mpg cyl disp hp drat wt qsec vs am gear carb
+## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
+## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
+## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
+## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
+## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
+## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
+# Mostra as 6 ultimas linhas
+tail(cars)
## mpg cyl disp hp drat wt qsec vs am gear carb
+## Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.7 0 1 5 2
+## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.9 1 1 5 2
+## Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.5 0 1 5 4
+## Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.5 0 1 5 6
+## Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.6 0 1 5 8
+## Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.6 1 1 4 2
+# data.frames possuem colnames e rownames
+colnames(dt)
## [1] "number" "letter" "logical"
+row.names(dt)
## [1] "1" "2" "3" "4" "5"
+# Podemos alterar colnames e rownames
+row.names(dt) <- c("a", "b", "c", "d", "e")
+
+# Alterando apenas a posição 2
+colnames(dt)[2] <- "letras"
+
+# Podemos alterar valores específicos de um data.frame
+dt[3,1] <- "10"
+dt$logical <- as.numeric(dt$logical)
+dt$letras <- NA
É possível verificar as ocorrencias de um data.frame em outro:
+biometria <- data.frame(nomes = c("Carlos", "Roberto", "Olivio", "Joel"),
+ altura = c(180, 187, 155, 168),
+ peso = c(80, 90, 98, 64))
+biometria
## nomes altura peso
+## 1 Carlos 180 80
+## 2 Roberto 187 90
+## 3 Olivio 155 98
+## 4 Joel 168 64
+esportes <- data.frame(nomes = c("Carlos", "Roberto", "Olivio", "Jomar"),
+ esportes = c("futebol", "remo", "sumo", "maratona"))
+esportes
## nomes esportes
+## 1 Carlos futebol
+## 2 Roberto remo
+## 3 Olivio sumo
+## 4 Jomar maratona
+# Retorna um vetor lógico
+biometria$nomes %in% esportes$nomes
## [1] TRUE TRUE TRUE FALSE
+# Pode ser usado como índice
+idx <- biometria$nomes %in% esportes$nomes
+x <- biometria[idx,]
+x
## nomes altura peso
+## 1 Carlos 180 80
+## 2 Roberto 187 90
+## 3 Olivio 155 98
+# Ordenando data.frames por uma coluna
+biometria <- biometria[with(biometria, order(altura)), ]
+biometria
## nomes altura peso
+## 3 Olivio 155 98
+## 4 Joel 168 64
+## 1 Carlos 180 80
+## 2 Roberto 187 90
+Unindo data.frames com a função merge()
:
# Independe da ordem dos data.frames, a busca é feita pelo nome, não pela ordem.
+# O resultado sempre virá em ordem alfabética
+unido <- merge(biometria, esportes, by = "nomes")
+unido
## nomes altura peso esportes
+## 1 Carlos 180 80 futebol
+## 2 Olivio 155 98 sumo
+## 3 Roberto 187 90 remo
+# As informações não disponíveis são preenchidas por NA
+# com todos os presentes no primeiro
+unido <- merge(biometria, esportes, by = "nomes", all.x = TRUE)
+unido
## nomes altura peso esportes
+## 1 Carlos 180 80 futebol
+## 2 Joel 168 64 <NA>
+## 3 Olivio 155 98 sumo
+## 4 Roberto 187 90 remo
+# Com todos os presentes no segundo
+unido <- merge(biometria, esportes, by = "nomes", all.y = TRUE)
+unido
## nomes altura peso esportes
+## 1 Carlos 180 80 futebol
+## 2 Jomar NA NA maratona
+## 3 Olivio 155 98 sumo
+## 4 Roberto 187 90 remo
+# Com todos presentes
+unido <- merge(biometria, esportes, by = "nomes", all = TRUE)
+unido
## nomes altura peso esportes
+## 1 Carlos 180 80 futebol
+## 2 Joel 168 64 <NA>
+## 3 Jomar NA NA maratona
+## 4 Olivio 155 98 sumo
+## 5 Roberto 187 90 remo
+Pacotes agrupam funções projetadas para atacar um problema
+específico. Como exemplo é possível citar o pacote
+data.table
do CRAN, que possui funções para
+manipulação de grandes quantidades de dados. Pacotes disponíveis em
+repositórios (como o CRAN, Bioconductor e Github) podem ser instalados
+por meio de poucas linhas de comando.
Verificando se um pacote já está instalado:
+require(affy)
O CRAN é o principal +repositório de pacotes do R e possui pacotes com finalidades variadas: +importação e exportação de dados, manipulação de grafos, desenvolvimento +de pacotes, plotagem, paralelismo.
+Instalando o pacote devtools
do CRAN:
install.packages("devtools")
O Bioconductor é o principal +repositório de pacotes voltados para a bioinformática. Todo pacote do +Bioconductor é classificado em uma de três categorias: pacote de +anotação (bancos de dados e dicionários), pacote de dados de +experimentos (datasets relacionados a um experimento) ou pacote de +software (funcionalidades).
+Instalando o pacote geneplast
do Bioconductor:
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
+BiocManager::install("geneplast")
O processo de alinhamento e/ou quantificação gera uma tabela de +contagens de features (genes, transcritos ou exons, por +exemplo) a qual será utilizada para a análise de expressão diferencial. +Nesta tabela, as features estão representadas nas linhas e as +amostras estão representadas nas colunas. Além disso, para a análise de +expressão diferencial, é necessário saber, pelo menos, a qual grupo +experimental cada amostra pertence.
+Para este curso, analisaremos a expressão de genes e transcritos de +amostras de adenocarcinoma de pulmão de estadiamento precoce (da sigla +em inglês esLUAD, early-stage lung adenocarcinoma). +Nosso objetivo é comparar a expressão gênica de amostras de +tumores metastáticos (invasive) e tumores não-invasivos +(indolent) por meio de metologias de expressão +diferencial amplamente difundidas. Além das informações sobre os grupos +de comparação, alguns dados clínicos adicionais dos pacientes que +forneceram as amostras podem ser obtidos, como idade, sexo e +características histológicas dos tumores. Para mais informações sobre +este conjunto de dados, o artigo está disponível neste link.
+Utilizaremos pacotes disponíveis no CRAN e no Bioconductor. A seguir, +listamos os pacotes de ambas as fontes e os passos para a +instalação.
+# Lista de pacotes do CRAN
+packs_cran <- c("dplyr", "ggplot2", "pheatmap", "UpSetR", "igraph", "classInt")
+
+# Lista de pacotes do Bioconductor
+packs_bioc <- c("tximport", "GenomicFeatures", "PCAtools", "DESeq2",
+ "DRIMSeq", "stageR", "biomaRt", "RedeR", "clusterProfiler",
+ "TreeAndLeaf", "GOSemSim", "org.Hs.eg.db", "CEMiTool",
+ "transcriptogramer", "edgeR")
+
+# Instalar pacotes CRAN
+lapply(packs_cran, function(i) {
+ if(!require(i, character.only = TRUE)) install.packages(i)
+})
+
+# Instalar pacotes Bioconductor
+if(!require("BiocManager")) install.packages("BiocManager")
+lapply(packs_bioc, function(i) {
+ if(!require(i, character.only = TRUE)) BiocManager::install(i)
+})
As reads foram quantificadas com o pseudoalinhador kallisto.
+Para cada amostra, o programa gera um diretório com um arquivo
+.tsv
, um arquivo .h5
e um arquivo de log
+run_info.json
, com informações sobre o processo de
+quantificação para aquela amostra. Utilizaremos apenas os arquivos
+.tsv
.
Cada arquivo .tsv
contém a quantificação dos
+transcritos. Este arquivo possui as seguintes informações para cada
+transcrito:
target_id
: identificador do transcrito (ex.:
+ENST00000631435.1);length
: tamanho do transcrito;eff_length
: tamanho efetivo dos transcritos (usado no
+cálculo do TPM);est_counts
: valores de contagem estimados;tpm (transcripts per million)
: abundância em
+transcritos por milhão (TPM);Utilizaremos o pacote tximport
+para fazer a importação dos valores de contagem dos arquivos do
+kallisto. O tximport
contém métodos para a
+sumarização dos valores de contagem de transcritos em genes, que serão
+utilizados para a expressão diferencial de genes. Além disso, o pacote
+dispõe de outros métodos de normalização para os dados de contagem, os
+quais são necessários nas abordagens em nível de transcrito. Além do
+kallisto, outros softwares de alinhamento/quantificação, como o Salmon,
+Alevin, RSEM e StringTie, também fornecem os valores de contagem de
+transcritos.
Para a sumarização das contagens de transcritos em genes, temos que +obter um dicionário relacionando cada um dos genes com seus respectivos +transcritos. É aconselhável que se use o mesmo arquivo de anotação +utilizado na etapa de alinhamento/quantificação. Aqui, será usado o +arquivo GTF do Ensembl/GENCODE (release 97) da etapa de quantificação do +kallisto. Neste exemplo, será criado um banco de dados do GTF e este +banco será usado para criarmos o dicionário.
+# Importar metadado
+meta <- read.csv("./dados_curso/inputs/metadata/metadata.csv")
+
+# O metadado a ser usado pelo DESeq2 precisa ter os rownames idênticos às colunas
+# da tabela de contagem
+rownames(meta) <- meta$run
+meta$run <- NULL
library(tximport)
+library(GenomicFeatures)
+
+# Caminho do diretório do GTF
+gtf <- "./dados_curso/inputs/annotation/Homo_sapiens.GRCh38.97.gtf.gz"
+
+# Caminho para o banco de dados a ser criado
+txdb.filename <- "./dados_curso/inputs/annotation/Homo_sapiens.GRCh38.97.gtf.sqlite"
+
+# Criar banco, caso não exista
+if(!("Homo_sapiens.GRCh38.97.gtf.sqlite" %in% list.files("./dados_curso/inputs/annotation/"))) {
+ txdb <- makeTxDbFromGFF(gtf, format = "gtf")
+ saveDb(txdb, txdb.filename)
+}
+
+# Carregar banco de dados
+txdb <- loadDb(txdb.filename)
+
+# Deste banco, selecionar os IDs de genes (GENEID) e transcritos (TXNAME)
+txdf <- AnnotationDbi::select(txdb, keys(txdb, "GENEID"), "TXNAME", "GENEID")
+tab <- table(txdf$GENEID)
+txdf$ntx <- tab[match(txdf$GENEID, names(tab))]
+
+# Criar dataframe com duas colunas: tx (transcritos) e gene (genes). Esta etapa é requerida
+# pelo tximport para fazer a sumarização de transcritos em genes.
+tx2gene <- data.frame(tx = txdf$TXNAME, gene = txdf$GENEID, stringsAsFactors = F)
+
+# Salvar dataframe para consultas posteriores
+save(tx2gene, file = "./dados_curso/outputs/annotation/dict_gene_tx.RData")
+
+# Obter os caminhos dos arquivos
+files <- list.files(path = "./dados_curso/inputs/kallisto", pattern = "tsv", recursive = TRUE, full.names = TRUE)
+files <- files[base::match(rownames(meta), sapply(strsplit(files, "\\/"), "[[", 5))]
+names(files) <- rownames(meta)
# Importar dados de contagem
+txi_gene <- tximport(files = files, type = "kallisto", tx2gene = tx2gene, ignoreTxVersion = T)
+
+save(txi_gene, file = "./dados_curso/outputs/tximport/txi_gene.RData")
O objeto txi
é uma lista contém três matrizes:
abundance
, matriz de abundâncias normalizadas (o
+kallisto fornece valores de abundância em TPM);counts
, matriz de valores de contagem estimados (como o
+kallisto estima as contagens, pode apresentar números
+fracionários);length
, matriz com os tamanhos efetivos dos
+transcritos.Para importar os valores de contagens de transcritos, para posterior
+análise de expressão diferencial de transcritos (DTE), basta usar a
+opção txOut = T
na função tximport
.
txi_tx <- tximport(files = files, type = "kallisto", txOut = T, ignoreTxVersion = T)
+save(txi_tx, file = "./dados_curso/outputs/tximport/txi_tx.RData")
Para a análise do uso diferencial de isoformas (DTU), os valores dos
+transcritos precisam estar normalizados de acordo com o tamanho da
+biblioteca. O tximport proporciona o argumento
+countsFromAbundance = scaledTPM
para a normalização das
+bibliotecas em TPM.
txi_tx_dtu <- tximport(files = files, type = "kallisto", txOut = T, countsFromAbundance = "scaledTPM", ignoreTxVersion = T)
+
+save(txi_tx_dtu, file = "./dados_curso/outputs/tximport/txi_tx_dtu.RData")
Antes de procedermos para a análise de expressão diferencial, é +necessário realizar uma análise exploratória do conjunto de dados. Aqui, +buscamos verificar a distribuição da expressão gênica nos grupos de +comparações (no nosso caso, o subtipo, se metastático ou não invasivo). +Para tal, precisamos de técnicas estatísticas robustas, as quais +consigam considerar a multidimensionalidade dos dados de expressão. De +uma forma simples, o caráter multidimensional dos dados de expressão é +decorrente do fato de que lidamos com milhares de genes ao mesmo tempo +e, portanto, precisamos “resumir” as informações destes genes. Por isso, +técnicas de redução de dimensionalidade, como a análise de componentes +principais (PCA) e técnicas de agrupamento, como a clusterização +hierárquica, são importantes para entendermos a estrutura do nosso dado +de expressão.
+A Análise de Componentes Principais (PCA) é uma técnica para análise +multivariada que visa sumarizar as informações provenientes de uma +amostra em novas características, ou componentes principais. Assim, para +nossos propósitos, a informação a ser captada é a variabilidade da +expressão dos genes de uma amostra. Por meio desta técnica, podemos +resumir uma amostra em um ponto do espaço com dimensões.
+A análise de componentes principais decompõe as informações dos genes
+em espaços de variabilidade. Quanto mais variabilidade há na expressão
+de um gene, mais este gene contribui para este espaço. Neste processo,
+definimos componentes que indicam em quais direções do espaço há maior
+variabilidade. A componente principal 1 (PC1) é direção onde há maior
+variabilidade, a componente principal 2 é direção onde há a segunda
+maior variabilidade, a componente principal 3 é direção onde há a
+terceira maior variabilidade, e assim por diante. Em geral, plotamos as
+componentes principais 1 e 2 e verificamos como as amostras se
+relacionam. O pacote PCAtools
+contém inúmeras funções customizadas para esta tarefa.
Antes de realizarmos análise, aplicaremos uma transformação dos dados +de expressão. Esta transformação é necessária para reduzir a dependência +da variância em relação à média, já que uma das principais +características de dados de expressão gênica é a heteroscedasticidade. +Ou seja, genes com expressão média pequena tem variância grande, +enquanto que genes com expressão média alta possuem menor variância. +Este passo de transformação dos dados não é necessário +para os testes de expressão diferencial, como veremos adiante. Alguns +algoritmos de transformação dos dados são o variance stabilizing +transformation (VST) e o regularized logarithm (rlog), +ambos implementados no DESeq2.
+library(PCAtools)
+library(DESeq2)
+
+# Carregar dados de expressão de genes
+load("dados_curso/outputs/tximport/txi_gene.RData")
+
+# Criar objeto DESeq2
+dds <- DESeqDataSetFromTximport(txi_gene, colData = meta, design = ~ subtype)
+
+# Transformar os dados por meio do VST
+dds_vst <- vst(dds)
+
+# Obter counts transformados
+counts_genes_vst <- assay(dds_vst)
+
+# Computar PCA
+pca_lung <- pca(counts_genes_vst, metadata = meta)
+
+# Plotar PC1 e PC2
+biplot(pca_lung, colby = "histology", shape = "subtype", legendPosition = "left")
# Proporção de variabilidade de cada componente
+screeplot(pca_lung)
Em resumo, pela PCA, podemos ver que as amostras se agrupam pelo +subtipo do tumor (não-invasivo ou metastático) e, em geral, pelo tipo de +histologia apresentada.
+Clusterizar implica em agrupar amostras de acordo com a similaridade
+de expressão gênica. Diferentes algoritmos realizam este tipo de tarefa,
+que em geral consiste em calcular distâncias entre quaisquer par de
+amostras. Amostras mais similares entre si são representadas com maior
+proximidade no dendrograma, criando, assim, uma estrutura hierárquica.
+Neste exemplo, calcularemos as distâncias euclidianas
+(method = "euclidian"
, da função dist
) entre
+as amostras e o método de agrupamento será o de ligação completa
+(method = "complete"
, da função hclust
).
# Transpor matriz de counts transformados
+counts_genes_vst_t <- t(counts_genes_vst)
+
+# Setar nomes dos grupos
+rownames(counts_genes_vst_t) <- meta$subtype
+
+# Clusterização hierárquica
+d <- dist(as.matrix(counts_genes_vst_t), method = "euclidean")
+clusters <- hclust(d, method = "complete")
+
+# Plotar
+par(mar = c(1, 1, 1, 1))
+plot(clusters)
Pelo dendrograma, verificamos que as amostras, em geral, se separam +por subtipo.
+A análise de expressão diferencial nada mais é que um modelo de +regressão, no qual a variável resposta é a expressão de cada um dos +genes e a variável explanatória é o grupo de comparação (no nosso caso, +os subtipos do tumor). Este modelo de expressão é aplicado a cada um dos +genes para verificarmos se há diferença na expressão quando consideramos +os diferentes níveis do grupo de comparação.
+Entretanto, outras variáveis, além do grupo de comparação, podem ser +importantes para explicar a expressão gênica. Como vimos na PCA, as +características histológicas são um importante fator que influencia a +expressão gênica. Além disso, o pesquisador pode julgar, com base em +evidências prévias, que uma variável é relevante para o modelo. Estas +variáveis adicionais podem ser acrescentadas ao modelo de expressão +gênica a fim de melhorar a resolução dos resultados.
+Uma forma de selecionar covariáveis é verificar quais se
+correlacionam significativamente com as componentes principais. A
+seguir, a função eigencorplot
fornece uma maneira fácil de
+visualizarmos as correlações significativas.
eigencorplot(pca_lung, metavars = c("subtype", "histology", "sex", "age"))
As escolha de covariáveis para o modelo é um processo subjetivo.
+Entretanto, sugerimos selecionar as covariáveis que se correlacionam com
+as primeiras componentes. Portanto, pelo gráfico de correlações acima,
+podemos adicionar a covariável histology
ao modelo de
+expressão de genes e transcritos.
Nesta etapa, realizaremos análise de expressão diferencial de genes +(DGE) e transcritos (DTE) e a análise de uso de isoforma (DTU).
+As análises DGE e DTE serão realizadas com o pacote DESeq2
.
+O DESeq2 implementa um modelo linear generalizado baseado na
+distribuição Binomial Negativa para a expressão de genes e transcritos.
+Esta metodologia é uma das mais utilizadas para análise de expressão
+diferencial, entretanto, outras metodologias podem ser utilizadas, como
+edgeR
+e o limma
.
A análise DTU será feita por meio do pacote DRIMSeq
,
+que implementa uma metodologia baseada na distribuição
+Dirichlet-multinomial para modelar a contribuição da expressão dos
+transcritos de cada gene. As três abordagens utilizadas aqui podem ser
+exemplicadas na figura abaixo. Na análise de expressão diferencial de
+genes (DGE), a expressão do gene é sumarizada na abundância coletiva dos
+transcritos anotados para aquele gene. Portanto, inferimos se aquele
+gene teve alteração ou não na sua expressão, dadas diferentes condições.
+No exemplo, temos um aumento significativo da expressão do gene A nas
+amostras do grupo tratado, mas não houve alteração do gene B.
A análise de expressão diferencial de transcritos (DTE) é similar à +de genes, porém compara-se a expressão de cada um dos transcritos de um +gene em diferentes condições. Na figura, podemos observar que houve +aumento da expressão dos transcritos A-2 e A-3 do gene A nas amostras do +grupo tratado. Ainda, houve diminuição da expressão do transcrito B-1 e +aumento da expressão do transcrito B-2 nas amostras TDM do gene B. A +diferença da expressão global do gene B não é significativa, mas esta +aparece quando consideramos a expressão individual de seus +transcritos.
+A análise diferencial de uso de isoforma (DTU) é um tipo de análise +diferencial em nível de transcrito, a qual leva em consideração a +expressão relativa de cada transcrito em relação à expressão total do +gene. Deste modo, compara-se a proporção da expressão de cada transcrito +em diferentes condições e verifica-se, por meio de testes estatísticos, +se houve diferenças significativas nas comparações. Na figura, por +exemplo, verifica-se que, para o gene A, houve uma mudança nas +proporções dos transcritos A-1 e A-2, na comparação entre amostras +controle e amostras do grupo tratado.
+É importante frisar que, independente da metodologia utilizada no
+alinhamento/quantificação das reads, o dado de entrada para o
+DESeq2
precisa ser uma matriz de valores de contagens
+brutos. Portanto, não utilize dados normalizados por
+TPM, R/FPKM ou CPM, por exemplo.
Primeiramente, precisamos criar o objeto DESeqDataSet
.
+Este objeto irá armazenar a matriz de contagens, o dataframe de metadado
+e o modelo da expressão diferencial.
O DESeq2
provém diferentes funções para criarmos este
+objeto a partir de uma matriz de contagens. A função
+DESeqDataSetFromMatrix
cria este objeto a partir de uma
+matriz de contagens. A função DESeqDataSetFromHTSeqCount
+cria este objeto a partir da tabela de contagens proveniente da
+quantificação do HTseq, quando alinhadores como o STAR e o HISAT2 são
+utilizados. Já a função DESeqDataSetFromTximport
cria o
+objeto a partir do dado importado pelo tximport
. Esta
+função é bastante útil para o caso da quantificação é realizada por
+pseudoalinhadores, como o Salmon e o kallisto, por exemplo.
Para verificar em detalhes os argumentos destas funções, ver
+?DESeqDataSet
. Aqui, utilizaremos a função
+DESeqDataSetFromTximport
para fazer a importação dos dados
+de contagem.
Além disso, cada uma das funções de importação possui o argumento
+design
, o qual especifica o modelo de expressão a ser
+utilizado. Na etapa de análise de covariáveis, verificamos que a
+variável histology
é importante para explicar a
+variabilidade da expressão gênica e, portanto, esta variável será
+adicionada ao modelo.
Outro detalhe importante quanto ao fluxograma do DESeq2
+é a necessidade de que haja correspondência entre os nomes das linhas do
+dataframe do metadado (rownames(meta)
) e as colunas de
+contagem (colnames(txi_gene$counts)
).
library(DESeq2)
+
+# Carregar dados de contagem
+load("./dados_curso/outputs/tximport/txi_gene.RData")
+
+# Importar metadado
+meta <- read.csv("./dados_curso/inputs/metadata/metadata.csv")
+rownames(meta) <- meta$run
+meta$run <- NULL
+
+# Checar se as linhas do metadado correspondem às colunas da tabela de contagens
+identical(rownames(meta), colnames(txi_gene$counts))
## [1] TRUE
+# Adicionar as covariáveis ao modelo de expressão. O argumento 'design' especifica o modelo com
+# as variáveis do metadado. Estas devem estar especificadas de acordo com os nomes das colunas
+# do metadado.
+dds <- DESeqDataSetFromTximport(txi_gene,
+ colData = meta,
+ design = ~ histology + subtype)
A função DESeq
congrega diferentes funções. Ela estima
+os fatores de escala de cada amostra, estima a dispersão por gene ou
+transcrito, aplica o modelo e testa a significância dos parâmetros.
dds <- DESeq(dds)
Em seguida, precisamos extrair os resultados do teste. A função
+resultNames
retorna os resultados estimados para todos os
+contrastes possíveis.
Definindo de forma coloquial, um contraste é uma comparação entre os
+grupos de uma variável categórica. Para o nosso exemplo, queremos
+verificar a diferença na expressão das amostras do tumor o metatático e
+o não-invasivo. Assim, nosso contraste de interesse é
+invasive-indolent
. A ordem com a qual o contraste é
+construído interfere na interpretação dos resultados. No caso descrito,
+a referência da nossa comparação é o grupo das amostras de tumor
+não-invasivo.
A função resultNames
retorna todos os contrastes
+(comparações) possíveis de acordo com as variáveis fornecidas no
+modelo:
resultsNames(dds)
## [1] "Intercept" "histology_LPA_vs_AIS"
+## [3] "histology_MP_vs_AIS" "histology_SOL_vs_AIS"
+## [5] "subtype_invasive_vs_indolent"
+Perceba que o DESeq2
retorna o contraste na seguinte
+estrutura:
+<nome da variável>_tratamento_vs_referência
. Para
+extrair os resultados da comparação de interesse, use a função
+results
e passe o contraste de interesse com o argumento
+name
:
DESeq2::results(dds, name = "subtype_invasive_vs_indolent")
## log2 fold change (MLE): subtype invasive vs indolent
+## Wald test p-value: subtype invasive vs indolent
+## DataFrame with 35628 rows and 6 columns
+## baseMean log2FoldChange lfcSE stat pvalue
+## <numeric> <numeric> <numeric> <numeric> <numeric>
+## ENSG00000000003 842.331003 1.507696 0.592114 2.546296 0.0108873
+## ENSG00000000005 0.815306 -0.668878 2.466598 -0.271174 0.7862571
+## ENSG00000000419 427.530483 -0.338197 0.166686 -2.028950 0.0424634
+## ENSG00000000457 408.875269 0.317035 0.341831 0.927461 0.3536872
+## ENSG00000000460 115.161825 0.674323 0.312095 2.160630 0.0307239
+## ... ... ... ... ... ...
+## ENSG00000287988 0.093669 -2.018418 5.132544 -0.393259 0.694128
+## ENSG00000288000 0.339414 -0.713262 5.167800 -0.138020 0.890224
+## ENSG00000288014 0.000000 NA NA NA NA
+## ENSG00000288031 9.571927 0.442393 0.679384 0.651167 0.514938
+## ENSG00000288053 0.383422 2.895889 5.089498 0.568993 0.569361
+## padj
+## <numeric>
+## ENSG00000000003 0.0941317
+## ENSG00000000005 NA
+## ENSG00000000419 0.2147291
+## ENSG00000000457 0.6509532
+## ENSG00000000460 0.1763767
+## ... ...
+## ENSG00000287988 NA
+## ENSG00000288000 NA
+## ENSG00000288014 NA
+## ENSG00000288031 0.769868
+## ENSG00000288053 NA
+Uma outra forma de extrair os contrastes é usando o argumento
+contrast
. Há diferentes formas de usar este argumento (ver
+descrição da função ?results
). Uma delas é passar um vetor
+de caracter de tamanho 3, onde o primeiro elemento é a variável de
+comparação, o segundo elemento é grupo correspondente ao tratamento e o
+terceiro elemento é o grupo correspondente à referência:
# Extrair resultados
+res_dge <- DESeq2::results(dds, contrast = c("subtype", "invasive", "indolent"))
# MA-plot
+plotMA(res_dge)
# Transformar em dataframe
+res_dge <- as.data.frame(res_dge)
O dataframe res_dge
contém os resultados da comparação
+de interesse. O dataframe contém as seguintes informações:
baseMean
: média dos valores de contagens normalizados
+para todas as amostras;log2FoldChange
: métrica que indica o aumento ou
+diminuição da expressão gênica na comparação de interesse;lfcSE
: desvio padrão do log2FC;pvalue
: p-valor do teste;padj
: p-valor do teste ajustado.Uma observação importante do resultado acima é que o
+DESeq2
coloca valores NA
para genes valores de
+contagens baixos. Pode-se considerar que o p-valor para estes genes é 1
+ou próximo de 1 e, por isso, estes não são significativos.
O fold change é uma medida que indica o quanto que um +determinado gene (ou transcrito) teve sua expressão alterada em uma +determinada comparação. Para entendermos esta medida, vamos considerar o +seguinte exemplo:
+Gene | +Expressão média normalizada nas amostras metastáticas | +Expressão média normalizada nas amostras não-invasivas | +Fold Change (comparação metastático VS não-invasivo) | +Log2FoldChange | +
---|---|---|---|---|
A | +10 | +20 | +0.5 | +-1 | +
B | +50 | +10 | +5 | +2.32 | +
Se o intuito da nossa análise é comparar a expressão gênica entre as
+amostras metastáticas e as não-invasivas, o fold change é a
+razão entre a expressão média das amostras metastáticas
+e a das amostras não-invasivas. A expressão média do gene A nas amostras
+metastáticas é a metade daquela das amostras não invasivas; enquanto que
+a expressão média do gene B é 5 vezes maior nas amostras metastáticas.
+Dizemos, portanto, que o gene A está downregulado e o gene B
+está upregulado nas nas amostras metastáticas. A fim de se
+obter uma escala simétrica para os valores de fold change, é
+comum transformá-los para a escala logaritmica de base 2. Perceba que,
+se a nossa comparação fosse indolent-invasive
, a referência
+para a comparação seriam as amostras metastáticas e nossa a
+interpretação seria o contrário: o gene A estaria upregulado e
+o gene B estaria downregulado nas amostras não-invasivas.
O critério de identificação para os genes significativos é variável.
+Entretanto, é imprescindível que a seleção dos genes diferencialmente
+expressos leve em consideração um teste estatístico, como o teste
+realizado pelo DESeq2
. Assim, é necessário considerar o
+p-valor corrigido do teste realizado, uma vez que realizamos milhares
+testes (um para cada gene testado). Aqui, utilizaremos o critério
+p-valor ajustado <= para
+idenficação dos genes diferencialmente expressos.
Podemos visualizar este resultado a partir de um volcano +plot. O volcano plot é um gráfico que relaciona o p-valor com o +log2 Fold Change. Com este gráfico, é fácil ter uma ideia geral sobre a +quantidade de genes downregulados e upregulados na +nossa análise.
+library(dplyr)
+library(ggplot2)
+
+# Tratar os NA's: se houver NA na coluna 'padj', colocar 1. Criar coluna signif, onde identificamos os genes significativos (upregulados ou downregulados) e os genes não significativos
+res_dge <- res_dge %>%
+ mutate(
+ padj = ifelse(is.na(padj), 1, padj),
+ signif =
+ case_when(
+ padj <= 0.05 & log2FoldChange < 0 ~ "Downregulated",
+ padj <= 0.05 & log2FoldChange > 0 ~ "Upregulated",
+ padj > 0.05 ~ "Not significant"
+ )
+ )
+
+ggplot(res_dge, aes(x = log2FoldChange, y = -log10(padj), color = signif)) +
+ geom_point() +
+ scale_color_manual(values = c("Downregulated" = "blue", "Not significant" = "grey", "Upregulated" = "red"))
Vamos plotar o resultado da expressão diferencial de genes em um
+heatmap. De forma análoga, podemos utilizar este gráfico para
+representarmos o resultados das outras análises de expressão. Iremos
+utilizar o pacote pheatmap
, que possui uma interface
+simples para a criação do gráfico. A vantagem deste pacote é
+proporcionar a clusterização das amostras e dos genes de acordo com a
+expressão.
library(pheatmap)
+
+load("./dados_curso/outputs/exp_diff/dds.RData")
+
+# Filtrar os genes diferencialmente expressos (padj <= 0.05)
+res_dge <- subset(res_dge, padj <= 0.05)
+
+# Selecionar os genes DGE
+dge <- rownames(res_dge)
+
+# Obter tabela de expressão
+exp_genes <- DESeq2::counts(dds, normalized = TRUE)
+
+# Desta tabela, verificar quais genes são DGE
+idx <- rownames(exp_genes) %in% dge
+exp_dge <- exp_genes[idx,]
+
+# Transformar valores de expressão em z-score
+exp_dge_z_score <- t(apply(exp_dge, 1, scale, center = T, scale = T))
+colnames(exp_dge_z_score) <- rownames(meta)
+
+# Criar coluna de anotação
+ann_col <- as.data.frame(colData(dds)[, c("subtype", "histology")])
+
+# Plotar gráfico
+pheatmap(exp_dge_z_score, show_rownames = F, annotation_col = ann_col, cluster_cols = F)
Pelo heatmap, podemos ver os genes diferencialmente expressos +segregam perfeitamente as amostras pelo subtipo do tumor.
+Podemos então selecionar apenas os genes diferencialmente expressos e +salvá-los em um dataframe para a análise funcional posterior.
+# Salvar
+save(res_dge, file = "./dados_curso/outputs/exp_diff/res_dge.RData")
Em resumo, identificamos 1447 genes diferencialmente expressos na +análise DGE.
+A análise de expressão de transcritos é similar à análise DGE.
+Utilizaremos os valores de contagens de transcritos, os quais foram
+importados pelo tximport
.
# Carregar dados de contagens de transcritos
+load("./dados_curso/outputs/tximport/txi_tx.RData")
+
+# Criar objeto DESeqDataSet
+dds_tx <- DESeqDataSetFromTximport(txi_tx,
+ colData = meta,
+ design = ~ histology + subtype)
+
+# Calcular fatores de normalização, estimar dispersão e aplicar modelo
+dds_tx <- DESeq(dds_tx)
# Obter resultados
+res_dte <- DESeq2::results(dds_tx, name = "subtype_invasive_vs_indolent")
# MA-plot
+plotMA(res_dte)
# Organizar resultado
+res_dte <- as.data.frame(res_dte)
+rownames(res_dte) <- gsub("\\.\\d+", "", rownames(res_dte))
+
+res_dte <- res_dte[complete.cases(res_dte),]
+
+# Salvar
+save(res_dte, file = "./dados_curso/outputs/exp_diff/res_dte.RData")
Um passo adicional é requerido na análise em nível de transcrito. É +necessário fazer a correção do p-valor de cada transcrito agrupando-os +por gene. Este passo será feito adiante.
+Assim como na análise DTE, para a análise DTU o nosso foco também
+será no transcrito. Utilizaremos o pacote DRIMSeq
, o qual
+implementa um modelo estatístico diferente daquele usado pelo
+DESeq2
. Além disso, assim como na análise DTE, precisaremos
+dar uma atenção especial à correção para múltiplos testes, processo
+importante para a identificação dos transcritos significativos quanto ao
+seu uso diferencial.
O fluxo de análise do DRIMSeq
é um pouco diferente
+daquele do DESeq2
. Primeiramente, precisamos organizar a
+tabela de counts, o dicionário de genes e transcritos.
library(DRIMSeq)
+
+# Carregar dados dicionário gene/transcrito
+load("./dados_curso/outputs/annotation/dict_gene_tx.RData")
+
+# Carregar valores de contagens normalizados
+load("./dados_curso/outputs/tximport/txi_tx_dtu.RData")
+
+# Organizar dataframe de contagens
+counts <- txi_tx_dtu$counts
+rownames(counts) <- gsub("\\.\\d+", "", rownames(counts))
+
+# Verificar se todos os transcritos da tabela de contagens estão representados no dicionário
+all(rownames(counts) %in% tx2gene$tx)
## [1] FALSE
+# Alguns não estão, filtrar a tabela de contagens
+counts <- counts[rownames(counts) %in% tx2gene$tx,]
+tx2gene <- tx2gene[match(rownames(counts), tx2gene$tx),]
+
+# Criar a coluna 'sample_id' no metadado (requerido pelo pacote)
+meta$sample_id <- rownames(meta)
Vamos juntar as informações do dicionário de genes e transcritos e a +tabela de counts em um único dataframe:
+# Preparar dataframe com os valores de contagens de transcritos
+counts_drimseq <- data.frame(gene_id = tx2gene$gene, feature_id = tx2gene$tx, counts, stringsAsFactors = F)
Em seguida, vamos construir o objeto do DRIMSeq:
+# Criar objeto DRIMSeq
+d <- dmDSdata(counts = counts_drimseq, samples = meta)
Nesta etapa, removemos os transcritos pouco expressos:
+# Filtrar transcritos
+d <- dmFilter(d,
+ min_samps_gene_expr = 20, min_samps_feature_expr = 20,
+ min_gene_expr = 10, min_feature_expr = 10)
Agora criamos a matriz modelo, com as variáveis a serem utilizadas na +regressão:
+# Matriz modelo
+design_full <- model.matrix(~ histology + subtype, data = DRIMSeq::samples(d))
+
+colnames(design_full)
## [1] "(Intercept)" "histologyLPA" "histologyMP" "histologySOL"
+## [5] "subtypeinvasive"
+Os próximos passos consistem em estimar a precisão (um parâmetro +utilizado para estimar a variabilidade das proporções de diferentes +transcritos de um gene), ajustar o modelo e testar a significância do +parâmetro para cada transcrito:
+# Calcular precisão
+d <- dmPrecision(d, design = design_full, verbose = 1)
+
+# Ajustar modelo
+d <- dmFit(d, design = design_full, verbose = 1)
+
+# Testar parâmetros
+d <- dmTest(d, coef = "subtypeinvasive")
# Obter as proporções de cada transcrito (isoforma) em cada amostra
+head(DRIMSeq::proportions(d))
# Obter os resultados em nível de gene
+res_dtu_gene <- DRIMSeq::results(d)
+
+# Obter os resultados em nível de transcrito
+res_dtu_tx <- DRIMSeq::results(d, level = "feature")
+
+# Salvar
+save(res_dtu_gene, res_dtu_tx, file = "./dados_curso/outputs/exp_diff/res_dtu.RData")
A correção de múltiplos testes pelos métodos convencionais, como a
+correção de Benjamini & Hochberg, pode ser utilizada para análises
+em nível de transcrito. Entretanto, estas não consideram o fato de que
+há uma quantidade bastante variável de transcritos por genes, o que faz
+com que haja uma penalização maior para genes com poucos transcritos.
+Por isso, a correção para múltiplos testes de transcritos deve ser
+realizada gene a gene. Para esta tarefa, dispomos da metodologia
+implementada no pacote stageR
,
+que realiza os seguintes passos:
library(stageR)
+
+# Verificar quantos transcritos não estão anotados para um gene
+sum(!(rownames(res_dte) %in% tx2gene$tx))
## [1] 4060
+# Remover os transcritos que não estão anotados para nenhum gene
+res_dte <- res_dte[rownames(res_dte) %in% tx2gene$tx,]
+
+# Etapa de screening:
+pScreen <- res_dte$padj
+names(pScreen) <- tx2gene$gene[match(rownames(res_dte), tx2gene$tx)]
+
+# Etapa de confirmação
+pConfirmationTx <- matrix(res_dte$pvalue, ncol = 1)
+rownames(pConfirmationTx) <- rownames(res_dte)
+
+# Criar objeto stageR
+stageRObj <- stageRTx(pScreen = pScreen,
+ pConfirmation = pConfirmationTx,
+ tx2gene = data.frame(transcripts = rownames(res_dte),
+ genes = tx2gene$gene[match(rownames(res_dte), tx2gene$tx)],
+ stringsAsFactors = F),
+ pScreenAdjusted = T)
+
+# Corrigir p-valores para DTE
+stageRObj <- stageWiseAdjustment(object = stageRObj, method = "dte", alpha = 0.05)
+
+# Obter resultados
+padj <- getAdjustedPValues(stageRObj, order = TRUE, onlySignificantGenes = T)
+res_dte_ajustado <- padj[!padj$transcript == 0,]
+
+save(res_dte_ajustado, file = "./dados_curso/outputs/exp_diff/res_dte_ajustado.RData")
# Verificar quantos transcritos não estão anotados para um gene
+sum(!(res_dtu_tx$feature_id %in% tx2gene$tx))
## [1] 0
+# Neste caso, todos estão anotados para um gene
+
+# Etapa de screening:
+pScreen <- res_dtu_tx$adj_pvalue
+names(pScreen) <- res_dtu_tx$gene_id
+
+# Etapa de confirmação
+pConfirmationTx <- matrix(res_dtu_tx$pvalue, ncol = 1)
+rownames(pConfirmationTx) <- res_dtu_tx$feature_id
+
+# Criar objeto stageR
+stageRObj <- stageRTx(pScreen = pScreen,
+ pConfirmation = pConfirmationTx,
+ tx2gene = data.frame(transcripts = res_dtu_tx$feature_id,
+ genes = res_dtu_tx$gene_id,
+ stringsAsFactors = F),
+ pScreenAdjusted = T)
+
+# Corrigir p-valores para DTE
+stageRObj <- stageWiseAdjustment(object = stageRObj, method = "dtu", alpha = 0.05)
+
+# Obter resultados
+padj <- getAdjustedPValues(stageRObj, order = TRUE, onlySignificantGenes = T)
+res_dtu_ajustado <- padj[!padj$transcript == 0,]
+
+save(res_dtu_ajustado, file = "./dados_curso/outputs/exp_diff/res_dtu_ajustado.RData")
Em resumo, as três análises (DGE, DTE e DTU), separadamente, +proporcionaram a seguinte quantidade de genes significativos:
+# DGE
+nrow(res_dge)
## [1] 1447
+# DTE
+length(unique(res_dte_ajustado$geneID))
## [1] 564
+# DTU
+length(unique(res_dtu_ajustado$geneID))
## [1] 28
+Ao todo, temos a seguinte quantidade de genes alterados (em pelo +menos uma das abordagens):
+length(unique(c(rownames(res_dge), res_dte_ajustado$geneID, res_dtu_ajustado$geneID)))
## [1] 1711
+Em comum às três análises, temos a seguinte quantidade de genes +alterados:
+library(UpSetR)
+
+genes_alterados <- list(DGE = unique(rownames(res_dge)),
+ DTE = unique(res_dte_ajustado$geneID),
+ DTU = unique(res_dtu_ajustado$geneID))
+
+upset(fromList(genes_alterados), order.by = "freq")
Em geral, a análise de expressão diferencial é seguida por outras +análises as quais têm o intuito de buscar interpretações funcionais para +o conjunto de genes alterados. Para simplificarmos nosso fluxograma, +iremos utilizar apenas os genes DGE (aqueles identificados como +significativos na análise de expressão diferencial de genes). +Entretanto, é possível obter uma lista de genes que combine os genes +alterados nas três abordagens utlizadas (DGE, DTE e DTU) e realizar as +análises funcionais abordadas neste tutorial.
+Na bioinformática, é comum utilizarmos diferentes bancos de dados
+para fazermos nossas análises. Com isto, frequentemente temos a
+necessidade de relacionarmos genes, transcritos, proteinas, metabólitos
+entre estas diferentes fontes. Vários pacotes no R tem a finalidade de
+fazer a interconversão entre os diferentes identificadores de bancos de
+dados. Neste exemplo, usaremos o pacote biomaRt
+a fim de buscarmos outras informações sobre os genes diferencialmente
+expressos.
Iremos buscar o HGNC Symbol para cada gene DGE:
+library(biomaRt)
+
+# Importar lista de genes diferencialmente expressos
+load("./dados_curso/outputs/exp_diff/res_dge.RData")
+genes_ensembl <- rownames(res_dge)
+
+# Use a função listMart() para ver quais canais estão disponíveis (use 'ENSEMBL_MART_ENSEMBL');
+listMarts()
## biomart version
+## 1 ENSEMBL_MART_ENSEMBL Ensembl Genes 107
+## 2 ENSEMBL_MART_MOUSE Mouse strains 107
+## 3 ENSEMBL_MART_SNP Ensembl Variation 107
+## 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 107
+# Use a função useMart() para setar o canal em um objeto (este servirá como conexão aos servidores do ENSEMBL)
+ensembl <- useMart("ENSEMBL_MART_ENSEMBL")
+
+# Com a função listDatasets(), verifique o identificador para o dataset de humanos;
+datasets <- listDatasets(ensembl)
+
+# Com a função useDataset(), estabeleça o organismo para a busca
+ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
+
+# Obtenha o atributo de interesse, ou seja, para qual(is) identificador(es) se deseja fazer o mapeamento. Use a função listAttributes().
+atributos <- listAttributes(ensembl)
+
+# Obtenha os filtros, ou seja, qual identificador está sendo usado para fazer o mapeamento. Neste caso, é o ensembl_gene_id.
+filtros <- listFilters(ensembl)
+
+# Use a função getBM() para obter o mapeamento. Podemos escolher quantos atributos (identificadores) desejarmos
+ids <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id", "ensembl_peptide_id", "entrezgene_id"),
+ filters = "ensembl_gene_id",
+ values = genes_ensembl,
+ mart = ensembl)
Quando fazemos a conversão de um banco para outro, é comum perdermos +alguma informação. Por exemplo, veja que para o gene “ENSG00000272410”, +não há um identificador HGNC Symbol correspondente. Vamos obter uma +lista contendo apenas HGNC Symbols e os ENTREZ ids válidos.
+# Selecionar hgnc symbol e remover observações faltantes
+ids$hgnc_symbol <- ifelse(ids$hgnc_symbol == "", NA, ids$hgnc_symbol)
+genes_hgnc <- unique(ids$hgnc_symbol)
+
+# Selecionar entrez ids e remover observações faltantes
+genes_entrez <- unique(ids$entrezgene_id)
+genes_entrez <- genes_entrez[!is.na(genes_entrez)]
O processo de montar uma rede de interação proteína-proteína consiste +de, partindo de uma lista de genes, por exemplo os genes +diferencialmente expressos que encontramos, buscar suas proteínas +correspondentes no banco de dados STRINGdb. Este banco cataloga e +classifica interações proteicas de acordo com diferentes metodologias, +desde ensaios bioquímicas à análises de co-ocorrência na literatura.
+Para obtermos as informações do STRING, iremos usar sua API. Esta +interface possui métodos para a obtenção dos dados disponíveis no banco +(Para saber mais sobre cada método, visite https://string-db.org/help/api/).
+O primeiro método a ser usado é get_string_ids, que +irá mapear uma lista de genes para os identificadores próprios do STRING +(Para saber sobre os parâmetros deste método, visite https://string-db.org/cgi/help.pl?subpage=api%23mapping-identifiers).
+# Para fazer este mapeamento, quando submetemos uma requisição à API do string programaticamente,
+# temos que concatenar os identificadores e separá-los com o símbolo "%0d".
+genes_hgnc_concatenado <- paste0(genes_hgnc, collapse = "%0d")
+
+# Fazer a solicitação a API do STRING
+req <- RCurl::postForm(
+ "https://string-db.org/api/tsv/get_string_ids",
+ identifiers = genes_hgnc_concatenado,
+ echo_query = "1",
+ species = "9606"
+)
+map_ids <- read.table(text = req, sep = "\t", header = T, quote = "")
Agora, de posse dos identificadores do STRINGdb, vamos usá-los para +obtermos a rede de interação entre estes genes. Usaremos o método network +para obter a rede de interação.
+#Concatenar os identificadores do string para a requisição
+genes_string_concatenado <- paste0(unique(map_ids$stringId), collapse = "%0d")
+
+# Requisição para o método 'network'
+req2 <- RCurl::postForm(
+ "https://string-db.org/api/tsv/network",
+ identifiers = genes_string_concatenado, # identificadores do stringID, obtidos na etapa anterior
+ required_core = "0", # score mínimo para cada interação
+ species = "9606" # espécie (H. sapiens)
+)
+int_network <- read.table(text = req2, sep = "\t", header = T)
+int_network <- unique(int_network)
O dataframe int_network
possui as informações sobre as
+interações entre as proteínas buscadas. A inferência das interações
+entre a proteína A e a proteína B é obtida de diferentes fontes. Da
+coluna nscore
à coluna tscore
estão as fontes
+e o escore para cada interação. Assim, temos:
stringId_A
: Identificador STRING para a proteína
+A;stringId_B
: Identificador STRING para a proteína
+B;preferredName_A
: Nome da proteína A;preferredName_B
: Nome da proteína B;ncbiTaxonId
: Identificador do NCBI para o taxon;score
: escore combinado;nscore
: escore para evidências de vizinhança
+gênica;fscore
: escore para fusão gênica;pscore
: escore para perfil filogenético;ascore
: escore para coexpressão;escore
: escore para evidências experimentais;dscore
: escore para evidências contidas em bancos de
+dados;tscore
: escore para evidências provenientes de
+textmining.A seguir temos uma função para selecionarmos as fontes de interação e +calcularmos o escore combinado entre elas:
+# Função para combinar os scores de acordo com o algoritmo usado pelo STRING
+combinescores <- function(dat, evidences = "all", confLevel = 0.4) {
+ if(evidences[1] == "all"){
+ edat<-dat[,-c(1,2,ncol(dat))]
+ } else {
+ if(!all(evidences%in%colnames(dat))){
+ stop("NOTE: one or more 'evidences' not listed in 'dat' colnames!")
+ }
+ edat<-dat[,evidences]
+ }
+ if (any(edat > 1)) {
+ edat <- edat/1000
+ }
+ edat<-1-edat
+ sc<- apply(X = edat, MARGIN = 1, FUN = function(x) 1-prod(x))
+ dat <- cbind(dat[,c(1,2)],combined_score = sc)
+ idx <- dat$combined_score >= confLevel
+ dat <-dat[idx,]
+ return(dat)
+}
Nossa rede de interação será baseada nas seguintes fontes: evidências +de co-expressão, evidências experimentais e evidências de bancos de +dados. Iremos combinar os escores destas 3 fontes e escolher as +interações com escore combinado maior que 0,9.
+# Escolher canais, combinar escore e filtrar escore combinado
+int_network <- combinescores(int_network, evidences = c("ascore", "escore", "dscore"), confLevel = 0.9)
A seguir, iremos preprocessar a rede criada e plotá-la com o pacote +RedeR:
+library(igraph)
+library(RedeR)
+
+# Remover o identificador de espécie em cada ENSP
+int_network$stringId_A <- substring(int_network$stringId_A, 6, 1000)
+int_network$stringId_B <- substring(int_network$stringId_B, 6, 1000)
+
+# Filtrar a rede, mantendo apenas os ENSP que estão presentes no nosso dataframe inicial
+idx1 <- int_network$stringId_A %in% ids$ensembl_peptide_id
+idx2 <- int_network$stringId_B %in% ids$ensembl_peptide_id
+int_network <- int_network[idx1 & idx2,]
+
+# Manter somente os nomes dos genes na tabela de interação
+int_network$gene_name_A <- ids$hgnc_symbol[match(int_network$stringId_A, ids$ensembl_peptide_id)]
+int_network$gene_name_B <- ids$hgnc_symbol[match(int_network$stringId_B, ids$ensembl_peptide_id)]
+int_network <- int_network[, c("gene_name_A", "gene_name_B", "combined_score")]
+
+# Criar objeto igraph. Redes PPI são não direcionadas, portanto 'directed = FALSE'.
+g <- graph_from_data_frame(int_network, directed = FALSE)
+
+# Montar dataframe de anotação para colorir os nós
+ann <- data.frame(ensembl_gene_id = rownames(res_dge),
+ log2FC = res_dge$log2FoldChange,
+ hgnc_symbol = ids$hgnc_symbol[match(rownames(res_dge), ids$ensembl_gene_id)])
+
+# Obter a tabela apenas para os genes contidos na rede
+ann <- ann[ann$hgnc_symbol %in% V(g)$name,]
+
+# Mapear dataframe de anotação para a rede
+g <- att.mapv(g, dat = ann, refcol = 3)
+
+# Atribuir valores de log2FC para as cores dos nós
+g <- att.setv(g, from = "log2FC", to = "nodeColor", breaks = seq(-3,3,0.4), pal = 2)
+
+# Definir outros atributos do grafo
+V(g)$nodeLineWidth <- 2
+V(g)$nodeFontSize <- 16
+
+# att.sete() imprime atributos validos do RedeR
+E(g)$edgeColor <- "grey50"
+E(g)$edgeWidth <- 2
# Abrir porta para o RedeR e adicionar o grafo
+rdp <- RedPort()
+calld(rdp)
+addGraph(rdp, g)
+
+# Adicionar legenda
+scl <- g$legNodeColor$scale
+leg <- g$legNodeColor$legend
+addLegend.color(rdp, colvec = scl, labvec = leg, title = "log2FC")
A análise de enriquecimento funcional tem como finalidade ajudar o +pesquisador a identificar vias biológicas chaves as quais podem estar +alteradas na condição de interesse. Tomando como base a nossa lista de +genes diferencialmente expressos nas amostras tumorais, a análise de +enriquecimento irá verificar se há um número muito maior de genes +alterados em uma determinada via biológica do que seria esperado ao +acaso.
+Com isso, é necessário tomar por base algum banco de dados que +contenha informações sobre vias biológicas, como o KEGG Pathways e o Gene Ontology.
+Para a análise de enriquecimento, iremos utilizar o pacote clusterProfiler
.
Primeiramente, faremos o enriquecimento funcional com base no KEGG
+Pathways. Utlizaremos a função enrichKEGG()
. Esta função
+requer como entrada uma lista de ENTREZ ids e o identificador do
+organismo no KEGG (ver este link).
+Escolheremos as vias significativas aquelas com qvalor < 0,05.
library(clusterProfiler)
+
+# Enriquecimento
+kegg_enrich <- enrichKEGG(
+ gene = genes_entrez,
+ organism = "hsa",
+ keyType = "ncbi-geneid",
+ pAdjustMethod = "BH",
+)
+
+# Obter dataframe com os resultados
+df_kegg <- kegg_enrich@result
+df_kegg <- df_kegg[df_kegg$qvalue < 0.05,]
Agora, façamos o enriquecimento funcional com base no Gene Ontology.
+A ontologia utilizada para o enriquecimento será a de processos
+biológicos (biological processes, argumento
+ont = "BP"
).
go_enrich <- enrichGO(
+ gene = genes_ensembl,
+ OrgDb = "org.Hs.eg.db",
+ ont = "BP",
+ pAdjustMethod = "BH",
+ keyType = 'ENSEMBL'
+)
+
+df_go <- go_enrich@result
+df_go <- df_go[df_go$qvalue <= 0.05,]
Além disso, o clusterProfiler
dispõe de alguns gráficos
+para visualização do enriquecimento.
# Barplot - KEGG
+barplot(kegg_enrich,
+ drop = TRUE,
+ showCategory = 10,
+ title = "KEGG Pathways")
# Barplot - GO
+barplot(go_enrich,
+ drop = TRUE,
+ showCategory = 10,
+ title = "GO Biological Pathways")
# Dotplot - KEGG
+dotplot(kegg_enrich)
# Dotplot - GO
+dotplot(go_enrich, showCategory = 20)
A análise com o KEGG resultou em apenas 4 vias enriquecidas: Cell +cycle, ECM-receptor interaction, Protein digestion and +absorption e Focal adhesion. Este enriquecimento sugere +que o tumor metastático possui alterações na interação do tumor com a +matriz extracelular e no ciclo celular. O resultado do Gene Ontology +apresenta uma quantidade maior de processos biológicos enriquecidos, +muitos deles relacionados também ao ciclo celular e à matriz +extracelular.
+É comum obtermos uma lista extensa de processos biológicos +enriquecidos. Para que possamos ter uma ideia mais precisa sobre quais +processos biológicos se sobressaem, agrupamos estes por similaridade +semântica no Gene Ontology e e verificamos quais tem uma proporção maior +de genes diferencialmente expressos.
+library(dplyr)
+library(clusterProfiler)
+library(TreeAndLeaf)
+library(RedeR)
+library(igraph)
+library(RColorBrewer)
+library(GOSemSim)
+library(classInt)
+library(org.Hs.eg.db)
+
+# Obter a lista de processos enriquecidos
+terms <- df_go$ID
+
+# Criar objeto para análise de similaridade de ontologias
+semData <- godata(ont = "BP")
+
+# Proporção de genes diferencialmente expressos em relação ao total de genes da via
+df_go$path_length <- as.integer(sapply(strsplit(df_go$GeneRatio, "/"), "[", 2))
+df_go$ratio <- df_go$Count / df_go$path_length
+
+# Calcular similaridade semântica entre os termos
+mSim <- mgoSim(terms, terms, semData = semData, measure = "Wang", combine = NULL)
+
+# Organizando tamanho dos nós
+size <- df_go$Count
+aux <- sort(unique(size))
+names(aux) <- as.character(1:length(aux))
+sizeRanks <- as.factor(names(aux[match(size, aux)]))
+sizeIntervals <- 7
+sizeMultiplier <- 5
+sizeBase <- 50
+df_go$size <- (sizeBase + (as.numeric(sizeRanks) * sizeMultiplier))
+
+# Clusterização dos termos
+hc <- hclust(dist(mSim), "average")
+
+# Criar objeto TreeAndLeaf
+tal <- treeAndLeaf(hc)
+
+# Setar atributos do grafo
+tal <- att.mapv(g = tal, dat = df_go, refcol = 1)
+pal <- brewer.pal(9, "OrRd")
+tal <- att.setv(g = tal, from = "Description", to = "nodeAlias")
+tal <- att.setv(
+ g = tal,
+ from = "ratio",
+ to = "nodeColor",
+ cols = pal,
+ nquant = 5
+)
+tal <- att.setv(
+ g = tal,
+ from = "size",
+ to = "nodeSize",
+ xlim = c((sizeBase + sizeMultiplier),
+ (
+ sizeBase + (sizeMultiplier * sizeIntervals)
+ ),
+ sizeMultiplier),
+ nquant = sizeIntervals
+)
+tal <-
+ att.addv(tal,
+ "nodeFontSize",
+ value = 15,
+ index = V(tal)$isLeaf)
+tal <- att.adde(tal, "edgeWidth", value = 3)
+tal <- att.adde(tal, "edgeColor", value = "gray80")
+tal <- att.addv(tal, "nodeLineColor", value = "gray80")
# Inicializar o RedeR
+rdp <- RedPort()
+calld(rdp)
+addGraph(obj = rdp, g = tal)
+
+addLegend.color(obj = rdp, tal, title = "Proportion",
+ position = "topright")
+addLegend.size(obj = rdp, tal, title = "Counts",
+ position = "bottomright")
As redes de co-expressão são redes biológicas que representam a +correlação da expressão de pares de genes. Diversos algoritmos podem ser +utilizados para calcular a co-expressão de uma lista de genes.
+O pacote do Bioconductor CEMiTool +(Co-Expression Modules identification Tool) fornece aos seus usuários +uma maneira simples de realizar análises de co-expressão, encontrando +módulos gênicos que representam a correlação da expressão entre os +genes.
+Além disso, o CEMiTool permite integrar outras informações, como vias +metabólicas do KEGG, para se realizar um enriquecimento funcional e +interações da rede proteína-proteína dos genes a serem analisados, +permitindo a construção de uma rede que ilustre os resultados da +co-expressão.
+Primeiro, iremos ler o dado de expressão dos genes diferencialmente +expressos e as informações da tabela de metadados, que inclui os grupos +ou contrastes.
+library(CEMiTool)
+library(DESeq2)
+library(dplyr)
+library(biomaRt)
+library(org.Hs.eg.db)
+
+# Carregar dado de expressão
+load("./dados_curso/outputs/exp_diff/dds.RData")
+
+# Carregar genes DGE
+load("./dados_curso/outputs/exp_diff/res_dge.RData")
+
+# Importar metadado
+meta <- read.csv("./dados_curso/inputs/metadata/metadata.csv")
+meta <- meta[, c("run", "subtype")]
+
+# Obter expressão normalizada e filtrar genes DGE
+vst <- vst(dds)
+counts <- assay(vst)
+counts <- as.data.frame(counts[rownames(counts) %in% rownames(res_dge),])
ensembl <- useMart("ENSEMBL_MART_ENSEMBL")
+ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
+ids <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id", "ensembl_peptide_id", "entrezgene_id"),
+ filters = "ensembl_gene_id",
+ values = rownames(counts),
+ mart = ensembl)
Para isso, usaremos o org.Hs.db, especificamente o org.Hs.egPATH2EG, +que possui uma lista de vias do KEGG e os Entrez gene IDs associados a +cada, que podemos então cruzar com nosso dado de expressão.
+kegg_to_entrez <- as.data.frame(org.Hs.egPATH2EG)
+
+ids$entrezgene_id <- as.character(ids$entrezgene_id)
+kegg_to_symbol <- merge(ids, kegg_to_entrez, by.x = "entrezgene_id", by.y = "gene_id")
+kegg_to_symbol <- unique(kegg_to_symbol)
+
+save(kegg_to_symbol, file="./dados_curso/outputs/cea/pathways_for_cea.RData")
# Para fazer este mapeamento, quando submetemos uma requisição à API do string programaticamente,
+# temos que concatenar os identificadores e separá-los com o símbolo "%0d".
+genes_hgnc_concatenado <- paste0(ids$hgnc_symbol, collapse = "%0d")
+
+# Fazer a solicitação a API do STRING
+req <- RCurl::postForm(
+ "https://string-db.org/api/tsv/get_string_ids",
+ identifiers = genes_hgnc_concatenado,
+ echo_query = "1",
+ species = "9606"
+)
+map_ids <- read.table(text = req, sep = "\t", header = T, quote = "")
+
+# Função para combinar os scores de acordo com o algoritmo usado pelo STRING
+combinescores <- function(dat, evidences = "all", confLevel = 0.4) {
+ if(evidences[1] == "all"){
+ edat<-dat[,-c(1,2,ncol(dat))]
+ } else {
+ if(!all(evidences%in%colnames(dat))){
+ stop("NOTE: one or more 'evidences' not listed in 'dat' colnames!")
+ }
+ edat<-dat[,evidences]
+ }
+ if (any(edat > 1)) {
+ edat <- edat/1000
+ }
+ edat<-1-edat
+ sc<- apply(X = edat, MARGIN = 1, FUN = function(x) 1-prod(x))
+ dat <- cbind(dat[,c(1,2)],combined_score = sc)
+ idx <- dat$combined_score >= confLevel
+ dat <-dat[idx,]
+ return(dat)
+}
+
+# Concatenar os identificadores do string para a requisição
+genes_string_concatenado <- paste0(unique(map_ids$stringId), collapse = "%0d")
+
+# Requisição para o método 'network'
+req2 <- RCurl::postForm(
+ "https://string-db.org/api/tsv/network",
+ identifiers = genes_string_concatenado, # identificadores do stringID, obtidos na etapa anterior
+ required_core = "0", # score mínimo para cada interação
+ species = "9606" # espécie (H. sapiens)
+)
+int_network <- read.table(text = req2, sep = "\t", header = T)
+int_network <- unique(int_network)
+int_network <- combinescores(int_network, evidences = c("ascore", "escore", "dscore"), confLevel = 0.9)
+
+# Remover o identificador de espécie em cada ENSP
+int_network$stringId_A <- substring(int_network$stringId_A, 6, 1000)
+int_network$stringId_B <- substring(int_network$stringId_B, 6, 1000)
+
+# Filtrar a rede, mantendo apenas os ENSP que estão presentes no nosso dataframe inicial
+idx1 <- int_network$stringId_A %in% ids$ensembl_peptide_id
+idx2 <- int_network$stringId_B %in% ids$ensembl_peptide_id
+int_network <- int_network[idx1 & idx2,]
+
+# Manter somente os nomes dos genes na tabela de interação
+int_network$ENSG_A <- ids$ensembl_gene_id[match(int_network$stringId_A, ids$ensembl_peptide_id)]
+int_network$ENSG_B <- ids$ensembl_gene_id[match(int_network$stringId_B, ids$ensembl_peptide_id)]
+int_network <- int_network[, c("ENSG_A", "ENSG_B")]
+
+save(int_network, file="./dados_curso/outputs/cea/network_for_cea.RData")
A função cemitool()
é a função mestra do pacote, é por
+meio dela que iremos obter o objeto com os módulos de co-expressão e as
+figuras. Para isso, adicionaremos, é claro, nosso dado de expressão, mas
+também nossos contrastes - especificando em qual coluna está o nome de
+cada amostra (Run) e o nome dos contrastes (subtype) -, nossas vias
+metabólicas, e nossas interações.
load("./dados_curso/outputs/cea/network_for_cea.RData")
+load("./dados_curso/outputs/cea/pathways_for_cea.RData")
+paths <- unique(kegg_to_symbol[, c("path_id", "ensembl_gene_id")])
+paths <- setNames(paths, c("term", "gene"))
+
+cem <-
+ cemitool(
+ counts,
+ annot = meta,
+ sample_name_column = "run",
+ class_column = "subtype",
+ paths,
+ interactions = int_network,
+ filter = TRUE,
+ plot = TRUE,
+ verbose = TRUE
+ )
## pickSoftThreshold: will use block size 135.
+## pickSoftThreshold: calculating connectivity for given powers...
+## ..working on genes 1 through 135 of 135
+## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k. Density
+## 1 1 0.828 1.330 0.808 57.600 60.600 80.60 0.43000
+## 2 2 0.195 0.213 0.172 31.900 31.800 55.90 0.23800
+## 3 3 0.506 -0.257 0.721 20.000 18.300 42.30 0.14900
+## 4 4 0.602 -0.607 0.618 13.500 11.200 33.80 0.10100
+## 5 5 0.743 -0.814 0.771 9.620 7.240 28.00 0.07180
+## 6 6 0.818 -0.882 0.818 7.140 4.710 23.70 0.05330
+## 7 7 0.846 -0.972 0.816 5.480 3.250 20.50 0.04090
+## 8 8 0.898 -0.993 0.875 4.310 2.330 17.90 0.03220
+## 9 9 0.821 -1.070 0.777 3.470 1.660 15.90 0.02590
+## 10 10 0.841 -1.100 0.798 2.850 1.230 14.20 0.02120
+## 11 12 0.869 -1.100 0.832 2.000 0.694 11.50 0.01490
+## 12 14 0.860 -1.120 0.821 1.470 0.397 9.59 0.01100
+## 13 16 0.837 -1.100 0.791 1.120 0.258 8.09 0.00839
+## 14 18 0.885 -1.130 0.853 0.884 0.162 7.11 0.00660
+## 15 20 0.942 -1.090 0.927 0.711 0.108 6.32 0.00531
+## Centralization Heterogeneity
+## 1 0.1750 0.261
+## 2 0.1820 0.423
+## 3 0.1690 0.552
+## 4 0.1540 0.665
+## 5 0.1390 0.769
+## 6 0.1260 0.865
+## 7 0.1140 0.955
+## 8 0.1030 1.040
+## 9 0.0938 1.120
+## 10 0.0857 1.200
+## 11 0.0722 1.340
+## 12 0.0615 1.470
+## 13 0.0528 1.600
+## 14 0.0472 1.710
+## 15 0.0425 1.820
+## ..connectivity..
+## ..matrix multiplication (system BLAS)..
+## ..normalization..
+## ..done.
+## ..cutHeight not given, setting it to 0.991 ===> 99% of the (truncated) height range in dendro.
+## ..done.
+## mergeCloseModules: Merging modules whose distance is less than 0.2
+## mergeCloseModules: less than two proper modules.
+## ..color levels are 0, 1
+## ..there is nothing to merge.
+## Calculating new MEs...
+Primeiro, podemos ver quantos módulos foram encontrados com a função
+nmodules()
, e também alguns dos genes que compoem cada
+módulo através da função module_genes()
.
nmodules(cem)
## [1] 2
+head(module_genes(cem))
## genes modules
+## 1 ENSG00000198786 M1
+## 2 ENSG00000168542 M1
+## 3 ENSG00000164692 M1
+## 4 ENSG00000103811 M1
+## 5 ENSG00000198695 M1
+## 6 ENSG00000212907 M1
+Podemos também ver quais genes possuem a maior conectividade.
+n <- 10
+hubs <- get_hubs(cem, n)
+hubs
## $M1
+## ENSG00000148848 ENSG00000151388 ENSG00000204262 ENSG00000168542 ENSG00000099994
+## 23.71241 21.85683 21.34914 20.49144 20.11550
+## ENSG00000133110 ENSG00000171848 ENSG00000122861 ENSG00000137573 ENSG00000164692
+## 20.09753 19.09620 18.87951 18.45229 18.09585
+##
+## $Not.Correlated
+## ENSG00000067048 ENSG00000129824 ENSG00000114374 ENSG00000275302 ENSG00000125538
+## 1.65134593 1.61257353 1.57913458 0.28989405 0.28487654
+## ENSG00000171346 ENSG00000215030 ENSG00000275385 ENSG00000184564 ENSG00000171557
+## 0.10376992 0.09397441 0.07870171 0.04399561 0.01522507
+Como fornecemos as informações dos contrastes, a função cemitool +também realizou um gene set enrichment analysis, usando a +função do pacote fgsea. A figura associada a esse resultado +possui um círculo em cada módulo e em cada contraste, com a intensidade +e o tamanho de cada círculo correspondendo ao Score de enriquecimento +normalizado (NES), que é o score de enriquecimento para o módulo em cada +classe normalizado pelo número de genes do módulo.
+# generate heatmap of gene set enrichment analysis
+show_plot(cem, "gsea")
## $enrichment_plot
+
+Também podemos visualizar o padrão de expressão de cada gene em cada +módulo.
+# plot gene expression within each module
+plots <- show_plot(cem, "profile")
+plots[1]
## $M1
+
+A função cemitool também determina se as vias metabólicas estão +associadas com os módulos através de uma over-representation +analysis.
+plots <- show_plot(cem, "ora")
+plots[1]
## $M1
+## $M1$pl
+
+##
+## $M1$numsig
+## [1] 2
+E, dado que fornecemos a informação das interações entre nossos +genes, também podemos gerar a imagem de uma rede que combina a +informação das interações junto aos resultados da análise de +co-expressão.
+plots <- show_plot(cem, "interaction")
+plots[1]
## $M1
+
+Todos os resultados encontrados pelo pacote CEMiTool na nossa +análise, além de todos os parâmetros que utilizamos, podem ser +facilmente salvos e compartilhados através de um relatório HTML, que se +gera com a função generate_report(). Tente ir ao diretório escolhido +abaixo e abrir o arquivo HTML com seu navegador!
+# Cria relatorio do CEMiTool em HTML no diretório abaixo
+# force = TRUE sobrescreve o reporte previamente criado
+generate_report(cem, directory = "./dados_curso/outputs/cea/cemitool_report/", force = TRUE)
o Transcriptograma é uma metodologia de Biologia de Sistemas que tem +como objetivo identificar grupos de genes alterados em conjunto. Nesta +metodologia, se projeta a expressão gênica sobre uma lista de proteínas +ordenadas de tal forma que a probabilidade de que duas delas participem +de uma mesma via biológica decresce com o quadrado da distância entre +elas.
+Esta metodologia está implementada no pacote transcriptogramer
.
+O pacote fornece diferentes ordenamentos de proteínas de humanos com
+base no interatoma do STRING.
Primeiramente, com o pacote biomaRt
, vamos obter um
+dicionário relacionando os Ensembl Gene ID com os STRING IDS.
library(biomaRt)
+library(DESeq2)
+library(edgeR)
+
+# Obter tabela de expressão normalizada
+load("./dados_curso/outputs/tximport/txi_gene.RData")
+
+# Filtrar genes com baixos valores de contagens
+counts <- txi_gene$counts
+cpm <- cpm(counts)
+keep <- rowSums(cpm > 5) >= 20
+counts <- counts[keep,]
+
+# LogCPM
+logCPM <- cpm(counts, log = T)
+
+# Obter dicionário
+ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
+dict <- getBM(attributes = c("ensembl_peptide_id", "ensembl_gene_id"), mart = ensembl)
+dict$ensembl_peptide_id <- ifelse(dict$ensembl_peptide_id == "", NA, dict$ensembl_peptide_id)
+dict <- dict[complete.cases(dict),]
+dict$ensembl_peptide_id <- paste("9606.", dict$ensembl_peptide_id, sep = "")
+
+# Organizar dicionário para que a primeira coluna corresponda aos string ids e
+# a segunda coluna corresponda aos ensembl gene ids
+names(dict) <- c("Protein", "Probe")
+
+# Selecionar ocorrências presentes na tabela de expressão
+logCPM <- logCPM[rownames(logCPM) %in% dict$Probe,]
library(transcriptogramer)
+
+# Obter metadado
+meta <- read.csv("./dados_curso/inputs/metadata/metadata.csv")
# Criar objeto do transcriptogramer com raio 80
+t <- transcriptogramPreprocess(association = association, ordering = Hs900, radius = 80)
+
+# Rodar o passo 1 do transcriptogramer
+t <- transcriptogramStep1(object = t, expression = logCPM, dictionary = dict, nCores = 2)
+
+# Rodar o passo 2 do transcriptogramer
+t <- transcriptogramStep2(object = t, nCores = 2)
+
+save(t, file = "./dados_curso/outputs/transcriptogramer/t.RData")
# Rodar a expressão diferencial entre os grupos com base ao raio = 80
+# TRUE = amostras controle (indolent, tumor não-invasivo)
+# FALSE = amostras tratamento (invasive, tumor metastático)
+levels <- ifelse(meta$subtype == "indolent", T, F)
+
+t <- differentiallyExpressed(object = t, levels = levels, pValue = 0.05,
+ species = "Homo sapiens", title = "radius 80")
save(t, file = "./dados_curso/outputs/transcriptogramer/t_diff.RData")
# Enriquecimento dos clusters
+t <- clusterEnrichment(object = t, species = HsBPTerms,
+ pValue = 0.05, nCores = 3)
+
+save(t, file = "./dados_curso/outputs/transcriptogramer/t_enrich.RData")
# Termos enriquecidos por cluster identificado
+termos <- Terms(t)
Prof. Dalmolin supervises in the post-graduate program in Bioinformatics at the Federal University of Rio Grande do Norte, which offers both Master and PhD courses. Come visit us or, better yet, join us!
+ + +Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. +Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat.
+ + +The network elaborated in our study contemplates protein classifications according to the AP detection algorithm,facilitating the visual identification of the most important proteins of a network.
+ + + + + + + + + + + +Transcriptograms are genome wide gene expression profiles that provide a global view for the cellular metabolism, while indicating gene sets whose expressions are altered.
+ + + + + + + + + + + +Geneplast is designed for evolutionary and plasticity analysis based on orthologous groups distribution in a given species tree.
+ + + + + + + + + + + +A pipeline for sensitive taxonomic classification and flexible functional annotation of shotgun metagenomic data.
+ + + + + + + + + + + +Dalmolin Group’s workflow for pre-processing, alignment and quantification of bulk RNA-seq data.
+ + + + + + + + + + + +