[RNA-Seq 분석 프로젝트 로그 #4] RNA-Seq 이상치 샘플, 포함할까 뺄까 : 판단 기준과 실제 사례

서론

Pseudomonas aeruginosa (PAO1) 처리에 따른 애기장대 RNA-Seq 분석 프로젝트에서, TopHat 실행 도중 R2 fastq 경로 문제로 재실행했던 4개 샘플(Root-PA01-01~04) 이야기를 지난 글에서 다뤘습니다.

재실행 자체는 성공했고 mapping rate도 91~98%로 합리적인 범위였지만, 그중 Root-PA01-01의 값이 유독 범위의 가장 아래쪽(91.4%)에 있었습니다.

 

이 정도 차이를 이상치로 볼 것인가, 그냥 넘어갈 것인가는 실무에서 항상 애매한 판단이기도 하며, 이번 분석에서도 크게 고민했던 부분입니다.

 

이번 글에서는 이 판단을 어떤 기준으로, 어떤 순서로 내렸는지 정리합니다.

결론부터 말하면 하나의 지표만으로는 판단하지 않았고, 네 가지 서로 다른 근거가 같은 방향을 가리키는지 확인한 뒤에 결정했습니다.

---

기준 1: 정렬 품질 — mapping rate와 CDS exon 비율

먼저 QC 전반을 확인했습니다.

16개 샘플 전체적으로 fastp 필터링 전후 Q30 rate는 약 95% 이상, Good Read Rate는 93~98%로 시퀀싱 품질 자체에는 문제가 없었습니다.

 

Duplication rate가 처리군(Treatment)이 아니라 조직(Tissue)에 따라 차이가 났던 것도(Root 평균 5.5% vs Shoot 평균 8.3%) PAO1 처리 효과가 아니라 조직별 RNA 구성 차이로 판단해 넘어갔습니다.

 

문제는 Mapping rate였습니다.

대부분의 샘플이 97~98%대였는데, Root-PA01-01만 91.4%로 낮게 나타났습니다.

이 원인을 확인하기 위해 RSeQC로 정렬된 리드가 유전자 영역(exon)에 얼마나 정확히 위치했는지 봤더니, 이 샘플의 CDS exon 비율이 64.1%로 나머지 샘플(83~90%)보다 뚜렷하게 낮았습니다.

mapping rate 저하와 일치하는 결과였습니다.

 

여기서 헷갈리지 않아야 할 부분이 있습니다.

같은 샘플의 fastp 결과(Q30 96.2%, GC 44.3%)는 다른 샘플과 유사한 정상 범위였습니다.

즉 전처리 단계는 문제가 없었고, 정렬(alignment) 단계에서 문제가 생겼다는 뜻이었습니다.

다만 이 수치들이 재정렬이 필요한 수준(70% 이하)에는 못 미쳐서, 이 시점에는 샘플을 제외하지 않고 다음 지표를 더 보기로 했습니다.

Mapping Rate by Sample 막대그래프

Group별 Boxplot

CDS Exon Ratio 그래프

---

기준 2: 발현량은 오히려 정상이었다 (헷갈리는 신호)

다음으로 확인한 건 Cuffnorm으로 정량한 발현량이었습니다.

그런데 여기서 예상과 다른 결과가 나왔습니다.

Root-PA01-01의 발현 유전자 수(21,893개)와 평균 발현량이 16개 샘플 중 가장 높게 나타난 것입니다.

낮은 발현량(FPKM 0~1) 유전자 수도 2,905개로 Root 8개 샘플 중 적은 편이라, 정렬 실패로 인한 발현량 왜곡이라고 보기도 어려웠습니다.

 

이 지점이 판단을 어렵게 만든 부분이었습니다.

정렬 지표는 이상 신호를 보이는데, 발현량 계산 결과는 정상이었기 때문입니다.

그래서 이 결과만으로 샘플을 정상이라고 단정하지 않고, 샘플 간 유사도를 보는 다음 단계로 넘어갔습니다.

---

기준 3: 샘플 간 유사도 — PCA와 Correlation

PCA 결과 PC1(49.8%)은 Root와 Shoot 조직을 뚜렷하게 구분했고, PC2(10.4%)는 상대적으로 설명력은 낮지만 Root 그룹 내 replicate들의 위치를 보는 데 활용했습니다.

Root-Control 그룹의 4개 replicate는 PC2 상에서 서로 가깝게 뭉쳐 재현성이 양호했던 반면, Root-PA01-01만 나머지 PA01 replicate(02, 03, 04)와 달리 PC2 상에서 크게 떨어진 위치에 나타났습니다.

PCA scatter plot

Sample Correlation heatmap에서도 같은 경향이 확인됐습니다.

Root-Control 그룹 내 상관계수는 0.98 이상으로 균일했고 Root-PA01 그룹의 02·03·04 사이도 0.97~0.99로 양호했지만, Root-PA01-01은 나머지 세 replicate와의 상관계수가 0.947~0.973으로 낮게 나타났습니다.

Sample Correlation heatmap

두 가지 독립적인 방법(PCA, Correlation)에서 같은 샘플이 같은 방향으로 이탈한다는 건, 이 이탈이 우연이 아니라는 뜻이었습니다.

다만 이 차이가 정렬 과정의 기술적 편향 때문인지, PAO1 처리에 대한 개체 반응 차이 때문인지는 이 시점에는 단정하지 않았습니다.

---

기준 4: 실제 결과에 미치는 영향을 확인한다

앞선 세 가지 근거(정렬 품질 이상, 발현량은 정상, 유사도 이탈)만으로는 최종 결정을 내리기 애매했습니다.

그래서 가장 직접적인 방법을 택했습니다.

이 샘플을 포함했을 때와 뺐을 때 DEG 분석 결과가 실제로 얼마나 달라지는지 직접 재실행해서 비교한 것입니다.

 

두 가지 기준으로 비교했습니다.

• q-value 기준만 적용했을 때: 8개 샘플(포함) 3,473개 → 7개 샘플(제외) 5,055개. 공통 DEG는 3,360개(96.7%)였습니다.
• 최종 DEG 선정 기준(q-value < 0.05 & |log2FC| ≥ 1)을 적용했을 때: 8개 샘플(포함) 2,175개 → 7개 샘플(제외) 2,840개로, 약 30.6% 차이가 났습니다.

어느 기준으로 봐도 결론은 같았습니다.

Root-PA01-01이 그룹 내 분산을 키워서 통계적 검정력을 낮추고 있었다는 것입니다.

참고로 fastp와 FastQC로 이 샘플을 교차 검증했을 때도 GC 비율이나 Per base sequence quality는 정상 범위로 일치했는데, 이는 문제가 전처리가 아니라 정렬 단계에 있었다는 앞선 판단과 다시 들어맞는 결과였습니다.

DEG Overlap 벤다이어그램(Include vs Exclude Root-PA01-01)

---

최종 판단: 네 가지 근거가 같은 방향을 가리켰다

근거	Root-PA01-01	나머지 7개
Mapping rate	91.4%	97.5~98%(PA01-04만 94.6%)
CDS exon 비율	64.1%	83~90%
발현 유전자 수·평균 발현량	오히려 가장 높음	—
PCA(PC2) / Correlation	이탈 / r=0.947~0.973	정상 군집 / r=0.97~0.99
DEG 포함 시 vs 제외 시	검정력 저하 확인	—

 

정렬 품질(기준 1), 샘플 간 유사도(기준 3), DEG 검정력(기준 4)까지 세 가지 독립적인 근거가 모두 같은 결론을 가리켰기 때문에, Root-PA01-01을 downstream DEG 분석에서 제외하기로 결정했습니다.

발현량 지표(기준 2)만 정상으로 나온 것은 결과를 뒤집을 근거가 아니라, "하나의 지표만 보고 판단하면 안 된다"는 걸 보여준 사례로 남았습니다.

 

최종적으로 DEG 분석에는 Root-Control 4개 + Root-PA01 3개(02, 03, 04), 총 7개 샘플을 사용했습니다.

---

다음 포스팅 예고

이렇게 확정된 7개 샘플로 DEG를 뽑기 시작하면, 다음 질문은 "어떤 기준으로 유의한 유전자를 골라낼 것인가"입니다.

다음 편에서는 q-value와 log2FC 컷오프를 어떻게 정했는지, 그리고 그 결과가 무엇을 보여주는지 다룹니다.

 

🔗 전체 분석 코드와 실행 로그는 GitHub 저장소에서 확인할 수 있습니다. github.com/ourkofe/rna_analysis_report