Analisamos amostras de Leucemia Mielóide do projeto ICGC Data Portal TCGA-LAML obtidas do banco de dados ICGC Data Potal https://dcc.icgc.org/. O objetivo foi processar e integrar a análise de “Copy Number Variations” e Sequênciamento de RNA (expressão de genes). Este workflow computacional criar novas funções, bem como utilizar diferentes ferramentas para trabalhar com dados de projetos TCGA e concatena os algoritimos computacionais em um pacote denominado CNVcExp.
Bibliotecas necessárias:
-dplyr
-biomaRt
-org.Hs.eg.db
-pbapply
cnvs <- read.delim("C:/Users/bruno/Desktop/working/LMA-CNV/Data/icgc/tcga lma us/copy_number_somatic_mutation.LAML-US.tsv")
exp_seq <- read.delim("C:/Users/bruno/Desktop/working/LMA-CNV/Data/icgc/tcga lma us/exp_seq.LAML-US.tsv")Inspeção Visual dos dataframes
cnvs[1:5, c(6,10,12,13,14)]## submitted_sample_id segment_mean chromosome chromosome_start
## 1 TCGA-AB-2905-03A-01D-0756-21 0.0100 7 143711169
## 2 TCGA-AB-2905-11A-01D-0756-21 1.1718 4 172374884
## 3 TCGA-AB-2905-11A-01D-0756-21 0.0053 X 2703633
## 4 TCGA-AB-2905-03A-01D-0756-21 -3.5068 7 22435472
## 5 TCGA-AB-2905-11A-01D-0756-21 0.0006 21 23673883
## chromosome_end
## 1 159127004
## 2 172377776
## 3 27579602
## 4 22436487
## 5 46891566 exp_seq[1:5,c(5,8,9,10,12)]## submitted_sample_id gene_id normalized_read_count raw_read_count
## 1 TCGA-AB-2977-03B-01T-0760-13 ASB5 5.039169e-08 3
## 2 TCGA-AB-2977-03B-01T-0760-13 ASB2 2.131995e-07 9
## 3 TCGA-AB-2977-03B-01T-0760-13 ASB3 3.984592e-05 1022
## 4 TCGA-AB-2977-03B-01T-0760-13 ASB4 0.000000e+00 0
## 5 TCGA-AB-2977-03B-01T-0760-13 ASB6 1.252928e-05 1206
## assembly_version
## 1 GRCh37
## 2 GRCh37
## 3 GRCh37
## 4 GRCh37
## 5 GRCh37Precisamos modificar os valores da coluna “submitted_sample_id” para
poder acessar amostras de uma forma que possamos integrar. Precisamos de
um ID de amostra genêrico e menos específico. Criamos uma função interna
chamada “mudar_identificador_tcga()”.
2.1 Substituir “-” por “_” nos identificadores tcga
Aqui precisamos substituir o separador interno dos valores da coluna “submitted_sample_id” para um separador mais seguro ( “_” ) invés de “-”.
2.1.1 Samples passa a ser a coluna “submitted_sample_id” do data frame cnvs
samples <- cnvs[["submitted_sample_id"]]2.1.2 Utilizamos a função gsub() para modificar “-” para “_” (um caracter seguro em R).
samples <- gsub(x = samples[1:100], "-", "_")2.2 Modificar os identificadores (ids): função criada para essa finalidade
Função:
mudar_identificador_tcga <-function(x, sep){
names=NULL
for (i in x) {
identificador<-i
partes <- strsplit(identificador, sep)[[1]]
novo_identificador <- paste(partes[1:4], collapse = sep)
names=c(names, novo_identificador)
}
return(names)
}Modificando o identificador a partir da coluna
"submitted_sample_id"
Modificamos apenas as 100 primeiras linhas
samples <- mudar_identificador_tcga(x = samples[1:100], sep = "_")Remover o sufixo A ou B: esse sufixo possue um significado de acordo ao TCGA Barcode.
A numeração ao final indica se é uma amostra de tumor ou de tecido normal. Utilizaremos mais a frente como referência e seleção de amostras.
samples <- gsub("3A", "3", x = samples)
samples <- gsub("11A", "11", x = samples)Adicionar samples ids ao objeto cnv
Estamos utilizando apenas as 100 primeiras linhas apenas para critério de demonstração.
cnvs <- cbind(samples, cnvs[1:100,])Uma boa e correta prática é verificar se cada novo identificador criado “samples” corresponde ao identificador na coluna “submitted_sample_id”.
Visualizar como ficou:
samples[1:4]## [1] "TCGA_AB_2905_03" "TCGA_AB_2905_11" "TCGA_AB_2905_11" "TCGA_AB_2905_03"Como era:
cnvs$submitted_sample_id[1:2]## [1] "TCGA-AB-2905-03A-01D-0756-21" "TCGA-AB-2905-11A-01D-0756-21"A Tabela de cnvs com a adição da nova coluna:
Apenas as primeiras 5 linhas e as 5 primeiras colunas
cnvs[1:5,1:5]## samples icgc_donor_id project_code icgc_specimen_id icgc_sample_id
## 1 TCGA_AB_2905_03 DO21295 LAML-US SP44989 SA496570
## 2 TCGA_AB_2905_11 DO21295 LAML-US SP44995 SA496583
## 3 TCGA_AB_2905_11 DO21295 LAML-US SP44995 SA496583
## 4 TCGA_AB_2905_03 DO21295 LAML-US SP44989 SA496570
## 5 TCGA_AB_2905_11 DO21295 LAML-US SP44995 SA496583Queremos estudar e analisar apenas as variantes de cópias (cnvs) relacionadas ao cancer. Além disso, nesse projeto, não foi realizado o sequenciamento de rna para as amostras de tecido normal, dessa forma, não será possível integrar variantes de cópias de tecidos normais.
keep = grep("-03A-", cnvs$submitted_sample_id)
keep = cnvs[keep,]
cnvs=keepAcima, mantemos apenas as amostras de tumor através da interpretação do padrão doTCGA Barcode, no qual nesse estudo as amostras tumorais são identificadas pelo sufixo “-03A-”. Verifique se a tabela “cnv” contém algum identificador com o sufixo “11A” ou “11B”, o qual nesse projeto é utilizado para identificar as amostras de tecido normal.
Selecionar colunas de interesse e renomea-las
cnvs <- dplyr::select(.data = cnvs, "samples", "icgc_donor_id", "chromosome", "chromosome_start","chromosome_end","segment_mean")renomear as colunas
colnames(cnvs) <- c("samples", "donor_id", "chr", "start", "end", "segment_mean")corrigir rownames
rownames(cnvs) <- NULLUtilizamos a função para modificar a nomenclatura da nomeação dos cromossomos e também para remover os cromossomos sexuais, as quais não desejamos em nossa análise.
cnvs <- cnv_chr_change(df = cnvs, col = "chr", remove_X = TRUE, remove_Y = TRUE)Isso é util a depender do genomas de referência que irá utilizar mais a frente, a depender da referência a nomenclatura dos cromossomos estará de uma forma ou de outra.
4.1 Podemos utilizar gsub() para remover o caracter “chr” da nomeação dos cromossomos.
cnvs$seqnames <- gsub(x = cnvs$seqnames, "chr", "")Utilizamos os valores em segment_mean para através de cut-offs definir ganho ou perda de copias. Para ganho/ duplicação utilizamos o cut-off de 0.2, para perda de cópias ou deleções utilizamos o cut-off de -0.2. Esses valores são os indicados para dados do tcga.
Função criada: cnv_change_df()
cnv_change_df <- function(df, col, upper= 0.5, lower=-0.1){
type <- NULL
for(i in df[[col]]){
if(i <= lower){CN <- "Del" }
else{
if (i >= upper) { CN <- "Dup" }
else{ CN <- "Neutral"}
}
type=c(type, CN)}
df <- cbind(df, type)
return(df)
}Utilizando a função acima.
cnv_new <- cnv_change_df(df = cnvs, "segment_mean", upper = 0.2, lower = -0.2)Utilizamos um genoma de referência para anotar genes que se encontram nas regiões genomicas as quais identificamos cnvs (ganhos ou perdas de segmentos de DNA).
O genoma de referêmcia “GRCh37” foi utilizado, de acordo ao genoma de referência utilizado para mapear as amostras deste estudo. Os genomas de referência podem ser obtidos no banco de dados UCSC genome.
cnvs_ann <- cnv_ann_par(input = cnvs, gtf_ref="Homo_sapiens.GRCh37.87.gtf", specie="Homo sapiens" )Foi adicionada a coluna gene_ids ao dataframe com os cnvs. Essa coluna contêm os ENSEMBL ids dos genes sobrepostos a cada região em formato de lista para cada grange. Algumas regiões não possuem genes, o que gera uma coluna vazia.
6.1 Remover linhas não identificadas (que não sobrepõem genes)
cnvs_ann_genes=NULL
for(i in 1:10){
if (S4Vectors::isEmpty(cnvs$gene_ids[i])==FALSE){
y=cnv[i,]
cnvs_ann_genes=rbind(cnvs_ann_genes,y)
}
}6.2 Expandir dataframe: separar genes em linhas individuais
# Use a função unnest para expandir a coluna gene_ids
cnvs_ann_genes_new <- cnvs_ann_genes %>%
tidyr::unnest(cols = gene_ids) %>%
filter(gene_ids != "") # Filtras as colunas vazias7.1 Conectar ao database biomaRt
library(biomaRt)
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")Objeto com os identificadores de genes
gene_ids <- cnvs_ann_genes_new$gene_ids7.2 Obter os simbolos genicos, ensemble e entrezid
gene_symbols <- getBM(attributes = c("ensembl_gene_id", "external_gene_name", "entrezgene_id"), filters = "ensembl_gene_id", values = gene_ids, mart = ensembl, filter)7.3 Remover duplicatas
colnames(gene_symbols) <- c("gene_ids", "symbols", "entrezgene_id")
gene_symbols <- gene_symbols %>% distinct(gene_ids, .keep_all = TRUE)
#gene_symbols <- gene_symbols[!duplicated(gene_symbols$gene_ids), ]7.4 Juntar ao dataframe originar com left join
Esta função manterá as colunas que não tiveram simbolos gênicos anotados.
cnvs_ann_genes_new <- cnvs_ann_genes_new %>%
left_join(gene_symbols, by = c("gene_ids" = "gene_ids"))Os dados com as informações que desejamos estão no objeto cnvs_ann_genes_new.
Observação: se o conjunto de dados for muito grande/ pesado, é possível realizar os processos de anotação de genes dividindo o dataframe. Além de outras maneiras computacionalmente mais eficazes através das utilização de funções como apply, lapply, pbapply, dentre outras. Basta modificar as funções acima para que elas aceitem apenas um valor como entrada. Temos essas outras maneiras de fazer isso, mas não está disponibilizados nessa demonstração.
Alternativa 1: Dividir os data frames e realizar a anotação em cada parte individualmente.
Esta função divide os dataframes de acordo a quantidade de partes que se definir. Como exemplo estaremos dividindo-o em 100 partes iguais. Cada parte ficará disponível no ambiente global da sessão a ser utilizada.
Função 1: percentiles_df_divide().
Divide o data frame.
percentiles_df_divide <- function(df, n_parts){
n_rows <- nrow(df)
perc <- seq(0, 100, length.out = n_parts + 1)
parts <- vector("list", n_parts)
for (i in 1:n_parts){
start_row <- floor(perc[i] * n_rows / 100) + 1
end_row <- floor(perc[i + 1] * n_rows / 100)
parts[[i]] <- df[start_row:end_row, ]
}
return(parts)
}Função 2: cnvr_to_genes_divided().
Realiza a anotação de genes que se sobrepõem as regiões de cnvs.
cnv_genes <- cnvr_to_genes_divided(df = cnv, n=100)Para esse passo acima, também se aplica o que está disposto na explicação abaixo.
Função 3: cnv_individual_cromossomes().
Optamos por utilizar essa função e realizar a anotação de genes em cada
cromossomo individualmente.
Esta função separa cada cromossomo para realizar a anotação individualmente. Falta a implementação para expandir o data frame de acordo a colunas gene_ids, adicionada durante o processo de sobreposição de genes, como realizado em etapas anteriores. A falta dessa etapa pode gerar problemas durante a anotação dos simbolos de genes e entrezid. Após aplicar essas função, check os resultados, verifique a coluna gene_ids e realize a etapa de expansão do data frame.
cnv_individual_chromossomes <- function(n=22,dataf=cnv,chr_colx="seqnames", genes_colz="gene_ids"){
pb <- txtProgressBar(min = 0, max = n,style = 3, width = 50, char = "=") 11
for (i in 1:n) {
assign(paste0("chr",i), cnv_extract_chromossomes(chr_df=dataf,
chr_col=chr_colx))
assign(paste0("chr",i), cnvr_to_genes(df = paste0("chr",i), col = genes_colz)) assign(paste0("chr",i), paste0("chr",i)[,-10])
saveRDS(object = paste0("chr",i),file = paste0("chr",i,".rds"))
write.table(x = paste0("chr",i), file = paste0("chr",i,"_cnvs.rds"),
quote = F, row.names = F)
}
}Mais a frente iremos realizar o cálculo da correlação da expressão de genes e as mudanças no número de cópias. Para isso, precisamos que a tabela com as expressões esteja em um formato especifíco, como veremos a seguir.
8.1 Criar a coluna samples, coluna samples+genes
Criar esse novo identificador torna possível individualizar/ relacionar cada gene corretamente a sua amostra.
* mudar o identificador do tcga para a tabela de expressão como feito para a tabela de cnvs em etapas anteriores*
samples <- chr_exp[["submitted_sample_id"]]
samples=gsub(x = samples, "-", "_")
samples <- pbapply(as.data.frame(samples), 1, FUN = mudar_identificador_tcga)
chr_exp <- cbind(samples, chr_exp)Vamos criar os identificadores samples+ genes: podemos fazer para ambas as tabelas, cnvs e expressão, utilizando uma das maneiras a seguir.
chr_cnv <- sample_genes_ids(x1=cnvs, col_sample="samples", col_gene="SYMBOL") chr_exp <- sample_genes_ids(x1=chr_exp, col_sample="samples", col_gene="gene_id")Ou assim, utilizado pbaply: mais rápida.
As funções internas são disponibilizadas mais abaixo do código de exemplificação.
library(pbapply)
chr1_cnv= dplyr::select(chr1_cnv, "samples", "SYMBOL", "segment_mean", "type" ) chr1_cnv_new <- pbapply::pbapply(chr1_cnv, 1, FUN = sample_genes_ids_new) chr1_exp_new <- pbapply::pbapply(chr1_exp, 1, FUN = sample_genes_ids_new)Observar a coluna que corresponde aos simbolos dos genes em cada objeto e imputar corretamente na função.
Funções:
sample_genes_ids()
Essa função junta amostras e simbolos de genes/nomes em única coluna.
sample_genes_ids <- function(x1=chr1_cnv, col_sample="samples", col_gene="SYMBOL"){
# Define a barra de progresso
pb <- txtProgressBar(min = 0, max = nrow(x1), style = 3, width = 50, char = "=")
# Loop para criar sample+genes ids
sample_genes <- NULL
for (rows in 1:nrow(x1) ) {
ids<-paste0(x1[[col_sample]][rows], "_",
x1[[col_gene]][rows])
sample_genes<-rbind(sample_genes, ids)
# Progress bar
setTxtProgressBar(pb, rows)
}
close(pb)
# Retorna a coluna contendo samples+genes
x2<-cbind(sample_genes, x1)
rownames(x2)<-NULL
return(x2)
}Alternativa a função acima. Essa é mais rápida.
sample_genes_ids_new().
Para ser utilizada com pbaply ou outras funções semelhantes. É a função acima simplificada.
sample_genes_ids_new <- function(x){
# OBS.: The function expect sample names in vector (column) 1 and genes names in vector (column) 2
# Loop to create the labels
sample_genes <- paste0(chr_cnv[1], "_", chr_cnv[2])
# return a data frame with the column containing samples+genes names
rownames(sample_genes)<-NULL
return(sample_genes)
}9.1 Criar a tabela para cálculo da correlação
cordata <- merge(x = chr_cnv, y = chr_exp, by = "sample_genes")Observe que utilizamos o identificador sample + genes que criamos.
9.2 Calcular a correlação entre a variação na expressão de acordo a variação do número de cópias.
Aqui aplicamos o método de “spearman”, visto que esse não exije que os dados sigam uma distribuição normal.
correlation_results_chr1 <- cnv_exp_correlation(df = cordata, methodx = "spearman")Função:cnv_exp_correlation()
Essa é uma função para testar a correlação entre as regiões de cnvs para cada gene individualmente anotado e a expressão de cada gene. Os métodos de correlação podem ser “pearson”, “kendall”, “spearman”.
df = dataframe com uma coluna com os genes e as contagens de segment e read count
gene_ids_col= coluna que contem os genes
cnv_exp_correlation <- function(df=z,
x_col="segment_mean",
y_col="normalized_read_count",
gene_ids_col="gene_id",
methodx="spearman"){
pb <- txtProgressBar(min = 0,
max = length(df[[gene_ids_col]]),
style = 3,
width = 50,
char = "=")
correlation_results<-NULL
for (gene in unique(df[[gene_ids_col]])){
setTxtProgressBar(pb, gene)
x1 <- df[df[[gene_ids_col]]==gene,]
cor <- cor.test(x=x1[[x_col]],
y=x1[[y_col]],
method = methodx)
statistic <- cor$statistic
estimate <- cor$estimate
p_value <- cor$p.value
method <- cor$method
res <- cbind(gene, statistic, estimate , p_value, method)
rownames(res) <- NULL
correlation_results <- rbind(correlation_results, res)
}
close(pb)
return(correlation_results)
} Quanto aos resultados obtidos pela função acima:
S é o valor do teste estatístico (S = 10.871)
p-value é o nível de significância do teste estatístico (p-value = 0.4397).
Hipotese alternativa é um caracter descrevendo a hipotese alternativa (VERDADEIRO rho não é igual a 0).
Estimativa da amostra é o coeficiente de correlação. Para correlação de Spearman o coeficiente de correlação é denominado rho (Cor.coeff=0.4564)
Avaliar a correlação através de um plot
Avaliar a correlação da tabela e plotar apenas genes únicos.
Em z teremos a tabela concatenada com apenas um gene, x é o dado normalizado da contagens das reads no experimento de expressão de rna e y= a media dos segmentos do experimento de cnvs.
Plot 1: Scatter Plots
library("ggpubr")
ggscatter(z,
x = "normalized_read_count",
y = "segment_mean",
add = "reg.line",
conf.int = TRUE,
cor.coef = TRUE,
cor.method = "pearson",
xlab = "CNVs segment mean",
ylab = "Normalized read counts")