Carregando Dados de csv.Rmd

---
title: "Carregando Dados Genótipos dos Arquivos `csv`"
author: "James Hunter"
date: "21 de agosto de 2014"
output: pdf_document
---

Localização dos Arquivos `csv`
==============================

No formato fornecido pelo Ministério de Saúde, todos os genótipos vivem numa pasta identificado pela identificação da amostra. Cada pasta deve conter um arquivo `csv` que resume as mutações para o genes RT (transcriptase reversa) e PR (protease) do vírus. Eu posso obter acesso para algumas dessas pastas diretamente de um arquivo no formato comprido `zip`, que tem como nome a identificação do genótipo. Entretanto, existem exemplos em que o arquivo `zip` está localizado numa outra pasta com a identificação como nome. Eu desenvolvi duas estratégias para obter esses arquivos `zip`, abrí-los, e extrair o arquivo `csv`. Eu coloco o arquivo numa pasta chamado "20xxWork" onde o "xx" representa o número do ano com 2 dígitos.

A estrutura dos arquivos vindo do Ministério coloca os genótipos individuais dentro das pastas para os UF's em que o laboratório está localizado. Estas pastas ficam dentro de uma pasta para cada ano começando com 2007. 

Estratégia para Extração dos Arquivos `csv`
-------------------------------------------

**Passo 1** Dentro de um ano individual, eu abro uma pasta de uma UF. Por exemplo, para 2009 no estado de Pernambuco, esta pasta teria o nome "02PE09". Se esta pasta contem somente arquivos `zip` (e outros arquivos de formato texto), eu seleciono todos os arquivos e instruo o sistema operacional de abrir esses arquivos (Finder|Open). Esta instrução abre todos os arquivos `zip` dento da pasta.

**Passo 1 Alternativo** Se a pasta "02PE09" contivesse uma lista dos outras pastas, eu uso um "script" do programa Automator que escrevi que examina a pasta da UF e pastas subordinadas e abre todos os arquivos `zip` que ele acha.

**Limpeza depois do Passo 1** Estes dois procedimentos deixam um número grande de janelas abertas na tela. Depois, fecho todas as janelas de Finder. 

**Passo 2** Escrevi um script de Automator para examinar a pasta UF e as subordinadas, identificar qualquer arquivo com o final `csv` e transferir esses arquivos para a pasta "20xxWork" para aquele ano.

**Passo 3** Depois, eu apago a pasta UF e todo o conteúdo dela para ele não possa ser reaberta e ter os arquivos `csv` extraídos uma segunda vez. Este passo também conserva memoria do computador. Eu tenho dois backups, um com os arquivos originais do Ministério e o segundo que atualizo todos os dias em que eu trabalho com os dados.

Estrutura do Arquivo `csv` 
-------------------------

O arquivo `csv` tem cinco colunas. Essas colunas todos tratem das mutações existente na genotipagem. Essas colunas:

1. ID do genotipagem (relacionada a coluna `amostraID` na base de dados dos pacientes)
2. PR (Protease) ou RT (Transcriptase Reversa) que refere ao gene que gere a enzima.
3. Número de Codon com a mutação
4. Aminoácido produzido
5. Variável logica que indica se o codon só tem a mutação ou uma mistura do tipo selvagem e mutação.

O programa seguinte **só precisa ser ativado uma vez** (Por isso, está apresentada como comentário aqui). Ele estabelece a data frame para as mutações e tem nome `mutations`.

```{r mutationssetup}
# library(dplyr)
# mutations <- tbl_df(data.frame(id = NA, 
#                               gene = NA,
#                               codon = NA,
#                               aa = NA,
#                               wtMut = FALSE))
``` 

Depois de adicionar o primeiro caso a data frame, eu tiro esse caso exemplar que estabeleceu a estrutura. De novo, este código só precisa ser ativado uma vez. Por isso, está tratado como comentários aqui.

```{r remove 1st record}
# mutations <- mutations[-1,]
```

Depois de cada sessão de carregar o arquivo das mutações, eu salvo eles para o disco rígido como `mut.RData`.


```{r save mutations file}
#save(mutations, file = "mut.RData")
```

Revisão Preliminar dos Dados de Mutação para 2008 - 2010
========================================================

Eu converti os dados dos arquivos `csv` e os coloquei na data frame das mutações através do programa seguinte (tratado como comentários aqui):


```{r read csv files}
# dir <- "2010/2010Work"  # modificado para cada ano
#filelist <- list.files(path = dir, pattern = "csv")
#numfiles = length(filelist)

# passar por todos os arquivos, ler o conteúdo de cada mutação e copiar para data frame # das mutações

#for (i in 1:numfiles) {
#   fullfile <- paste(dir, "/", filelist[i], sep = "")

   # ler o arquivo na memoria como variável temp
#   temp <- read.csv(fullfile, header = FALSE, stringsAsFactors = FALSE)

   #contar o número de mutações (fileiras) no arquivo
#   rows <- nrow(temp)
   
   # processar cada fileira (mutação) em sequência
#   for (j in 1:rows) {
      # V1 - campo ID 
#      ID <- substr(temp$V1[j], 1, 10)
      # Si tiver carateres adicionais, tirar eles

      # V2 - Gene (ou "RT" ou "PR")
#      gene <- temp$V2[j]
      # V3 - codon
#      codon <- temp$V3[j]
      # V4 - aminoácidos (aa)
#      aa <- temp$V4[j]
      # v5 - wt (logical with TRUE = wild type + mutation, FALSE = mutation only)
#      wt <- as.logical(temp$V5[j])
      # colocar a mutação na data frame mutations
#      mutations <- rbind(mutations, data.frame(id = ID, gene = gene, codon = codon, 
#                                               aa = aa, wtMut = wt))
#   } # j processar linhas
#} # i processar arquivos
```


Este programa acima vai parar quando encontra um arquivo que não tem o formato esperado. Fiz isso de propósito para eu possa captar todos os arquivos que estão apresentando problemas e procurar alternativas para captar as mutações.

```{r echo=FALSE}
load(file = "mut.RData")
```

Problemas nos Arquivos `csv`
----------------------------

Os arquivos `csv` mostraram um número de problemas na gravação. Do total de 3.425 arquivos de 2008, somente 2.862 (83,6%) teve um formato que traduziu as mutações corretamente. Os arquivos vindo de Minas Gerais representaram a maior parte do problema: dos 520 arquivos de MG, 213 (41,0%) tiveram formatos corruptos e não podem ser incluídos na base de mutações. 

Também para 2008, houve um grupo de 10 arquivos `csv` que não teve nenhum esquema compreensível de nome que correspondeu ao campo `amostraid` na data frame `pacientes`. É provável que esses problemas apareceram porque os laboratórios estavam começando sua participação no programa de genotipagem. 

A preocupação central sobre esses arquivos inúteis é que eles ocorrem numa forma sistemática que pode comprometer a representividade dos dados. Se eles ficam distribuídos aleatoriamente entre as UF's, os resultados da análise não devem ser comprometidos. Complicando o problema com os arquivos `csv` de 2008 é que só 79 arquivos "BR" existem para aquele ano. Esses vem dos estados seguintes: 21PR(50), 13SP(8), 09RJ(1) e 05RS(20).

Para 2009, a taxa de erros é muito menor, 9,2%. Dos 4.864 arquivos `csv` daquele ano, 4.418 forneceram mutações úteis para a data frame `mutations` (90.8%). Para 2010, do total de 3.373 arquivos `csv`, somente 83 (2,5%) estavam corruptos. Desses, 26 (31,3%) vieram do grupo 25RS. 

Suplementando os Arquivos `csv` com os Arquivos BR
--------------------------------------------------

Para suplementar a busca de genotipagens vindo dos arquivos `csv`, eu quis usar os arquivos "xBR.pdf" que muitos das pastas de genotipagem têm. Criei um programa que constrói uma lista dos identificadores `amostraid` na data frame `mutations` para um ano e corresponder eles para a lista dos arquivos BR daquele ano. Quando existe um arquivo BR disponível para um `amostraid` que não existe ainda na data frame `mutations`, o programa criará uma lista desses arquivos para que os dados das mutações podem ser transferidos para `mutations`. 

Selecionando os Arquivos BR
============================

```{r select files}

# load libraries
library(dplyr)
library(tidyr)
library(limma)
load(file = "mut.RData")
options(warn = -1)
# Prepare list of files to to scrape within the BR template directories
year <- "10"

cases <- unique(mutations$id[substr(mutations$id, 5, 6) == year])

# Make list of BR files available for that year

dir <- paste("20", year, "/BR files", sep = "")
filelist <- substr(list.files(path = dir, pattern = "BR.txt"), 1, 10)
# set up vector to show matches; files we want will be FALSE

idmatch <- substr(filelist, 1, 10) %in% substr(cases, 1, 10) # substr to deal with 12 char id's

# identify files to be scraped

filelist <- filelist[idmatch == FALSE]

numfiles = length(filelist)
```

Programa para Processar Arquivos em `filelist`
==============================================

Função para Analisar os Codons nas Mutações
-------------------------------------------

Uma função em R é um pequeno bloco de código que o programa maior vai utilizar várias vezes para completar uma tarifa especial. Neste caso, a função `codaa` analisa todos os codons que encontra e separar eles em número de codon (localização no gene).

```{r function}
# função para testar o número de carateres no codon e aminoácido sequências
codaaReturn <- c(NA, NA, NA, FALSE) # codon, aa, wt, ok
codaa <- function(txt, chars) {
   ok <- TRUE
   if (chars == 3) {
      cod <- substr(txt, 1, 2)
      AA <- substr(txt, 3, 3)
      wt <- FALSE} else 
         if (chars == 4) {
            cod <- substr(txt, 1, 3)
            AA <- substr(txt, 4, 4)
            wt <- FALSE} else 
               if (chars == 5) {
                  cod <- substr(txt, 1, 2)
                  AA <- substr(txt, 3, 5)
                  wt <- TRUE} else 
                     if (chars == 6) {
                        cod <- substr(txt, 1, 3)
                        AA <- substr(txt, 4, 6)
                        wt <- TRUE} else {
                           cat("problem with", amostraid, tipo, txt, sep = " ")
                           ok <- FALSE
                        }
   codaaReturn <- c(cod, AA, wt, ok)
   return (codaaReturn)
} # termina codaa
```

Programa para processar os arquivos BR para um ano 
--------------------------------------------------

```{r file loop}
# Passar pelos arquivos para captar os codons
#for (i in 1:numfiles) {
#   fullfile <- paste(dir, "/", filelist[i], "BR.txt", sep = "")
#   filetext <- readLines(con = fullfile, skipNul = TRUE, encoding = "UTF8")

# Procurar linha que tem amostraid
#   x <- grep("Paciente:", c(filetext[2], filetext[3], filetext[4]))
#   idline <- x + 1 
#   amostraid <- substring(filetext[idline], 32, 41)

# Processar RT codons
#   gene <- "RT"

# achar linhas começando com texto 'IRTN'; aumentar por 1
#   tipo <- "ITRN"
#   itrnLine <- idline + 2
#   itrnText <- filetext[itrnLine]

# Criar um vetor das mutações ITRN
#   itrn <- as.vector(strsplit2(itrnText, split = " "))

# Processar todos os codons na lista das mutações
#   for (i in 1:length(itrn)){
#      txt <- itrn[i]
#      chars <- nchar(txt)
   # Analisar o campo de codon aminoácido baseado no # de carateres
#      codaaReturn <- c(NA, NA, NA, FALSE) # codon, aa, wt, ok  
#      codaaReturn <- codaa(txt, chars)  
    # termina if
    # anexar esta mutação para base mutations
#      if (codaaReturn[4]) testframe <- rbind(testframe, data.frame(id = amostraid, 
#                                 gene = gene, codon = codaaReturn[1], aa = codaaReturn[2], 
#                                 wtMut = codaaReturn[3])) 
#   }
# Processar codons ITRNN, se tem
#   tipo <- "ITRNN"
# Pular 2 linhas para ver se começar com a palavra "Outros"; se for, avançar para 
# Protease
# senão, processar os codons 'ITRNN'
#   if(substr(filetext[itrnLine + 2], 1, 5) != "Outro") {
#      itrnnText <- filetext[itrnLine + 2]
   # Criar um vetor das mutações ITRNN 
#      itrnn <- as.vector(strsplit2(itrnnText, split = " "))
   # Processar todos os codons na lista das mutações
#      for (i in 1:length(itrnn)){
#         txt <- itrnn[i]
#         chars <- nchar(txt)
      # Analisar o campo de codon aminoácido baseado no # de carateres
#      codaaReturn <- c(NA, NA, NA, FALSE) # codon, aa, wt, ok  
#      codaaReturn <- codaa(txt, chars)  
      # termina if
      # anexar esta mutação para base mutations
#      if (codaaReturn[4]) testframe <- rbind(testframe, data.frame(id = amostraid, 
#                               gene = gene, codon = codaaReturn[1], aa = codaaReturn[2], 
#                               wtMut = codaaReturn[3])) 
#   }

# Processar Inhibidor de Protease codons
#   gene <- "PR"
#   tipo <- "PR"

# Achar linha com a palavra "Protease""; mutações começam +1 linha
#   x <- grep("Protease", c(filetext[10], filetext[11], filetext[12], filetext[13]))
#   prLine <- x[1] + 10 
# 10 é a primeira linha possível; mutações ficam na próxima
# Se prLine começa com "Outros"; se, for, avançar para final do programa
# Senão, processar os codons como anteriormente
#   if(substr(filetext[prLine], 1, 5) != "Outro") {
#      prText <- filetext[prLine]

# dividir a linha nos campos componentes
#   pr <- as.vector(strsplit2(prText, split = " "))

# Processar as mutações PR e colocá-las na base mutations
#   for (i in 1:length(pr)){
#      txt <- pr[i]
#      chars <- nchar(txt)
   # Analisar o campo de codon aminoácido baseado no # de carateres
#   codaaReturn <- c(NA, NA, NA, FALSE) # codon, aa, wt, ok  
#   codaaReturn <- codaa(txt, chars)  
    # termina if
    # anexar esta mutação para base mutations
#   if (codaaReturn[4]) testframe <- rbind(testframe, data.frame(id = amostraid, 
#                                    gene = gene, codon = codaaReturn[1], aa = codaaReturn[2], 
#                                    wtMut = codaaReturn[3])) 
#}}}} # termina loop dos arquivos
# reiniciar advertências para 0
#options(warn = 0)
```

Problema com 25RS
-----------------

Para os arquivos vindo de UF 25RS, todos os arquivos para 2009 tiveram 0 bytes, ou seja, estavam vázios e por isso inútil. Para 2010, esta situação foi revertido. Para este ano, `r length(mutations$id[substr(mutations$id, 1, 6) == "25RS10"])` arquivos produziram mutações úteis.


```{r echo=FALSE}
# number of useful samples included in mutations from 2008 MG
load(file = "mut.RData")
(MG08cases <- length(unique(mutations$id[substr(mutations$id, 3, 6) == "MG08"])))
```

Número de Genotipagens por Ano
==============================

Utilizando a informação no campo `id` (caráteres 5 e 6), os arquivos `csv` têm o número de genotipagens seguintes para cada ano.

```{r count cases}
load(file = "pacientes.RData")
cases08 <- unique(mutations$id[substr(mutations$id, 5, 6) == "08"])
cases09 <- unique(mutations$id[substr(mutations$id, 5, 6) == "09"])
cases10 <- unique(mutations$id[substr(mutations$id, 5, 6) == "10"])
length(cases08) 
length(cases09)
length(cases10)
# number of useful samples included in mutations from 2008 MG
(MG08cases <- length(unique(mutations$id[substr(mutations$id, 3, 6) == "MG08"])))
```

Concordância com os `amostraid`'s da Data Frame dos Pacientes
-------------------------------------------------------------

Eu tirei amostras da cabeça e da cauda (300 cada amostra) da data frame `mutations`. Também tirei uma amostra aleatória (1.000 mutações) para ver quantos correspondem aos dados dos pacientes válidos. (`amostraid`). 

```{r test against patients amostraid}
load(file = "mut.RData")
x <- unique(mutations$id[1:300])
y <- unique(tail(mutations$id, 300))
z <- unique(sample(mutations$id, 1000))
length(x) ; length(y) ; length(z)
x
y
xtest <- which(substr(x, 1, 10) %in% substr(pacientes$amostraid, 1, 10)) 
# testar com os identificadores dos pacientes
ytest <- which(substr(y, 1, 10) %in% substr(pacientes$amostraid, 1, 10))
ztest <- which(substr(z, 1, 10) %in% substr(pacientes$amostraid, 1, 10))
pacientes[xtest, c(1:3,13:14)]
pacientes[ytest, c(1:3,13:14)]
length(ztest)
        
```

Esses testes são preliminares. As primeiras duas amostras incluem um número grande de duplicações porque todos os registros na data frame são mutações individuais e alguns casos tiveram até 30 mutações e polimorfismos. Como resultado, a primeira amostra teve  `r length(x)` casos únicos e a amostra de cauda `r length(y)`. A amostra aleatória teve muito menos duplicação porque o software não tirou registros contíguos. Assim, esta amostra teve `r length(z)` casos únicos para testar com a data frame dos pacientes.

O identificador do paciente na data frame `pacientes` (`codepac`) liga o resultado da genotipagem (através da variável `amostraid`) para a carga viral e a contagem de CD4+, esta correspondência deve ser exata. Infelizmente, os dados não são assim. `r length(xtest)` casos na amostra de cabeça correspondem às `amostraid`'s na data frame `pacientes` (`r sprintf("%.1f",100*length(xtest)/length(x))`%). Para a amostra de cauda, `r length(ytest)` casos corresponderam às `amostraid`'s na data frame `pacientes`.

Os números para a amostra aleatória são aproximadamente iguais a porcentagem dos sucessos para a amostra de cabeça. Dos `r length(z)` casos ´únicos, `r length(ztest)` corresponderam às `amostraid`'s na data frame `pacientes` (`r sprintf("%.1f",100*length(ztest)/length(z))`%). Na amostra aleatória, os casos que não corresponderam aos registros na data frame `pacientes` eram os seguintes:

```{r non-matches}
ztf <- substr(z, 1, 10) %in% substr(pacientes$amostraid, 1, 10)
# print which cases in sample are false
zfalse <- which(ztf == FALSE)
z[zfalse]
```

Acho que a taxa de correspondência de aproximadamente 90% nos primeiros três anos de um cadastro abrangente de genotipagens é aceitável. Espero que esta taxa aumentará a medida que transferir os dados dos arquivos `csv` para 2011 e adiante porque os operadores do sistema, os laboratórios autorizados, entendem melhor o emprego dos sistemas do Ministério.

Fora dos problemas com 2008 em Minas Gerais, não houve outros padrões persistentes que podem prejudicar a utilidade dos dados.

Anos 2011 - 2014
================

Eu comecei de processar as pastas desses anos usando os programas e procedimentos descritos acima. Vou completar esse trabalho em volta de 25 de agosto.