/
hapest.caco.R
executable file
·78 lines (53 loc) · 2.22 KB
/
hapest.caco.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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
hapest.caco <-
function ( geno , trait , lim=0.05 ) {
ca <- trait==1
co <- trait==0
nca <- sum(ca)
nco <- sum(co)
ck <- 2
# infer haplotypes in pooled sample (cases and controls together)
hpool <- itegeppXXR ( geno , des = 0 , lim = lim )
# failed inferring haplotypes ?
if ( all( is.na( hpool$hap )) ) {
return ( list ( haplotypes=NA , desres=NA ) )
}
hapcaco.pool <- colSums ( hpool$desres ) / (ck*dim(hpool$desres)[1])
hapcaco.ca <- colSums ( hpool$desres [ca,] ) / (ck*nca)
hapcaco.co <- colSums ( hpool$desres [co,] ) / (ck*nco)
hapcaco.pool <- hapcaco.pool[names(hapcaco.pool)!="R",drop=F]
hapcaco.ca <- hapcaco.ca [names(hapcaco.ca)!="R" ,drop=F]
hapcaco.co <- hapcaco.co [names(hapcaco.co)!="R" ,drop=F]
hapcaco <- data.frame (
as.character(names(hapcaco.pool)),as.numeric(hapcaco.pool),
as.numeric(hapcaco.ca),
as.numeric(hapcaco.co) , stringsAsFactors=F)
desres <- hpool$desres
hapcaco[,1] <- as.character(hapcaco[,1])
hapcaco[,2] <- as.numeric(as.character(hapcaco[,2]))
hapcaco[,3] <- as.numeric(as.character(hapcaco[,3]))
hapcaco[,4] <- as.numeric(as.character(hapcaco[,4]))
colnames (hapcaco) <- c("Hap","Pool","Case","Control")
# Threshold for protecting haplotypes
bin <- ( hapcaco$Pool >= lim )
hapcaco <- hapcaco[bin,,drop=F]
if ( dim(hapcaco)[1]==0 ) {
print ( paste("All inferred haplotypes or haplotype pairs less than ",
"lim = ",lim,sep="" ) )
return ( list ( haplotypes=NA , desres=NA ) )
}
hapsi <- which(!is.na(match(colnames(desres),hapcaco[,1])))
desres <- as.matrix(desres[ , hapsi, drop=FALSE ])
# desres <- desres/2
if ( !any(colnames(desres)=="R") ) {
desres <- cbind(desres, R=2-rowSums(desres))
}
desres <- as.matrix(desres)
# round to 6 decimal places
hapcaco[,2] <- round ( hapcaco[,2] , 6 )
hapcaco[,3] <- round ( hapcaco[,3] , 6 )
hapcaco[,4] <- round ( hapcaco[,4] , 6 )
hi <- match ( hapcaco[,1] , colnames(desres) )
hapcaco <- hapcaco [ !is.na(hi) , ,drop=FALSE ]
rownames (hapcaco) <- NULL
return ( list ( haplotypes=hapcaco , desres=desres ) )
}