Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data

2019, BMC Bioinformatics, 247citation (3.327 impact factor)

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-019-2599-6

Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data - BMC Bioinformatics

Abstract
(background)
scRNAseq 데이터의 분석은 내재적 및 외재적인 세포 과정을 이해하는 데 중요한 역할을 한다.
scRNAseq 데이터는 매우 다양하며 많은 수의 제로 카운트(zero counts)를 가지고 있어 DE 유전자의 탐지에 어려움이 있다.
(zero counts 의 경우 발현량을 의미하는 데, bulk RNA-seq과 다르게 scRNA-seq은  single cell level로 발현량을 봄으로써, 보다 세부적이며, 데이터가 sparse 되는 경향이 있다. 예를 들어 gene 1 에 대한 cell 별 발현량은 50 0 0 0 0 이었다면 bulk RNA-seq에서는 50/5 10 (평균 발현)으로 표현된다. 즉 보다 많은 0을 확인하게 된다. 따라서 새로운 계산법이 필요하다)
이를 고려한 11가지 scRNAseq 데이터 분석 방법의 성능을 평가하고 비교한다.
(Result)
시뮬레이션 및 실제 데이터를 사용하여 accuracy & precision 성능 평가
- 시뮬레이션 데이터 : 샘플 크기가 도구의 탐지 정확도에 미치는 영향을 조사
- 실제 데이터를 사용하여, DE 유전자를 식별하는 도구 간의 일치성, 도구의 실행 시간, 검출된 DE 유전자의 생물학적 관련성을 검토
(Conclusion)
- 도구간 DE 유전자 호출의 일치 정도는 높지 않다.
- 정밀도가 높으면  적은 수의 DE 유전자를 식별함으로써 낮은 참 양성률을 보인다.(precision과 TP는 trde-off 관계이다.)
- scRNAseq 데이터용으로 설계된 방법들은 일반적으로 bulk RNAseq 데이터를 위해 설계된 방법들에 비해 더 나은 성능을 보이지 않는 것으로 관찰

 

Introduction

background

일반적으로 RNAseq 데이터 분석에서 하나의 샘플의 유전자의 발현값은 모든 세포 집단의 발현값의 평균을 나타냅니다.

순환종양세포 및 줄기세포와 같은 다른 생물학적 연구에서는 single cell type 별 발현분석이 필요합니다.

 

기존의 RNA seq 데이터 DEG 분석용 도구인 DESeq 및 edgeR 으로도

평균 발현 차이의 정의를 이용하여 Single-cell differential expression (SCDE) 분석이 가능하다.

드롭아웃(drop-out)한 세포에서는 높은 발현을 보이는 유전자가 다른 세포에서는 누락될 수 있습니다.

 

the challenges in scRNAseq

scRNAseq 데이터에서 단일 세포의 RNA 분자 수가 적고 캡처 효율이 낮다.

드롭아웃(drop-out) : 한 세포에선 높은 발현을 보이는 유전자가 다른 세포에서 누락될 수 있다.

세포 주기나 세포 상태의 차이, 세포 간 상호작용, 외부 환경 요인 등으로 동일한 뇌조직에서 얻은 세포여도, 세포마다 큰 이질성을 보인다. 

multimodal expression values : 같은 조직내 세포간의 이질성때문에 한 유전자의 발현이 서로 다른 분포 형태를 갖게된다.

 

scRNA-seq 전용 분석 툴
Monocle2 : 단일 세포 실험의 변동성을 제거하고 카운트를 더 잘 정규화하기 위해 정규화된 전사체 카운트 대신 인구 조사 카운트를 사용
scDD : 유전자 발현 값 분포 내에서 생물학적 조건 내에서의 네 가지 다른 모드 상황을 고려
DEsingle : 제로 팽창 음이항 (ZINB) 회귀 모델을 사용하여 실제 및 드롭아웃 제로의 비율을 추정하고, 차이가 있는 유전자를 세 가지 범주로 분류
SigEMD [37], EMDomics [31], D3E : 비모수적인 방법, a distance metric between the distributions of genes in two conditions.
목적 : differential gene expression analysis of heterogeneous data(다양한 실험 조건의 데이터 등)
(Without modeling the distributions of gene expression values and estimating their parameters)

 

scRNAseq 데이터에 대한 기존의 벤치마킹 논문
1. Jaakkola et al. [40]은 scRNAseq 데이터를 위해 개발된 다섯 가지 통계 분석 방법을 비교
- SCDE, MAST / DESeq, Limma / 관련 기능이 없는 방법ROST

- Methods developed originally for bulk RNA-seq, DESeq and Limma, were not suitable for analyzing scRNA-seq data.

2. Miao et al. [41]

 - DEGSeq는 일반적으로 더 좋은 성능을 보이는 경향이 있습니다. 하지만 데이터의 특성에 따라 다른 방법들이 장점을 가질 수도 있습니다.
- SCDE, monocle, D3E
- BPSC, DESeq, edgeR,bySeq, NGPSeq, Cuffdiff, DEGseq, TSPM, limma, ballgown, SMAseq

Table 1. Information of gene differential expression analysis methods used.

3. 한 연구 [42]

- scRNAseq 데이터에 대해 개발된 네 가지 differential expression 분석 도구와 bulk RNAseq에 사용되는 두 가지 도구를 비교하였습니다. 
- MAST, SCDE, Monocle, D3E / DESeq, edgeR

- 두 연구에서 모두 Islam 데이터셋을 사용했지만, 데이터 처리 방식이 다른 경우에도 SCDE가 DESeq와 MAST보다 우수한 성능을 보였다는 점을 강조

- 데이터 시뮬레이션에서는 드롭아웃 구성 요소가 고려되지 않았기 때문에, 벌크 방법의 성능이 단일 셀 도구의 성능과 크게 다르지 않은 이유가 될 수 있음

 

Fig1

- ES는 Embryonic Stem cell (배아 줄기 세포)

- MEF는 Mouse Embryonic Fibroblast (마우스 배아 섬유아세포)

- DU: Dual mode (이중 모양)을 의미하며, 분포가 두 개의 주요 모드 또는 피크를 가지고 있는 상태를 나타냅니다.
- DP: Peak mode (피크 모양)을 의미하며, 분포가 하나의 주요 모드 또는 피크를 가지고 있는 상태를 나타냅니다.
- DM: Multi-mode (다중 모드)을 의미하며, 분포가 여러 개의 주요 모드 또는 피크를 가지고 있는 상태를 나타냅니다.
- DB: Broad mode (넓은 모양)을 의미하며, 분포가 넓게 퍼져있는 상태를 나타냅니다.

왜 scRNAseq 데이터는 매우 다양하며 많은 수의 제로 카운트(zero counts) 를 갖는가?

(자신의 의견을 작성합니다.) 

 

Methods

Differential gene expression analysis methods for scRNAseq data

 

Table 1 Software tools for identifying DE genes using scRNAseq data

method의 원리 및 수식을 이해하고 정리합니다.

1.SCDE

gene expression values에 대한 혼합 확률 모델(포아송 분포 + BM)

drop out 이벤트는 read count의 발현 값이 예상보다 작기 때문에, 드물게 발생하는 사건에 대한 모델링에 자주사용하는 포아송 분포를 사용

 

 

2. MAST

3. scDD

4. EMDomics

5. Monocle2

6. D3E

7. SINCERA

8. edgeR

9. DESeq2

10. DEsingle

11. SigEMD

 

(이해를 바탕으로 작성합니다.) 

Results

결과를 요약하고, 그 결과가 시사하는 바를 정리합니다.

(결과를 비판적으로 생각하며 작성합니다.) 

Accuracy of identification of DE genes

시뮬레이션 된 데이터 사용 (시뮬레이션된 데이터에 다중 형식과 0 카운트가 포함)

- TP를 진짜 DE 유전자로 정의
- FP를 유의미하지만 실제로 DE 유전자가 아닌 유전자
- TN은 실제로 DE되지 않은 유전자이면서 유의미하지 않게 호출된 유전자
- FN은 실제로 DE인 유전자이지만 유의미하게 호출되지 않은 유전자

 

시뮬레이션 데이터를 사용하여 DE 유전자를 탐지하는 도구들의 진짜 감도와 정확도를 계산
Number of deteccted DE genes : Monocle2>EMDomics >DESeq2 >D3E >DEsingle >SigEMD >edgeR >scDD >MAST >SINCERA >SCDE 


TPs 중에서 제대로 감지된 비율로, DE 유전자 중에서 실제로 감지된 DE 유전자의 비율을 나타냅니다. 
높은 감도는 DE 유전자를 더 잘 식별할 수 있다는 것을 의미합니다. 
Monocle2이 가장 높은 감도를 보이고, SCDE이 가장 낮은 감도를 보입니다.

False positive rate (거짓 양성율): SCDE < MAST <SINCERA  < MAST <DEsingle <edgeR<SigEMD < scDD <D3E <DESeq2 <EMDomics  <Monocle2 
DE가 아닌데 DE로 잘못 분류된 비율을 나타냅니다.

Precision (정밀도): SCDE > MAST > SINCERA > DEsingle > edgeR > SigEMD > scDD > > DESeq2 >D3E > EMDomics > Monocle2
TPs 중에서 실제로 DE인 비율로, DE로 식별된 유전자 중에서 실제로 DE인 비율을 나타냅니다.
높은 정밀도는 DE 유전자 식별의 정확성을 나타냅니다. 
SCDE 가장 높은 정밀도를 보이고, Monocle2 가 가장 낮은 정밀도를 보입니다.

Accuracy (정확도): DEsingle > SigEMD > edgeR > MAST < scDD > SCDE > DESeq2 > SINGCERA > D3E > EMDomics > Monocle2
전체 유전자 중에서 제대로 분류된 비율로, DE와 non-DE를 포함한 모든 유전자에 대한 분류의 정확성을 나타냅니다. 
높은 정확도는 전체적으로 잘 분류되었음을 나타냅니다. 
DEsingle 가장 높은 정확도를 보이고, Monocle2 가장 낮은 정확도를 보입니다.

F1 score: DEsingle > SigEMD > DESeq2> edgeR > scDD > MAST> D3E > EMDomics > SCDE > SINCERA > Monocle2
정밀도와 감도의 조화 평균으로, DE 유전자 식별의 균형을 나타내는 지표입니다.
 높은 F1 점수는 정밀도와 감도가 모두 높은 모델을 나타냅니다. 
DEsingle이 가장 높은 F1 점수를 보이고, Monocle2 가 가장 낮은 F1 점수를 보입니다.

 

ROC curves for the eleven differential gene expression analysis tools using simulated data

-  Monocle2는 가장 많은 수의 실제 DE 유전자를 식별했지만, 동시에 가장 많은 거짓 DE 유전자도 도입하여 식별 정확도가 낮게 나타났습니다(0.824). 
- 비모수적인 방법인 EMDomics와 D3E는 모수적 방법보다 더 많은 수의 실제 DE 유전자를 식별했습니다(각각 2465.8개와 1683.4개). 그러나 이들은 또한 많은 거짓 양성율을 도입하여 정확도가 낮아졌습니다(각각 0.910과 0.929). 
- 반면에 정밀도가 0.9보다 큰 값인 MAST, SCDE, edgeR 및 SINCERA와 같은 도구들은 적은 수의 거짓 양성율을 도입하지만, 작은 수의 TP를 식별했습니다. 
- 흥미롭게도, DESeq2와 edgeR과 같은 전통적인 bulk RNA-seq 데이터를 위한 도구들은 scRNA-seq 데이터를 위해 설계된 도구들과 비교했을 때 성능이 나쁘지 않았습니다. 
- DEsingle과 SigEMD는 정확도와 F1 점수 측면에서 가장 우수한 성능을 보여주었으며, 높은 TP를 식별하면서도 많은 수의 FP를 도입하지 않았습니다.

 

우리는 1000개의 검증된 유전자를 골드 표준 유전자 집합으로 사용했습니다. 
우리는 DE 유전자 중 도구에 의해 호출되고 1000개의 골드 표준 DE 유전자 중에 있는 유전자를 참으로 감지된 DE 유전자로 정의했습니다. 
FDR 또는 조정된 p-value 값이 0.05일 때, 각 도구의 검출된 DE 유전자 수와 1000개의 골드 표준 유전자 대비 실제로 감지된 DE 유전자 수(민감도)
Monocle2, EMDomics, SINCERA, D3E, DEsingle은 0.7 이상의 민감도를 가지고 첫 번째 수준에 위치합니다. 
edgeR, DESeq2, SigEMD는 0.4에서 0.7 사이의 민감도를 가지고 두 번째 수준에 위치합니다. 
SCDE, scDD, MAST는 0.4 미만의 민감도를 가지고 세 번째 수준에 위치합니다. 

그림 4에서 파란색 막대는 골드 표준 유전자와 방법에 의해 호출된 DE 유전자의 교집합(실제로 감지된 DE 유전자)을 나타내고, 노란색 막대는 골드 표준 유전자 중에 포함되지 않은 유의한 DE 유전자의 수를 나타냅니다.

그러나 더 높은 민감도를 보이는 방법들은 또한 7000개 이상의 유전자를 유의한 DE 유전자로 호출했습니다. 그림 4에서 파란색 막대는 골드 표준 유전자와 방법에 의해 호출된 DE 유전자의 교집합(실제로 감지된 DE 유전자)을 나타내고, 노란색 막대는 골드 표준 유전자 중에 포함되지 않은 유의한 DE 유전자의 수를 나타냅니다.

 

Agreement among the methods in identifying DE genes

우리는 도구들이 DE 유전자를 식별하는 데 얼마나 일치하는지 알아보기 위해 각 도구 쌍에서 공통적으로 식별된 DE 유전자의 수를 조사했습니다. 
먼저, 유전자를 조정된 p-value 또는 FDR로 정렬한 다음 상위 1000개의 DE 유전자를 선택했습니다. 
우리는 도구 쌍에 의해 공통적으로 식별된 DE 유전자의 수를 쌍별 일치로 정의했습니다. 도
구들 간의 공통적인 DE 유전자의 수는 시뮬레이션 데이터에서 770에서 1753개 사이이며 (추가 파일 1: 도표 S7), 실제 데이터에서는 142에서 856개 사이입니다 (도표 5). 
우리는 시뮬레이션 데이터와 실제 데이터 모두에서 도구들이 높은 쌍별 일치를 가지고 있지 않음을 관찰했습니다.

Effect of sample size

우리는 시뮬레이션 데이터를 사용하여 샘플 크기가 DE 유전자 감지에 미치는 영향을 조사했습니다. TPR, FPR, 정밀도 및 정확도 관점에서의 영향을 조사하였습니다. 정밀도는 TP/(TP+FP)로 정의되며, 정확도는 (TP+TN)/(TP+TN+FP+FN)로 정의되었습니다. 각 조건마다 10개, 30개, 50개, 75개, 100개, 200개, 300개 및 400개의 셀로 구성된 8가지 케이스를 생성하였습니다. 우리는 모든 도구에 대해 기본 FDR 또는 조정된 p-value (<0.05) 하에서 식별된 DE 유전자의 수와 탐지율(TPR)이 10에서 400까지의 샘플 크기 증가와 함께 증가하는 경향을 관찰했습니다 (도표 7).

Enrichment analysis of real data

생물학적 과정과 관련된 DE 유전자가 의미 있는지 확인하기 위해,

웹 기반 GSEA 도구의 Investigate Gene Sets 기능을 통해 유전자 집합 풍부도 분석을 수행했다.

비모수적인 방법(EMDomics 와 D3E)으로 식별된 상위 n=300개 DE 유전자들이 다른 방법들에 비해 더 많은 KEGG 경로에 풍부하게 분포되어 있는 것을 관찰할 수 있다.

 

우리는 각 도구에 의해 식별된 사우이 n=300개의 DE유전자를 대상으로 Gene Onthology Process 풍부도를 분석하기 위해 DAVID를 사용했습니다. 

FDR 임계값이 0.05인 경우 Gene Onthology 용어의 수를 표 5에서 보여줍니다.

EMDomics, D3E, Monocle, DESeq2 에서 식별된 상위 DE 유전자들은 다른 도구들보다 더 많은 KEGG 경로 및/또는 GO 용어에 풍부하게 분포되어 있는 것을 관찰 할 수 있습니다.

 

Conclusions

평가 기준

- DE 유전자 식별하는 정확도

- 도구들 간의 DE유전자 검출 일치도

- 걸리는 시간

- 실제 데이터를 통하  pathway & GO analysis 결과

Agreement among the methods

더 높은 민감도를 보이는 도구들은 보다 낮은 정밀도를 가지는 경향이 있습니다. 
scRNAseq를 위해 설계된 DEsingle과 SigEMD는 TPR과 정밀도 사이의 균형을 더 잘 유지하는 경향이 있습니다.

 

(scRNA 챌리지한 조건이 없을 때에는 모두 잘 동작함, 넘 당연한 말임;;)
모든 도구는 다중모드가 없거나 낮은 수준의 다중모드인 경우에는 모두 잘 작동합니다. 
또한, 희소성(0 값)이 적을 때 성능이 더 좋습니다. 

결론적으로 다중모드 수준이 낮을 경우에는 SCDE, MAST, 그리고 edgeR이 보다 높은 정밀도를 제공할 수 있습니다.

(다중모드 수준이 높은 데이터에서)

요약 : 결국 다중모드 수준이 높은 데이터에서는DESeq2, EMDomics, Monocle2, DEsingle, 그리고 SigEMD 이 잘 작동하고, 그러한 이유는 다중 모드를 잡아내는 통계모델을 사용하고 있기 때문이다.

 

DESeq2, EMDomics, Monocle2, DEsingle, 그리고 SigEMD와 같이 각 유전자의 행동을 고려하는 방법들이 더 높은 TPR을 보입니다. 
EMDomics와 SigEMD는 두 분포 사이의 거리를 비모수적인 방법으로 계산하여 다중모드를 포착할 수 있습니다. 
DEsingle은 발현값에서 실제 값과 드롭아웃 값의 비율을 추정하기 위해 zero inflated negative 모델을 사용하여 드롭아웃 이벤트를 잘 모델링합니다. 
Monocle2는 TPM 값과 같은 정규화된 read count 대신 각 유전자의 상대적인 전사체 수를 추정하기 위해 census 알고리즘을 사용합니다. 
DESeq2는 read count에 대한 음이항 모델을 적합시키기 위해 유전자별로 슈링크 추정을 사용하여 적산 효과 추정을 수행합니다. 

Sample size effect

모든 도구들은 각 조건별로 샘플이 많을수록 성능이 향상됩니다. 
샘플 크기가 10에서 75로 증가함에 따라 TPR이 크게 향상되지만, 샘플 크기가 100 이상이 되면 향상 속도가 둔해지며, 300 이상의 샘플 크기에서는 TPR과 FPR에 거의 변화가 없습니다. 
Monocle2, EMDomics, DESeq2, DEsingle 및 SigEMD는 샘플 크기를 증가시킴으로써 TPR을 거의 100%에 가깝게 달성할 수 있습니다. 

 

DEsingle은 많은 수의 zero counts 또는 적은 수의 샘플에 대해 잘 작동합니다. 
zero counts의 수가 적고 샘플 수가 많을 때는 해당 모델이 dropout 이벤트를 잘 포착하지 못할 수 있습니다.

Enrichment analysis

줄기세포 생물학과 관련된 GO(Gene Ontology) 용어로 풍부하게 어휘화되어 있음을 확인

 

scDD와 SCDE는 높은 유의성에서 줄기세포와 관련된 용어를 복구하지 못했습니다.

대신, 그들은 세포 기능과 관련된 용어를 얻었습니다.

이 결과는 다중모드를 고려하지 않는 모델 기반의 단일세포 DE 분석 방법은 데이터에서 생물학적으로 관련 있는 유전자 집합을 잘 추출하지 못하는 것으로 나타났습니다.

- scRNAseq 데이터에 특화된 도구들은 주로 zero counts나 다중모드를 처리하는 데 초점을 맞추고 있지만, 둘 다 처리하는 도구는 없습니다. 
- 일반적으로, 다중모드를 포착할 수 있는 방법(비모수적 방법)이 zero counts를 처리하는 모델 기반 방법보다 성능이 우수
- 그러나 drop-out 이벤트를 잘 모델링할 수 있는 모델 기반 방법은 true positive과 false positive 측면에서 더 좋은 성능을 발휘할 수 있습니다.

 scRNAseq 데이터에 특화된 방법들은 일반적으로 bulk RNAseq 데이터를 대상으로 설계된 방법들과 비교하여 크게 우수한 성능을 보이지 않았으며, 

이 도구들 간 DE gene의 식별에서의 일치 부족

실제 DE gene 및 생물학적으로 관련성 있는 유전자 세트를 감지하는 한계

 

scRNAseq 데이터의 DE 분석을 위해 더 정확한 방법의 개발이 필요함을 나타냅니다. 

 

다중모드, 이질성 및 희소성(zero counts가 많은)은 새로운 방법을 개발할 때 고려해야 할 scRNAseq 데이터의 주요 특성입니다.