-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathExtractGenoFromCross.R
More file actions
65 lines (46 loc) · 1.52 KB
/
ExtractGenoFromCross.R
File metadata and controls
65 lines (46 loc) · 1.52 KB
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
# code written by Paul Tanger, Brook Moyers or John Lovell. Please cite appropriately.
ExtractGenoFromCross = function(qtlobject) {
# extract genotype data from object to use for joinmap, MSTmap, or our CS tool
genos = as.data.frame(qtlobject$geno[[1]]$data)
i = 1
for (i in 2:length(qtlobject$geno)) {
genos = cbind(as.matrix(genos), qtlobject$geno[[i]]$data)
}
# what about pull.geno? or write.cross
# change back to A and B
genos[genos == 1] = "A"
genos[genos == 2] = "B"
genos[is.na(genos)] = "-"
genos = as.data.frame(genos)
# add line names back
finalgenos = cbind(qtlobject$pheno[1], genos)
geno_joinmap = finalgenos
# transpose
geno.t = t(geno_joinmap)
# add col names
colnames(geno.t) <- geno_joinmap[,1]
# remove first row
geno.t = geno.t[-1,]
# make row.names into column
geno.t2 = as.data.frame(cbind(line = rownames(geno.t), geno.t))
# get map
themap = pull.map(qtlobject, as.table=T)
# merge map distances and geno data
genomap = merge(themap, geno.t2, by="row.names")
# sort
genomap2 = with(genomap, genomap[order(chr, pos), ])
# transpose
genomap = t(genomap2)
# add line col as "ID"
genomap <- cbind(ID = rownames(genomap), genomap)
# rename cols
colnames(genomap) <- as.character(genomap[1,])
# fix ID col
colnames(genomap)[1] = "ID"
# remove old row with col names
genomap = genomap[-c(1,4),]
# remove chrom and marker and cM from ID col
genomap[1,1] = ""
genomap[2,1] = ""
return(genomap)
}