-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathhapConten.R
77 lines (67 loc) · 3.4 KB
/
hapConten.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
hapConten <- function(data=NULL, min.freq=50, max.freq=2000, pop.list=1, snpSites = NULL, ...) {
if (!is.null(snpSites) && length(snpSites)>=1) {
data <- data[rownames(data) %in% snpSites, ]
}
dat.mat <- data
dat.mat[is.na(dat.mat)] <- "-"
dat.mat <- t(dat.mat)
dat.mat.bin <- ape::as.DNAbin(dat.mat)
h <- pegas::haplotype(dat.mat.bin)
h <- pegas:::sort.haplotype(h)
h.sub <- pegas:::subset.haplotype(h, minfreq = min.freq, maxfreq = max.freq)
dimnames(h.sub)[[1]] <- as.character(as.roman(1:nrow(h.sub)))
h.sub.seq <- sapply(1:nrow(h.sub), function(x){
paste0(h.sub[x, ], sep="", collapse = "")
})
h.sub.seq <- toupper(h.sub.seq)
names(h.sub.seq) <- rownames(h.sub)
if (pop.list == 3) {
acc.info$Ecotype[acc.info$Ecotype %in% c("Ind_Int", "IndI", "indica", "IndII")] <- "Ind"
acc.info$Ecotype[acc.info$Ecotype %in% c("Jap_Int", "TeJ", "TrJ")] <- "Jap"
acc.info$Ecotype[acc.info$Ecotype %in% c("Int", "VI")] <- "Other"
acc.info$Ecotype[acc.info$Ecotype %in% c("Or-I", "Or-II", "Or-III")] <- "Wild"
} else if (pop.list == 2) {
acc.info$Ecotype[acc.info$Ecotype %in% c("Ind_Int", "IndI", "indica", "IndII")] <- "Ind"
acc.info$Ecotype[acc.info$Ecotype == "Jap_Int"] <- "Jap"
acc.info$Ecotype[acc.info$Ecotype %in% c("Int", "VI")] <- "Other"
acc.info$Ecotype[acc.info$Ecotype %in% c("Or-I", "Or-II", "Or-III")] <- "Wild"
} else if (pop.list == 1) {
acc.info$Ecotype[acc.info$Ecotype %in% c("Ind_Int", "IndI", "indica", "IndII")] <- "Ind"
acc.info$Ecotype[acc.info$Ecotype == "Jap_Int"] <- "Jap"
acc.info$Ecotype[acc.info$Ecotype %in% c("Int", "VI")] <- "Other"
}
net.pie <- lapply(1:length(attr(h.sub, "index")), function(x){
x.acc <- rownames(dat.mat)[attr(h.sub, "index")[[x]]]
x.acc.info <- acc.info$Ecotype[acc.info$ID %in% x.acc]
if (pop.list == 3) {
return(c(length(which(x.acc.info=="Ind")), length(which(x.acc.info=="Jap")),
length(which(x.acc.info=="Aus")), length(which(x.acc.info=="Wild")),
length(which(x.acc.info=="Other"))))
} else if (pop.list == 2) {
return(c(length(which(x.acc.info=="Ind")), length(which(x.acc.info=="Jap")),
length(which(x.acc.info=="TeJ")), length(which(x.acc.info=="TrJ")),
length(which(x.acc.info=="Aus")), length(which(x.acc.info=="Wild")),
length(which(x.acc.info=="Other"))))
} else if (pop.list == 1) {
return(c(length(which(x.acc.info=="Ind")), length(which(x.acc.info=="Jap")),
length(which(x.acc.info=="TeJ")), length(which(x.acc.info=="TrJ")),
length(which(x.acc.info=="Aus")), length(which(x.acc.info=="Or-I")),
length(which(x.acc.info=="Or-II")), length(which(x.acc.info=="Or-III")),
length(which(x.acc.info=="Other"))))
}
})
net.pie.df <- do.call(rbind, net.pie)
if (pop.list == 3) {
colnames(net.pie.df) <- c("Ind", "Jap", "Aus", "wild", "Other")
} else if (pop.list == 2) {
colnames(net.pie.df) <- c("Ind", "Jap", "TeJ", "TrJ", "Aus", "wild", "Other")
} else if (pop.list == 1) {
colnames(net.pie.df) <- c("Ind", "Jap", "TeJ", "TrJ", "Aus", "Or-I", "Or-II", "Or-III", "Other")
}
h.sub.name <- lapply(1:length(attr(h.sub, "index")), function(x){
x.acc <- rownames(dat.mat)[attr(h.sub, "index")[[x]]]
return(x.acc)
})
names(h.sub.name) <- names(h.sub.seq)
return(list(h.sub.seq, h.sub.name, net.pie.df))
}