edgeR & DESeq2

Intro

  `edgeR`과 `DESeq2`는 RNAseq downstream analysis를 하다 보면 한번쯤은 보게되는 논문입니다. 두 가지 모두 `biomarker selection`에 대한 방법론과 해석을 제시합니다. 먼저 `edgeR`(2010)이 나오고 `DESeq2`(2014)가 나왔으며 두 논문 citation수가 (2022/10/12기준) 각각 27575회, 42124회 입니다. 엄청난 파워가 있는 논문들로 아직까지 커뮤니티에 질문과 답변이 활발하게 올라오고 있습니다. 

edgeR

DESeq2

  이 논문들은 `biomarker selection`을 위한 `Differentially Expressed Gene(DEG)` 분석 논문입니다. DEG: control(normal healthy group)과 case(patient group)간의 차이나는 feature(transcriptome)를 분석하는 방법론, t-test를 생각하면 쉽다. 해당 포스트에서는 `edgeR`과 `DESeq2`가 무엇인지, 시용된 모형 부분만 간단하게 살펴보고 정리할 예정입니다. 두가지 방법 모두 Bayesian 부분이 중요하고 추정하는데 사용되지만... motivation이나 수식들을 drive 할 능력은 없으므로 깊게 다루지 않겠습니다. 그 외의 수학적인 사실들은 본문의 수식을 보면 최대한 이해 가능하도록 글로 설명 하겠습니다. 본 논문의 수식과 같이 보세요!

edgeR & DESeq2

  `DEG`, 조금 더 자세히는, 비교군과 대조군 같은 condition에 따라 transcript, exom의 수가 통계적으로 유의미하게 다른지 확인하기 위한 workflow를 제시합니다. 여기서 사용되는 transcript, exom 자료는 `RNAseq` 기술로 count수로 정량화된 tabular자료로 생각하면 됩니다. Count수가 많을 수록 해당 gene의 발현이 많이 됬다는 의미. 이 table의 row와 column은 sample과 feature이며 일반적으로 n<<p 양상을 띄고 있어, 단순한 biomarker selection 뿐 아니라, condition의 유의미한 prediction을 위한 feature selection의 일환으로 이 방법들을 적용 할 수 있습니다.

 

  일반적으로 RNAseq 기술을 사용하여 정량화된 자료는 생물학적, 기술적인 노이즈가 끼어있습니다. 이로 인해 발생하는 RNAseq 자료의 고질적인 문제는 크게 `overdispersion`과 `sparseness` 두 가지가 있습니다.

Overdispersion 

  우선 RNAseq 자료는 count data이다 보니 Poisson 분포를 고려하는데, 여기서 분산이 평균보다 (큰 차이로) 커지는 문제가 바로 `overdispersion` 입니다. 참고로 Poisson분포는 이론적으로 분산과 평균이 같아야 합니다. 때문에 이 논문들은 이를 고려하여 modeling시  negative binomial(NB) 분포의 모수를 활용한 `overdispersion Posisson model`을 사용합니다. 구체적으로 이야기 하자면, 자료에 대하여 NB 분포의 모수를 구한 뒤, NB 분포와 Poisson 분포의 모수들간의 관계성을 이용하여 Poisson으로 모델링 한다는 의미입니다. 실제로 NB 분포는 count 자료에 적용 가능하며, 분산이 평균보다 커도 문제되지 않고, dispersion parameter가 0일 때 Poisson 분포와 같아집니다. 분포들의 특징을 잘 이해하여 적용된, 굉장히 스마트한 케이스로 생각됩니다.

Sparseness

  그리고 `sparseness`는 다른 말로 `excess zeros`라고도 하는데, 이는 기술적인 문제로 인해 자료에 0이 너무 많아지는 0 과잉 현상입니다. 때문에 이런 특성이 있는 자료를 모델링 할때에는 정말로 발현이 되지않은 0인지 (True 0) 아니면 기술적인 노이즈로 인한 가짜 0인지 고려해야 하기 때문에, 0 부분과 그 이외의 부분을 따로 modeling하는 `mixed model`을 많이 사용합니다. 이 논문들에서는 해당 부분에 대하여는 거의 언급을 하지 않으며, 나중에 다룰 다른 여러 방법론에서 어떠한 방식으로 이를 해소하는지 확인 할 예정입니다.

etc.

  그 외에 edgeR에서는 empirical Bayes method를 사용하여 feature 전반에 걸친 overdispersion 정도를 조절하여 추론의 신뢰성을 향상시킨다거나, replicate를 최소 수준으로 사용하여도 적용 가능하다는 실험적인 이야기, count data 형태를 띄는 다른 종류의 시퀀싱 자료에도 적용 가능하다는 이야기, 두 그룹 이상에서의 검정에 대한 약간의 언급과 R에서 사용 가능하다는 사실이 추가적으로 나옵니다.

 

  `DESeq2`는 `edgeR`과 거의 대부분 같은 방식이며 크게 `normalization`과 `dispersion estimation` 두 가지에서 차이를 보입니다. Normalization의 경우 두 방법 모두 `대부분의 gene은 그룹간의 차이가 없다` 라는 전제 가설로 부터 시작되며, normalization factor를 곱함으로 normalization을 합니다. 다만 두 방법이 취하는 방식이 `edgeR`은 `weighted mean log ratio (TMM)`를 사용하고, `DESeq2`는 `geometric concept (RLE)`를 사용함으로 달라진 것 뿐입니다.

 

  이러니 저러니 해도 결국 핵심은 아래와 같습니다.

 

"edgeR과 DESeq2는 RNAseq data의 DEG 분석을 위해 overdispersion Poisson (혹은 negative binomial)모델을 사용한다."

---

Volcano plot

  `DEG`분석에서 가장 중요한 output plot을 하나를 뽑으라면 `volcano plot`이라고 말하고 싶습니다. 위 두 방식을 포함한 모든 `DEG` 분석은 `p-value`와 `Fold Change(FC)`를 결과물로 냅니다. P-value는 두 그룹간의 발현량 차이가 없다는 귀무가설로 부터 나오는 값으로, 다중검정문제를 교정한 뒤의 p-value가 0.05 이하라면 두 그룹간 유의미한 차이가 있다는 결론을 내립니다. FC는 두 그룹간의 평균이 몇 배나 차이 나는지 직관적으로 보게하는 값입니다. 예를들어, 그룹 1에 대한 그룹 2의 FC값은 다음과 같이 계산됩니다. 그룹 1의 평균값이 50이고, 그룹 2의 평균값이 100인 경우 FC = 100/50 = 2 가 됩니다. 단순 FC는 계산이 용이하지만, 1을 기준으로 비대칭적인 FC값이 많이 등장하여(예시 FC의 역수를 생각해보라) log를 취하여 문제를 해소합니다. 이렇게 되면 비대칭 문제를 포함하여, log FC의 값이 양수면 "더 크구나", 음수면 "더 적구나" 하는 직관적인 해석이 가능해집니다. 다만 log를 취할 때 각 그룹 평균이 0인 경우가 종종 있어 pseudo count를 더하여 분석을 진행할 때도 많습니다. edgeR과 DESeq2 사용시 모든 포멧을 맞추었는데도 estimation이 안된다면 높은 확률로 pseudo count 문제다.

  

  `Volcano plot`은 log FC와 -log p-value를 각각 x, y축으로 사용한 scatter plot입니다. 통계적으로 얼마나 유의미한 결과가 나왔으며, 각 그룹간에 얼마나 큰 차이가 나는지 확인이 가능한 plot입니다. 

 

Volcano plot, from Maria Doyle (Galaxy project)

  이 plot의 경우 0.05 이상의 p-value는 유의미한 차이가 없으니 회색으로 표시했으며, log FC가 음수면 적다는 파란색으로, 양수면 크다는 빨간색으로 표시했습니다. 더 많이 차이 날수록 통계적으로 더 유의미하다는 상식에 맞게, 일반적으로 위와같은 모양의 plot이 나오게 됩니다. 이 모양이 화산이 분출하는 분화구같이 보입니다. 만약 이 plot이 이런 상식적인 모양에서 벗어났다면 무언가 잘못된 것은 없는지 코드를 전반 다시 확인해 볼 필요가 있습니다.

Outro 1

  사실 `DESeq2`와 `edgeR`은 regression에 몇 가지 Bayesian concept를 사용한 것으로 그렇게 대단한 수학이 사용 된 것은 아닙니다. 그런데 아직까지도 이런 엄청난 파워를 내며 높은 citation을 갖는 이유는 크게 두 가지로 생각됩니다. 첫째로 해석가능한 conventional statistical method를 사용했기 때문입니다. `Bioinformatics`에서 분석은 결과에 대한 `해석`이 매우 중요하기 때문에, prediction은 SOTA를 적용 하더라도 해석을 위해 conventional statistical method를 사용하기도 합니다. 둘째로 RNAseq 뿐 아니라 한참 인기 있던 `microbiome`과 computing power가 증가하며 가능해진 `single cell`분석에도 이 방법이 많이 사용되어 citation수가 계속 증가하는 것 같습니다. 

 

  오죽 이쪽 분야에서 입김이 강하면 아래와 같은 제목을 갖는 포스터가 2022년도에도 나왔습니다. 비록 태클 걸고싶은 부분이 많은 포스터지만, 충분히 생각해볼 만 하고 읽어볼만 한 내용입니다. 

https://towardsdatascience.com/deseq2-and-edger-should-no-longer-be-the-default-choice-for-large-sample-differential-gene-8fdf008deae9

DESeq2 and edgeR should no longer be the default choices for large-sample differential gene…

 

참고로 이 포스터에서는 Wilcoxon rank sum test에서 confounding factor를 보정 못한다고 하지만, 사실은 가능하다. 여러가지 방법이 있겠지만, 내가 가장 즐겨 사용하는 방식은 residual을 사용하는 것이다. 나중에 이에 관한 내용도 포스팅으로 찾아오겠다. 아래에 내용을 정리해두었다. 아래는 그에 관한 게시물의 첫번째 이다.

https://hello-world-jhyu95.tistory.com/entry/Permutation-of-Regressor-Residual-1https://jaehong-data.tistory.com/39

Permutation of Regressor Residual

 

Outro 2

  누군가는 그냥 더 좋은 방법을 사용하면 되는것이 아니냐는 질문을 합니다. 우선 어떠한 모형도 현실을 완벽히 설명 불가능합니다. 부분적으로 만족되는 현실의 현상들로부터 가설과 분포가정을 setting 할 수 있으며, 각각의 setting들은 만족되는 부분들을 설명합니다. `edgeR`과 `DESeq2`역시 마찬가지 입니다. 두 가지 모두 같은 현실에 대해 설명하지만, 각기 다른 부분을 (때로는 비슷한 부분을) 비출 수 밖에 없으므로 다른 결과가 나옵니다. 때문에 PCA의 PC1 마냥 가장 많은 정보를 보여주는 방법을 상황에 맞게 사용하던가, comparative하게 여러가지 방법을 사용 한 뒤 보수적으로 공통된 marker를 선택하던가 실험의 목적에 맞게 설계하는것이 좋아보입니다. 다만, 여러 방법 중 실험에 유리하도록 부분적으로 사용하는 manipulate은 꼭 피하도록 합니다. 개인적으로는 cheating이라 생각된다...!!

 

  결론적으로는 이런 방법들을 포함하여 여러가지 방법들을 이해하며, 현실을 어떠한 방향에서 비추고 있는지, 무슨 말을 하고싶은지에 대한 이해가 매우 중요해 진것같습니다.

 

+ Tip

  만약 `DESeq2`를 처음 사용한다면 `bioconductor`에 올라온 tutorial을 꼭 한번쯤 진행하고 넘어가길 바랍니다. (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html) 정말 친절하게 잘 설명합니다.

 

  Comparative research를 위해 이 두 가지 `R library`를 같이 사용 할 때가 있습니다. 이때 `calcNormFactors()` function을 조심합시다. 두 라이브러리 모두 같은 이름의 함수를 가지고 있어 충돌이 일어나기 때문이죠. 혹시라도 자동화 할 때 저 함수에서 같은 에러가 발생하면 각 방법 코드 앞에 라이브러리 이름과 ::을 붙여주면 됩니다. (예시) `edgeR::calcNormFactors()`

 

  `DESeq2`는 커뮤니티에서 활발하게 질문을 주고 받습니다. 이중 종종 답이 이상하거나 없는 경우가 있어 시간 낭비할 때가 많은데, `Michael I. Love`(저자)의 답변은 정답이고 많이도 달아두어서, 그 사람의 답만 보아도 대부분의 경우 커버가 가능합니다.

 

  개인적으로는 robust 결과를 위한 permutation 방식을 사용하는 버젼도 있으면 좋겠고, 그에 따라 Wald test보다 Score test를 사용 하여 코드가 도는 시간을 줄이면 하는 아쉬움이 있습니다. 하긴 nonparametric 하게 하려면 그냥 다른 방법을 쓸듯...  Permutation test과 test statistic 각각의 걸리는 시간에 대하여도 나중에 수식과 함께 포스터로 찾아오겠다.

 

출처

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 - Genome Biology

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2796818/

edgeR: a Bioconductor package for differential expression analysis of digital gene expression data - PMC

https://towardsdatascience.com/deseq2-and-edger-should-no-longer-be-the-default-choice-for-large-sample-differential-gene-8fdf008deae9

DESeq2 and edgeR should no longer be the default choices for large-sample differential gene…

https://www.biostars.org/p/284775/

How do I explain the difference between edgeR, LIMMA, DESeq etc. to experimental Biologist/non-bioinformatician

http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

Analyzing RNA-seq data with DESeq2

https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/rna-seq-viz-with-volcanoplot/tutorial.html

Transcriptomics / Visualization of RNA-Seq results with Volcano Plot / Hands-on: Visualization of RNA-Seq results with Volcano P