source: trunk/brs2010p35/src/jag/distribution_plot.R @ 2

Last change on this file since 2 was 2, checked in by tim.te.beek@…, 8 years ago

Add revision 338 snapshot of previous SVN instance, omitting PLINK binaries

File size: 2.2 KB
Line 
1dist_plot<-function(emp_file,outfile,prefix,pheno){
2  #################################### DISTRIBUTION PLOTS  #####################################
3  ## Met de onderstaande code kan je voor iedere groep die je hebt getest de distributie van de permutaties plotten.
4  ## Ik heb deze code aangepast op de JAG output files
5  ##  Ik heb de code getest op 100 permutaties. Input data staat in 'input_output_plots'.
6
7  ## READ IN FILE WITH EMP PVALUES AND GET ORDER FROM LOWEST Pemp to HIGHEST Pemp
8  distr.sumlogs <- read.table(emp_file, sep="\t", header=TRUE, colClasses='character')
9  ordered.sumlogs <- order(as.numeric(distr.sumlogs[,4]))
10
11  ## GET PERMUTED SUMLOGS (NOTE: de naamgeving van deze file is nu variabel!)
12  permuted.sumlogs <- read.delim(outfile, sep = "\t", colClasses ='character')
13
14  nr.groups <- ncol(permuted.sumlogs)-1
15  pdf.title <- paste(prefix,"distribution_sumlogs.",pheno,".pdf", sep = '')
16  pdf(file=pdf.title)
17 
18  for(i in 1:length(ordered.sumlogs))
19  {     
20          # OPEN PNG > note that name of file is not specific now!
21
22          print(paste("plotting group ", distr.sumlogs[ordered.sumlogs[i],1], sep='')) 
23         
24          # Generate 100 (or another number) equal bins and take the min and max of permuted data as extremes of bins 
25          seq.data <- seq(from = min(as.numeric(permuted.sumlogs[,ordered.sumlogs[i]])),to = max(as.numeric(permuted.sumlogs[,ordered.sumlogs[i]])), length.out = 100)
26         
27          values <- c(as.numeric(permuted.sumlogs[,ordered.sumlogs[i]]),distr.sumlogs[ordered.sumlogs[i],2])
28          min <- min(as.numeric(values))/100*95
29          max <- max(as.numeric(values))/100*105
30         
31          # plot histogram of permuted distribution
32          hist(as.numeric(permuted.sumlogs[,ordered.sumlogs[i]]), seq.data, main = distr.sumlogs[ordered.sumlogs[i],1], xlab = "sumlog values", xlim = c(min, max), cex.main=1.8, cex.lab=1.5, col="lightgrey")
33          abline(v =distr.sumlogs[ordered.sumlogs[i],2], col="red", lty=2, cex=2, lwd=2) 
34  }
35  dev.off()
36}
37
38#read all arguments given by python
39#DO NOT REMOVE!
40args=(commandArgs(TRUE))
41if(length(args)==0){
42    print("No arguments supplied.")
43}else{
44    for(i in 1:length(args)){
45         eval(parse(text=args[[i]]))
46    }
47}
48dist_plot(emp_file,out_file,prefix,pheno)
Note: See TracBrowser for help on using the repository browser.