RNAseq: Fisher-exact test with Benjamini-Hochberg multiple-testing correction using R

Question

[Ankit]

Hi everyone,

I would like to use R for performing statistical analysis on RNA-seq data. I have Deseq2 normalized count matrix. I want to perform Fisher-exact test with Benjamini-Hochberg multiple-testing correction to calculate p-values for each genes with respect to different conditions. In that way I want to to check if each genes have significant expression across two different condition WT vs KO.

For examples this table:

Gene    WT_allele1      WT_allele2      KO_allele1  KO_allele2  **Fisher_Exact_Test(two-sided)** 

Xkr4    1      0    0   0   **1**

Gm1992  0   1   0   0   **1**

Rp1 0   0   1   0   **1**

Sox17   1   1   0   2   **1**

Mrpl15  93  94  127 147 **0.506745441**

Lypla1  92  48  87  59  **0.328408209**

Tcea1   290 158 225 149 **0.192618917**

Gm6104  9   21  15  36  **1**

Rgs20   2   1   3   0   **1**

.   .   .   .   .   .

.   .   .   .   .   .

.   .   .   .   .   .

so on

I know tools like edgeR, DEseq2, NOIseq and others perform differential expression and provide statistical values. But is it possible if I use fisher exact test and Benjamini–Hochberg multiple correction to get statistical values for them using R manually? Am I thinking correctly or not? I appreciate any suggestions.

Thanks

[fr]

Ankit, may I ask why do you want to use Fisher's instead of those tools?

[Ankit]

It is mainly because previously we have performed Fisher for our old RNAseq data (not using R). I just wanted to understand that if I can derive same values by using Rscripts/functions.

I would also like to implement the same for my new data. It is mostly for the comparison purposes.

It will be very helpful if I will get the bit of guidance on this.

[fr]

I understand. However, DESeq, edgeR, or limma should give you more robust results than Fisher's exact test from my understanding. My suggestion is: Perform the Fisher's exact test > correct its p values as @Jean-Karim said below; Compare your significant (FDR < 0.05?) genes with those that com up from DESeq

[Ankit]

Thanks. Fisher.test gives only one p-value for contingency table. How will I get gene-wise p-values (as show in the table above or as obtained from DESeq2 normalisation)?

Answer 1

[Jean-Karim Heriche]

I think the function you're looking for is p.adjust().

[Ankit]

Thanks for the response.

This function will work on p-values as input.

As in manual:

Given a set of p-values, returns p-values adjusted using one of several methods.

I have DEseq2 normalized count matrix. How should I proceed in that case?

[Jean-Karim Heriche]

It seems you can start the DESeq2 analysis pipeline from a count matrix.

EDIT: You may want to check if this works with normalized counts.

[Ankit]

Hi,

Thank you for the suggestion.

Using p.value provided by Deseq2, I recalculated padj (using p.adjust(p.val, method ="BH")).It matched with p.adj value provided by Deseq2.

But

Can I apply fisher.test function on raw/normalized count matrix to get p.value gene-wise?

Later I can apply the function p.adjust() with BH method to get my adjusted p.value.

Please suggest