cgi-bin/phenome/cvterm_qtl.r

   1
   2 ##SNOPSIS
   3
   4 ##A banch of r and r/qtl commands for reading input data,
   5 ##running qtl analysis, and writing output data. Input data are feed to R
   6 ##from and output data from R feed to a perl script called
   7 ##../phenome/population_indls.
   8 ##QTL mapping visualization is done by the perl script mentioned above.
   9 ## r/qtl has functions for data visualization but the perl script
  10 ## allows more control for publishing on the browser,
  11 ##cross-referencing to other SGN datasets and CONSISTENT formatting/appearance
  12
  13
  14 ##AUTHOR
  15 ## Isaak Y Tecle (iyt2@cornell.edu)
  16
  17
  18 options(echo = FALSE)
  19
  20 library(qtl)
  21
  22 allargs<-commandArgs()
  23 warning()
  24
  25 infile <- grep("infile_list", allargs, ignore.case=TRUE, perl=TRUE, value=TRUE)
  26
  27 outfile<-grep("outfile_list", allargs, ignore.case=TRUE, perl=TRUE, value=TRUE)
  28
  29 statfile<-grep("stat", allargs, ignore.case=TRUE, perl=TRUE, value=TRUE)
  30
  31
  32 ##### stat files
  33 statfiles<-scan(statfile, what="character")
  34 #print(statfiles)
  35
  36 ###### QTL mapping method ############
  37 qtlmethodfile<-grep("stat_qtl_method", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
  38 qtlmethod<-scan(qtlmethodfile, what="character", sep="\n")
  39 #print(qtlmethod)
  40 if (qtlmethod == "Maximum Likelihood") {
  41   qtlmethod<-c("em")
  42   #print(qtlmethod)
  43
  44 } else
  45
  46 if (qtlmethod == "Haley-knott Regression") {
  47   qtlmethod<-c("hk")
  48
  49 } else
  50
  51 if (qtlmethod == "Multiple Imputation") {
  52   qtlmethod<-c("imp")
  53
  54 }
  55
  56 ###### QTL model ############
  57 qtlmodelfile<-grep("stat_qtl_model", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
  58 #print(qtlmodelfile)
  59 qtlmodel<-scan(qtlmodelfile, what="character", sep="\n")
  60 #print(qtlmodel)
  61 if (qtlmodel == "Single-QTL Scan") {
  62   qtlmodel<-c("scanone")
  63
  64 } else
  65 if  (qtlmodel == "Two-QTL Scan") {
  66   qtlmodel<-c("scantwo")
  67 }
  68
  69 ###### permutation############
  70 userpermufile<-grep("stat_permu_test", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
  71 #print(userpermufile)
  72 userpermuvalue<-scan(userpermufile, what="numeric", dec = ".", sep="\n")
  73 #print(userpermuvalue)
  74 if (userpermuvalue == "None") {
  75   userpermuvalue<-c(0)
  76 }
  77 userpermuvalue<-as.numeric(userpermuvalue)
  78
  79 #####for test only
  80 #userpermuvalue<-c(0)
  81
  82 ######genome step size############
  83 stepsizefile<-grep("stat_step_size", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
  84 stepsize<-scan(stepsizefile, what="numeric", dec = ".", sep="\n")
  85 stepsize<-as.numeric(stepsize)
  86
  87
  88 ######genotype calculation method############
  89 genoprobmethodfile<-grep("stat_prob_method", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
  90 genoprobmethod<-scan(genoprobmethodfile, what="character", dec = ".", sep="\n")
  91
  92
  93 ########No. of draws for sim.geno method###########
  94 drawsnofile<-c()
  95 if (is.logical(grep("stat_no_draws", statfiles))==TRUE) {
  96   drawsnofile<-(grep("stat_no_drwas", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE))
  97 }
  98
  99 if (is.logical(drawsnofile) ==TRUE) {
 100   drawsno<-scan(drawsnofile, what="numeric", dec = ".", sep="\n")
 101   drawsno<-as.numeric(drawsno)
 102 }
 103 ########significance level for genotype probablity calculation###########
 104 genoproblevelfile<-grep("stat_prob_level", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
 105 genoproblevel<-scan(genoproblevelfile, what="numeric", dec = ".", sep="\n")
 106 genoproblevel<-as.numeric(genoproblevel)
 107
 108
 109 ########significance level for permutation test###########
 110 permuproblevelfile<-grep("stat_permu_level", statfiles, ignore.case=TRUE, fixed = FALSE, value=TRUE)
 111 permuproblevel<-scan(permuproblevelfile, what="numeric", dec = ".", sep="\n")
 112 permuproblevel<-as.numeric(permuproblevel)
 113
 114
 115 #########
 116 infile<-scan(file=infile,  what="character")#container for the ff
 117
 118 cvtermfile<-infile[1]#file that contains the cvtername
 119 popid<-infile[2]#population dataset identifier
 120 genodata<-infile[3] #file name for genotype dataset
 121 phenodata<-infile[4] #file name for phenotype dataset
 122 permufile<-infile[5]
 123 crossfile<-infile[6]
 124
 125 print(crossfile)
 126 cross<-scan(crossfile, what="character", sep="\n")
 127 #print(cross)
 128
 129 popdata<-c()
 130 if (cross == "f2") {
 131   popdata<- read.cross("csvs", genfile=genodata, phefile=phenodata, na.strings=c("NA"), genotypes=c("1", "2", "3", "4", "5"), estimate.map=TRUE, convertXdata=TRUE)
 132   popdata<-jittermap(popdata)
 133 } else
 134 if (cross == "bc") {
 135   popdata<- read.cross("csvs", genfile=genodata, phefile=phenodata, na.strings=c("NA"), genotypes=c("1", "2"), estimate.map=TRUE, convertXdata=TRUE)
 136   popdata<-jittermap(popdata)
 137 }
 138
 139
 140
 141 if (genoprobmethod == "Calculate") {
 142   popdata<-calc.genoprob(popdata, step=stepsize, error.prob=genoproblevel)
 143   #calculates the qtl genotype probablity at the specififed step size and probability level
 144 } else
 145 if (genoprobmethod == "Simulate") {
 146   popdata<-sim.genoprob(popdata, n.draws= drawsno, step=stepsize, error.prob=genoproblevel, stepwidth="fixed")
 147   #calculates the qtl genotype probablity at the specififed step size
 148 }
 149
 150
 151
 152 cvterm<-scan(file=cvtermfile, what="character") #reads the cvterm
 153 cv<-find.pheno(popdata, cvterm)#returns the col no. of the cvterm
 154
 155
 156 permuvalues<-scan(file=permufile, what="character")
 157
 158 permuvalue1<-permuvalues[1]
 159 permuvalue2<-permuvalues[2]
 160 if ((is.logical(permuvalue1) == FALSE)) {
 161   if (qtlmodel == "scanone") {
 162     if (userpermuvalue == 0 ) {
 163       popdataperm<-scanone(popdata, pheno.col=cv, model="normal",  method=qtlmethod)
 164       #permu<-summary(popdataperm, alpha=c(0.05))
 165     } else
 166     if (userpermuvalue != 0) {
 167       popdataperm<-scanone(popdata, pheno.col=cv, model="normal", n.perm = userpermuvalue, method=qtlmethod)
 168       permu<-summary(popdataperm, alpha=permuproblevel)
 169       #print(permu)
 170     }
 171   }else
 172   if (qtlmodel == "scantwo") {
 173     if (userpermuvalue == 0 ) {
 174       popdataperm<-scantwo(popdata, pheno.col=cv, model="normal", method=qtlmethod)
 175       #permu<-summary(popdataperm)
 176     } else
 177     if (userpermuvalue != 0) {
 178       popdataperm<-scantwo(popdata, pheno.col=cv, model="normal", n.perm = userpermuvalue, method=qtlmethod)
 179       permu<-summary(popdataperm, alpha=permuproblevel)
 180       #print(permu)
 181
 182     }
 183   }
 184 }
 185
 186
 187 chrlist<-c("chr1")
 188
 189 for (no in 2:12) {
 190   chr<-paste("chr", no, sep="")
 191   chrlist<-append(chrlist, chr)
 192 }
 193
 194 chrdata<-paste(cvterm, popid, "chr1", sep="_")
 195 chrtest<-c("chr1")
 196 for (ch in chrlist) {
 197   if (ch=="chr1"){
 198     chrdata<-paste(cvterm, popid, ch, sep="_");
 199   }
 200
 201   else {
 202     n<-paste(cvterm, popid, ch, sep="_");
 203     chrdata<-append(chrdata, n)
 204
 205   }
 206 }
 207
 208 chrno<-1
 209
 210 datasummary<-c()
 211 confidenceints<-c()
 212
 213
 214 for (i in chrdata){
 215
 216
 217   filedata<-paste(cvterm, popid, chrno, sep="_");
 218   filedata<-paste(filedata,"txt", sep=".");
 219
 220   i<-scanone(popdata, chr=chrno, pheno.col=cv)
 221   position<-max(i, chr=chrno)
 222
 223   p<-position[2]
 224   p<-p[1, ]
 225   peakmarker<-find.marker(popdata, chr=chrno, pos=p)
 226   confidenceint<-bayesint(i, chr=chrno, prob=0.95, expandtomarkers=TRUE)
 227   confidenceint<-rownames(confidenceint)
 228   print (confidenceint)
 229   confidenceint<-c(chrno, confidenceint, peakmarker)
 230   print(confidenceint)
 231   if (chrno==1) {
 232     datasummary<-i
 233     confidenceints<-confidenceint
 234 }
 235   if (chrno > 1 ) {
 236     datasummary<-rbind(datasummary, i)
 237     confidenceints<-rbind(confidenceints, confidenceint)
 238   }
 239
 240 chrno<-chrno + 1;
 241
 242 }
 243 #print("cis")
 244 #print(confidenceints)
 245 #print("data summary")
 246 #print(datasummary)
 247
 248 outfiles<-scan(file=outfile,  what="character")
 249 qtlfile<-outfiles[1]
 250 confidenceintfile<-outfiles[2]
 251
 252 write.table(datasummary, file=qtlfile, sep="\t", col.names=NA, quote=FALSE, append=FALSE)
 253 write.table(confidenceints, file=confidenceintfile, sep="\t", col.names=NA, quote=FALSE, append=FALSE)
 254
 255 if (userpermuvalue != 0) {
 256   if ((is.logical(permuvalue1) == FALSE)) {
 257     write.table(permu, file=permufile, sep="\t", col.names=NA, quote=FALSE, append=FALSE)
 258   }
 259 }
 260
 261
 262 q(runLast = FALSE)