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".

Vocabulary for reading DNA methylation papers

 If you are new to the field of DNA methylation, here is some introductory vocabulary.

• 27k array - short for Infinium HumanMethylation27 BeadChip (approx 27,000 methylation sites). An early DNA methylation arrays that samples CpG sites across the human genome. Not used anymore. 
• 450k array - short for Illumina HumanMethylation450K BeadChip (>450,000 methylation sites). This used to be the gold standard for array-based DNA methylation measurements for several years, though now the EPIC array is available with more methylation sites.
• Beta value and M value are two terms used to describe the methylation measurement.
• Beta values range from 0 to 1 and describe the proportion of methylation for a specific site in the sample, from completely unmethylated to completely methylated. 
• M values can be negative or positive and are a result of data normalization. M values are more useful than beta values as input for statistical models. If you want to describe differentially methylated probes (DMPs) then you need M values for your analysis.
• Beta values or sometimes "delta beta", annoyingly, can also refer to the model coefficient result from a generalized linear model. In that case, they would describe if a site is more or less methylated in group1 versus group2. Positive and negative values refer to the direction of higher methylation. Consult the manuscript to figure out the direction.
• Bisulfite conversion - this is a chemical change done by the lab to make it easier to distinguish unmethylated versus methylated sites on the DNA. It is done before the DNA methylation measurement step, whether that be bisulfite sequencing or bisulfite PCR or a methylation array. Bisulfite conversion does NOT mean bisulfite sequencing.
• Bisulfite sequencing - this is a method of measuring DNA methylation in a sample, using bisulfite conversion followed by DNA sequencing. 
• Bonferroni correction - this is a way of correcting p-values when you have a lot of measurements and thus a higher risk of false positives. The Bonferroni method is much stricter than the FDR method.
• P<0.05 is "nominally" significant, meaning it seems significant for a single site or gene but if you're looking at many sites or genes then you have a higher chance that your data has some false positives, and therefore P<0.05 isn't good enough for big data. You need to adjust for multiple comparisons to cut down the risk of false positives.
• FDR<0.05 is actually significant for big data.
• Bonferroni<0.05 is actually significant for big data and much stricter than FDR<0.05
• Bonferroni<0.05 data points are also going to meet criteria for FDR<0.05 and P<0.05, since Bonferroni<0.05 is stricter than both
• CpG - dinucleotide grouping of a cytosine base followed by a guanine base on the same DNA strand. It is not the same as C and G binding across DNA strands. CpG methylation (with a methyl group on the C base) is what people usually mean when they say DNA methylation, since it's the most common form of DNA methylation.

Why do we adjust p-values for multiple comparisons in big data analysis?

For RNA-sequencing, DNA methylation arrays, GWAS, and other big data projects, we run a statistical test to compare two groups for each row in the data matrix (each gene, each methylation site, each SNP, etc). That gives us a p-value for each reading. 

However, with big data, p<0.05 is only suggestive, not sufficiently strict. When working with thousands (miRNA-seq), tens of thousands (mRNA-seq or total RNA-seq), or >450,000 sites (DNA methylation), the occurence of false positives is much greater. This webcomic shows why we don't rely only on simple p-values to make conclusions: "Significant" (xkcd #882)

To reduce risk of false positives and narrow down to significant values, we adjust the original p-values for multiple comparisons. These two methods are popular:

• Benjamini-Hochberg's False Discovery Rate (FDR)
• Most common p-value adjustment for RNA-seq data.
• Sometimes used for DNA methylation data.
• Citation: Benjamini, Y., & Hochberg, Y. (1995). "Controlling the false discovery rate: a practical and powerful approach to multiple testing." Journal of the Royal Statistical Society. Series B (Methodological), 57(1):289-300.
• In R, make a new "FDR" column in your complete dataset using the full list of p-values, before adding any filters for expression or fold-change or otherwise subsetting the data:
  data$FDR = p.adjust(data$pvalues, method="BH")
  or, synonymously,
  data$FDR = p.adjust(data$pvalues, method="fdr")
• The R package DESeq2 creates a "padj" column using the p.adjust() R function as above. 
  
  
• Bonferroni correction
• Much stricter (more conservative) than the Benjamini-Hochberg FDR method.
• Leaves less significant values overall. Drastically reduces false positives (type I error), but increases false negatives (type II error).
• Sometimes used for DNA methylation data if other quality control steps suggest the data has high inflation (higher risk of false positives).
• Not typically used for RNA-seq.
• In R, make a new "bonferroni" column in your complete dataset using the full list of p-values, before adding any filters for expression or fold-change or otherwise subsetting the data:
  data$bonferroni = p.adjust(data$pvalues, method="bonferroni")
  
  

Additional reading:

• Mohieddin Jafari & Naser Ansari-Pour. (2019) "Why, When and How to Adjust Your P Values?" Cell J. 20(4): 604–607. doi: 10.22074/cellj.2019.5992
• R's p.adjust() function documentation

How to look up gene synonyms quickly with db2db (and fix the calendar gene issue)

If you have a list of gene symbols or Ensembl IDs and need to identify gene synonyms, try the Database to Database Conversions tool from bioDBnet: https://biodbnet-abcc.ncifcrf.gov/db/db2db.php

This is also a great way to correct calendar gene name conversion errors when opening a csv file in Excel. 

Quick tutorial for db2db

DNA methylation: Illumina methylation array annotation sources

Here I list of available resources and annotation files for Illumina methylation arrays to measure cytosine nucleotide DNA methylation.

There are three different probeType column values: "cg", "ch", "rs"

• "cg" - CpG sites, the most common DNA methylation sites
• "ch" - CpH sites (where the C base is followed by something other than a G base, H=A,T,C). Rare methylation sites. They are slightly more common in stem cells, but still not as common as CpG sites.
• "rs" - methylation sites that overlap single nucleotide polymorphisms (SNPs). These should be excluded from most analysis.
• If you don't have the probeType column in your data, you can identify these sites by searching the probeID prefixes, for example: grepl("^cg", df$probeID) is a regular expression for R programming that will return a column of TRUE/FALSE depending on the text values in column df$probeID. The caret symbol (^) in "^cg" means that the text match needs to be at the beginning to be TRUE. An example probeID that would return TRUE is "cg07881041". 
  
  

Official Illumina manifest files for methylation arrays

A "manifest" is a list of all items (sometimes samples, in this case methylation probes) and information about them. Illumina provides manifests for each of their methylation arrays with information such as probeID, probe type, probe chemistry details, chromosome position, and associated genes.  

Infinium MethylationEPIC BeadChip (>860k sites)

• Annotation files (Infinium MethylationEPIC Product Files, get the "manifest" files)
• Infinium MethylationEPIC v1.0 B5 Manifest File (CSV Format), Mar 13, 2020
• GRCh36, GRCh37, and GRCh38 human reference genomes
• See annotation website for pdf explaining each annotation column

Illumina HumanMethylation450K BeadChip (450k sites)

• Annotation files (Infinium HumanMethylation450K v1.2 Product Files, get the "manifest" files)
• HumanMethylation450K v1.2 Manifest File (CSV Format), May 23, 2013
• GRCh36 and GRCh37 human reference genomes

Infinium HumanMethylation27 BeadChip (27k sites)

• Annotation spreadsheet (Excel spreadsheet), Aug 27, 2013
• Uses GRCh36 human genome reference