학습 목표
-
카운트 데이터 변환 방법의 중요성 이해
-
PCA(주성분 분석) 이해 -
PCA와 계층적 군집화를 사용하여 샘플 품질을 평가하는 방법 이해 -
품질 관리
DESeq2 워크플로우의 다음 단계는 QC이며, 이는 샘플과 유전자 수준에서 카운트 데이터에 대한 QC 검사를 수행하여 샘플 또는 반복이 양호하게 보이도록 보장하는 것을 포함합니다.
- 샘플 QC
RNA-seq 분석에서 유용한 초기 단계는 일반적으로 샘플 간의 전반적인 유사성을 평가하는 것입니다:
- 어떤 샘플이 서로 유사하고 어떤 샘플이 다른가요?
- 이것이 실험 설계의 예상과 일치합니까?
- 데이터 세트의 주요 변동 원인은 무엇입니까?
샘플의 유사성을 탐색하기 위해 우리는 주성분 분석(PCA)과 계층적 군집화 방법을 사용하여 샘플 수준의 QC를 수행합니다. 이러한 방법 또는 도구는 반복들이 서로 유사한 정도(군집)를 확인하고 실험 조건이 데이터 변동의 주요 원인인지 확인할 수 있게 해줍니다. 샘플 수준의 QC는 또한 잠재적인 이상값을 보이는 어떤 샘플을 식별하는 데 도움이 될 수 있습니다; 우리는 DE 분석 전에 제거해야 할지 확인하기 위해 잠재적인 이상값을 추가로 탐색할 수 있습니다.
이러한 비지도 군집화 방법은 log2 변환된 정규화 카운트로 실행됩니다. log2 변환은 시각화 거리를 개선합니다. 우리는 일반적인 log2 변환 대신 저 카운트 유전자로 인한 편향을 피하기 위해 정규화 로그 변환(rlog)을 사용할 것입니다;
- 왜 데이터 변환이 필요한가요?
다차원 데이터 탐색적 분석에 사용되는 많은 통계 방법, 특히 군집화 및 정렬 방법(예: 주성분 분석 등)은 최소한 근사적으로 동분산 데이터에 가장 적합합니다; 이는 관측량의 분산(여기서는 유전자 발현 값)이 평균에 의존하지 않음을 의미합니다. 그러나 RNA-seq 데이터에서 분산은 평균이 증가함에 따라 증가합니다. 예를 들어, 정규화된 읽기 카운트 행렬에 직접 PCA를 실행하면 결과는 일반적으로 가장 큰 절대적 차이를 보이는 소수의 고발현 유전자에만 의존합니다. 이를 피하는 간단하고 자주 사용되는 전략은 정규화 카운트 값에 작은 의사 카운트를 더한 로그를 취하는 것입니다; 그러나 이제 저 카운트를 가진 유전자들은 작은 카운트 값 내재의 강한 포아송 노이즈로 인해 샘플 간에 가장 강한 상대적 차이를 보이기 때문에 결과를 지배하는 경향이 있습니다.
rlog사용의 장점:
따라서 DESeq2는 정규화 로그 변환, 또는 간단히 rlog를 제공합니다. 고 카운트 유전자의 경우, rlog 변환은 일반적인 log2 변환과 크게 다르지 않습니다. 그러나 저 카운트 유전자의 경우, 이 값들은 모든 샘플에서 유전자의 평균으로 축소됩니다. 이렇게 하는 이유는 rlog 변환된 데이터가 근사적으로 동분산이 되도록 하기 위함입니다.
DESeq2는
rlog함수가 오랜 시간이 걸릴 수 있기 때문에 대규모 데이터 세트(100개 샘플)에는 카운트 변환을 위해vst(분안 안정 변환)를 사용할 것을 권장합니다.vst()함수는 유사한 상황에서 더 빠릅니다.
- PCA
주성분 분석(PCA)은 변동을 강조하고 데이터 세트의 차원을 축소하는 기술입니다. 이는 품질 관리와 Bulk RNA-seq 및 단일 세포 RNA-seq 데이터 분석에 매우 중요한 기술입니다.
3.1. PCA 플롯
본질적으로 두 샘플의 유전자 발현 수준이 유사하면, 이러한 유전자가 주어진 PC(주성분)가 나타내는 변동에 상당한 기여를 하므로, 이 샘플들은 해당 PC를 나타내는 축에 가까이 그려집니다. 따라서 우리는 생물학적 반복들이 유사한 점수를 가질 것으로 기대합니다(우리의 기대는 동일한 유전자들이 변화하고 있다는 것입니다) 그리고 함께 군집을 형성합니다. 이를 이해하기 몇 가지 예제 PCA 플롯을 시각화하는 것이 가장 쉽습니다.
아래에 예제 데이터 세트와 관련된 PCA 플롯이 있습니다. 실험의 메타데이터는 다음과 같이 표시됩니다. 관심의 주요 조건은 처리입니다.
PC1과 PC2에서 시각화할 때 처리로 인해 샘플이 분리되지 않는 것을 보았으므로, 데이터에 존재하는 다른 변동 원인을 탐색하기로 결정했습니다. 우리가 메타데이터 테이블에 이미 모든 가능한 알려진 변동 원인을 포함했고 이러한 요소들을 사용하여 PCA 플롯에 색을 칠할 수 있다고 가정합니다.
cage 요소부터 시작하지만 cage 요소는 PC1 또는 PC2의 변동을 설명하지 못하는 것 같습니다.
그런 다음 sex 요소로 색을 칠하는데, 이는 PC2에서 샘플을 분리하는 것 같습니다. 이는 우리가 모델에서 sex로 인한 변동을 설명하고 회귀시키기 위해 다운스트림에서 사용할 수 있으므로 주의해야 할 좋은 정보입니다.
다음으로 strain 요소를 탐색하고 이것이 PC1의 변동을 설명할 수 있음을 발견합니다.
PC1과 PC2의 변동 원인을 식별할 수 있어서 기쁩니다. 이를 모델에 고려함으로써 우리는 처리로 인해 차이 발현하는 더 많은 유전자를 감지할 수 있을 것입니다.
우리는 두 샘플이 올바른 strain으로 군집화되지 않은 것을 보고 있습니다. 이는 샘플 교환이 있을 수 있음을 나타내며, 이 샘플들이 실제로 표시된 strain인지 확인하기 위해 조사해야 합니다. 교환이 있다는 것을 발견하면 메타데이터에서 샘플을 교환할 수 있습니다. 그러나 올바르게 표시되었다고 생각하거나 확실하지 않다면 데이터 세트에서 샘플을 제거할 수 있습니다.
우리는 여전히 처리가 strain과 sex 후 변동의 주요 원인인지 발견하지 못했습니다. 따라서 우리는 PC3와 PC4를 탐색하여 처리가 이 두 PC 중 어느 하나가 나타내는 변동을 주도하는지 확인합니다.
우리는 처리로 인해 샘플이 PC3에서 분리되는 것을 발견했고, 우리가 관심 있는 조건인 처리가 PC3에서 분리되고 PC1과 PC2의 변동을 회귀시킬 수 있기 때문에 우리의 DE 분석에 대해 낙관적입니다.
처음 몇 개의 주성분이 변동의 얼마를 설명하는지에 따라 더 많은 것을 탐색하고 싶을 수 있습니다(즉, 더 많은 성분을 고려하고 쌍 조합으로 플롯). 샘플이 실험 변수에 의해 명확하게 분리되지 않더라도 여전히 DE 분석에서 생물학적으로 관련된 결과를 얻을 수 있습니다. 효과 크기가 매우 작을 것으로 예상한다면 신호가 무관한 변동 원원에 의해 묻힐 수 있습니다. 이러한 원인을 식별할 수 있는 경우 모델에서 이를 고려하는 것이 중요합니다. 왜냐하면 이는 DE 유전자를 감지하는 도구에 더 많은 기능을 제공하기 때문입니다.
- 계층적 군집화
PCA와 유사하게, 계층적 군집화는 데이터 세트의 패턴과 잠재적 이상값을 식별하는 또 다른 보완적인 방법입니다. 히트맵은 데이터 세트의 모든 쌍 샘플 조합의 유전자 발현 상관관계를 보여줍니다. 대부분의 유전자는 차이 발현하지 않기 때문에 샘플 간에는 일반적으로 높은 상관관계(0.80 이상의 값)가 있습니다. 0.80 미만의 샘플은 데이터 및/또는 샘플 오염에 이상값이 있을 수 있음을 나타냅니다.
축을 따르는 계층적 나무는 어떤 샘플이 서로 더 유사한지, 즉 군집을 형성하는지를 나타냅니다. 상단의 색상 블록은 데이터의 하위 구조를 나타내며, 우리는 각 샘플 그룹으로서 반복들이 함께 군집을 형성하기를 원합니다. 우리의 기대는 PCA 플롯에서 관찰한 군집과 유사한 방식으로 샘플이 군집을 형성하는 것입니다.
아래 그림에서, Wt_3와 KO_3 샘플은 다른 반복과 함께 군집을 형성하지 않습니다. 우리가 동일한 샘플 군집을 보는지 확인하기 위해 PCA를 탐색하고 싶습니다.
- Mov10 QC
이제 우리가 RNA-seq에 일반적으로 사용되는 QC 단계를 잘 이해했으므로, Mov10 데이터 세트에 대해 QC를 수행해 보겠습니다.
5.1. 데이터 변환
- 정규화 카운트 변환
PCA 및 계층적 군집화 시각화 방법의 거리 또는 군집을 촉진하기 위해 정규화 카운트에 rlog 변환을 적용하여 평균의 분산을 조절해야 합니다.
정규화 카운트의 rlog 변환은 이러한 시각화 방법에 대한 품질 평가 기간 동안에만 필요합니다. 우리는 이 변환된 카운트를 사용하여 차이 발현을 결정하지 않습니다.
# rlog 변환
rld <- rlog(dds, blind=TRUE)
blind=TRUE 매개변수는 rlog() 함수가 샘플 그룹을 고려하지 않도록 하는 것입니다—즉, 편향 없는 방식으로 변환을 수행합니다. 품질 평가를 수행할 때 이 옵션을 포함하는 것이 중요합니다.
rlog() 함수는 DESeqTransform 객체를 반환합니다. 이는 DESeq에 특정된 또 다른 객체 유형입니다. 단순히 변환된 값 행렬을 얻는 이유는 rlog 변환을 계산하는 데 사용되는 모든 매개변수(즉, 크기 인자)가 해당 객체에 저장되기 때문입니다. 우리는 품질 평가를 위해 PCA 및 계층적 군집화 플롯을 그리기 위해 이 객체를 사용합니다.
5.2. PCA
이제 우리는 QC 단계를 준비했으므로 PCA부터 시작해 보겠습니다!
DESeq2에는 ggplot2를 백그라운드에서 사용하여 PCA 플롯을 생성하는 내장 함수가 있습니다. 이는 우리가 코드 줄을 입력할 필요가 없고 다양한 ggplot2 레이어를 다룰 필요가 없게 해주기 때문에 훌륭합니다. 또한 rlog 객체를 직접 입력으로 받으므로 관련 정보를 추출하는 번거로움을 덜어줍니다.
plotPCA()는 두 개의 매개변수를 입력으로 필요로 합니다: DESeqTransform 객체와 intgroup, 즉 메타데이터에 실험 샘플 그룹 정보가 포함된 열의 이름입니다.
# PCA 플롯
plotPCA(rld, intgroup="sampletype")
기본적으로 plotPCA()는 가장 변이성이 높은 500개의 유전자를 사용합니다. ntop= 매개변수를 추가하고 함수가 고려할 유전자 수를 지정하여 이 설정을 변경할 수 있습니다.
plotPCA() 함수는 PC1과 PC2의 값만 반환합니다. 데이터의 다른 PC를 탐색하거나 PC에 기여하는 유전자를 식별하려면 prcomp() 함수를 사용할 수 있습니다. 예를 들어, 모든 PC를 플롯하려면 다음 코드를 실행할 수 있습니다:
# 입력은 로그 변환된 값의 행렬입니다
rld <- rlog(dds, blind=T)
rld_mat <- assay(rld)
pca <- prcomp(t(rld_mat))
# ggplot에 입력할 메타데이터와 PC3 및 PC4 값이 있는 데이터 프레임 생성
df <- cbind(meta, pca$x)
ggplot(df) + geom_point(aes(x=PC3, y=PC4, color = sampletype))
5.3. 계층적 군집화
MOV10데이터 세트 계층적 군집화
DESeq2에는 모든 샘플 간의 쌍 상관관계와 계층적 군집화 정보를 보여주는 히트맵을 그리는 내장 함수가 없습니다; 우리는 pheatmap 패키지의 pheatmap() 함수를 사용할 것입니다. 이 함수는 DESeqTransform 객체를 입력으로 사용할 수 없지만 행렬이나 데이터 프레임이 필요합니다. 따라서 해야 할 첫 번째 작업은 assay()라는 함수를 사용하여 rld 객체에서 이 정보를 검색하는 것입니다. 이 함수는 DESeqTransform 객체의 데이터를 단순한 2차원 데이터 구조로 변환합니다.
# 객체에서 rlog 행렬 추출
rld_mat <- assay(rld)
다음으로 모든 샘플의 쌍 상관 값을 계산해야 합니다. cor() 함수를 사용하여 이를 수행할 수 있습니다:
# 쌍 상관 값 계산
rld_cor <- cor(rld_mat)
상관 행렬의 열 이름과 행 이름을 살펴보겠습니다.
head(rld_cor)
head(meta)
샘플에 대해 제공한 이름이 시작 시 메타데이터 데이터 프레임에서 사용한 이름과 일치하는 것을 알 수 있습니다. 이는 아래 주석 매개변수를 사용하여 상단에 색상 블록을 그릴 수 있도록 중요합니다. 이 블록은 계층적 군집화의 시각화를 쉽게 구현할 수 있습니다.
# pheatmap 패키지 로드
library(pheatmap)
# 상관 행렬과 메타데이터 객체를 사용하여 히트맵 플롯
pheatmap(rld_cor, annotation = meta)
pheatmap()로 플롯을 그릴 때 계층적 군집화 정보는 유사한 샘플을 함께 배치하는 데 사용되며, 이 정보는 축을 따르는 나무 구조로 표시됩니다. 주석 매개변수는 입력으로 데이터 프레임을 받으며, 우리의 경우에는 메타데이터 프레임입니다.
전반적으로 우리는 이상 샘플이 없음을 나타내는 높은 상관관계(> 0.999)를 관찰했습니다. 또한 PCA 플롯과 유사하게 샘플이 샘플 그룹으로 군집을 형성하는 것을 볼 것입니다. 결론적으로 이러한 플롯은 데이터 품질이 양호하다는 것을 보여주며 우리는 차이 발현 분석을 진행할 수 있는 자신감을 가지고 있습니다.