MahShaaban / pcr

Quality assessing, analyzing and testing the statistical significance of real-time quantitative PCR data
https://CRAN.R-project.org/package=pcr
GNU General Public License v3.0
27 stars 8 forks source link

genorm calculations #7

Open MahShaaban opened 6 years ago

MahShaaban commented 6 years ago

Related #4

Here, I tried to quickly reproduce the normalization factor and the relative fold change using an example of the geNorm method.

The following code shows the data in the attached file and compares the calculations from pcr_nf and pcr_genorm to the original calculations from the example (a clean version of the file is attached). clean_example.xlsx

# load required libraries
library(tidyverse)
library(readxl)
library(FSA)
#> ## FSA v0.8.20. See citation('FSA') if used in publication.
#> ## Run fishR() for related website and fishR('IFAR') for related book.
library(pcr)

# load and show raw data from the example
raw_data <- read_excel('./clean_example.xlsx', sheet = 1)
print(raw_data)
#> # A tibble: 16 x 8
#>    Target Sample    `Cq value` `Mean Cq` `SD Cq`      RQ  `SD RQ`  `SE RQ`
#>    <chr>  <chr>          <dbl>     <dbl>   <dbl>   <dbl>    <dbl>    <dbl>
#>  1 REF1   Control         24.5      24.6  0.212   1       0.147    0.104  
#>  2 <NA>   Control         24.8      NA   NA      NA      NA       NA      
#>  3 <NA>   Treatment       25.4      25.6  0.283   0.518   0.101    0.0718 
#>  4 <NA>   Treatment       25.8      NA   NA      NA      NA       NA      
#>  5 REF2   Control         18.3      18.4  0.0707  1       0.0490   0.0347 
#>  6 <NA>   Control         18.4      NA   NA      NA      NA       NA      
#>  7 <NA>   Treatment       18.8      18.8  0.0707  0.758   0.0371   0.0263 
#>  8 <NA>   Treatment       18.7      NA   NA      NA      NA       NA      
#>  9 REF3   Control         32.3      32.3  0.0707  1       0.0490   0.0347 
#> 10 <NA>   Control         32.4      NA   NA      NA      NA       NA      
#> 11 <NA>   Treatment       33        33.2  0.354   0.536   0.131    0.0929 
#> 12 <NA>   Treatment       33.5      NA   NA      NA      NA       NA      
#> 13 GOI1   Control         22.3      22.2  0.141   1       0.0980   0.0693 
#> 14 <NA>   Control         22.1      NA   NA      NA      NA       NA      
#> 15 <NA>   Treatment       26.8      26.7  0.141   0.0442  0.00433  0.00306
#> 16 <NA>   Treatment       26.6      NA   NA      NA      NA       NA

# clean data
ct <- raw_data %>%
  fill(Target) %>%
  select(Target, Sample, `Cq value`) %>%
  setNames(c('gene', 'group', 'ct')) %>%
  group_by(gene, group) %>%
  mutate(row_id = row_number()) %>%
  spread(gene, ct) %>%
  ungroup() %>%
  select(-row_id, -group)

group <- rep(c('control', 'treatment'), each = 2)

# calculate a normalization factor
pcr_nf(ct,
       group_var = group,
       reference_group = 'control',
       reference_gene = paste0('REF', 1:3))
#> Joining, by = c("group", "gene")
#> # A tibble: 2 x 4
#>   group      mean     sd factor
#>   <chr>     <dbl>  <dbl> <chr> 
#> 1 control   1     0.0542 NF    
#> 2 treatment 0.595 0.0630 NF

# load and show normalization factor from the example
read_excel('clean_example.xlsx', sheet = 2)
#> # A tibble: 2 x 4
#>   `Normalization Factor (NF)`    NF `NF SD` `NF SE`
#>   <chr>                       <dbl>   <dbl>   <dbl>
#> 1 NF Control                  1      0.0542  0.0383
#> 2 NF Treatment                0.595  0.0630  0.0445

# calculate relative expression change
pcr_genorm(ct,
           group_var = group,
           reference_group = 'control',
           reference_gene = paste0('REF', 1:3))
#> Joining, by = c("group", "gene")
#> Joining, by = c("group", "gene")
#> Joining, by = "group"
#> # A tibble: 2 x 5
#> # Groups:   group, gene [?]
#>   group     gene  factor change  error
#>   <chr>     <chr> <chr>   <dbl>  <dbl>
#> 1 control   GOI1  NF     1      0.112 
#> 2 treatment GOI1  NF     0.0743 0.0107

# load and show relative expression from the example
read_excel('clean_example.xlsx', sheet = 3)
#> # A tibble: 2 x 4
#>   `Normalized Expression Level (NRQ)`    NRQ `NRQ SD` `NRQ SE`
#>   <chr>                                <dbl>    <dbl>    <dbl>
#> 1 GOInorm Control                     1        0.112   0.0792 
#> 2 GOInorm Treatment                   0.0743   0.0107  0.00758

Created on 2018-10-17 by the reprex package (v0.2.1)