[Genabel-commits] r2048 - tutorials/GenABEL_general

[noreply at r-forge.r-project.org]

Author: lckarssen
Date: 2016-02-03 14:49:45 +0100 (Wed, 03 Feb 2016)
New Revision: 2048

Modified:
   tutorials/GenABEL_general/GenABEL-tutorial.Rnw
   tutorials/GenABEL_general/ImputedDataAnalysis.Rnw
Log:
Some updates to the Imputed Data Analysis chapter of the GenABEL tutorial

Based on feedback during a course I was teaching.


Modified: tutorials/GenABEL_general/GenABEL-tutorial.Rnw
===================================================================
--- tutorials/GenABEL_general/GenABEL-tutorial.Rnw	2016-02-03 13:24:02 UTC (rev 2047)
+++ tutorials/GenABEL_general/GenABEL-tutorial.Rnw	2016-02-03 13:49:45 UTC (rev 2048)
@@ -81,7 +81,7 @@
 \newcommand{\GA}{\texttt{GenABEL-package}}
 \newcommand{\MA}{\texttt{MetABEL-package}}
 \newcommand{\MxA}{\texttt{MixABEL-package}}
-\newcommand{\PA}{\texttt{ProbABEL-package}}
+\newcommand{\PA}{\texttt{ProbABEL}}
 \newcommand{\DA}{\texttt{DatABEL-package}}
 
 % Commonly used abbreviations

Modified: tutorials/GenABEL_general/ImputedDataAnalysis.Rnw
===================================================================
--- tutorials/GenABEL_general/ImputedDataAnalysis.Rnw	2016-02-03 13:24:02 UTC (rev 2047)
+++ tutorials/GenABEL_general/ImputedDataAnalysis.Rnw	2016-02-03 13:49:45 UTC (rev 2048)
@@ -116,24 +116,32 @@
 write.table(data.frame(id=idnames(df500), rcT=rcT),
        file="rcT.PHE", quote=FALSE, row.names=FALSE)
 @
-next, try the command '\texttt{system("head rcT.PHE")}' to check the few first
-lines of the file.
+next, try the command \texttt{system("head rcT.PHE")} to check the
+few first lines of the file\footnote{The R command
+  \texttt{system("some command")} runs \texttt{some command} on the
+  Linux (or Windows or Apple OS X) command line. So in this case we
+  run the Linux command \texttt{|head} on the file \texttt{rcT.PHE}
+  that you just created.}.
 
-At this moment, leave \texttt{R} (or, rather, start new console!),
-copy the mach-files to the working directory with
+Now that the file with phenotype data is prepared, open a new terminal
+or SSH connection to the server, go to the directory where you started
+your R session and run the following commands on the Linux command
+line\footnote{The \texttt{yourname at server>} indicates the prompt of
+  the Linux command line. Only type the commands that follow it.}.
+
 \begin{verbatim}
 yourname at server> cp RData/mach1* .
 \end{verbatim}
 and run \PA{} analysis with
 \begin{verbatim}
-yourname at server> palinear --pheno rcT.PHE --info mach1.out.mlinfo /
+yourname at server>  palinear --pheno rcT.PHE --info mach1.out.mlinfo /
                  --dose mach1.mldose.fvi
 \end{verbatim}
 Do not forget to check that you start the analysis in right directory, i.e.
 all files (\texttt{rcT.PHE}, \texttt{mach1.out.mlinfo},
 \texttt{mach1.mldose.fvi}, \texttt{mach1.mldose.fvd})
 are present in the working directory (use the '\texttt{ls}' command
-from the console).
+from the Linux command line).
 
 <<echo=false,hide=true>>=
 system("cp RData/mach1* .")
@@ -151,13 +159,13 @@
 we need to compute it:
 <<>>=
 qtsPal$Chisq <- (qtsPal$beta/qtsPal$sebeta)^2
-qtsPal$"P-value" <- pchisq(qtsPal$Chisq,1,low=F)
+qtsPal$"P-value" <- pchisq(qtsPal$Chisq, 1, low=F)
 qtsPal[1:5,]
 @
 
 Let us have a look at the top 20 associated SNPs:
 <<>>=
-qtsPal[order(qtsPal$Chisq,decreasing=T)[1:20],]
+qtsPal[order(qtsPal$Chisq, decreasing=T)[1:20],]
 @
 and produce the plot of the results with
 <<>>=
@@ -169,13 +177,13 @@
 To compare with the results obtained using 500 directly typed
 SNPs only, we can add the points with
 <<>>=
-points(map(qts500),-log10(qts500[,"P1df"]),col="red",pch=19,cex=0.5)
+points(map(qts500), -log10(qts500[,"P1df"]), col="red", pch=19, cex=0.5)
 @
 \begin{figure}[t]
 \centering
 <<fig=true,echo=false,results=hide>>=
 plot(map5k,-log10(qtsPal[,"P-value"]))
-points(map(qts500),-log10(qts500[,"P1df"]),col="red",pch=19,cex=0.5)
+points(map(qts500), -log10(qts500[,"P1df"]), col="red", pch=19, cex=0.5)
 abline(h=-log10(5e-8))
 @
 \caption{
@@ -188,9 +196,9 @@
 Now, we can also check different characteristics of
 the SNPs, which are claimed to be significant in our analysis:
 <<>>=
-gwsSnpsImp <- qtsPal$name[which(qtsPal[,"P-value"]<=5e-8)]
+gwsSnpsImp <- qtsPal$name[which(qtsPal[,"P-value"] <= 5e-8)]
 gwsSnpsImp
-mlinfo <- read.table("mach1.out.mlinfo",head=T,strings=F)
+mlinfo <- read.table("mach1.out.mlinfo", head=T, strings=F)
 mlinfo[mlinfo$SNP %in% gwsSnpsImp,]
 @
 
@@ -206,15 +214,15 @@
 \begin{Answer}
 Here is the sequence of commands leading you to the answer:
 <<>>=
-qts5k <- mlreg(rcT~1,df5k)
-bestHits5k <- descriptives.scan(qts5k,top=20)
+qts5k <- mlreg(rcT~1, df5k)
+bestHits5k <- descriptives.scan(qts5k, top=20)
 bestHits5k
 plot(qts5k)
 abline(h=-log10(5e-8))
-gwsSnps5k <- rownames(bestHits5k)[bestHits5k$P1df<=6e-8]
+gwsSnps5k <- rownames(bestHits5k)[bestHits5k$P1df <= 6e-8]
 summary(gtdata(df5k[,gwsSnps5k ]))
 
-plot(-log10(qts5k[,"P1df"]),-log10(qtsPal[,"P-value"]))
+plot(-log10(qts5k[,"P1df"]), -log10(qtsPal[,"P-value"]))
 abline(h=-log10(5e-8))
 abline(v=-log10(5e-8))
 
@@ -238,7 +246,7 @@
 \begin{figure}[t]
 \centering
 <<fig=true,echo=false,results=hide>>=
-plot(-log10(qts5k[,"P1df"]),-log10(qtsPal[,"P-value"]))
+plot(-log10(qts5k[,"P1df"]), -log10(qtsPal[,"P-value"]))
 abline(h=-log10(5e-8))
 abline(v=-log10(5e-8))
 @