sstar_scripts/MigrationEffectOnAdmixture_TPR_FDR.Rmd at master · abwolf/sstar_scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
---
title: "Untitled"
author: "Aaron Wolf"
date: "10/3/2018"
output: html_document
---

```{r}
library(data.table)
library(ggplot2)
library(dplyr)
library(tidyr)
library(grid)
library(gridExtra)
library(stringr)
```


```{r}

mdls <- list.files('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Tenn/', pattern="n1_0.05")

TPR_FDR <- data.table(NULL)

for (mdl in mdls){
    print(mdl)

    infile <- paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Tenn/',mdl)

    chr_list <- fread(paste0(infile,'/Tenn.chr_list'), col.names = c('chr'), colClasses = c('numeric'))

    dt.TPR <- fread(paste0('cat ',infile,'/TPR_FDR.txt | grep TPR'), select=c(1,4), col.names=c('chr','TPR'))
    dt.FDR <- fread(paste0('cat ',infile,'/TPR_FDR.txt | grep FDR'), select=c(1,4), col.names=c('chr','FDR'))


    dt.TPR %<>% rowwise() %>% mutate(chr=as.numeric(strsplit(chr,split = ':')[[1]][[2]])) %>% as.data.table()
    dt.TPR %<>% rowwise() %>% mutate(TPR=as.numeric(strsplit(TPR,split = ':')[[1]][[2]])) %>% as.data.table()

    dt.FDR %<>% rowwise() %>% mutate(chr=as.numeric(strsplit(chr,split = ':')[[1]][[2]])) %>% as.data.table()
    dt.FDR %<>% rowwise() %>% mutate(FDR=as.numeric(strsplit(FDR,split = ':')[[1]][[2]])) %>% as.data.table()

    dt <- left_join(chr_list,dt.TPR, by='chr') %>% left_join(dt.FDR, by='chr') %>% as.data.table()

    mean.TPR <- round(t.test(dt$TPR,conf.level = 0.95)$estimate[[1]],digits = 4)
    TPR.CI.lwr <- round(t.test(dt$TPR,conf.level = 0.95)$conf.int[[1]],4)
    TPR.CI.upr <- round(t.test(dt$TPR,conf.level = 0.95)$conf.int[[2]],4)

    mean.FDR <- round(t.test(dt$FDR,conf.level = 0.95)$estimate[[1]],4)
    FDR.CI.lwr <- round(t.test(dt$FDR,conf.level = 0.95)$conf.int[[1]],4)
    FDR.CI.upr <- round(t.test(dt$FDR,conf.level = 0.95)$conf.int[[2]],4)
############
############
    dt.TPR.30kb <- fread(paste0('cat ',infile,'/TPR_FDR.30kb.txt | grep TPR'), select=c(1,4), col.names=c('chr','TPR'))
    dt.FDR.30kb <- fread(paste0('cat ',infile,'/TPR_FDR.30kb.txt | grep FDR'), select=c(1,4), col.names=c('chr','FDR'))


    dt.TPR.30kb %<>% rowwise() %>% mutate(chr=as.numeric(strsplit(chr,split = ':')[[1]][[2]])) %>% as.data.table()
    dt.TPR.30kb %<>% rowwise() %>% mutate(TPR=as.numeric(strsplit(TPR,split = ':')[[1]][[2]])) %>% as.data.table()

    dt.FDR.30kb %<>% rowwise() %>% mutate(chr=as.numeric(strsplit(chr,split = ':')[[1]][[2]])) %>% as.data.table()
    dt.FDR.30kb %<>% rowwise() %>% mutate(FDR=as.numeric(strsplit(FDR,split = ':')[[1]][[2]])) %>% as.data.table()

    dt <- left_join(chr_list,dt.TPR.30kb, by='chr') %>% left_join(dt.FDR.30kb, by='chr') %>% as.data.table()

    mean.TPR.30kb <- round(t.test(dt$TPR,conf.level = 0.95)$estimate[[1]],digits = 4)
    TPR.CI.lwr.30kb <- round(t.test(dt$TPR,conf.level = 0.95)$conf.int[[1]],4)
    TPR.CI.upr.30kb <- round(t.test(dt$TPR,conf.level = 0.95)$conf.int[[2]],4)

    mean.FDR.30kb <- round(t.test(dt$FDR,conf.level = 0.95)$estimate[[1]],4)
    FDR.CI.lwr.30kb <- round(t.test(dt$FDR,conf.level = 0.95)$conf.int[[1]],4)
    FDR.CI.upr.30kb <- round(t.test(dt$FDR,conf.level = 0.95)$conf.int[[2]],4)
###########
###########

    row <- data.table(mdl, mean.TPR, mean.TPR.30kb, TPR.CI.lwr, TPR.CI.upr, TPR.CI.lwr.30kb, TPR.CI.upr.30kb, mean.FDR, mean.FDR.30kb, FDR.CI.lwr, FDR.CI.upr, FDR.CI.lwr.30kb, FDR.CI.upr.30kb)
    TPR_FDR <- rbind(TPR_FDR, row)
}

setnames(TPR_FDR, c('mdl', 'mean.TPR', 'mean.TPR.30kb', 'TPR.CI.lwr', 'TPR.CI.upr', 'TPR.CI.lwr.30kb', 'TPR.CI.upr.30kb','mean.FDR', 'mean.FDR.30kb', 'FDR.CI.lwr', 'FDR.CI.upr', 'FDR.CI.lwr.30kb', 'FDR.CI.upr.30kb'))

```


```{r}
ggplot(data=TPR_FDR) + theme_bw() +
  geom_line(aes(x=as.factor(mdl), y=as.numeric(mean.TPR), color='mean.TPR',group=1)) +
  geom_line(aes(x=as.factor(mdl), y=as.numeric(mean.TPR.30kb), color='mean.TPR.30',group=1))
```


```{r}

write.table(x = TPR_FDR, file = "~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Tenn/TPR_FDR.txt",quote = FALSE, sep = "\t",row.names = FALSE, col.names = TRUE)

```