Body: <p>I have done the limma association analysis using below model design but not sure how to  evaluate Potential systematic biases using Q-Q plots and the genomic inflation factor (k). Could you please help with this?</p>
  <pre class="lang-r prettyprint-override"><code>#model matrix
  var&lt;-model.matrix(~BMI + as.factor(Gender) + Age +CD8T +CD4T +NK + Bcell +Mono ,data=targets2)
  fit&lt;-lmFit(mval,var)
  fit2&lt;-eBayes(fit,trend=TRUE, robust=TRUE)
  probe&lt;-topTable(fit2,adjust=&quot;BH&quot;,coef=2,num=Inf)
  
  write.table(probe, file = &quot;limma.output.txt&quot;)
  data &lt;- read.delim(&quot;limma.output.txt&quot;)
  head(data)
  
                 logFC       AveExpr.    t         P.Value  adj.P.Val         B
  cg03546163  0.41018776  -0.8013218  7.85104192  1.13E-12    8.65E-07    12.3051291
  cg17944885  -0.4884859  2.31974279  -5.7319645  6.33E-08    0.02431349  5.15237109
  </code></pre>
  <pre><code>#For Genomic inflation factor
  
  # For p-values, calculate chi-squared statistic
  chisq &lt;- qchisq(1-data$P.Value,1)
  
  #Calculate lambda gc (λgc)
  median(chisq)/qchisq(0.5,1)
  
  #QQ PLOT
  
  library(gaston)
  pdf(&quot;qq.plot.pdf&quot;, width = 12, height = 12)
  qqplot.pvalues(data$P.Value)
  dev.off()
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-18T18:27:35.283Z
Comments.Text: What have you already tried? see e.g. [here](https://rdrr.io/cran/GenABEL/man/estlambda.html) and [here](https://stats.stackexchange.com/questions/230226/how-to-manually-make-a-qq-plot-of-p-value-against-uniform-distribution) that i found by google quickly.
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-20T08:19:20.377Z
Comments.Text: thank you so much @MaximilianPress, I am not sure do we need to filter limma output based on adj.P.val before running QQ plot or not? and second I am not sure how to proceed after this step.
Comments.UserDisplayName: bioinfonext
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-20T15:22:35.287Z
Comments.Text: `adj.P.val` computes q-values (or applies some other correction), so its output are not really p-values. You will need to plot the p-values themselves, not a transformation.
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-20T15:45:30.113Z
Comments.Text: thanks@MaximilianPress, I updated code above for QQ plot and genomic inflation factor. I hope its correct?
Comments.UserDisplayName: bioinfonext
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-20T16:21:18.083Z
Comments.Text: I don't see any obvious issues at this point but I'll let other possibly more expert people comment.
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-29T16:01:28.750Z
Comments.Text: You do not need to filter the limma output based on adj.P.Val before generating the Q-Q plot / getting the inflation factor. Your code looks good to me. Is there anything unusual about the Q-Q plot or inflation factor that makes you think something went wrong?
Comments.UserDisplayName: BigMistake
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-11-03T09:43:30.313Z
Comments.Text: Thanks for all help.
Comments.UserDisplayName: bioinfonext
Tags: <r><limma>
Title: Q-Q plots and the genomic inflation factor (k)
ThreadText: # Q-Q plots and the genomic inflation factor (k)
  
  (Asked by: username_0 on 2023-10-13)
  
  I have done the limma association analysis using below model design but not sure how to evaluate Potential systematic biases using Q-Q plots and the genomic inflation factor (k). Could you please help with this?
  
  ```r
  #model matrix
  var<-model.matrix(~BMI + as.factor(Gender) + Age +CD8T +CD4T +NK + Bcell +Mono ,data=targets2)
  fit<-lmFit(mval,var)
  fit2<-eBayes(fit,trend=TRUE, robust=TRUE)
  probe<-topTable(fit2,adjust="BH",coef=2,num=Inf)
  
  write.table(probe, file = "limma.output.txt")
  data <- read.delim("limma.output.txt")
  head(data)
  
                 logFC       AveExpr.    t         P.Value  adj.P.Val         B
  cg03546163  0.41018776  -0.8013218  7.85104192  1.13E-12    8.65E-07    12.3051291
  cg17944885  -0.4884859  2.31974279  -5.7319645  6.33E-08    0.02431349  5.15237109
  
  ```
  
  ```
  #For Genomic inflation factor
  
  # For p-values, calculate chi-squared statistic
  chisq <- qchisq(1-data$P.Value,1)
  
  #Calculate lambda gc (λgc)
  median(chisq)/qchisq(0.5,1)
  
  #QQ PLOT
  
  library(gaston)
  pdf("qq.plot.pdf", width = 12, height = 12)
  qqplot.pvalues(data$P.Value)
  dev.off()
  
  ```
  
  **Tags:** r, limma
  
  > **Comments:**
  > **username_1**: What have you already tried? see e.g. [here](https://rdrr.io/cran/GenABEL/man/estlambda.html) and [here](https://stats.stackexchange.com/questions/230226/how-to-manually-make-a-qq-plot-of-p-value-against-uniform-distribution) that i found by google quickly.  
  > **username_0**: thank you so much @username_1, I am not sure do we need to filter limma output based on adj.P.val before running QQ plot or not? and second I am not sure how to proceed after this step.  
  > **username_1**: `adj.P.val` computes q-values (or applies some other correction), so its output are not really p-values. You will need to plot the p-values themselves, not a transformation.  
  > **username_0**: thanks@username_1, I updated code above for QQ plot and genomic inflation factor. I hope its correct?  
  > **username_1**: I don't see any obvious issues at this point but I'll let other possibly more expert people comment.  
  > **username_2**: You do not need to filter the limma output based on adj.P.Val before generating the Q-Q plot / getting the inflation factor. Your code looks good to me. Is there anything unusual about the Q-Q plot or inflation factor that makes you think something went wrong?  
  > **username_0**: Thanks for all help.  
  
  ---
  
  ## 0 answers
  
  