Tutorial - PCA

Introdução

Esse documento serve como estudo próprio e também como uma rápida demonstração da Análise de Componentes Principais a partir de dados genéticos. Fui altamente inspirado pelo Gábor Mészáros do Genomics Boot Camp no youtube. Passamos brevemente pela obtenção dos dados, pelo controle de qualidade e finalizamos com o PCA.

Obtenção dos dados

Um lugar onde podemos obter datasets disponíveis publicamente é o datadryad.org
Esse site possui diversos datasets, em seus diferentes formatos, prontos para serem utilizados.
Para esse tutorial, escolhi utilizar os dados do artigo Aedes aegypti in North America (Microsatellite and SNP array)
Basta clicarmos no link de download e obtemos o dataset em formato plink.

Que tal darmos uma olhadinha nesses dados, só para irmos nos acostumando?

O arquivo .bim é aquele que contém informações sobre os SNPs presentes no dado. Cada linha corresponde a um SNP diferente.

bim <- read.table("./data/doi_10.5061_dryad.5x69p8d5j__v2/DatasetS2_NorthAmericaCombinedSnps_reorder.bim")
head(bim)

##     V1          V2 V3     V4 V5 V6
## 1 3003 AX-93233745  0  44357  T  G
## 2 3003 AX-93233746  0  44795  A  G
## 3 3003 AX-93233749  0  56930  G  A
## 4 3003 AX-93233752  0  66000  G  T
## 5 3003 AX-93233758  0 118514  A  T
## 6 3003 AX-93233759  0 131149  T  C

nrow(bim)

## [1] 20003

O arquivo .fam é aquele que contém informações sobre os indivíduos. Cada linha contém um indivíduo diferente.

fam <- read.table("./data/doi_10.5061_dryad.5x69p8d5j__v2/DatasetS2_NorthAmericaCombinedSnps_reorder.fam")
head(fam)

##   V1 V2 V3 V4 V5 V6
## 1 SM  3  0  0  0 -9
## 2 SM  8  0  0  0 -9
## 3 SM 11  0  0  0 -9
## 4 SM 12  0  0  0 -9
## 5 SM 13  0  0  0 -9
## 6 SM 14  0  0  0 -9

nrow(fam)

## [1] 381

O arquivo .bed é um arquivo binário contendo a informação genética desses indivíduos.

Estrutura dos Dados

Essa olhada geral nos dados nos revelou que estamos lidando com um dataset com 20003 SNPs e 381 indivíduos. Além disso, os nomes dos cromossomos estão desformatados (não vai impedir em nada nossa análise), e não temos informações de parentais nem de sexo dos indivíduos.

Controle de Qualidade

Após baixarmos os dados, vamos passá-los pelo controle de qualidade.
A etapa de controle de qualidade envolve passarmos os dados por uma série de filtros comumente utilizados para se assegurar uma alta taxa de genotipagem, ou seja, um dado de alta qualidade. Isso é feito no próprio plink, e os filtros são:

Missingness por SNP, ou seja, SNPs que possuem uma taxa de missing data maior que 10% serão removidos: –geno 0.1
Missingness por Indivíduo, ou seja, indivíduos com mais que 10% de missing data também serão removidos: –mind 0.1
Minor Allele Frequency, filtro para alelos raros: –maf 0.05

system("./programs/plink.exe --bfile ./data/doi_10.5061_dryad.5x69p8d5j__v2/DatasetS2_NorthAmericaCombinedSnps_reorder --geno 0.1 --mind 0.1 --maf 0.05 --allow-extra-chr --make-bed --out ./data/aedes_QC")

## [1] 0

Aqui vemos que após os filtros de qualidade, de 20003 SNPs e 381 indivíduos, com taxa de genotipagem de 0.958638, ficamos com 19806 SNPs e 40 indivíduos removidos. Além disso, a taxa de genotipagem subiu para 0.987982.

Calculando as distâncias genéticas

Após passado pelo controle de qualidade, o dado precisa ser tratado de algumas outras formas antes de colocarmos no gráfico do PCA. Para construir esse gráfico, nós vamos precisar calcular as distâncias genéticas entre cada indivíduo. Isso é feito no plink, com o comando:

# Calculating genetic distances
system("./programs/plink.exe --bfile ./data/aedes_QC --distance-matrix --allow-extra-chr --allow-no-sex --nonfounders --out ./data/dist_ma_aedes")

## [1] 0

Classic Multidimensional Scaling

Nessa etapa, com a nossa matriz de distância entre os indivíduos em mãos, vamos trabalhar em cima desses dados.
Por enquanto, vamos carregar a matriz de distância:

dist_ma <- read.table("./data/dist_ma_aedes.mdist", header = F)

Aqui estamos criando dois data frames, que serão utilizados posteriormente

# Family ID Column of .fam file
family <- tibble(family = read.table("./data/dist_ma_aedes.mdist.id")[,1])

# Individual ID column of .fam file
ind <- tibble(IID = read.table("./data/dist_ma_aedes.mdist")[,2])

Aplicamos a função cmdscale() na nossa matriz de distância genética. Essa função representa o multidimensional scaling clássico de uma matriz de dados.

mds <- cmdscale(dist_ma, eig = T, 5)

Finalizamos essa seção computando os autovetores e os autovalores. Estamos interessados apenas nos dois primeiros componentes principais:

eigenvec <- cbind(family, ind, mds$points)

eigen_percent <- round(((mds$eig)/sum(mds$eig)) * 100, 2)

Ajustando as diferentes populações nas regiões que ocorrem

Queremos fazer com que o gráfico do PCA seja útil para nós vermos se as populações irão agrupar de acordo com as regiões em que elas ocorrem. Para isso, eu tive que coletar, no meta-dado, a tabela que conecta as abreviações com as regiões.
Essa parte foi um pouco trabalhosa. Consegui essa tabela neste link.
Com a tabela em mãos, tive que arranjar alguma forma de tirar ela do formato .docx. Nesse caso, eu copiei e colei num bloco de notas, e depois salvei em .tsv.
Vamos trabalhar em cima do dado agora:

aedes_pops <- read_tsv("./data/aedes_pops.tsv")

## Rows: 70 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): Pop. number, Population, Abbrev., Region name
## dbl (3): Year, N1, N2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# 77 abbrv
aedes_tidy <- aedes_pops |>
  select(Abbrev., `Region name`) |> 
  rename("abbrv" = Abbrev.,
         "region" = `Region name`) |> 
  separate(abbrv, into = c('a', 'b'), sep = " and ") |> 
  mutate(rows = row_number()) |>
  pivot_longer(c(a, b), names_to = "key", values_to = "values") |> 
  drop_na() |> 
  select(values, region)

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 63 rows [1, 2, 3,
## 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].

eigenvec <- cbind(family, ind, mds$points) |>
  left_join(aedes_tidy, by = c("family" = "values")) 

eigenvec <- eigenvec |>
  mutate(region = case_when(
    family == "MAD13" ~ "Northern CA",
    family == "MAD15" ~ "Northern CA",
    family == "Alamagordo" ~ "Central",
    family == "SanDiego" ~ "Southern CA",
    family == "Tucs" ~ "Southwest",
    family == "LC" ~ "Central",
    family == "Lubbock" ~ "Central",
    family == "Bexar" ~ "Central",
    family == "NO11" ~ "Southeast",
    family == "Mia" ~ "Southeast",
    TRUE ~ as.character(region)))

eigen_percent <- round(((mds$eig)/sum(mds$eig)) * 100, 2)

Principal Component Analysis (PCA)

Finalmente, com os autovalores e autovetores calculados, vamos plotar nossos indivíduos no gráfico usando ggplot

eigenvec |>
  ggplot() +
  geom_point(aes(`1`, `2`, color = region), size = 1.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    x = paste0("Componente Principal 1 (", eigen_percent[1], "%)"),
    y = paste0("Componente Principal 2 (", eigen_percent[2], "%)")) +
  theme_bw()