Operadores

• Aritméticos

# 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

• Relacionais

# 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

Variáveis

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

Funções

Chamando funções:

sum(1, 3, 5)

## [1] 9

a <- rep("Aluno", times = 3)
a

## [1] "Aluno" "Aluno" "Aluno"

Acessando a documentação

Estas funções buscam e exibem a documentação de funções:

help(sum)
?sd
??plot

Diretório de trabalho

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()

Salvando e carregando objetos

# 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

Vetores

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"

Hierarquia de classes

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 1², 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 ou impossíveis

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

Atributos de objetos

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"

Fator

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

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

Arrays

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

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

Trabalhando com listas

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

Data.frames

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"

Trabalhando com data.frames

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

Importação de valores de contagens de genes

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.

Importação dos valores de contagens de 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")

Importação de valores de contagens de transcritos normalizados

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")

Análise de Componentes Principais (PCA)

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 $n$ 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.

Clusterização hierárquica

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.

Representações gráficas do resultado da expressão diferencial

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.

Correção para DTE

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")

Correção para DTU

# 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")

CEMiTool

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.

Ler dados de 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),])

Traduzir identificadores

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)

Adquirir vias do KEGG

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")

Adquirir interações PPI

# 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")

Rodar o CEMiTool

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...

Investigação inicial do resultado

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

Enriquecimento dos módulos

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

Visualizar padrões de expressão nos módulos

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

Over-representation analysis

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

Interações

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

Gerando relatório

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)