Results of Scantwo (with permutations) of different phenotypes is the Same

[Robert Mattera]

Hey Dr. Broman,

I am working on this RIL mapping population with 66 different phenotypes. I am at the step of looking at scantwo results and I came across this odd result and I am not sure if it is something innately wrong with my script that is causing this result. Basically scantwo works fine but when I run scantwo with operms and print my summary results into a .csv file all the phenotypes have the same two chromosomes and position interactions. 

Firstly, I ran scantwo with the following command:

> out2.em.all.ALL <- scantwo(qtl.finalmap, method="em", pheno.col = phenotypes)

I have attached a document that shows an snippet of the summary results from my scantwo in the first tab. 

I used the following script to export this into a .csv file:

> for (x in 1:66) {

phenotype.name <- phenames(qtl.finalmap)[x]  ## This allows me to pull the name of the phenotype and print it atop the summary results

write.table(phenotype.name, "summary.out2.em.all.ALL.csv", append = TRUE, sep = ",") ## Prints name of phenotype for results below to .csv file

temp <- summary(out2.em.all.ALL, threshold=c(6.08, 4.91, 4.71, 4.67, 2.84), lodcolumn = x) ## temp value serves as something that can be written over each time this loops 

write.table(temp, "summary.out2.em.all.ALL.csv", append = TRUE, sep = ",")} ## Pulls the current temp value and prints it to the csv file under the correct name 

It worked really well in exporting my results. The thresholds were established from a previous run on a smaller subset of phenotypes but I wanted to continue on and use operms2 to determine qtl interactions.

I used the following script to generate my operm2 results with the tictoc package to track how long it was taking. As you can imagine running permutations with the "em" method took a while so it was helpful to track how long it was taking while doing my for loops in smaller chunks:

library(tictoc) ## Load package
tic.clearlog() ## Clear tictoc log
tic() ## start timer
for (x in 52:66) {
  tic(x) ## start timer
  phenotype_name <- phenames(qtl.finalmap)[x]
  print(phenotype_name) # prints pheno name to visually track what stage analysis is
  operm2_em_temp <- scantwo(qtl.finalmap, method = "em", n.perm = 1000, pheno.col = x, n.cluster = 15) ## assigns scantwo results to operm2_em_temp
  temp <- summary(out2.em.all.ALL, perms=operm2_em_temp, alphas=c(0.05, 0.05, 0, 0.05, 0.05), pvalues=TRUE) ## Creates a temporary data value to hold the summary.scantwo results with the operm2 values
  write.table(phenotype_name, "summary.operm2.em.All.csv", append = TRUE, sep = ",") ## Prints name of phenotype for results below to .csv file
  write.table(temp, "summary.operm2.em.All.csv", append = TRUE, sep = ",") ## Prints summary results to summary .csv file
  toc(log = TRUE, quiet = TRUE) ## stop timer for that iteration
  }
log.txt <- tic.log(format = TRUE) ## Fetch results of toc()
log.lst <- tic.log(format = FALSE) ## Extract the list
timings <- unlist(lapply(log.lst, function(x) x$toc - x$tic)) ## Convert list elements to timings
mean(timings) ## average time for each iteration
writeLines(unlist(log.txt)) ## Print text output

In the attached file the second tab is the same select phenotypes but the results of the scantwo summary with the operm instead. I might be interpreting this incorrectly but from what i can see all of these have interactions between chr 2 and 9 and all have the same chromosome position as well. This doesn't seem probable (as this pattern continues with all the rest of the phenotypes). Like I mentioned maybe my for loop is wrong and acting on the incorrect phenotype and I am just missing something. 

Any help would be appreciated. This whole process took me quite some time and I am just confused as to how this happened. 

Thanks

Rob

[Karl Broman]

It seems like in this line within the loop over 52:66:

    temp <- summary(out2.em.all.ALL, perms=operm2_em_temp …

you are summarizing the same set of scantwo results, just with different permutations to establish the significance thresholds.

You either need to re-run scantwo without the permutations within that loop, as well, or you need to save the results from above, so that you can re-use them here.

I would recommend saving the full results and not just the summaries to files, for example using saveRDS().

That way you can analyze them at some leisure and not need to re-run the long-running calculations.

karl

[Robert Mattera]

Thanks for your quick response! Okay - so I think I now realize. I am missing the lodcolumn variable in that line of code which is so obvious -- I do not know how I missed it. 

What would saving the full results file look like? saveRDS(out2.em.all.ALL) ? This would generate a file containing all the scantwo results for all the phenotypes? And then using saveRDS(operm2_em_temp) to save the permutation results for each phenotype? Also in theory could I save both out2.em.all.ALL and operm2_em_temp to .csv files so I can just append each new iteration that I run to the .csv?

Thanks for your help

Rob

--
You received this message because you are subscribed to the Google Groups "R/qtl discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rqtl-disc+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/rqtl-disc/84da2fd2-165f-458e-a4ff-40e5734b4b11n%40googlegroups.com.

--

Robert Mattera III, MS

Pronouns: he/him

Research Project Coordinator and Ph.D Candidate

 Department of Plant Biology & Pathology, Foran Hall

 School of Environmental and Biological Sciences

 Rutgers, The State University of New Jersey

 59 Dudley Road

 New Brunswick, NJ 08901-8520

CoFounder and Head Mixologist

 3BR Distillery, LLC

 7 Main Street

 Keyport, NJ 07735

 www.3brdistillery.com

[Karl Broman]

You can use saveRDS() and readRDS() to save/read a single R object to/from an .rds file.

You can use save() and load() to save/load a set of R objects to/from an .RData file.

karl