raw-data-process.Rmd

---
title: "raw data process"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## growth curve data

```{r}
## read qPCR raw data
result_file <- list.files(path="./data-raw/qPCR_growth_curve", 
                          pattern = "Results_ViiA7_export.txt",
                          full.names = T)
result1 <- read_delim(result_file[[1]],"\t",skip=43,na="Undetermined")
result2 <- read_delim(result_file[[2]],"\t",skip=43,na="Undetermined")
result2 <- filter(result2,Task!="STANDARD")
results <- rbind(result1,result2)

## read qPCR meta data
meta_file <- list.files(path="./data-raw/qPCR_growth_curve", 
                        pattern = "meta.txt",
                        full.names = T)
meta <- read_delim(meta_file,"\t")

## combine results with meta data
results <- results %>% 
  select(Well,`Sample Name`,`Target Name`,CT) %>% 
  dplyr::rename(sample=`Sample Name`,species=`Target Name`) %>%
  left_join(meta,by="sample")

## remove unspecific control
results <- results %>%
  filter(!(condition=="all" & species=="PP")) %>%
  filter(!(condition=="none" & species=="EC"))

## quantify
cal_quantity <- function(ct){10^(-.2922*ct+12.077)}
results <- results %>% mutate(quantity=cal_quantity(CT)*dilution) %>%
  filter(quantity > 1e5)


qPCR_data <- results %>% 
  select(sample, condition, time, species, quantity) %>%
  arrange(as.numeric(str_extract(sample,"[0-9]+")), condition)

write.csv(qPCR_data, file = "data/qPCR-data.csv",row.names = FALSE)
```


## RNA-seq

### Identifing differentially expressed genes

Start from ht-seq counts.

```{r deseq}
ht_counts <- readRDS("data/rna/ht_counts.rds")
library(DESeq2)

myDESeqMatrix <- function(ht_counts=ht_counts, org=NULL) {
  require(dplyr)
  require(tidyr)
  require(tibble)

  if (!is.null(org)){
    ht_counts <- ht_counts %>% filter(organism==org) %>%
      filter(ratio0 != ifelse(org=="EC","none","all"))
  }
  
  count_data <- ht_counts %>% dplyr::select(gene,count,sample_id) %>%
    group_by(sample_id) %>%
    spread(sample_id,count) %>% 
    as.data.frame() %>%
    column_to_rownames(var = "gene")
  col_data <- ht_counts %>% dplyr::select(sample_id,ratio0,time,group) %>% 
    arrange(sample_id) %>% distinct()
  return(list("count_data"=count_data,"col_data"=col_data))
}

mat.EC <- myDESeqMatrix(ht_counts,org = "EC")
mat.PP <- myDESeqMatrix(ht_counts,org = "PP")

dds.EC <- DESeqDataSetFromMatrix(countData = mat.EC$count_data, 
                                 colData   = mat.EC$col_data,
                                 design    = ~ group)
dds.PP <- DESeqDataSetFromMatrix(countData = mat.PP$count_data, 
                                 colData   = mat.PP$col_data,
                                 design    = ~ group)

# run DESeq()
dds.EC <- DESeq(dds.EC)  # this step is time consuming, normalization
dds.PP <- DESeq(dds.PP)

# get deg
comparisons.EC <- list(c("less_0h","all_0h"),
                       c("equal_0h","all_0h"),
                       c("more_0h","all_0h"),
                       c("less_4h","all_4h"),
                       c("equal_4h","all_4h"),
                       c("more_4h","all_4h"),
                       c("less_8h","all_8h"),
                       c("equal_8h","all_8h"),
                       c("more_8h","all_8h"),
                       c("less_24h","all_24h"),
                       c("equal_24h","all_24h"),
                       c("more_24h","all_24h"))

comparisons.PP <- list(c("less_0h","none_0h"),
                       c("equal_0h","none_0h"),
                       c("more_0h","none_0h"),
                       c("less_4h","none_4h"),
                       c("equal_4h","none_4h"),
                       c("more_4h","none_4h"),
                       c("less_8h","none_8h"),
                       c("equal_8h","none_8h"),
                       c("more_8h","none_8h"),
                       c("less_24h","none_24h"),
                       c("equal_24h","none_24h"),
                       c("more_24h","none_24h"))

DEG_results.EC <- myDEG_Results(dds=dds.EC, comparison = comparisons.EC)
DEG_results.PP <- myDEG_Results(dds=dds.PP, comparison = comparisons.PP)
gene_expression.EC <- myDEG_Results(dds.EC,comparison = comparisons.EC,filtered = F)
gene_expression.PP <- myDEG_Results(dds.PP,comparison = comparisons.PP,filtered = F)


saveRDS(dds.EC,"data/dds.EC.2.rds")
saveRDS(dds.PP,"data/dds.PP.2.rds")
saveRDS(DEG_results.EC,"data/DEG_results.EC.rds")
saveRDS(DEG_results.PP,"data/DEG_results.PP.rds")
saveRDS(gene_expression.EC, "data/gene_expression.EC.rds")
saveRDS(gene_expression.PP, "data/gene_expression.PP.rds")
```