-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path.Rhistory
52 lines (52 loc) · 1.79 KB
/
.Rhistory
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
library(methylKit)
library(reshape2)
colnames(methpos.chick)
firstC<-which(colnames(methylation.Pos)=="numCs1$" )
firstC
length(columnsC )
methylation.Pos<-methylation[,c(1,3,columnsC,columnsT)]
rm(list=ls())
rm(list=ls())
?group_split
colnames(methylationData)
function_glmm <- function(D_analysis) {tryCatch({
meth_temp_final <- as.data.frame(methylationData)
CHR_POS <- as.character(meth_temp_final[1,2])
null <- glmer(cbind(numCs, numTs) ~ 1 + (1|indNumber), meth_temp_final,
family="binomial")
model <- glmer(cbind(numCs, numTs) ~ 1 + site + (1|indNumber), meth_temp_final,
family="binomial")
lrt <- anova(null, model)
pvalue=lrt$"Pr(>Chisq)"[2]
deviance <- lrt$deviance[2]
Df <- lrt$Df[2]
Chisq <- lrt$Chisq[2]
message<-data.frame(model@optinfo$conv$lme4)
mess1<-str_c(" ",message$messages[1])
mess2<-str_c(" ",message$messages[2])
mess3<-str_c(" ",message$messages[3])
mess1.1<-str_c(" ",message[1,1])
mess2.1<-str_c(" ",message[1,2])
mess3.1<-str_c(" ",message[1,3])
mess4.1<-str_c(" ",message[1,4])
message1.2 <- data.frame(model@optinfo$message)
mess1.2 <-str_c(" ",message1.2[1,1])
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
ratio <- Pearson.chisq/rdf
pval_disp <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
return(data.frame(CHR_POS=CHR_POS, df=Df, Chisq=Chisq, pval_lrt=pvalue, df_residual=rdf,
disp_ratio=ratio message1=mess1,
message2=mess2, message3=mess3, message1.1=mess1.1, message2.1=mess2.1,
message3.1=mess3.1, message4.1=mess4.1, message1.2=mess1.2))
}, error = function(e){cat("ERROR :", conditionMessage(e), "\n");print(CHR_POS)})
}
?split
coefsOut<-data.frame(matrix(ncol=6))
colnames(randomOut)<-paste(random[1,],random[4,],sep="_")
nullModel<-"1 + site + (1|Run)"
65/3
indNames
write.table(nLociData,"results/nLociPerInd.tsv")
max(nLociData)/100