Mestranda no programa de Endocrinologia da Faculdade de Medicina da USP - FMUSP. Faz pesquisa na área de biologia molecular e bioinformática.

exec

Trabalho Final - Propostas A e B

PROPOSTAS PARA O TRABALHO FINAL

Proposta A

Função plot.coverage

Função irá criar um gráfico detalhando a cobertura do sequenciamento por base de uma ou várias amostras sequenciadas por exoma. O eixo X terá os números de pares de bases, enquanto o eixo Y tem a cobertura base a base, e uma linha pelo gráfico indicando os pontos da cobertura. A importância de se ter uma cobertura boa de cada base é a confiança do que se está vendo é real ou não, ou seja, que a presença de variante, indels ou até mesmo a falta de variantes na sequencia esteja realmente no material genético do paciente. Muitos laboratórios padronizam o mínimo de cobertura desejável como 20 reads por base, e pedaços do genoma que não alcançam essa cobertura acabam por ter de ser refeitos, o que gera um grande custo para o laboratório. A ideia desse gráfico, então, é de facilitar a ilustração e, consequentemente, a visualização da quantidade de reads por base em um gene, ou seja, a cobertura do sequenciamento. Será tomado como exemplo os gráficos de cobertura do site do ExAC (http://exac.broadinstitute.org/). Todas as informações necessárias para a criação do gráfico serão retiradas de um arquivo .BED provindo de um .BAM que contenha as informações de cobertura. Este arquivo possui formato de tabela, universal para o tipo de arquivo, e poderá ser lido como um data frame. Cada uma das colunas possui uma informação importante para a estimativa dessa cobertura por base, sendo elas o cromossomo, início do transcrito, final do transcrito, nome do gene, tamanho do fragmento, base do fragmento e cobertura na base. O gene a ser plotado também poderá ser escolhido pelo usuário como argumento da função. A função também terá como opção fazer o gráfico de apenas uma amostra, que no caso só iria plotar as informações presentes no arquivo, dependendo do gene escolhido pela pessoa (que deve estar no arquivo de origem, se não a função retornará erro), ou então fazer o gráfico de cobertura de diversas amostras (que entrarão em um for), podendo serem feitas a média (mean()), mediana (median()) ou variancia (var()) por base, dependendo da escolha do usuário entre cada uma ou todas as três, que serão então plotadas em gráficos diferentes.

Proposta B

Função plot.saturation

Função irá criar um gráfico de saturação de cobertura por quantidade total de reads na amostra, para responder a questão de quantos reads são necessários para se atingir uma média/mediana de 100 vezes de cobertura. O gráfico então terá o eixo X com a quantidade de reads em milhões, e o eixo Y com a média/mediana de 0 a 100. Para isso será necessário como entrada uma tabela com a informação de cobertura com a quantidade de reads desejados para várias amostras, a fim de criar o gráfico. Cada uma das amostras vira um ponto nos gráficos, que serão unidos por uma linha. A ideia do gráfico é ajudar o pesquisador a decidir a quantidade de reads necessários para o sequenciamento chegar a uma cobertura ideal, para que informações importantes como por exemplo variantes pontuais, não sejam perdidas na análise.

Apesar das funções propostas não serem funções genéricas, são função que poderiam ser adaptadas para outros usos, uma vez que se utilizam de uma entrada que possui um formato simples e padronizado, e tem como finalidade fazer medidas estatísticas de variáveis, afim de criar gráficos. Além disso, são funções interessantes para a área de sequenciamentos de larga escala, como aquelas presentes em pacotes do Bioconductor, que possui 1383 pacotes de R para análises e compreensão de dados genômicos. Portanto, são funções que poderiam ser utilizadas por todos aqueles que trabalham com sequenciamento de larga escala, uma vez que a cobertura é muito importante para a análise final e validação das variantes encontradas.

Proposta A

plot.baseCoverage(x, geneID, med = c(“mean”, “median”, “var”))

Estrutura do arquivo de dados (universal, explicado por extenso no FAQ do UCSC), que possui:

Sendo que as duas últimas colunas são específicas do OUTPUT do coverageBed com a opção -d, utilizado para verificar a cobertura de regiões de interesse do pesquisador. Também são OUTPUTs deste programa o default e o -hist, que são como seguem:

Default:

-hist:

E, no final do arquivo, um histograma resumindo a cobertura de todas as regiões.

Entretanto, para esses dois outros OUTPUTs não seria possível fazer um gráfico xy, e sim um histograma. Por esse motivo, é mais interessante que essas funções fossem desenvolvidas separadamente. Dessa forma, optou-se pelo desenvolvimento da primeira função aqui descrita na proposta A, dado que envolve mais passos do que apenas a criação de um gráfico.

Antes de iniciar a função:

• Ler o arquivo de dados e salvar em um objeto
• OU ler todos os arquivos de uma pasta e salvar em um objeto

Na função:

• SE x == 1:
• gene ← tabela aonde aparece o gene passado como argumento
• xmin ← (primeiro valor da coluna 2) + 5
• xmax ← (ultimo valor da coluna 3) + 5
• cov ← ultima coluna da tabela gene
• eixo.x ← (coluna 2) + (coluna 7)
• Plotar um gráfico com eixo.x e cov, utilizando como limites do eixo x os objetos xmin e xmax (plot(type=”l”)).
• SE x > 1:
• eixo.x, xmin, xmax, covs ← NULL
• Entrar em um FOR (i in files)
• x ← ler(arquivo)
• gene ← tabela do gene dado como argumento
• eixo.x ← (coluna 2) + (coluna 7)
• xmin ← (primeiro valor da coluna 2) + 5
• xmax ← (ultimo valor da coluna 3) + 5
• covs ← cbind(covs, gene$ncol(gene))
• SE o argumento em med for “mean”
• mea ← apply(covs, 1, mean)
• Plotar gráfico eixo.x e mean
• SE o argumento em med for “median”
• med ← apply(covs, 1, median)
• Plotar gráfico eixo.x e med.
• SE o argumento em med for “var”
• vari ← apply(covs, 1, var)
• Plotar gráfico eixo.x e vari.
• SE não houver argumento, fazer os três.

Proposta B

plot.saturation(x)

Fará um gráfico xy de variáveis quantitativas que se encontrem em mais de um arquivo, onde o dado se localize na mesma coluna. O arquivo de dados utilizados pode ser uma tabela simples em que os dados se encontram.

Antes da função: folders ← nomes(pastas em um diretório) para colocar na função

Quando na função:

• Criar objetos nulos
• FOR (i in 1:(ultimo vetor em folders))
• Ler arquivos presentes em folders[i]
• Salvar em diferentes objetos para cada folder
• FOR (arquivos no objeto pasta 1..ultima)
• ler(arquivo) → tabela
• apply(tabela, 2, mean) → mean
• rbind(objeto, mean) → objeto
• plot(medias)

Downloads

Página Help

Código

plot.baseCoverage <- function(x, geneID, measure){ ##função e seus argumentos padrões
  if(missing(geneID)) { ##se não for colocado algo no argumento do geneID
    stop("A gene identification is necessary for this function") ##dar um aviso que ele é um argumento necessário para que a função rode
  }
  if(class(x) == "data.frame"){ ##testa se o primeiro objeto inserido com um argumento é um dataframe, que é o objeto criado pelo read.table
    genes.ids <- unique(x$V4) ##cria um objeto com todos os nomes de genes como aparecem no arquivo
    if (any(geneID == genes.ids)){ ##se o argumento dado for encontrado no arquivo, continuar
      gene <- x[x$V4 == geneID,] ##salva apenas a tabela do gene de interesse no objeto, dado como argumento no geneID
      gene.min <- gene[1,2] ##determina a primeira base do gene
      gene.max <- gene[nrow(gene),3] ##determina a ultima base do gene
      gene$exact.base <- gene$V2 + gene$V7 ##cria uma coluna na tabela do gene correspondente a cada uma das bases que estão na tabela original
      dataf1 <- cbind(gene$exact.base, gene$V8) ##cria um dataframe com a cobertura em cada base e sua base correspondente
      gene.extension <- seq(gene.min, gene.max) ##cria um vetor com cada base no gene inteiro (incluindo introns, se não estiverem na tabela original)
      fragments <- c(unique(gene$V2), unique(gene$V3)) ##cria um objeto com o começo e o final de todos os fragmentos, ou seja, inicio e fim de todos os exons
      fragments <- sort(fragments) ##ordenar os fragmentos em ordem numérica
      dataf2 <- NULL ##cria um objeto nulo para:
      dataf2 <- cbind(dataf2, gene.extension) ##cria um dataframe com todas as bases do gene inteiro. a ideia é que se tiverem bases que não tiveram leitura, não se tenha um ponto para elas. por isso:
      dataf3 <- merge(dataf2, dataf1, by.x = "gene.extension", by.y = "V1", all = T) ##junta os dois dataframes anteriores criados em um só pela base, colocando NA nas bases que não possuem leitura
      par(las = 1, bty = "l") ##especifica para o gráfico o tipo de caixa, números na horizontal
      plot(dataf3, type = "l", xaxt = "n", ##plota o gráfico final, com uma linha ligando todos os pontos, sem marcas no eixo X
         ylim = range(0:max(gene$V8)), ##padroniza o eixo Y para que seja possível se ver quais bases chegam até 100,
         xlab = geneID, ylab = "Coverage", main = "Sample") ##nomeia os eixos X e Y e o gráfico
      axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
    } else { ##caso o gene dado no argumento não seja encontrado no arquivo
      stop("The geneID was not found in the file! Please make sure it is the same as appears on your file or it is not misspeled!") ##retornar mensagem de erro
    }
  } else { ##caso o primeiro if retorne negativo, entrar na opção de vários arquivos
    ##criação de objetos que serão utilizados após o for
    coverages <- NULL ##para salvar as coberturas
    gene.extension <- NULL ##para salvar a extensão do gene
    exact.base <- NULL ##para salvar a base exata de onde se tem cobertura
  
    for (i in x){ ##para cada um dos arquivos cujo nomes estão em x
      table <- read.table(i, header = F, sep = "\t", as.is = T) ##leitura dos arquivos como dataframes
      genes.ids <- unique(table$V4) ##cria um objeto com todos os nomes de genes como aparecem no arquivo
     if (any(geneID == genes.ids)){ ##se o argumento dado for encontrado no arquivo, continuar
        gene <- table[table$V4 == geneID,] ##salva apenas a tabela do gene de interesse dado no argumento geneID no objeto
        gene.min <- gene[1,2] ##determina a primeira base do gene
        gene.max <- gene[nrow(gene),3] ##determina a ultima base do gene
        exact.base <- gene$V2 + gene$V7 ##cria um vetor que possua cada uma das bases que possuem cobertura
        fragments <- c(unique(gene$V2), unique(gene$V3)) ##cria um objeto com todos os inícios e fins dos fragmentos, ou seja, o inicio e o fim dos exons, que são a parte que terão cobertura no caso do exoma
        fragments <- sort(fragments) ##ordena esse objeto em ordem numérica
        coverages <- cbind(coverages, gene$V8) ##cria dataframe apenas com as coberturas de todos os arquivos
        gene.extension <- seq(gene.min, gene.max) ##cria um vetor com cada base no gene inteiro (incluindo introns, se eles não estiverem na tabela original)
    } else { ##caso o gene dado no argumento não seja encontrado no arquivo
        stop("The geneID was not found in the file! Please make sure it is the same as appears on your file or it is not misspeled!") ##retornar mensagem de erro
      }
        }
    dataf2 <- NULL ##cria um objeto nulo para:
    dataf2 <- cbind(dataf2, gene.extension) ##transforma o gene.extension criado anteriormente em dataframe
    par(las = 1, bty = "l") ##padroniza para todos os gráficos feitos o tipo de caixa e os números nas horizontais nos eixos
    if(missing(measure)){ ##se não tiver o argumento med, retornar os 3 gráficos
        media <- apply(coverages, 1, mean) ##aplicar a função média escolhida no argumento
        dataf1a <- cbind(exact.base, media) #juntar as médias com as bases aonde se encontram em um dataframe
        dataf3a <- merge(dataf2, dataf1a, by.x = "gene.extension", by.y = "exact.base", all = T) ##juntar os dois dataframes criados em um só pela base, afim de se ter um dataframe no qual não se tem nada nas bases que não tiveram cobertura nenhuma
        plot(dataf3a, type = "l", xaxt = "n", ##plota o gráfico com a média das bases e retira as marcas e números no eixo X
             ylim = range(0:max(media)), ##padroniza o eixo Y
             xlab = "geneID", ylab = "Coverage", main = "Mean Coverage") ##nomeia o gráfico e seus eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
        mediana <- apply(coverages, 1, median) ##aplicar a função média escolhida no argumento
        dataf1b <- cbind(exact.base, mediana) #juntar as médias com as bases aonde se encontram em um dataframe
        dataf3b <- merge(dataf2, dataf1b, by.x = "gene.extension", by.y = "exact.base", all = T) ##juntar os dois dataframes criados em um só pela base, afim de se ter um dataframe no qual não se tem nada nas bases que não tiveram cobertura nenhuma
        par(ask=T) ##pergunta antes de plotar os gráficos, pedindo para que o enter seja apertado
        plot(dataf3b, type = "l", xaxt = "n", ##plota o gráfico com a mediana e retira os números do eixo X
             ylim = range(0:max(mediana)), ##padroniza o eixo Y
             xlab = geneID, ylab = "Coverage", main = "Median Coverage") ##nomeia o gráfico e seus eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
        dataf1c <- cbind(exact.base, coverages) #juntar as médias com as bases aonde se encontram em um dataframe
        plot(NULL, ##plota o gráfico vazio
             ylim = range(0:max(coverages)), xlim = range(gene.min, gene.max), ##especifica os limites dos eixos X e Y
             xaxt = "n", ##retira as marcas do eixo X
             xlab = geneID, ylab = "Coverage", main = "Variance") ##nomeia o gráfico e os dois eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
        for(i in 2:ncol(dataf1c)){ ##entra em um for para plotar todos os pontos de cobertura encontrados no dataframe
          points(dataf1c[,1], dataf1c[,i], pch = 20) #plota cada um dos pontos no gráfico
        }
        lines(dataf3a, col = "red") ##plota a linha de média de referência
      }
      else if(measure == "mean"){ ##ou, caso o argumento dado em med for "mean", ou seja, a média
        media <- apply(coverages, 1, mean) ##aplica a função média a cada linha do dataframe, para determinar a média de cobertura de cada base 
        dataf1 <- cbind(exact.base, media) #junta as médias com suas bases em um dataframe
        dataf3 <- merge(dataf2, dataf1, by.x = "gene.extension", by.y = "exact.base", all = T) ##juntar os dois dataframes criados em um só pela base, afim de se ter um dataframe no qual as bases que não tiveram cobertura são NAs
        plot(dataf3, type = "l", xaxt = "n", ##plota o gráfico com a média das bases e retira os números do eixo X
           ylim = range(0:max(media)), ##padroniza o eixo Y
           xlab = geneID, ylab = "Coverage", main = "Mean Coverage") ##nomeia o gráfico e seus eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
      }
      else if(measure == "median"){ ##ou, caso o argumento dado em med for "median", ou seja, mediana
        mediana <- apply(coverages, 1, median) ##aplica a função mediana em cada linha do dataframe, e guarda em um objeto com apenas essas medianas
        dataf1 <- cbind(exact.base, mediana) #junta as medianas e as bases correspondentes em um dataframe
        dataf3 <- merge(dataf2, dataf1, by.x = "gene.extension", by.y = "exact.base", all = T) ##junta os dois dataframes criados em um só pela base, afim de se ter um dataframe no qual as bases que não tiveram coberturas apresentam NAs
        plot(dataf3, type = "l", xaxt = "n", ##plota o gráfico com a mediana  e retira os números do eixo X
           ylim = range(0:max(mediana)), ##padroniza o eixo Y
           xlab = geneID, ylab = "Coverage", main = "Median Coverage") ##nomeia o gráfico e seus eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
      }
      else if(measure == "variance" || measure == "var"){ ##ou, caso o argumento dado em med for "var", ou seja, variancia
        dataf1 <- cbind(exact.base, coverages) #juntar as coberturas com as bases
        media <- apply(coverages, 1, mean) ##aplica a função média a cada linha do dataframe, para determinar a média de cobertura de cada base 
        datafm <- cbind(exact.base, media) #junta as médias com suas bases em um dataframe
        datafm <- merge(dataf2, datafm, by.x = "gene.extension", by.y = "exact.base", all = T) ##juntar os dois dataframes criados em um só pela base, afim de se ter um dataframe no qual as bases que não tiveram cobertura são NAs
        plot(NULL, ##plota o gráfico vazio
           ylim = range(0:max(coverages)), xlim = range(gene.min, gene.max), ##especifica os limites dos eixos X e Y
           xaxt = "n", ##retira as marcas do eixo X
           xlab = geneID, ylab = "Coverage", main = "Variance") ##nomeia o gráfico e os dois eixos
        axis(1, at = fragments, labels = F) ##coloca marcas em cada começo e fim dos fragmentos no eixo X do gráfico
        for(i in 2:ncol(dataf1)){ ##entra em um for para plotar todos os pontos de cobertura encontrados no dataframe
          points(dataf1[,1], dataf1[,i], pch = 20) #plota cada um dos pontos no gráfico
        }
        lines(datafm, col = "red") ##plota a linha de média de referência
      }
  }
}