-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathload-frs.R
58 lines (49 loc) · 1.9 KB
/
load-frs.R
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
##
## load an FR file along with its parameters
## [nodes] size support case ctrl OR p pri
## [7,15,33,48,93] 159 354 49 305 0.161 8.59E-48 794
## for plotting and analysis
library(dplyr)
prefix = readline(prompt="FR file prefix (ex. HTT.400-0.5-1): ")
frFilename = paste(prefix, ".frs.txt", sep="")
## read the FRs table
## [nodes] size support case ctrl OR p pri
frs = read.table(frFilename, header=TRUE, stringsAsFactors=FALSE)
## remove extra header lines when multiple files have been catted
frs = frs[frs$nodes!="nodes",]
## convert to numeric
frs$size = as.numeric(frs$size)
frs$support = as.numeric(frs$support)
frs$case = as.numeric(frs$case)
frs$ctrl = as.numeric(frs$ctrl)
frs$OR = as.numeric(frs$OR)
frs$p = as.numeric(frs$p)
frs$pri = as.numeric(frs$pri)
## remove dupes
frs = distinct(frs)
## rownames(frs) = frs$nodes
## frs$nodes = NULL
## frs$size = as.numeric(frs$size)
## frs$support = as.numeric(frs$support)
## frs$case = as.numeric(frs$case)
## frs$ctrl = as.numeric(frs$ctrl)
## frs$OR = as.numeric(frs$OR)
## frs$p = as.numeric(frs$p)
## frs$pri = as.numeric(frs$pri)
prefix.parts = strsplit(prefix, "-", fixed=TRUE)
graphPrefix = prefix.parts[[1]][1]
alpha = as.numeric(prefix.parts[[1]][2])
kappa = as.numeric(prefix.parts[[1]][3])
## load the parameters from the params.txt file
source("load-params.R")
## label counts (if exists)
labelFile = paste(graphPrefix,".labelcounts.txt",sep="")
labelsExist = file.exists(labelFile)
if (labelsExist) {
labelCounts = read.delim(file=labelFile, header=FALSE, stringsAsFactors=FALSE, row.names=1)
colnames(labelCounts) = c("count")
}
## nodes
## 1 rs114039523 6 29910286 29910286 T/T 0.6643768400392541
## 4 rs114039523 6 29910286 29910286 ./. 8.92140244446427E-5
nodes = read.table(file=paste(graphName,"nodes","txt",sep="."), row.names=1, col.names=c("node","rs","chr","start","end","genotype","p"))