before_after.R 2.25 KB
Newer Older
mariefbourdon's avatar
mariefbourdon committed
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
##################################################################
########### before_after.R
########### R script for stuart article figures 
##################################################################

## Goal and data
###########################

# The goal of this script is to calculate est.map and permutations of cross_before et cross_after
# 
# It imports the csv files "cross_before.csv" and "cross_after.csv"

## Import packages
library(qtl)

##############################
## Before
##############################

## read cross
cross_before <- read.cross(format="csv",file="./cross_before.csv",
                           genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

## calc geno prob
cross_before <- calc.genoprob(cross_before, step=2.0, off.end=0.0,
                             error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")

## calculate est map
newmap_before <- est.map(cross=cross_before,error.prob=0.01)
save(newmap_before,file="./newmap_before.rda")


## calculate 1000p
before_1000p <- scanone(cross=cross_before, 
                       chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), 
                       pheno.col="Pheno", model="normal", method="em", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE) 
save(before_1000p,file="./before_1000p.rda")

##############################
## After
##############################

## read cross
cross_after <- read.cross(format="csv",file="./cross_after.csv",
                          genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

## calc geno prob
cross_after <- calc.genoprob(cross_after, step=2.0, off.end=0.0, 
                              error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")

## calculate est map
newmap_after <- est.map(cross=cross_after,error.prob=0.01)
save(newmap_after,file="./newmap_after.rda")


## calculate 1000p
after_1000p <- scanone(cross=cross_after, 
                        chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), 
                        pheno.col="Pheno", model="normal", method="em", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE) 
Marie Bourdon's avatar
Marie Bourdon committed
60
save(after_1000p,file="./after_1000p.rda")