correlation plot, FDR value

Question

[Rob]

I am using the following code to create a correlation plot between genes and a continuous variable. It gives me the correlation coefficient (R) and p-value on the plot. I want this to:

1. give me FDR values
2. How can I get R values and FDR values in a csv or excel file, Not just written on the image of plot?

code:

library(ggplot2)
library(plyr)
library(reshape2)
dat <- read.csv("Gene_variable_file.csv", check.name = FALSE, stringsAsFactors = FALSE, header = T)
df <- as.data.frame(dat)
# Basic scatter plot
ex <- melt(df, id.vars="continousVariable")

colnames(ex) <- c("continousVariable", "gene", "exprs")
#######################
  ####
  ggscatter(ex, x = "continousVariable", y = "exprs",
            add = "reg.line",                         
            add.params = list(color = "blue", fill = "lightgray"),
            color = "black", palette = "jco", fill = "lightgray",          
            #shape = "cyl",                            
            fullrange = TRUE,                       
            rug = TRUE, facet.by = "gene", cor.coef = T,
            title = "Sarco-Specific DEGs-SMA Correlation-FEMALE",
            conf.int = TRUE, 
            cor.coeff.args = list(),
            cor.method = "spearman",
            cor.coef.coord = c(NULL, NULL),
            cor.coef.size = 4,                               
  )+
    geom_vline(xintercept = 102, colour="red", linetype = "longdash")

This is how my data look like:

Answer 1

[i.sudbery]

This is the sort of thing the broom packages was created for. Look at grouping your DF by gene, and then using broom to run a linear model on each gene, and extract the coefficients and pvalues of the per gene fits. You can then use p.adjust(pvalues, method="BH") to calculate the FDR values from the p.values.

[Rob]

Thank you. It tried broom from this link: https://cran.r-project.org/web/packages/broom/vignettes/broom.html

1. But which method I should use?ttest?glm?lm? I used lm and I got all my p-values 1,1,1,1... which cannot be correct.

2. which one from the following code will be my data in broom code?"dat" or "ex"

 dat <- read.csv("Gene_variable_file.csv", check.name = FALSE, stringsAsFactors = FALSE, header = T)
    df <- as.data.frame(dat)


ex <- melt(df, id.vars="continousVariable")

colnames(ex) <- c("continousVariable", "gene", "exprs")
#######################
  ####
  ggscatter(ex, x = "continousVariable", y = "exprs",
            add = "reg.line",                         
            add.params = list(color = "blue", fill = "lightgray"),
            color = "black", palette = "jco", fill = "lightgray",          
            #shape = "cyl",                            
            fullrange = TRUE,                       
            rug = TRUE, facet.by = "gene", cor.coef = T,
            title = "Sarco-Specific DEGs-SMA Correlation-FEMALE",
            conf.int = TRUE, 
            cor.coeff.args = list(),
            cor.method = "spearman",
            cor.coef.coord = c(NULL, NULL),
            cor.coef.size = 4,                               
  )+
    geom_vline(xintercept = 102, colour="red", linetype = "longdash")

I used ex as my data in broom

[i.sudbery]

By the looks of things, it would be ex. You'd need to group_by gene in some way. I've always used group_by, but this vignette suggests nesting https://cran.r-project.org/web/packages/broom/vignettes/broom_and_dplyr.html.

If you want to do pearson correlation, then lm should work fine, as the p-value on a linear model is identical to the p-value of a pearson's correlation. I'm not sure you can do spearman's that way though. You can use cor.test(method="spearman") though.

I agree that it is unlikely that all your p-values or 1. However, all your FDRs/adjusted p-values being 1 is entirely possible.