May 11, 2024

R code to identify the FDR=0.05 equivalent P-value (for plotting purposes)

For Manhattan plots and other plots where you're plotting -log10(P-values) or untransformed P-values, you sometimes want to draw line to identify when FDR=0.05 is reached. But where do you draw the line? Below is R code that identifies the nearest P-value to FDR=0.05. You can then use variable FDR5.equiv.P to create a line on your Manhattan plot.

Code assumes you have a data frame with two columns, "Pval" and "FDR".


## identify FDR equivalent in Pvalue-------------

## sort Pvalues smaller-to-larger

  data = data[order(data$Pval), ]

  head(data$Pval)

  

## calculate absolute difference between FDR and 0.05

  data$fdr_diff = abs(data$FDR - 0.05)

  

## find row index with the smallest absolute value

  nearest_row = which.min(data$fdr_diff)

  

## get the equivalent Pvalue for FDR=0.05

  data[c(nearest_row+1, nearest_row, nearest_row-1) , c("Pval", "FDR")]

  FDR5.equiv.P = mean(data[c(nearest_row+1, nearest_row), "Pval"]); FDR5.equiv.P


NOTE: this method uses an average of the two Pvalues closest to FDR=0.05, which is fine for drawing a line. If you're using this to create a threshold for point formatting, then you should instead get the first Pvalue match that is below FDR=0.05 by comparing nearest_row-1, nearest_row, and nearest_row+1.

No comments:

Post a Comment

Bookmarks: single cell RNA-seq tutorials and tools

These are my bookmarks for single cell transcriptomics resources and tutorials. scRNA-seq introductions How to make R obj...