시계열 예측에서 기존의 LSTM, GRU 기반 모델들은 단일 시점 예측에는 강한 성능을 가지지만, 미래의 여러 시점 (Multi-Horizon)을 동시에 예측하는데에는 한계가 있다.

현재 Forecasting분야의 SOTA(State of the art) 모델인 Temporal Fusion Transformer모델은 Transformer의 Attention 메커니즘을 도입하여 장기 의존성(Long-term dependency)을 포착하고, 다양한 형태의 변수를 처리할 수 있다. 특히, 딥러닝 모델의 고질적인 문제인 블랙박스 문제를 해결하고 해석 가능성을 제공한다는 점에서 유통, 헬스 케어, 금융 등 신뢰성이 중요한 분야에서 큰 잠재력을 가지고 있다.

Idea

Multi-Horizon Forecasting을 위한 TFT모델의 Input으로 다음 변수들이 필요하다.

Observed Input : 관측을 통해 알 수 있는 값이지만, 현 시점에서 미래의 값은 알 수 없는 변수. (Target) Time Vary Known Input : 시간에 따라 달라지지만, 현재에도 혹은 미래에도 그 값을 알 수 있는 변수 (Week, Holiday 등) Static Covariates : 시간과 관계없이 변하지 않는 정적 공변량 (성별, 나이 등)

기존 다른 모델들은 이 Static 데이터들을 제대로 활용하기 힘들었지만, TFT 는 이 데이터를 모델 전체에 Context로 뿌려 활용한다는 특징이 있다 .

Architecture

Gating Mechanism(GLU) 데이터의 흐름을 제어하여 현 시점에서 불필요한 정보는 억제한다.

Variable Selection Network(VSN) 수많은 Feature 중에서, 어떤 변수가 예측에 중요한지 스스로 선택한다.

Static Covariate Encoder 정적 변수를 인코딩 하여 시간적 변수와 결합한다.

Temporal Fusion Decoder LSTM의 지역적인 패턴과 Multi-Head Attention을 결합하여 최종적인 예측을 한다.

Interpretability

TFT의 장점은 예측의 이유를 시각화 할 수 있다는 것이다. 주로 사용하는 차트는 Variable Importance(변수 중요도)와 Attention Weights이다.

Attention Weights 에서, 모델은 예측을 위해 어느 시점의 Attention에 가중치를 두엇는지 측정할 수 있으며, 여러 Variables 중 어느 변수에 더 가중치를 두엇는지 확인할 수 있다.

또, TFT는 단순히 하나의 값을 예측하는 것이 아니라, Quantile Loss를 사용하여 예측의 범위(신뢰구간)을 제공한다. 이를 통해 예측의 불확실성 까지 파악할 수 있어 리스크 관리에 용이하다.

Implementation

import torch
from pytorch_forecasting import TemporalFusionTransformer, TimeSeriesDataSet

# 1. 데이터셋 정의 
training = TimeSeriesDataSet(
    data[lambda x: x.time_idx <= training_cutoff],
    time_idx="time_idx",
    target="vital_sign", # 예측하고자 하는 변수 
    group_ids=["patient_id"],

    ...

    static_categoricals=["gender"], # 정적 변수  
    time_varying_known_reals=["time_of_day"], # 미래에도 아는 변수
    time_varying_unknown_reals=["heart_rate", "bp_systolic"], # 관측 변수

    ...
)

# 2. 모델 설정 
tft = TemporalFusionTransformer.from_dataset(
    training,
    learning_rate=0.03,
    hidden_size=16,
    attention_head_size=1,
    dropout=0.1,
    output_size=7, 
    loss=QuantileLoss(),
)

# 3. 학습 ...

Molecular Niche란 특정 세포의 생존, 분화, 그리고 기능을 조절하는 국소적인 분자적 환경을 의미한다. 단순히 물리적인 위치를 넘어, 세포가 노출된 화학적 신호, 대사 물질, 물리적 지지체 등의 총체적인 집합이다.

2. Components of Molecular Niche

수용성 인자

• Cytokines & Chemokines : 면역 반응과 세포의 Migration을 유도한다.

• Growth Factors(성장 인자) : 세포의 증식을 촉진한다. (ex. VEGF, EGF)

• Metabolites(대사 산물) : 세포 대사에 직접적인 영향을 주는 저분자 화합물이다.

  세포 간 상호작용 (Cell to Cell Interaction, CCI)

  이웃한 세포가 직접적인 접촉이나 근거리 신호 전달을 통해 형성하는 Niche이다.

• Juxtacrine Signaling : 세포막에 붙어 있는 단백질(Ligand)와 옆 세포의 수용체(Receptor)가 직접 결합하는 방식이다.

• Paracrine Signaling : 가까운 거리의 세포에게 신호 물질을 분비하는 방식이다.

  세포외 기질 (Extracellular Matrix, ECM)

  세포가 부착되어 있는 구조적 지지체이다.

• Composition : Collagen, Fibronectin, Laminin 등으로 구성된다.

• Mechanotransduction(기계적 신호 전달) : ECM의 강도나 탄성 자체가 세포에게 물리적 신호를 주어 유전자 발현을 바꿀 수 있다.

• Heart Failure에서는 ECM이 과도하게 쌓이는 Fibrosis가 일어나는데, 이때 변한 경도 자체가 Niche가 되어 심장 기능을 약화시킨다.

  3. Bioinfomatics에서의 Niche

  1. Cell to Cell Interaction

  scRNAseq 데이터는 세포의 위치정보를 소실하지만, 이를 CellChat, NicheNet등의 툴을 사용하여 주고받는 신호를 복원할 때 사용한다.

2. 공간 전사체 (Spatial Transcriptomics Analysis)

최근 유망한 분야로, 세포의 위치 정보(Spatial Coordinates)가 있을 때 Niche는 각 구역을 구분하는 기준이 된다.

3. 세포의 분화 궤적 (Trajectory Inference)

외부의 신호에 의해 다른 산물을 만드는 경우, 해당 신호를 분석하기 위해 Niche를 사용한다.

KNN 기본 원리

새로운 데이터가 주어졌을 때, 기존 데이터 중에서 가장 가까운 $k$ 개의 이웃을 찾는다. 그 이웃의 label을 확인하여 새로운 데이터의 label을 결정한다.

• Classification : $k$ 개의 이웃 중 가장 많이 등장하는 label로 분류한다.
• Regression : $k$개의 이웃들이 가진 값의 평균을 구한다.

KNN Algorithm

1. $k$ 설정 : 몇 개의 이웃을 살펴볼지 결정한다
2. 거리 계산 : 새로운 데이터와 기존의 모든 데이터 사이의 거리를 주로 계산한다.
   • 변수가 연속형 변수일 경우, 주로 유클리드 거리(Euclidean Distance) 또는 맨하튼 거리(Manhattan Distance) 가 많이 사용된다.
     • Euclidian distance : 가장 가까운 직선 거리 $$d(p,q) = \sqrt{\sum_{i = 1}^{n}(q_i - p_i)^2}$$
     • Manhattan : 길이 존재하는 최단 루트 거리 $$d(p,q) = \sum_{i=1}^{n}{|q_i - p_i|^2}$$
       • 변수가 범주형 변수일 경우, Hamming distance를 주로 사용한다. $$d(p,q) = \sum_{i=1}^{n}{I(p_i \neq q_i)}$$
3. 이웃 선택 : 계산된 거리를 기준으로 가장 가까운 $k$ 개의 데이터를 선택한다.
4. 결과 도출 : 분류 혹은 회귀에 따라 결과를 도출한다.

$k$ 값의 중요성 (Hyperparameter)

$k$ 값을 어떻게 설정하느냐에 따라 성능이 크게 달라진다.

• $k$ 가 너무 작을 때
  • 가장 가까운 데이터에만 의존하므로, Noise & Outlier에 민감해진다.
  • Overfitting 될 위험이 크다.
• $k$ 가 너무 클 때
  • 데이터의 특성을 잃어버릴 수 있다.
  • Underfitting 될 위험이 크다.

Anndata(Annotated Data) 는 유전자 발현량과 같은 행렬 데이터와 그에 대한 주석(Annotation)을 하나로 묶어 관리하기 위해 설계된 데이터 구조이다.

아래에서 설명될 내용이지만, 만약 obs에서 특정 조건에 의해 필터링하면, 다른 계층에 있는 (X, obsm, layers 등) 모든 연관 데이터가 자동으로 함께 잘려나간다.

2. Anndata 구조

1. X - Main matrix

• 실제 데이터(발현량 등)가 저장되는 영역
• 보통 Cells X Genes 의 크기를 가지며, 메모리 절약을 위해 0이 많은 경우 Sparse Matrix형태로 저장된다.

2. obs - Observations

• 행에 관한 정보인 Cell에 대한 annotated가 담겨있다.
• 보통 세포ID, 샘플 그룹, 클러스터 번호, 세포 종류가 기록된다.

3. var - Variables

• 열에 대한 정보인 Gene에 대한 annotated가 담겨있다.
• 보통 유전자 이름, 미토콘드리아 유전자 여부 등이 기록되어있다.

4. obsm, varm - Multi-dimensional Annotation

• 2차원 이상의 다차원 데이터를 저장한다.
• obsm의 경우 PCA, UMAP, t-SNE 좌표값처럼 각 세포가 여러 값을 갖는 데이터가 들어간다.
• varm은 각 유전자가 가지는 다차원 정보를 가진다.

5. uns - Unstructured

• 행렬 형태가 아닌 비정형 데이터를 저장한다.
• 예시로 분석에 사용한 파라미터 정보, 클러스터별 지정 색상 등

6. layers - Data의 여러 버전 관리

• 동일한 차원의 데이터를 여러 버전으로 저장할 수 있다.

7. obsp - Pairwise Annotation

• 세포와 세포사이의 관계를 저장하는 행렬이다.

3. 살펴보기

다음은 본인이 실제로 Covid19 - Single cell RNA seq 를 진행하던 anndata를 살펴보았다.

print(adata)

AnnData object with n_obs × n_vars = 16716 × 15002
    obs: 'n_genes', 'samples', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'doublet_score', 'predicted_doublet', 'leiden05', 'cell_type', 'condition'
    var: 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'mean', 'std'
    uns: 'cell_type_colors', 'dendrogram_leiden05', 'hvg', 'leiden05', 'leiden05_colors', 'leiden05_rank', 'log1p', 'neighbors', 'pca', 'scrublet', 'umap'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'counts', 'lognorm'
    obsp: 'connectivities', 'distances'

1. library & data load

library(dplyr) # 데이터 집계 및 처리용 
library(Seurat) # scRNAseq
library(patchwork) # 그래프

# 2700개의 Peripheral Blood MononuclearCells(PBMC) data Load
pbmc.data <- Read10X(data.dir= "/Users/juno/Documents/GitHub/seurat_prac/filtered_gene_bc_matrices/hg19")

# Seurat Object 지정 
pbmc <- CreateSeuratObject(
  counts = pbmc.data, # row data 지정
  project = "pbmc3k", # 프로젝트 이름
  min.cells = 3, #최소 3개 이상의 세포에서 발견된 유전자만 남김
  min.features = 200 #유전자가 최소 200개 이상 검출된 세포만 남김
)

• dplyr : 데이터 가공
• Seurat : scRNA-seq
• patchwork : 그래프 배치 그리고 CreateSeuratObject를 활용해서 Seurat object를 생성하는데, min.cells와 min.features를 이용하여 너무 적은 유전자만 가진 세포나, 너무 적은 세포에서만 발현된 유전자를 1차로 필터링 한다.

2. Quality Control

데이터의 노이즈를 필터링하여 분적의 정확도를 높인다. 여기서, 노이즈는 다음 경우들로 간주한다.

1. 죽은 세포 : 미토콘드리아 유전자(^MT-) 비율이 비정상적으로 높은 세포

• 싱글 셀 실험과정 중에 세포막이 깨지면(Lysis), 세포질 RNA는 유실되지만 미토콘드리아는 막으로 싸여있어 RNA가 보존되므로 죽은 세포 이거나 상태가 안좋은 세포의 경우 상대적으로 미토콘드리아 유전자 비율이 높게 잡힌다.

2. Empty Droplets : 유전자 수가 너무 적은 경우
3. Doublets(두 세포가 뭉친 경우) : 유전자 수가 비정상적으로 많은 경우

   # 미토콘드리아 비율 계산
   pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
   

QC 지표 시각화

VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

필터링 (유전자 수 200~2500 사이, 미토콘드리아 < 5%)

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

![](https://velog.velcdn.com/images/door_jono/post/9c3a8a38-c3b7-42ac-a9ad-5ded9d5d45f7/image.png)
- `nFeature_RNA` : 세포 하나에서 검출된 유전자의 가짓수
- `nCount_RNA` : 세포 하나에서 검출된 RNA 분자의 총 개수(UMI count)
- `percent.mt` : 전체 유전자 중 **미토콘드리아 유전자**가 차지하는 비율

# 3. Normalization & Feature Selection
세포마다 읽힌 RNA의 양이 다르기 때문에 정규화를 진행한 후 분석에 유의미한 **Feature gene**을 찾는다.

Normalizing the data

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

위의 파라미터들은 기본 값으로, pbmc <- NormalizeData(pbmc) 로 적용 가능

highly variable features 식별

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

top10 <- head(VariableFeatures(pbmc), 10)

plot1 <- VariableFeaturePlot(pbmc) plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) plot1 + plot2

![](https://velog.velcdn.com/images/door_jono/post/d79d5fb7-eb4b-4e8f-9dd9-e4b0dea439e0/image.png)

# 4. Scaling & PCA
데이터를 표준화(Scaling) 한 후, 주성분 분석(PCA)을 통해 데이터의 복잡도를 줄인다.
- `ScaleData()` : 유전자 발현량의 평균을 0, 분산을 1로 조절한다. 발현량이 매우 높은 특징 유전자가 분석을 방해하는 것을 방지한다.
- `RunPCA()` : 수천 개의 유전자 정보를 50개 정도의 주성분으로 압축한다. 이를 통해 데이터의 주요 패턴을 요약 할 수 있다. 

PCA를 할수있도록 데이터의 차원을 축소한다.

mean == 0, variance == 1

all.genes <- rownames(pbmc) pbmc <- ScaleData(pbmc, features = all.genes)

scaling이 끝난 데이터를 가지고 PCA (주성분 분석) 시작

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 5. Determine Dimensionality
PCA를 통해 압축된 주 성분 중 몇 개를 분석에 사용할 지 결정한다.
전통적으로 아래의 `JackStraw` 를 사용했지만 요즘은 `ElbowPlot`을 사용한다.

우리의 데이터는 여러 유전자의 발현도를 다양한 각도로 분석한 데이터를 가지고 있기 때문에

이것을 그래프로 표현하기 위해 차원을 축소해나가야 한다.

JackStraw 를 활용한 data resampling (구버전)

pbmc <- JackStraw(pbmc, num.replicate = 100) pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

총 데이터에서 랜덤으로 1% 묶어 PCA를 돌려 차원을 축소하고, 그걸 이용해 null distribution을 만든다.

시간이 오래 걸린다 .

JackStrawPlot(pbmc, dims = 1:15)

ElbowPlot(pbmc)

![](https://velog.velcdn.com/images/door_jono/post/f03aa337-dc3f-47a8-b7fd-adc7fc7e9fa9/image.png)
선형적인 추세를 가지는 PC가 7~12 사이 임을 알수있다. 그러므로 분석에 유의미한 PC는 처음부터 7~12 까지이다.
![](https://velog.velcdn.com/images/door_jono/post/b803d050-87f4-4d80-9b2f-a98f4e39bc00/image.png)
아래의 `ElbowPlot`을 보았을때, 감소 추세가 줄어드는 지점이 역시 PC 7~10 부근이므로, 역시 분석에 유의미한 PC는 0에서 부터 7~10까지이다.

# 6. Clustering 
이제 우리는 PC를 10으로 기준점 삼아 유사한 유전자 발현 패턴을 가진 세포들끼리 Clustering을 할 것이다.

위에서 결정한 첫 10개의 PC만을 사용한다 .

pbmc <- FindNeighbors(pbmc, dims = 1:10) # PCA공간에서의 유클리드 거리 기반 KNN 그래프 구성

pbmc <- FindClusters(pbmc, resolution = 0.5) # resolution 은 세분화 정도를 지정.

- `FindNeighbors` : PCA공간에서 세포 간의 거리를 계산하여 그래프를 만든다.
- `FindClusters` : Louvain 알고리즘을 이용하여 세포들을 그룹화 한다.

# 7. 시각화 (UMAP / tSNE)
고차원 데이터를 우리가 볼 수 있는 2차원 평면에 투사한다.
- **UMAP** : 클러스터링 된 세포들이 2차원 상에서 끼리끼리 모여 있는지 시각적으로 확인한다. 가장 대중적인 시각화 방법이다.

비선형 차원 축소 (UMAP and tSNE)

클러스터링 결과를 시각화하고 데이터의 기본 구조 탐색

pbmc <- RunUMAP(pbmc, dims = 1:10) DimPlot(pbmc, reduction = "umap")

![](https://velog.velcdn.com/images/door_jono/post/cd5c4e1e-ac5d-4267-a6f3-85ec58890b5c/image.png)

# 8. Finding Biomarker
각 클러스터가 생물학적으로 어떤 **Cell Type** 인지 알아내기 위해, 해당 클러스터에만 특이적으로 많이 발현한 유전자를 찾는다.

모든 클러스터를 돌며 바이오마커를 발견한다.

pbmc.markers <- FindAllMarkers(pbmc,only.pos = TRUE)

발현량이 증가한 유전자만 찾는다.

pbmc.markers %>% group_by(cluster) %>% dplyr::filter(avg_log2FC > 1)

특정 마커의 Violin Plot

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

you can plot raw counts as well

VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE) FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP","CD8A"))

Heatmap

pbmc.markers %>% group_by(cluster) %>% dplyr::filter(avg_log2FC > 1) %>% slice_head(n = 10) %>% ungroup() -> top10 DoHeatmap(pbmc, features = top10$gene) + NoLegend()

![](https://velog.velcdn.com/images/door_jono/post/19e29230-069f-4de0-ab1e-0c13b888c9ae/image.png)
- `MS4A1 & CD79A` 는 면역 세포 중 **B세포** 에서만 특이적으로 발현되는 마커이다.
- `cluster 3`에서 특징적으로 분포함을 알 수 있어, `cluster 3`는 B세포임을 알 수 있다.

![](https://velog.velcdn.com/images/door_jono/post/73500517-4aca-4d1d-aa66-929df47ffc95/image.png)
- `NKG7`은 세포를 죽이는 독성 물질 유전자로, **NK세포** 혹은 **CD8+ T세포**에서 주로 나온다.

    - 주로 4,5,6번 클러스터에서 높게 나타나는데, 그중 `cluster 6`이 아주 높고 밀도가 높아 **NK 세포**임을 알 수 있으며, `cluster 4`가 CD8 T세포일 가능성이 높으나 정확한 구분을 위해 추가적인 마커 확인이 필요하다.

- `PF4`는 혈소판 인자로, `cluster 8`에서 높은 분포를 가지므로, 혈소판임을 알 수 있다.

![](https://velog.velcdn.com/images/door_jono/post/8469cd6f-2f74-42ae-a9c1-b1c758a3c16e/image.png)
이는 마커들이 위에서 봣던 UMAP상에서 어디에 위치하는가를 표현하는 **Feature Plot**이다.
- `CD3E`는 T Cell의 마커인데, plot을 보았을때, 우측 상단의 영역에 많이 분포함을 알 수 있다.

이런식으로 생물학적 특징이 기하학적 위치와 어떻게 매칭되는지 알 수 있다.
![](https://velog.velcdn.com/images/door_jono/post/f7f42251-e4e9-4579-a619-20909b63823c/image.png)
이는 **Heatmap**으로 모든 분석의 최종 기록이라고 볼 수 있다.
- 가로축 : Cell이며, 위쪽의 숫자가 클러스터 번호이다.
- 세로축 : Gene이며, 각 클러스터를 대표하는 Top10 마커 유전자들이다.


- 노란색 : 유전자 발현량이 높음
- 보라색/검은색 : 유전자 발현량이 낮음 혹은 없음


- 읽는법 
    - `cluster0,2` 의 경우 최상단의 `LDHB`,`CD3D`,`CD3E` 등 T세포의 마커들의 유전자 발현량이 높다. 전반적으로 비슷한 경향성을 가지고 있지만, `CCR7` 마커의 강도차이 등으로 인해 알고리즘이 다른 클러스터로 분리해둠을 알 수 있다. (`cluster0`이 Naive T, `cluster 2`가 Memory T로 분리된다.)
    - `cluster 1`의 경우 전형적인 **단핵구(Monocyte)**의 패턴이다.
    - `Cluster 3`의 경우 확실한 **B세포**이다.
- 위 히트맵과 같이 노란 영역이 대각선에 집중된다면 각 클러스터가 서로 겹치지 않는 고유한 특징을 가지고 있다는 뜻이다.


# 9. Annotation
마지막으로, 찾아낸 마커 유전자를 마탕으로 클러스터에 이름표를 붙여준다.

cell type identity to cluster

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet") names(new.cluster.ids) <- levels(pbmc)

숫자로 된 cluster ID를 실제 세포 유형 이름으로 변환

pbmc <- RenameIdents(pbmc, new.cluster.ids) DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

``` RenameIdents함수를 이용하여 숫자 ID를 실제 이름으로 바꿔줌을 알 수 있다.

우리는 전 게시글에서 ERR458493 의 Saccharomyces cerevisiae 샘플을 정렬하였다.

정렬의 결과물 중 _ReadsPerGene.out.tab 파일은 각 유전자에 몇 개의 시퀀싱 Read가 Mapping 되었는지 그 Count를 정리한 데이터 인데, 해당 데이터로 DEG 분석을 진행해보고자 한다.

우선적으로, ERR458493 는 Wild Type(WT) 샘플이며, 추가적으로 1개의 WT 샘플과 2개의 snf2돌연변이(MU) 샘플을 준비해 총 4개의 RNA-seq 샘플을 대상으로 분석을 진행하였다.

• 분석 대상
  • 대조군(WT) : ERR458493 , ERR459175
  • 실험군(MU) : ERR458506 , ERR458515

해당 파일들을 각각 4회 STAR 정렬하여 , 4개의 ...ReadsPerGene.out.tab 파일이 있다는 전제로 진행하겠다.

2. Gene Count Matrix 생성

STAR의 --quantMode GeneCounts 옵션은 각 샘플에 대해 독립적인 유전자 발현량(Read Count) 파일을 생성한다. DEG 분석을 수행하기 위해 이 4개의 데이터들을 하나의 통합 행렬로 병합하여야한다. 이 행렬은 유전자(행) X 샘플(열)을 구조로 가진다.

2-1. ReadsPerGene.out.tab 살펴보기

N_unmapped    44500    44500    44500
N_multimapping    111531    111531    111531
N_noFeature    35757    445692    474155
N_ambiguous    60129    3329    3260
YDL246C    0    0    0
YDL243C    28    14    14
YDR387C    59    27    32
YDL094C    2    3    2
YDR438W    27    15    12
YDR523C    2    1    1
...

다음은 ERR458493_ReadsPerGene.out.tab의 내용인데, 최 상단 4줄은 통계에 대한 요약 행이고, 그 아래로 GeneID, Unstranded Count, Starnded(Forward) Count, Strand(Reverse) Count 를 열로 가지는 데이터가 있다.

2-2. Count Matrix 생성

import pandas as pd
import os

BASE_DIR = "./alignment_output"

samples = {
    "WT_1": "ERR458493_ReadsPerGene.out.tab",
    "WT_2": "ERR459175_ReadsPerGene.out.tab",
    "MU_1": "ERR458506_ReadsPerGene.out.tab",
    "MU_2": "ERR458515_ReadsPerGene.out.tab"
}

counts_list = []

print("Create Matrix ...")

for sample , filename in samples.items():
    file_path = os.path.join(BASE_DIR, filename)

    if not os.path.exists(file_path):
        print(f"ERROR: {file_path} not found")
        continue
    # 상단 4줄은 통계 요약이니 제외
    df = pd.read_csv(file_path, sep='\t', header=None, skiprows=4)

    # GeneID, Count만 추출 후 rename
    df = df.iloc[:, [0,1]]
    df.columns = ['GeneID', sample]
    df.set_index('GeneID', inplace=True)

    counts_list.append(df)
    print(f"-{sample} load complete.")

if counts_list :
    count_matrix = pd.concat(counts_list, axis = 1)
    count_matrix.to_csv("count_matrix.csv")
    print("COMPLETE : save [count_matrix.csv]")
    print(count_matrix.head())
else : 
    print("ERROR : data not founds")

해당 코드의 결과로 4개의 카운트가 하나로 병합된 count_matrix.csv가 생성된다.

GeneID,WT_1,WT_2,MU_1,MU_2
YDL246C,0,0,0,0
YDL243C,28,31,52,33
YDR387C,59,56,131,103
YDL094C,2,0,6,2
YDR438W,27,23,67,46
YDR523C,2,4,10,3
YDR542W,1,0,0,1
YDR492W,57,70,159,64
YDR018C,5,6,39,31
...

이제, 이 파일을 가지고 DEG 분석을 실행 할 수 있다.\

3. 차등 발현 유전자 분석

3.1 Workflow

우리는 DEG분석을 위해 PyDESeq2 를 사용할 것인데, 이는 RNA-seqCount 데이터의 특성을 고려하여 통계적으로 유의미한 DEG를 식별하는 표준 도구이다. 정규화 - 분산 추정 - 가설 검정 - 다중 검정 보정 의 과정을 거치며 해당 내용은 추후에 게시글을 업로드 하겠다.

3.2 in Python code

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

print("Load data...")

counts_df = pd.read_csv("count_matrix.csv", index_col = 0)

# PyDESeq2 는 sample X gene 형태를 요구한다.
counts_df = counts_df.T
counts_df = counts_df.fillna(0)
counts_df = counts_df.astype(int)

# 샘플 합계가 10 미만인 노이즈 유전자 제거
counts_df = counts_df.loc[:, counts_df.sum(axis=0) >= 10]
print(f"-> analsis target : {counts_df.shape[0]} sample, {counts_df.shape[1]} gene")

# Create Metadata
metadata = pd.DataFrame(index=counts_df.index)
metadata['condition'] = ['WT' if 'WT' in sample else 'MU' for sample in metadata.index]
print("-> Metadata:")
print(metadata)

# DESeq2
print("analsis DESeq2...")

dds = DeseqDataSet(
    counts=counts_df,
    metadata=metadata,
    design_factors='condition', # 기준 열
    ref_level=["condition","WT"], # Reference is WT
    n_cpus=4
)

dds.deseq2() # 정규화 및 분산 추정

# 통계분석
stat_res = DeseqStats(dds, contrast=["condition","MU","WT"], n_cpus=4)
stat_res.summary()

res = stat_res.results_df
res.to_csv("DESeq2_results.csv")
print("save to [DESeq2_results.csv]")

# 시각화
print("\n3. 볼케이노 플롯 생성 중...")

# 시각화를 위한 데이터 준비 (p-value 없는 것 제거)
plot_data = res.dropna(subset=['padj', 'log2FoldChange']).copy()

# 색상 기준 설정
# Red: 유의미하게 증가 (Up) / Blue: 유의미하게 감소 (Down) / Grey: 의미 없음
conditions = [
    (plot_data['log2FoldChange'] > 1) & (plot_data['padj'] < 0.05),
    (plot_data['log2FoldChange'] < -1) & (plot_data['padj'] < 0.05)
]
choices = ['red', 'blue']
plot_data['color'] = np.select(conditions, choices, default='lightgrey')

# 그래프 그리기
plt.figure(figsize=(10, 8))

# 모든 점 찍기
plt.scatter(
    plot_data['log2FoldChange'], 
    -np.log10(plot_data['padj']), 
    c=plot_data['color'], 
    alpha=0.6, 
    s=15
)

# 기준선 그리기
plt.axvline(x=1, color='black', linestyle='--', linewidth=0.8)  # 2배 증가
plt.axvline(x=-1, color='black', linestyle='--', linewidth=0.8) # 2배 감소
plt.axhline(y=-np.log10(0.05), color='black', linestyle='--', linewidth=0.8) # p-value 0.05

plt.title("Volcano Plot: Mutant vs WT (PyDESeq2)", fontsize=15)
plt.xlabel("log2 Fold Change", fontsize=12)
plt.ylabel("-log10 Adjusted P-value", fontsize=12)

top_genes = plot_data.sort_values('padj').head(5)
for gene_name, row in top_genes.iterrows():
    plt.text(
        row['log2FoldChange'], 
        -np.log10(row['padj']) + 0.2, 
        gene_name, 
        fontsize=9, 
        ha='center',
        color='black',
        fontweight='bold'
    )

# 저장 및 출력
plt.tight_layout()
plt.savefig("Volcano_Plot_Python.png", dpi=300)
print("[완료] 그래프가 'Volcano_Plot_Python.png'에 저장되었습니다.")
plt.show()

4. 결론

위 단계를 통해 우리는 DEG 결과인 DESeq2_results.csv와 Volcano_plot을 얻었다.

4.1 .csv

GeneID,baseMean,log2FoldChange,lfcSE,stat,pvalue,padj
YDL243C,35.23241178557316,-0.46323511752066077,0.4113790443156044,-1.1260542410256393,0.26014256106731215,0.5698109731025939
YDR387C,81.22576980096602,0.04681292104056832,0.2963043850078386,0.15798929549871463,0.8744652337332883,0.9598669278117801
...

해당 파일에서 가장 중요한 열은 GeneID, log2FoldChange,padj 이다. 해당 파일에서 padj < 0.05 && |log2FoldChange| > 1 인 유전자를 필터링 하여 유의미한 최종 DEG 목록을 확인할 수 있다.

4.2 Volcano_plot

이 이미지에서 파란색 그룹은 DOWN그룹으로 유의미하게 발현이 감소한 그룹을 의미하고 빨간색 그룹은 UP 그룹으로 유의미하게 발현이 증가한 그룹을 의미한다. 이 최종 DEG 목록(UP 그룹, DOWN 그룹)을 유전자 기능 강화 분석 (GO Enrichment Analysis) 또는 경로 분석 (KEGG Pathway Analysis)에 사용하여 생물학적 의미를 해석하는 후속 연구로 이어질 수 있다.

우리는 효모(Saccharomyces cerevisiae)의 유전자 정보들을 STAR를 이용해 정렬하기 위해 다음과 같은 파일이 필요하다.

1. 참조 유전체 파일 (.fasta) : 효모의 전체 DNA 구성이 적힌 파일으로서 STAR는 해당 파일을 기준으로 정렬해나간다.
   • Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa
2. 유전자 주석 파일 (GTF) : 참조 유전체 파일에서 실제 유전자 (Exon, Intron)의 위치가 어디인지 좌표를 기록한 파일이다.
   • Saccharomyces_cerevisiae.R64-1-1.111.gtf
3. RNA-Seq Read (.Fastq) : 우리가 분석할 실제 시퀀싱 데이터이다.
   • ERR458493.fastq

2. 인덱싱 단계

STAR \
--runMode genomeGenerate \
--runThreadN 2 \
--genomeDir ./genome_index \
--genomeFastaFiles Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa \
--sjdbGTFfile Saccharomyces_cerevisiae.R64-1-1.111.gtf \
--sjdbOverhang 99 \
--genomeSAindexNbases 10

• --runMode genomeGenerate : STAR의 실행 모드는 genomeGenerate(유전자 인덱스 생성) 과 alignReads(정렬) 로 구분된다.

• --genomeDir : 인덱스 파일들의 출력 디렉토리를 지정한다.

• --genomeFastaFiles : 인덱싱의 기본 입력 1 파일로써, 참조 유전체 파일을 지정한다.

• --sjdbGTFfile : 인덱싱의 기본 입력 2 파일로써, GTF파일을 지정한다. STAR는 이. 파일을 읽어 Splice Junctions(SJ) 목록을 미리 만들어 둔다.

• --sjdbOVerhang : SJ를 만들때, Exon의 경계에서 양쪽으로 몇 개의 염기를 더 읽어서 인덱스에 포함할 것인가?이다.

  • 일반적으로 리드 길이 - 1로 설정하는 것이 표준이다.

• --genomeSAindexNbase : SA인덱스의 정밀도를 조절하는 메모리 관련 옵션이다.

3. 정렬 단계

STAR \
--runMode alignReads \
--runThreadN 2 \
--genomeDir ./genome_index \
--readFilesIn ERR458493.fastq \
--outFileNamePrefix ./alignment_output/ERR458493_SE_ \
--outSAMtype BAM Unsorted \
--quantMode GeneCounts

• --runThreadN : 정렬 작업을 수행할 때 사용할 CPU 스레드 수를 지정한다.

• --genomeDir ./genome_index : 1단계에서 출력했던 디렉토리를 이번에는 입력으로 사용해야한다.

• --readFilesIn : 정렬한 시퀀싱 데이터 파일이다.

• --outFileNamePrefix : 결과 파일이 저장될때, 공통으로 붙을 접두사를 지정한다.

• --outSAMtype BAM Unsorted

  -BAM : 결과를 BAM(표준 압축 바이너리 형식)으로 저장한다. (SAM의 압축버전 )

  • Unsorted : 좌표순 정렬은 하지말고, 찾은 순서대로 저장한다.

• --quantMode GeneCounts : 정렬과 동시에 GTF를 기준으로 유전자 별로 리드가 몇 개인지 카운트하는 명령이다.

4. 결과 파일

Log.final.out

해당 파일은 STAR가 생성한 요약 리포트이다.

• Number of input reads : 총 몇개의 리드를 읽었는지,
• Uniquely mapped reads % : 몇 %가 유전체에 성공적으로 정렬 되었는지,
• Number of splices : 스플라이싱을 몇개나 발견했는지 등을 수치로 내세워 정렬 분석의 품질을 설명한다.

ReadsPerGene.out.tab

각 유전자에 리드가 몇 개씩 카운트 되엇나에 관한 텍스트 파일으로, 다음 통계적 분석 단계의 입력 자료가 된다.

.Bam & .bai

정렬 결과 파일으로, 모든 리드가 어디에 어떻게 정렬되어있는지에 대한 모든 정보가 담겨있다. 해당 파일들로 IGV에서 시각화가 가능하며, 역시 다음 분석 단계의 입력 파일이 된다.

STAR는 Spliced Transcripts Alignment to a Reference의 약자로 스플라이싱이 일어난 RNA 리드를 유전체에 정렬하는 빠른 프로그램이다.

STAR는 인덱싱을 하여 접미사 배열(Suffix Array)을 만든후, 시드 탐색을 통해 최대 매핑 가능 접두사(Maximal Mappable Prefix, MMP) 를 찾는다. 이후 시드 클러스터링 및 스플라이싱을 한 후 해당 정렬본을 평가하고 출력하는 매커니즘을 가지고 있다.

1. 인덱싱

STAR는 정렬을 시작하기 전에, 참조 유전체(Reference Genome, .fasta)를 가지고 접미사 배열(Suffix Array) 이라는 자료구조를 사용하여 인덱스 파일을 만든다.

접미사 배열이란?

유전체 서열에서 가능한 모든 부분 문자열(Suffix)을 가져와서 사전순으로 정렬한 리스트이다. 예를 들어, BANANA라는 문자열이 있다면 A,ANA,BANANA,NA,NANA라는 조각들이 나올텐데, ANA라는 조각을 찾고 싶다면 이 배열에서 A로 시작하는 부분을 보면 되므로 전체에서 ANA를 탐색하는것 보다 빠르다. 시간 복잡도 = $O(M \log N)$

다음과 같은 접미사 배열을 유전체에 적용하여 참조가능한 색인(index) 를 만들고, 이 정렬된 목록에 '위치 값'을 기록해 둔다.

ex : GATTACA sorted : A / ACA / ATTACA / CA / GATTACA / TACA / TTACA position A: 6 / 4 / 1 / 5 / 0 / 3 / 2

STAR는 전체 접미사 배열을 메모리에 올리면 너무 과도한 용량을 요구하기에, genomeParameters.txt에 정의된 특정 길이를 기반으로 압축되고 희소한 Suffix Array를 생성하여 메모리 사용량을 조절한다. ( --genomeSAindexNbits)

2. 시드 탐색

1. MMP ( Maximal Mappable Prefix ) 탐색

STAR는 주어진 리드에서 최대 매핑 가능 접두사(Maximal Mappable Prefix, MMP) 를 찾는것을 최우선으로 한다. MMP 는 리드의 맨 앞에서부터 시작해서, 인덱싱된 참조 유전체와 100% 일치하는 가장 긴 조각 을 의미한다.

예를 들어, 1번~28번 염기(28bp)가 100% 일치하고, 29번 염기에서 mismatch가 발생했다면, MMP1은 1번 염색체 1000 ~ 1027 위치에 매핑된다.

2. 연속 탐색

mismatch가 발생한 29번 염기로부터 리드의 끝을 새로운 리드로 간주하고, 두번째 MMP를 탐색한다.

결론적으로, 한 리드가 여러개의 시드(seed) 조각으로 분해되고, 각 시드가 유전체 어디에 매핑되는지에 대한 위치 목록이 나온다.

• Seed 1 (28bp) : chr1:1000-1027
• Seed 2 (27bp) : chr1:1028-1054
• Seed 3 (40bp) : chr1:50000-50039

다음을 보면, Seed1과 Seed2의 위치는 가까운 반면 Seed3과는 굉장히 멀리 있음을 알 수 있다. 이런 분해된 시드들을 어떻게 연결하고 판단하는지에 대한 과정이 다음 단계인 시드 체이닝 및 스코어링 이다.

3. 시드 체이닝(Chaining) 및 스코어링(Scoring)

동적 프로그래밍을 사용하여 분해된 seed 사이의 최고 점수 경로를 찾는다.

• Seed 1 + Seed 2 연결 :
  • Seed 1의 끝과 Seed 2의 시작이 정확히 이어진다.
  • Gap = 0
  • 패널티 없음
• Seed 2 + Seed 3 연결 :
  • Seed 2의 끝과 Seed 3의 시작 사이에는 거대한 Gap이 있다.
  • STAR는 이 Gap이 이미 알려진 것인지 스플라이스 정션 데이터베이스에서 조회를 한다.
  • Case 1 : 이미 알려진 경우 SJDB에 이미 있는 경우, 해당 구간은 Intron구간이므로 패널티 없음.
  • Case 2 : DB에 없지만, 표준 스플라이싱 신호와 일치 DB에 없지만, 틈의 경계가 GT-AG등 표준 스플라이싱 신호(Canoical motif)와 일치하다면 낮은 패널티.(--alignIntroMin에 기반)
  • Case 3 : DB에도 없고, 경계도 이상함.. 매우 높은 패널티.
  • Case 4 : 틈이 매우 작을때 (스플라이싱이 아님) Indel 패널티(삽입,삭제 패널티)가 적용된다.

4. 최종 정렬

동적 프로그래밍을 통해 최고 점수를 받은 시드 체인은 이 리드의 최종 정렬 후보가 된다. 시드와 시드 사이의 작은 틈이나 리드의 양 끝부분은 로컬 정렬(Local Alignment) 알고리즘을 사용하여 정리한 후, SAM/BAM 포맷으로 출력된다.

Q-Learing은 Reinforcement의 대표적인 가치 기반(Value-based) 알고리즘이다.

• 핵심 : "어떤 상태(state)에서 어떤 행동(Action)을 하는 것이 얼마나 좋은가?" 를 Q값으로 학습하는 것이다.
• 목표 : Agent가 Q-Table를 작성하게 하고, 완성된 Q-Table을 활용하여 어떤 상태에 있든 Q-값이 가장 높은 행동(High Value)을 하게 하는 방식이다.
• 학습 방식 : 탐험과 Q값을 활용한 행동을 반복한다.
  • 에이전트가 상태 S에서 행동 A를 한다.
  • 환경으로부터 보상 R과 다음 상태 S'를 받는다.
  • (S , A, R, S') 을 바탕으로 Q-table을 업데이트 한다.

2. Bellman Equation (벨만 방정식)

벨만 방정식은 Q-Table을 업데이트하는 공식이다. $$ 새로운 Q값 = 현재 Q값 + 학습률 * (목표 Q값 - 현재 Q값) $$ 그리고, 목표 Q값은 다음 과 같이 계산된다. $$ 목표 Q값 = 받은 보상(R) + 할인률 * (다음상태에서 얻을 수 있는 최대 Q값) $$

3. Q-Table

Q-Table 은 Q-Leaning의 가치 저장소 인 행렬이다.

• ROW : 환경의 모든 가능한 상태
• COL : 할 수 있는 모든 행동
• Value($Q(s,a)$) : s상태에서 a행동을 했을 때,미래까지 포함하여 받을 것으로 기대되는 총 보상의 합

다음은 아래에서 설명할 FrozenLake를 학습한 에이전트의 Q-Table이다.

• 첫번째 요소를 보면, 상황0에서 L/D/R/U 행동을 했을 때 얻을수 있는 기댓값이 표시되어있다.
• 에이전트는 다음 값을 참고하여 가장 높은 행동인 Right 방향으로 이동할 것이다.

  4. FrozenLake Game을 활용한 예시

• FrozenLake* 게임은 S에서 시작하여 G에 도착하면 이기는 게임으로, F는 지나갈수있으나, H에 도착하면 지는 게임이다.
• Agent 입장에서, F,G는 숨겨져있다.
• 44 게임 기준으로, 상태는 16가지, 행동은 4가지(Left,Down,Right,Up)가 있으므로 Q-Table의 크기는 164 행렬이다.

학습중에는, Epsilon-Greedy방식을 채택하였다.

Epsilon-Greedy 정책

이 정책은 학습중에는 탐험과 활용을 적절히 섞어야한다는 것이다.

• epsilon : 탐험을 할 확률

  if random.uniform(0,1) < epsilon :
    탐험 행동을 한다 !
  else :
    활용 행동을 한다 !

  학습 초기에는 epsilon을 1.0으로 두엇다가, episode가 진행될수록 줄여나간다.

다시 돌아와서..

아무튼 .. 이런 방식으로 학습을 하고 완성된 164 사이즈의 Q-table을 보고, 우리는 최대요소들만을 골라 최적 경로를 찾을 수 있다. 이를 *최적 정책(Optimal Policy)** 라고 한다. 다음과 같이 해석할 수 있다.

• (0,0) 상태에서, 최적의 행동은 2(right)이다.

5. 코드 구현

import gymnasium as gym
import numpy as np
import random
import time

# 1. Environment Setting
env = gym.make("FrozenLake", is_slippery=False)
# Q-Table init
q_table = np.zeros((env.observation_space.n, env.action_space.n))

# 2. Hyperparameter
num_episodes = 20000
alpha = 0.1     # Leaning Rate
gamma = 0.99    # Discount Factor

# Epsilon - Greedy
epsilon = 1.0
max_epsilon = 1.0
min_epsilon = 0.01
decay_rate = 0.0005     # Epsilon 감소율

# 3. Q-Learning
print("Training Started ...")
for episode in range(num_episodes):
    state, info = env.reset()
    done = False

    while not done:
        if random.uniform(0, 1) < epsilon:
            action = env.action_space.sample()  # 탐험 : 무작위 행동
        else:
            action = np.argmax(q_table[state, :])

        new_state, reward, terminated, truncated, info = env.step(action)

        done = terminated or truncated

        # Q값 업데이트 - 벨만 방정식 사용
        q_table[state, action] = q_table[state, action] + alpha * \
            (reward + gamma *
             np.max(q_table[new_state, :]) - q_table[state, action])

        state = new_state

    epsilon = min_epsilon + (max_epsilon - min_epsilon) * \
        np.exp(-decay_rate * episode)

print("Training Finished")
print("final Q-Table")
print(np.argmax(q_table, axis=1).reshape(4, 4))

env.close()

# 4. 학습된 에이전트 실행 (검증)
print("\n--- Evaluating Trained Agent ---")
env_eval = gym.make("FrozenLake", is_slippery=False, render_mode="rgb_array")
state, info = env_eval.reset()
done = False

while not done:
    # 검증 시에는 Epsilon-Greedy를 쓰지 않고, 오직 Q-Table에서 가장 좋은 행동만 선택
    action = np.argmax(q_table[state, :])

    new_state, reward, terminated, truncated, info = env_eval.step(action)
    done = terminated or truncated
    state = new_state


print("Evaluation finished.")
print(f"Final reward: {reward}")
env_eval.close()

Alignment(정렬) 란, NGS 장비가 읽어낸 수억 개의 짧은 DNA 염기서열 조각(Read)를 참조 유전체(Reference Genome) 와 대조하여 Read가 유전체의 어디에 위치하는지 찾아내는 과정이다. 앞서 Qulity Control 하여 가공된 데이터를 정렬 시키는 과정이 선행되어야만 비로소 유전 변이(Variant), 유전자 발현량(Gene Expression) 등 의미있는 생물학적 분석을 시작할 수 있다.

2. BWA & Samtools

BWA (Burrows-Wheeler Aligner )

BWA 는 Alignment를 수행하는 가장 대표적인 소프트웨어이다. Burrows-Wheeler 는 매우 빠른 문자열 검색을 하는 알고리즘이며, 해당 알고리즘을 사용하여 BWA가 수많은 Read를 거대한 유전체에 신속하게 정렬시킨다.

1. bwa index : Alignment를 시작하기전, 참조 유전체의 'index'를 만드는 명령어이다.
2. bwa mem : 실제 Alignment를 수행하는 명령어로, index를 기반으로 FASTQ파일의 Read들을 참조 유전체에 정렬하고 결과를 SAM형식으로 출력한다.

   Samtools

   BWA가 정렬을 마친 결과물인 sam파일은 거의 사람이 읽을 수 있는 텍스트 형식의 모든 정보를 담은 대용량 파일이다. 이를 컴퓨터가 처리하기 쉽도록 binary 파일로 변환하고 가공하는 툴이다.
3. samtools view : SAM파일을 BAM 파일로 변환한다.
4. samtools sort : BAM 파일 내의 Read들을 좌표 순서대로 재배열한다.
5. samtools index : 정렬된 BAM파일의 index를 만들어 .bai파일을 만든다. 이 파일은 IGV같은 시각화 도구가 특정 유전 영역을 지연 없이 보여주도록 도와준다.

3. Alignment 실습해보기

# 실습용 데이터 다운받기 (PhiX 바이러스 유전체 및 시퀀싱 데이터)
wget -O phix_ref.fna.gz ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/819/615/GCF_000819615.1_ViralProj14015/GCF_000819615.1_ViralProj14015_genomic.fna.gz
wget -O phix_reads_1.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR852/006/SRR8527756/SRR8527756_1.fastq.gz
wget -O phix_reads_2.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR852/006/SRR8527756/SRR8527756_2.fastq.gz
gunzip phix_ref.fna.gz

# BWA로 인덱싱 및 정렬 수행
bwa index phix_ref.fna
bwa mem -t 4 phix_ref.fna phix_reads_1.fastq.gz phix_reads_2.fastq.gz > aligned.sam

# Samtools로 후처리
samtools view -bS aligned.sam > aligned.bam
samtools sort aligned.bam -o aligned.sorted.bam
samtools index aligned.sorted.bam

# 작동 검증 ( 0이 아닌 큰 숫자가 나와야한다. )
samtools view -c -F 4 aligned.sorted.bam

4. IGV로 시각화 하기

IGV는 정렬된 유전체 정보를 시각화 할 수 있는 프로그램이다. Genomes -> Load Genome From File -> phix_ref.fna 를 선택하고, File->Load from File -> aligned.sorted.bam 을 선택해 파일을 시각화할 수 있다.

차세대 염기 서열 분석 (Next-Generation Sequencing, NGS) 기술의 발전으로 방대한 양의 유전체 데이터를 얻게 되었는데, 이 데이터를 저장하는 방식이 FASTQ이다.

@SEQ_ID                (1. 시퀀스 고유 ID)
GATTTGGGGTTCAAAGCAG... (2. 염기 서열 (A, T, C, G))
+                      (3. 구분자)
!''*((((***+))%%%++... (4. 염기 서열의 품질 점수 (Phred Score))

Phred 품질 점수

FASTQ 파일의 네 번째 줄의 기호들은 염기 서열의 각 자리가 얼마나 정확한지에 대한 신뢰도 점수이다. 이 점수는 해당 염기가 잘못 분석되었을 확률(P, error probability)를 바탕으로 계산된다. $$ Q = -10 \log_{10}{(P)} $$

2. QC 과정

Quality Control의 핵심 원리는 기준에서 벗어나는 데이터를 찾아내고 원인을 파악하는 것이다.

1. 품질 지표 확인 (Quality Metrics)

FastQC와 같은 툴을 사용해서 데이터의 다양한 품질 지표를 시각적으로 확인한다.

• 전반적인 염기 서열 품질은 괜찮은가? (Per base sequence quality)
• GC 함량 분포는 정상적인가? (Per sequence GC content)
• 중복된 Read가 비정상적으로 많지는 않은가? (Sequence Duplication Levels)
• 시퀀싱 과정에서 사용된 어댑터(Adapter) 서열이 남아있지는 않은가? (Adapter Content)

2. 문제점 확인

• 품질이 낮은 염기 서열이 있다면 -> Trimming(잘라내기)
• 어댑터 서열 발견 -> Adapter Trimming(어댑터 제거)
• 오염 의심 -> 원인 파악 및 오염 서열 제거

3. 데이터 정제 (Data Cleaning)

Trimmomatic,fastq와 같은 툴을 사용하여 실제 문제점을 해결한다.

4. 정제 후 QC (Post-Cleaning QC)

정제된 데이터를 다시 FastQC를 실행하여 문제가 해결되었는지 확인한다

해당 QC과정을 거쳐 신뢰할 수 있는 데이터를 가지고 다음 단계인 Alignment(정렬) 분석을 실행한다.

3. macOS 터미널에서 FastQC 실습하기

# FastQC 설치
brew install fastqc

# 데이터 다운로드
prefetch SRR1553606
fastq-dump --split-files SRR1553606

# fastQC 실행
fastqc SRR1553606_1.fastq 

SRR1553606은 효모(Saccharomyces cerevisiae)의 RNA - Seq 데이터이다. 해당 데이터를 가지고 fastQC를 실행하면 html파일이 생기는데, 해당 파일을 열어 분석 결과를 볼 수 있다.

4. FastQC Report

우선 각 항목의 지표들이 무슨 의미를 가지는지 알아보자.

Basic Statistics

• Total Sequences :Read의 총 개수(FASTQ 파일의 4줄짜리 한 세트가 Read 1개 이다.)
• Total Bases : 총 염기의 수 (약 3,810만개의 염기가 존재)
• Sequences flagged as poor quality : 품질 저하 서열 수 (0은 품질이 심각하게 나쁜 Read는 없다는 의미로, 좋은 결과이다.)
• Sequence length : Read 하나의 염기 서열 길이 (101로 표시된 것은 모든 Read의 길이가 101bp로 동일하다는 의미이다.)
• %GC : 전체 염기중 구아닌(G)와 사이토신(C)이 차지하는 평균 비율. 이 값은 생물 종의 유전체 특성을 반영한다.

Per base sequence quality ❌

해당 지표는 QC에서 가장 중요한 부분 중 하나이며, 각 위치에서의 Phred 품질 점수의 분포를 박스 플롯으로 설명하고 있다.

결과 분석

• 초반부 ✅ (1 ~ 15bp) : 품질이 약간 흔들리지만, 전반적으로 좋음. 이는 시퀀싱 초반에 나타날 수 있는 자연스러운 현상이다.
• 중반부 ✅ ( 15 ~ 65bp) : 품질이 매우 높고 안정적이다. 대부분의 Read들이 Phread 36 이상인 최상위 품질을 가진다.
• 후반부 ❌ ( 65bp ~ ) : 품질이 급격하게 저하하여 대부분의 염기들을 신뢰할 수 없다.

데이터의 후반부의 품질 저하가 심하므로 Trimming 작업이 필요함을 짐작할 수 있다.

Per sequence quality scores ✅

해당 지표는 개별 Read의 평균 품질 점수 분포를 보여준다. 데이터의 전반적으로 품질이 좋음을 알 수 있다.

Per base sequence content ❌

이 그래프는 Read의 각 자리마다 네 종류의 염기가 얼마나 차지하는지 보여준다. 이상적으로는, 4개의 염기가 모든 위치에서 거의 동일한 비율을 차지하여 4개의 선이 평행하게 그려지는 것이다.

결과 분석

초반부의 심한 편향(1 ~ 15bp) 에 의해 해당 그래프의 초반부의 요동이 심한 것을 볼 수 있다. 이는 시퀀싱 초반 염기 서열이 특정 염기들로 쏠려있다는 염기 구성 편향(sequence bias) 가 있다는 강력한 증거이다. 이러한 현상은 RNA-Seq 실험에서 매우 흔하게 나타나는데, 대표적인 원인은 라이브러리 제작 과정에서 사용하는 Random Hexamers(무작위 6개 염기 조합 프라이머) 때문이다.

Per sequence GC content ⚠️

이 그래프는 Read들의 평균 GC함량 분포를 나타낸다. 빨간 선이 실제 데이터 분포이고, 파란 선은 컴퓨터가 계산한 이상적인 정규 분포 모양이다.

결과 분석

파란 색 선이 빨간색 선과 완전히 일치하지 않고, 모양이 약간 비정상적이기 때문에 ⚠️ 표시가 되었다. 주 봉우리 오른쪽(60-65%)근처에 작은 봉우리가 하나 보이는 것이 문제임을 알 수 있다. 이러한 비정상적인 분포는 오염이거나 생물학적 특성임을 시사한다.

Per base N content ✅

시퀀싱 장비가 특정 위치의 염기를 특정할 수 없을때 그 자리를 N으로 표시하는데, 그러한 N의 비율을 나타낸 그래프이다.

해당 데이터는 N이 완벽하게 없으므로, 어떠한 위치에서도 염기를 판독하지 못한 경우가 없었다라는 의미이다.

Sequence Length Distribution ✅

서열 길이 분포를 나타내며, 해당 데이터의 모든 Read의 길이가 101bp로 동일하다는 좋은 결과를 나타낸다. 대부분의 illumina 시퀀싱 장비는 길이가 고정된 데이터를 생성하기 때문에, 하나의 길이만이 나타나는 그래프가 가장 이상적이다.

Sequence Duplication Levels ❌

이 항목은 ❌ 가 표시되었지만 RNA-Seq 데이터에서는 지극히 정상적인 결과일 수 있기 때문에 해석에 주의가 필요하다. 이 그래프는 완전한 동일한 서열이 데이터에 얼마나 많이 중복되었는지 나타낸다. RNA-Seq는 유전자 발현량을 측정하는 실험인데, 세포 내에는 발현량이 매우 높은 유전자(Housekeeping gene)와 발현량이 낮은 유전자가 공존한다. 따라서 발현량이 매우 높은 유전자에서 유래한 RNA 서열은 당연하게도 데이터에 수만 번 이상 중복되어 나타난다.

Overrepresented sequences ❌

이 표는 데이터에 비정상적으로 많이 나타내는 서열들의 목록이다. 'Possible Source' 항목을 보면 FastQC가 이 서열들의 정체를 RNA PCR Primer또는 Illumina PCR Primer 등으로 정확히 찾아낸 것을 볼 수 있다. 이들은 모두 라이브러리 제작에 사용되는 일종의 어댑터 서열이다.

Adapter Content ❌

이 그래프는 Read의 위치에 따라 Adapter 서열이 얼마나 포함되었는지 나타내는데, 그래프를 보면 Read의 뒷부분(약 55bp 이후)으로 갈수록 어댑터 서열의 비율이 급격하게 증가하는 것을 볼 수 있다. 이는 실제 분석 대상인 RNA 조각의 길이가 Read길이(101bp)보다 짧아 시퀀싱 장비가 RNA 조각을 다 읽고 그 이후의 Adapter 서열을 읽었음을 알 수 있다. 그러므도 해당 서열을 Trimming을 통해 반드시 제거 해야함을 알 수 있다.

정리를 해보자면 ..

첫번째 QC 결과를 따르면, 우리 데이터에는 여러 문제들이 있었다.

1. Adatper Sequences
2. 후반부의 Low Quality bases 다음 문제들을 위해 Trimming을 해보도록 하자.

5. Trimming

# Trimming 을 위한 라이브러리 설치
brew install fastq

fastp \
  -i SRR1553606_1.fastq \
  -I SRR1553606_2.fastq \
  -o trimmed.SRR1553606_1.fastq \
  -O trimmed.SRR1553606_2.fastq \
  -h trimmed.report.html

  # Trimming 후 post - QC 실행
 fastqc trimmed.SRR1553606_1.fastq

• -i,-I : 입력 파일
• -o,-O : 출력 파일
• -h : Trimming 결과에 대한 리포트 생성

다음 과정을 걸치고 나면 trimmed. 접두사가 붙은 fastq 파일과 trimmed.report.html파일이 생김을 알 수 있다.

6. 결과

아까의 Per base sequence quality 그래프와 비교해보면 후반부의 저품질의 데이터가 삭제되어 좋은 품질의 데이터만이 남음을 알 수 있다. 즉, 분석에 방해가 되던 서열들이 제거 되며 데이터의 품질이 극적으로 향상되었음을 알 수 있다.

상태가 바뀐 항목들

• Per base sequence quality (염기 서열 품질): (❌ Fail → ✅ Pass) Read 뒷부분의 품질이 좋지 않은 염기들이 성공적으로 잘려나가 전체적인 품질이 향상 되었다.
• Adapter Content (어댑터 포함): (❌ Fail → ✅ Pass) 분석을 방해하던 어댑터 서열이 제거되었다.
• Overrepresented sequences (과대 발현 서열): (❌ Fail → ✅ Pass) 어댑터 서열이 사라지면서 비정상적으로 많던 서열 문제가 해결되었다.
• Sequence Length Distribution (서열 길이 분포): (✅ Pass → ⚠️ Warn) Trimming 전에는 모든 Read의 길이가 101bp로 일정했지만, 각 Read의 품질에 따라 서로 다른 길이로 잘려 나갔기 때문에 분포가 다양해 졌다. 이는 Trimming의 정상적인 결과이다.
• Per sequence GC content (GC 함량 분포): (⚠️ Warn → ❌ Fail) Trimming의 결과로 전체 Read들의 GC분포가 달라졌다. 데이터가 나빠진 것이 아니라 구성이 바뀐것이므로 신경쓰지 않아도 된다.
• Sequence Duplication Levels (중복 서열 수준): (❌ Fail → ❌ Fail) 이는 RNA-Seq의 유전자 발현량 차이를 의미하므로 Trimming으로 해결되지 않는 문제이다.

이제 Trimming의 결과로 분석을 방해하던 낮은 품질의 데이터, 어댑터문제는 해결되었다. 이 신뢰성 있는 데이터를 가지고 다음 단계인 Alignment(유전체에 서열을 정렬하는 작업) 을 진행하면 된다.

Volcano Plot이란 수많은 변수들 중 통계적으로 유의미한 변화를 보이는 데이터를 직관적으로 식별하기 위해 사용하는 산점도(Scatter plot)의 일종이다. 이름처럼 화산이 폭발하는 듯한 모양을 띄고있다.

주로 두 그룹(ex. 실험군 vs 대조군)간의 유전자 발현량, 단백질 수치 등의 차이를 분석할 때 사용된다. 수천, 수만 개에 달하는 데이터 포인트 중에서 어떤 유전자가 의미 있는 변화를 보였는지 한눈에 알려주는 아주 강력한 도구이다.

• X축 : Log2FoldChange

  • X축에는Log2FoldChange를 사용하는데, 데이터의 분포를 대칭으로 만들어주고, 변화량이 큰 값들의 차이를 시각적으로 잘 표현할 수 있다.
  • 결과적으로 X축의 절댓값이 클수록 두 그룹 간의 변화량 차이가 크다는 것을 의미한다.

• Y축 : -Log10P-value

  • 일반적으로 P-value는 0에 가까운 아주 작은 값을 가지기 때문에, -log10을 취해 값의 분포를 넓혀준다.
  • 결과적으로 Y축의 값이 클수록 두 그룹 간의 차이가 통계적으로 유의미하다는 것을 의미한다.

2. R로 구현해보자

library(DESeq2)
library(airway)
library(ggplot2)
library(ggrepel)
library(dplyr)


# 2. 데이터 준비 및 차등 발현 분석 (Analysis)
data(airway)

# DESeq2 분석을 위한 데이터 객체를 생성
dds <- DESeqDataSet(airway, design = ~ cell + dex)

# 차등 발현 분석을 실행
dds <- DESeq(dds)

# 분석 결과를 'trt'(처리군) vs 'untrt'(대조군) 기준으로 추출
res <- results(dds, contrast = c("dex", "trt", "untrt"))


# 3. 결과 데이터프레임 가공 (Data Wrangling)

res_df <- as.data.frame(res)
res_df <- na.omit(res_df)

res_df$gene <- rownames(res_df)

# 기준: |log2FoldChange| > 1 그리고 padj < 0.05
res_df <- res_df %>%
  mutate(change = case_when(
    log2FoldChange > 1 & padj < 0.05  ~ "Upregulated",
    log2FoldChange < -1 & padj < 0.05 ~ "Downregulated",
    TRUE                             ~ "Not Significant"
  ))


# 4. 최종 시각화 (Visualization)
# 플롯에 라벨링할 상위 10개 유의 유전자를 선택
top_genes <- res_df %>%
  filter(change != "Not Significant") %>%
  arrange(padj) %>%
  head(10)

# ggplot2를 사용하여 Volcano Plot을 생성
ggplot(data = res_df, aes(x = log2FoldChange, y = -log10(padj), color = change)) +
  # 점 그래프 (Scatter plot)
  geom_point(alpha = 0.6, size = 2) +

  # 색상 직접 지정
  scale_color_manual(name = "Expression Change",
                     values = c("Upregulated" = "#E64B35FF", "Downregulated" = "#0073C2FF", "Not Significant" = "grey")) +

  # 유의성 기준선 추가 (점선)
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +

  # ggrepel을 이용해 겹치지 않게 유전자 라벨 추가
  geom_text_repel(data = top_genes,
                  aes(label = gene),
                  box.padding = 0.5,
                  max.overlaps = Inf) +

  # 전체적인 테마 및 제목, 축 이름 설정
  theme_bw(base_size = 14) +
  labs(title = "Differential Gene Expression",
       subtitle = "Dexamethasone Treatment vs. Control in Airway Cells",
       x = expression(paste("Log"[2], " Fold Change")),
       y = expression(paste("-Log"[10], " Adjusted P-value")))

3. 결과 분석

• 빨간점 : 통계적으로 유의미하게 발현이 증가한 유전자들. 이들은 약물 처리의 직접적인 타겟이거나 그로 인해 활성화된 중요한 바이어 마커 후보가 될 수 있다.
• 파란점 : 통계적으로 유의미하게 발현이 감소한 유전자들. 약물에 의해 억제된 유전자들로, 질병과 관련된 경로를 차단하는 역할을 할 수 있다.
• 회색점 : 발현량 변화가 통계적으로 유의미하지 않거나, 변화 자체가 미미한 유전자들이다.

DEG(Differentially Expressed Genes) 란 차등 발현 유전자라고 하며 서로 다른 조건의 샘플 그룹 간에 발현량이 통계적으로 유의미한 차이를 보이는 유전자를 의미한다. 예를 들어, 암세포와 정상 세포 간에 유전자 발현 차이를 비교하여 암세포에서 높은 수준으로 발현되는 유전자들을 찾는 경우가 있다.

DEG를 쓰는 이유와 장점

• 핵심 원인 파악 : 수만 개의 전체 유전자 중에서 특정 조건(질병, 약물 반응 등)과 직접적으로 관련된 핵심 유전자를 소수로 압축하여 연구의 효율성을 높인다.
• 생물학적 기전 이해 : DEG를 분석하면 특정 생물학적 경로(Pathway)니 기능이 두 그룹 간에 어떻게 다른지 파악하여 현상의 메커니즘을 이해하는 단서를 얻을 수 있다.
• 객관적인 바이오마커 획득 : 통계적 기준을 바탕으로 유전자를 선별하므로, 질병 진단이나 치료 반응 예측을 위한 객관적인 바이오마커 후보를 발굴하는 데 매우 유용하다.

2. DEG Analysis Workflow

1. 데이터 전처리 (Data-Pre-processing)

RNA-Seq 장비에서 생성된 원본 데이터(Raw data)의 품질을 확인(Quality Control)하고, 불필요한 서열들을 제거(Trimming)한다.

2. 서열 정렬(Alignment)

전처리된 데이터를 기준 유전체(Reference Genome)에 매핑하여 각 유전자가 얼마나 읽혔는지 확인한다.

3. 유전자 발현 정량화 (Quantification)

각 샘플에서 유전자별로 읽힌 서열의 수(Read Count)를 측정하여 Count Matrix를 생성한다.

4. 정규화 (Normalization)

각 샘플의 총 RNA 양(Sequencing Depth) 차이 등을 보정하여 샘플 간 발현량을 공정하게 비교할 수 있도록 한다.

5. 통계 분석 (Statistical Analysis)

DESeq2,edgeR 과 같은 전문 통계 패키지를 사용하여 그룹 간 발현량 차이의 통계적 유의성(P-value)과 변화량의 크기(Log2 Fold Change) 를 계산한다.

6. 결과 해석 및 시각화

통계적으로 유의미한 DEG목록을 최종적으로 선별하고, Volcano Plot, Heatmap 등으로 시각화하여 결과를 해석한다.

3. R을 활용한 DEG Analysis

우리는 GSE153494데이터를 사용해 보도록 하겠다. 이 데이터는 사람 폐암 세포주(A549)에 특정 약물(Cisplatin)을 처리한 그룹과 대조군을 비교한 RNA-Seq 데이터이다. 항암제 반응과 관련된 유전자 변화를 볼 수 있다. NCBI GEO데 접속하여 GSE1534947_raw_counts_GRCh38.p13_NCBI.tsv.gz파일을 다운로드 받자.

library(DESeq2)

count_matrix <- read.table("GSE183947_raw_counts_GRCh38.p13_NCBI.tsv", header = TRUE, row.names = 1, sep = "\t")

# GEO 사이트에서 확인한 실제 샘플 정보 (12개로 수정)
sample_info <- data.frame(
  condition = c(
    "Cisplatin", "Cisplatin", "Cisplatin", "Cisplatin", # GSM5574715 to 718
    "Control", "Control",                              # GSM5574723, 725
    "Cisplatin", "Cisplatin", "Cisplatin",             # GSM5574728 to 730
    "Control", "Control", "Control"                    # GSM5574733, 739, 740
  )
)
# count_matrix의 열 이름과 순서를 정확히 일치시킴
rownames(sample_info) <- colnames(count_matrix)
meta_data <- sample_info

-
print("--- 최종 Metadata ---")
print(meta_data)
cat("\n메타데이터 행 개수:", nrow(meta_data), "\n")
cat("카운트 매트릭스 열 개수:", ncol(count_matrix), "\n\n")

stopifnot(all(colnames(count_matrix) == rownames(meta_data)))

dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = meta_data,
    design = ~ condition
)

# 분석 실행
dds <- DESeq(dds)

# 결과 추출
res <- results(dds, contrast = c("condition", "Cisplatin", "Control"))
# padj값 오름차순 정렬하기
res_df <- as.data.frame(res)
res_df <- na.omit(res_df)

res_sorted <- res_df[order(res_df$padj),]

# 정렬 결과 확인
cat("--- padj 값으로 정렬된 결과 (상위 6개) ---\n")
print(head(res_sorted))

significant_degs <- subset(res_sorted, padj < 0.05 & abs(log2FoldChange) > 1)

# 유의미한 DEG 필터링
cat("\n--- 최종 선별된 DEG (상위 6개) ---\n")
print(head(significant_degs))
cat("\n총", nrow(significant_degs), "개의 유의미한 DEG를 찾았습니다.\n")

--- padj 값으로 정렬된 결과 (상위 6개) ---
            baseMean log2FoldChange     lfcSE       stat    pvalue      padj
100287102  28.600427     -0.4801389 1.0355784 -0.4636433 0.6429033 0.9999907
653635    129.621323      0.6706098 0.5537453  1.2110438 0.2258786 0.9999907
102466751   7.901483      0.3402459 1.3465738  0.2526753 0.8005192 0.9999907
107985730  14.035188      0.8373989 1.7231867  0.4859595 0.6269959 0.9999907
100302278   1.591026      0.6523731 2.8006700  0.2329347 0.8158121 0.9999907
645520     13.850072     -0.3686650 1.6929921 -0.2177594 0.8276166 0.9999907

--- 최종 선별된 DEG (상위 6개) ---
[1] baseMean       log2FoldChange lfcSE          stat           pvalue        
[6] padj          

• 이 표에서 가장 중요한 두가지 정보는 얼마나 변했는가(log2FoldChange) 와 변화가 통계적으로 유의한가(padj)이다.
  • baseMean : 모든 샘플에 걸친 유전자의 평균 발현량이다.
  • log2FoldChange : Control 그룹 대비 Treatment 그룹에서 발현량이 얼마나 변했는 지를 로그2 스케일로 나타낸다.
    • 양수이면 Treatment 그룹에서 발현량이 증가했다는 의미이다.
  • pvalue : 이 발현량의 차이가 우연히 발생했을 확률이다. 값이 낮을수록 의미 있는 차이 임을 나타낸다.
  • padj : adjusted p-value이며, 수천 개의 유전자를 동시에 검정하면서 발생하는 통계적인 오류를 보정한 값이다. 이 값이 0.05보다 작은지를 기준으로 통계적 유의성을 판단한다.
• 그래서, 일반적으로 다음 두가지 기준을 동시에 만족하는 유전자를 DEG라고 한다.
  1. 통계적 유의성 : padj < 0.05
  2. 생물학적 유의성 : log2FoldChange 의 절댓값 > 1
  • 절댓값이 1보다 크다는 것은 2배 이상 증가했거나 작아졌다는 의미이다.

단백질 3D 구조(Protein 3D Structure) 란 아미노산 서열이 3차원 공간상에서 접혀서(folding) 만들어내는 고유한 입체구조를 의미한다. 단백질은 이. 입체 구조를 형성해야 비로소 생명 활동에 필요한 기능을 수행할 수 있으므로, 그 구조를 이해하는 것은 매우 중요하다.

단백질 구조 분석 도구

• AlphaFold : 구글 딥마인드가 개발한 인공지능으로, 아미노산 서열만으로 단백질의 3차원 구조를 매우 높은 정확도로 예측한다. AlphaFold Protein Structure Database를 통해 누구나 예측된 구조 데이터를 쉽게 얻을 수 있다.
• py3Dmol : 파이썬 환경, 특히 Jupyter Notebook에서 단백질과 같은 분자의 3D 구조를 시각화할 수 있도록 도와주는 라이브러리이다.

구조 분석의 중요성

• 기능 이해 : 단백질의 구조는 그 기능을 결정한다. 특정 분자와 결합하는 부위(활성 부위, active site) 등을 파악할 수 있다.
• 질병 연구 : 많은 질병이 단백질의 구조 이상과 관련이 있으며, 단백질 구조 분석은 질병의 원인을 밝히는데 결정적인 단서를 제공한다.
• 신약 개발 : 질병 관련 단백질의 3D 구조를 타겟으로 하여 약물 분자를 설계하는 데 활용된다.

2. py3Dmol을 활용한 시각화

본 라이브러리는 Jupyter Notebook과 같은 ipynb 환경에서 작동하므로 미리 준비해두자.

데이터 준비하기

https://alphafold.ebi.ac.uk/ 우선 AlphaFold Protein Structure DB 사이트에 접근하여 원하는 데이터를 가져오자. 나는 인간의 p53 단백질을 가져왓고, 단백질 id(Uniprot)는 P04637이다.

저 화면에서 눈여겨 보아야할 것은 Average pLDDT인데, 예측한 구조의 평균 신뢰도를 의미한다.

시각화하기

import py3Dmol

# 뷰어 객체 생성
view = py3Dmol.view(width=800, height=600)

# 파일 열기
with open('AF-P04637-F1-model_v4.pdb','r') as f:
    pdb_data = f.read()
view.addModel(pdb_data,'pdb')

view.setStyle({'cartoon':{'color':'spectrum'}})

view.zoomTo()

view.show()

추가적으로 view.setStyle()의 속성을 변경해가며 다른 기준으로 3D 모델링을 할 수 있다.

AlphaFold 신뢰도 기준

# 1. AlphaFold 신뢰도(pLDDT)에 따라 색칠 (bFactor 값 활용)
    # 파란색(신뢰도 높음) ~ 빨간색(신뢰도 낮음)
    view.setStyle({'cartoon': {'color': 'bFactor'}})

    view.zoomTo()
    view.show()

• 파란색 부분(안정적 코어) : 파란 부분은 예측 신뢰도가 매우 높은 DNA 결합 도메인(DNA-binding domain, DBD) 이다. DNA에 직접 결합하여 유전자 발현을 조절하는 p53의 핵심 기능 부위이다. 대부분의 암 관련 돌연변이가 이 영역에 집중된다.
• 주황색/빨간색 부분(유연한 말단) : 신뢰도가 낮다고 예측 실패가 아닌, 다른 단백질과 상호작용하기 위해 의도적으로 구조가 정해져있지않은 유연한 영역(Flexible Region) 임을 의미한다.

아미노산 종류(Residue Type) 기준

단백질을 구성하는 20가지 아미노산은 각기 다른 화학적 특성을 가지기때문에 , 단백질 표면의 화학적 성질 분포를 시각적으로 이해하는데 도움이 된다.

view.setStyle({'cartoon': {'color': 'spectrum'}})

2차 구조(Secondary Structure) 기준

단백질은 Alpha-helix, Beta-sheet, loop 등과 같은 기본적인 2차 구조를 가지고있다.

view.setStyle({'cartoon': {'color': 'ss'}})

특정 기능 부위 강조

이번에는 전체를 단색으로 칠하고, 암과 관련된 돌연변이가 자주 발생하는 '핫스팟'부위만 다른 색과 스타일로 지정해보겠다.

# 1. 전체를 연한 회색 카툰으로 초기화
view.setStyle({'cartoon': {'color': 'lightgrey'}})

# 2. 암에서 자주 발견되는 돌연변이 '핫스팟' 잔기(175, 248, 273번)를 공 모양(sphere)으로 강조
hotspot_residues = [175, 248, 273]
view.addStyle({'resi': hotspot_residues}, 
              {'sphere': {'color': 'red'}})

view.show()

여러 데이터를 가지고 직접 시각화 하면 좋을거같다 :)

계통수(Phylogenetic Tree) 란 생물학에서 여러 종(또는 개체군, 유전자 등)들이 진화 과정에서 어떻게 갈라져 나왔는지, 그들 사이의 진화적 관계와 유연관계를 Tree 형태로 나타낸 그림이다.

계통수의 구성요소

• Branch : 선으로 표현되며, 각 종이나 그룹 간의 진화적 경로를 나타낸다. branch의 길이는 진화적 시간이나 유전적 변화의 양을 의미하기도 한다.
• Node : 가지가 갈라지는 지점으로, 공통 조상을 의미한다. 해당 지점에서 새로운 종으로 분화가 일어났음을 의미한다.
• Leaf or Terminal : 현재 살아있는 종이나 분석 대상이 되는 특정 그룹을 나타낸다.
• Root : 계통수의 가장 시작점이 되는 부분으로, 모든 생명체의 가장 오래된 공통 조상을 의미한다. 뿌리가 있는 계통수를 유근 계통수(Rooted tree) 라고 한다.

계통수에 쓰이는 데이터

• 형태학적 데이터 : 생물의 겉모습, 해부학적 구조, 화석 기록 등 눈으로 관찰 가능한 특징들을 비교하여 분석한다. 전통적으로 많이 사용된 방식이다.
• 분자생물학적 데이터 : DNA, RNA 서열이나 단백질 아미노산 서열과 같은 유전 정보를 비교 분석한다. 현대 계통분류학에서는 이 데이터가 매우 중요하게 사용된다.

계통수의 중요성

• 생물 분류 : 생물을 진화적 관계에 따라 분류하는 역할을 한다.
• 진화 연구 : 특정 형질이 어떻게 진화했는지, 종 분화가 언제 어떻게 일어났는지 등을 연구하는 데 사용된다.
• 질병 연구 : 바이러스나 세균의 확산 경로를 추적하고 변이 과정을 이해하는데 활용된다.
• 생태학 및 보전 생물학 : 종의 다양성을 파악하고 보전 우선순위를 정하는데 도움을 준다.

2. Neignbor-Joining(NJ) 알고리즘

NJ 알고리즘(Neignbor-Joining) 은. ㅕ러 종 간의 유전적 또는 형태적 거리 를 바탕으로 계통수를 구축하는 거리 기반(Distance-based) 방법이다. 이 알고리즘의 핵심 아이디어는 전체 게통수의 길이를 최소화하는 방향으로 가장 가까운 이웃(Neighbor)을 순차적으로 묶어 나가는 것이다. 단순히 거리가 가장 가까운 쌍을 묶는 것이 아니라, 전체 Tree의 균형을 고려하여 최적의 쌍을 찾는다는 점에서 다른 군집화 방법(ex. UPGMA)과 차이가 있다.

NJ 알고리즘의 단계별 과정

1. 거리 행렬(Distance Matrix) 계산

분석하고자 하는 모든 종들 간의 거리(유전적 차이 등)를 계산하여 행렬 형태로 만든다.

2. Q - 행렬 (Q-Matrix) 계산

단순히 두 종 사이의 거리만 보는 것이 아니라, 다른 모든 종들과의 평균 거리를 보정하여 '진짜 이웃'을 찾기 위한 Q-value를 계산한다. $$ Q(i,j) = (N-2) \times d(i,j) - \sum_{k=1}^{N}{d(i,k)} - \sum_{k=1}^{N}{d(j,k)} $$

• $N$ : 현재 남아있는 종의 수
• $d(i,j)$ : 종 $i$와 종 $j$ 사이의 거리
• $\sum d(i,k)$ : 종 $i$에서 다른 모든 종까지의 거리의 합

3. 가장 가까운 이웃 찾기

계산된 Q-matrix에서 가장 작은 값을 갖는 쌍(pair)을 찾는다. 이 쌍은 현재 단계에서 가장 가까운 '이웃'으로 간주된다.

4. 새로운 마디(Node) 생성 및 가지 길이 계산

3단계에서 찾은 이웃을 묶어 새로운 내부 마디(Internal node)를 생성한다. 그리고 각 이웃에서 새로운 마디까지의 Branch 길이를 계산한다.

5. 업데이트

이제 묶인 두 종은 하나의 새로운 마디로 취급된다. 기존 거리 행렬에서 A와B를 제거하고, 새로운 마디를 추가한다. 남아있는 다른 종들과 새로운 마디 사이의 거리를 다시 계산하여 행렬을 업데이트한다.

6. 반복

종의 개수가 3개가 남을 때 까지 반복을 하며, 마지막 남은 3개의 마디는 하나의 뿌리없는 트리(Unrooted Tree)로 연결되어 계통수가 최종적으로 완성된다.

NJ 알고리즘의 장단점

장점

• 빠른 속도
• 분자 시계 불필요 : 진화 속도가 종마다 달라도 적용할 수 있어, 모든 가지의 길이가 같아야한다는 '분자 시계(molecular clock)'가정을 필요로 하지 않는다.
• 대부분의 경우에서 합리적이고 정확한 계통수를 생성한다.

단점

• 하나의 결과만 생성한다.
• Long-Branch Attraction(LBA) : 진화적으로 멀리 떨어져있지만 진화 속도가 매우 빨라 유전적 변화가 많이 축적된 두 종을 실수로 가깝게 묶는 현상이 발생할 수 있다.
• Unrooted Tree 생성 : 종들 간의 상대적인 관계만을 보여줄뿐, 어디가 가장 오래된 공통 조상(Root)인지는 알 수 없다.

3. Python으로 계통수 만들기

우선 먼저 데이터를 가져와보자. 예제를 위해 나는 가상의 데이터를 가져왔다. 실제 데이터가 아님에 주의바란다. 아래 내용은 사람, 침팬지, 쥐, 닭의 미토콘드리아 유전자(Cytochrome B)일부의 데이터이다.

sequences.fasta 중 
>Human
ATGGCCCCAAATCTCATCATCATCAACAACTTCCTGATCTGCTCCGCCACCCAA
GACACCTACTTCGCCTTCATCATAGCCGCTATCTCCACCGCAACAGGA
>Chimpanzee
ATGGCCCCAAATCTCGTCATCATCAACAACTTCCTGATCTGCTCCGCCACCCAA
GACACCTACTTCGCCTTCATCATAGCCGCTATCTCCACCGCAACAGGA
>Mouse
ATGGCCCCAAACCTAGTCGTCTTCAACAATTTTCTGATCTGCTCCGCCACCCAA
GGCTCCTACTTCGCCTTCATCATCGCCGCCATTTCCACCGCCACCGGC
>Chicken
ATGGCCTCAAACCTAGTCGTCATCAACAATTTCTTGATTTGCTCCGTCATCCAA
GGCTCCCCCCACGCCTTCGTCGTCGCCGCCATCTCCACCACCACCGGG

그리고 다음 sequences.fasta파일을 clustalo를 활용하여 MSA처리 하겠다.

import subprocess

clustalw_exe = "/opt/homebrew/bin/clustalo"
input_file = "sequences.fasta"

cmd = [clustalw_exe, f"-INFILE={input_file}"]

try:
    subprocess.run(cmd, check=True)
    print(f"ClustalW 실행 완료! {input_file}.aln 파일과 {input_file}.dnd 파일이 생성되었을 것입니다.")
except FileNotFoundError:
    print(f"오류: '{clustalw_exe}' 파일을 찾을 수 없습니다. 경로를 다시 확인해주세요.")
except subprocess.CalledProcessError as e:
    print(f"ClustalW 실행 중 오류가 발생했습니다: {e}")

clustal.aln 중...

CLUSTAL O(1.2.4) multiple sequence alignment


Human           ATGGCCCCAAATCTCATCATCATCAACAACTTCCTGATCTGCTCCGCCACCCAAGACACC
Chimpanzee      ATGGCCCCAAATCTCGTCATCATCAACAACTTCCTGATCTGCTCCGCCACCCAAGACACC
Mouse           ATGGCCCCAAACCTAGTCGTCTTCAACAATTTTCTGATCTGCTCCGCCACCCAAGGCTCC
Chicken         ATGGCCTCAAACCTAGTCGTCATCAACAATTTCTTGATTTGCTCCGTCATCCAAGGCTCC
                ****** **** **  ** ** ******* **  **** ******* ** ***** * **

Human           TACTTCGCCTTCATCATAGCCGCTATCTCCACCGCAACAGGA
Chimpanzee      TACTTCGCCTTCATCATAGCCGCTATCTCCACCGCAACAGGA
Mouse           TACTTCGCCTTCATCATCGCCGCCATTTCCACCGCCACCGGC
Chicken         CCCCACGCCTTCGTCGTCGCCGCCATCTCCACCACCACCGGG
                  *  ******* ** * ***** ** ****** * ** ** 

그리고, 이 정렬된 데이터를 가지고 계통수를 그려보자!

import matplotlib.pyplot as plt
from Bio import AlignIO, Phylo
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

alignment_file = "clustal.aln" 

print(f"'{alignment_file}' 파일을 읽어 계통수를 생성합니다...")

# 1. 정렬 파일 읽기
try:
    aln = AlignIO.read(alignment_file, 'clustal')
except FileNotFoundError:
    print(f"오류: '{alignment_file}' 파일을 찾을 수 없습니다. 파일 이름을 확인해주세요.")
    exit()

# 2. 거리 행렬(Distance Matrix) 계산
# 'identity' 모델은 서열이 얼마나 다른지를 기준으로 거리를 계산합니다.
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(aln)

# 3. 계통수 구축 (Neighbor-Joining 알고리즘 사용)
constructor = DistanceTreeConstructor(calculator)
tree = constructor.nj(dm) # nj() 메서드가 Neighbor-Joining을 수행합니다.

# 뿌리(root)를 가장 긴 가지의 중간 지점으로 설정하여 보기 좋게 만듭니다.
tree.root_at_midpoint()

# 4. 계통수 시각화
print("계통수 생성이 완료되었습니다. 이제 그림으로 표시합니다.")

# Matplotlib
fig = plt.figure(figsize=(10, 6), dpi=100)
axes = fig.add_subplot(1, 1, 1)
Phylo.draw(tree, axes=axes, do_show=False) # do_show=False로 설정 후

# 보기 좋게 레이아웃 조정 및 제목 추가
plt.title("Phylogenetic Tree (Neighbor-Joining)")
plt.xlabel("Evolutionary Distance")
plt.ylabel("Taxa")
plt.show()

4. 데이터 분석하기

사람과 침팬지가 inner1이라는 가장 가까운 공통 조상에서부터 갈라져 나옴을 알 수 있으며, 두 종으로 뻗어 나가는 branch 길이가 매우 짧은 것을 볼 수 있는데, 이는 두 종의 유전적 차이가 아주 적다는 것을 의미한다.

그 다음으론 쥐, 닭 순서대로 유전적으로 가까움을 알 수 있다.

다중 서열 정렬(MSA) 는 셋 이상의 DNA 또는 단백질 서열을 가져와, 진화적으로 동일한 위치(상동 부위)를 찾아 나라히 정렬하는 분석 기법이다.

다중 서열 정렬의 주요 목적

• 진화적 관계 추정 : 여러 종에 걸쳐 보존된 서열 영역을 분석하여 계통 발생 관계(Phylogenetic relationship)을 추정하고 공통 조상을 찾는 데 사용된다.
• 기능적 부위 예측 : 여러 서열에서 공통으로 보존되는 단백질은 중요한 기능을 담당할 가능성이 높다.
• *단백질 3차 구조 예측 : 상동성 모델링(homology modeling)에서 알려진 구조의 단백질 서열과 구조를 모르는 단백질 서열을 정렬하여 구조를 예측하는 데 중요한 기초 자료로 사용된다.
• PCR 프라이머 설계 : 다양한 종의 특정 유전자를 증폭시키기 위해 공통적으로 보존된 서열 영역에 결합할 수 있는 축퇴 프라이머(Degenerate primer)를 설계할 수 있다.

2. 다중 서열 정렬 알고리즘

정확한 최적의 MSA를 찾는 것은 NP-Hard 문제 이기 때문에, 주로 경험적이고 근사적인 방법을 사용한다.

점진적 정렬(Progressive Alignment)

가장 널리 사용되는 방식으로, 모든 서열 쌍의 유사도를 계산하여 가이드 트리(Guide tree)를 만들고, 이 트리에 따라 가장 가까운 서열부터 차례대로 정렬을 추가해 나간다. 초기 정렬의 오류가 최종 결과에 계속 영향을 끼친다는 단점이 있다.

• 대표 도구 : ClustalW / Clustal Omega, T-Coffee

반복적 정렬 (iterative Alignment)

초기 정렬을 만든 후, 이를 반복적으로 수정하고 재정렬하여 점진적 방법의 단점을 보완하고 더 나은 정렬 결과를 찾아가는 방식이다.

• 대표 도구 : MUSCLE, MAFFT

3. Python으로 TP53 단백질 정렬하기

최신 트렌트 : MAFFT와 대화형 시각화

• MAFFT : MUSCLE이나 ClustalW보다 더 빠른 속도와 높은 정확도를 보여주는 MSA 프로그램으로, 현재 연구실에서 매우 선호된다.
• 대화형 시각화 : 정적인 이미지 대신, 마우스로 직접 탐색하고 확대/축소하며 서열을 확인할 수 있는 대화형 시각화 라이브러리(plotly)를 사용한다.

3-1. 분석 데이터 준비

우선 사용 데이터는 인간, 고릴라, 오랑우탄의 TP53 단백질 서열을 사용할 것이며, 이를 NCBI에서 FASTA 형식으로 다운로드 하겠다.

from Bio import Entrez

Entrez.email = "bublman3375@knu.ac.kr"

# 가져올 단백질 서열의 NCBI id들
# 인간 , 고릴라, 오랑우탄
ids = ["NP_000537.3","XP_063556659.1","XP_002827020.3"]

output_filename = "tp53_primates.fasta"

print(f"NCBI에서 다음 ID의 서열을 다운로드합니다:\n{ids}")

try :
    handle = Entrez.efetch(db="protein",
                           id=",".join(ids),
                           rettype="fasta",
                           retmode="text")
    fasta_data = handle.read()
    handle.close()

    with open(output_filename,"w") as f :
        f.write(fasta_data)

    print(f"\n성공. '{output_filename}' 파일이 생성되었습니다.")

except Exception as e :
    print(f"오류발생 : {e}")

NCBI에서 다음 ID의 서열을 다운로드합니다:
['NP_000537.3', 'XP_063556659.1', 'XP_002827020.3']

성공. 'tp53_primates.fasta' 파일이 생성되었습니다.

그리고 tp53_primates.fasta파일을 열어보면 데이터가 잘 받아와진 것을 볼 수 있다.

3-2. MAFFT와 Plotly를 사용하여 MSA 하기

우선적으로 plotly는 pip install plotly로 터미널에서 다운로드 받을수 있으며, MAFFT는 사이트에 접속하여 직접 다운로드 하여야한다. https://mafft.cbrc.jp/alignment/software/

이제 필요한 라이브러리와 프로그램, 데이터가 준비되었으니 MSA를 해보도록하자.

from Bio.Align.Applications import MafftCommandline
from Bio import AlignIO
import plotly.express as px
import pandas as pd
import os

# --- 설정 부분 ---
input_file = "tp53_primates.fasta"
output_file = "tp53_aligned.fasta"
mafft_exe = "mafft"
# -----------------

try:
    if not os.path.exists(input_file):
        raise FileNotFoundError(f"입력 파일 '{input_file}'을 찾을 수 없습니다.")

    print(f"MSA - MAFFT를 사용하여 정렬을 시작합니다...")

    # --- 1. MAFFT 실행 ---
    mafft_cline = MafftCommandline(mafft_exe, input=input_file)
    stdout, stderr = mafft_cline()

    with open(output_file, "w") as f:
        f.write(stdout)
    print(f"정렬 완료. 결과가 '{output_file}' 파일에 저장되었습니다.")

    # --- 2. 정렬 결과 읽기 ---
    alignment = AlignIO.read(output_file, "fasta")

    aligned_sequences = [list(record.seq) for record in alignment]
    sequence_ids = [record.id for record in alignment]
    df = pd.DataFrame(aligned_sequences, index=sequence_ids)

    # --- 3. [수정] 시각화를 위한 데이터 전처리 ---
    # 각 아미노산 문자('M', 'E', 'P'...)를 고유한 숫자(0, 1, 2...)로 변환
    print("\n시각화를 위해 아미노산 문자를 숫자로 변환합니다...")

    # 모든 아미노산 종류를 찾아서 숫자와 매핑하는 딕셔너리 생성
    amino_acids = sorted(list(set(df.values.flatten())))
    aa_to_int = {aa: i for i, aa in enumerate(amino_acids)}

    # 데이터프레임의 모든 문자를 숫자로 변환
    df_numeric = df.replace(aa_to_int)

    # --- 4. Plotly를 활용한 대화형 시각화 ---
    print("Plotly를 활용하여 시각화를 시작합니다...")
    fig = px.imshow(df_numeric, # !! 숫자 데이터프레임을 사용 !!
                    labels=dict(x="아미노산 위치", y="종(Species)", color="아미노산 ID"),
                    title="<b>TP53 단백질 다중 서열 정렬(MSA) 시각화</b>",
                    color_continuous_scale="Viridis")

    # 마우스를 올렸을 때 원래 아미노산 문자가 보이도록 설정
    fig.update_traces(
        customdata=df.values, # 원본 문자 데이터를 customdata로 전달
        hovertemplate="<b>위치</b>: %{x}<br><b>종</b>: %{y}<br><b>아미노산</b>: %{customdata}<extra></extra>"
    )

    fig.show()

except FileNotFoundError as e:
    if "mafft" in str(e):
        print("[에러] MAFFT 프로그램을 찾을 수 없습니다.")
    else:
        print(f"[파일 없음 에러] {e}")
except Exception as e:
    print(f"알 수 없는 에러가 발생했습니다: {e}")

4. 결과 분석

결과적으로 다음 그래프가 나오게 된다.

• X축 : 단백질을 구성하는 393개의 아미노산의 위치를 나타낸다.
• Y축 : 3종의 TP53 단백질 서열을 나타낸다.

그래프의 모든 열이 위아래로 완전히 동일한 색상 패턴을 보이는것은 "인간, 고릴라, 오랑우탄의 TP53 단백질은 아미노산 서열이 100% 동일하다" 라고 결론 지을 수 있다. 수백만 년에 걸쳐 인간, 고릴라, 오랑우탄이 공동 조상으로부터 현재까지의 진화 과정속에서 이 단백질의 서열에는 어떤 변화도 허용되지 않았음을 의미한다. 또, 해당 유전자에 돌연변이가 발생했다면, 그 개체는 암 발생 억제 능력이 떨어져 생존에 불리하여 도태되었다는 음성 선택(Negative Selection) 가설을 강력하게 뒷받침 한다.

TP53 유전자는 DNA 손상을 복구하거나, 복구가 불가능할 경우 세포 사멸을 유도하여 암을 억제하는 '게놈의 수호자'라고 불리는 중요한 유전자 이다. 이 유전자에 돌연변이가 발생하면 DNA 손상을 입은 세포가 계속 증식하여 다양한 종류의 암 이 발생할 수 있으며 , 리-프라우메니 증후군(Li-fraumeni like syndrome)과 같은 유전성 암 발생 위험을 높이기도 한다.

특히 폐암, 유방암, 대장암 등을 포함한 많은 종류의 암에서 50 ~ 60% 정도로 높은 빈도의 변이를 보이는데, 이는 인체 종양에서 발견되는 발암 또는 종양 억제 유전자들의 돌연변이 빈도중에서도 가장 높다.

1-1. TP53(p53 단백질)의 주요 기능

• 세포 주기 정지 : DNA가 손상된 것을 감지하면, p53 단백질은 세포 분열 주기를 일시적으로 멈춘다.
• DNA 복구 : p53 단백질은 손상된 DNA를 수리하는 다른 유전자들을 활성화시켜 DNA를 원상태로 복구하도록 돕는다.
• 세포 사멸 유도(Apoptosis) : 만약 DNA 손상이 너무 심각해 복구가 불가능하다면, p53 단백질은 세포 자살 프로그램 을 가동한다.

1-2. 임상적 중요성

• 진단 및 예후 예측 : 암 조직에서의 TP53 돌연변이의 유무와 종류를 분석하여 암의 악성도를 추적할 수 있다. 일반적으로 TP53 돌연변이가 있는 경우 암은 공격적인 성향을 보인다.
• 치료 반응 예측 : 특정 항암 치료에 대한 반응성을 예측하는 지표로도 사용된다. 일부 항암제는 p53의 정상적인 기능에 의존하여 암세포를 사멸시키므로, TP53에 돌연변이가 있는 경우 치료 효과가 떨어질 수 있다.

1-3. 생화학 경로

정상 세포에서 p53은 음성조절자 mdm2 에 의해 불활성화 된다. 그러나 DNA손상이나 다른 스트레스가 있을시 다양한 경로들이 p53과 mdm2를 분리시켜 활성화 된다. *p53은 세포주기를 멈추게 하고 손상을 고치거나, 세포를 없앨 것을 결정하는데 어떤 원리로 p53이 결정을 내리는지는 아직 알 수 없다. * 이는 현재 분자생물학에서 활발히 연구되고 있는 중요한 질문 중 하나이다.

2. FASTA 형식과 BLAST

2-1. FASTA

FASTA 형식은 유전자를 구성하는 염기서열이나 단백질을 구성하는 아미노산 서열을 텍스트로 저장하는 표준 방식이다.

• Description Line

  • 반드시 '>'로 시작한다.
  • '>' 뒤에는 해당 서열의 이름, 아이디, 출처 등을 식별할 수 있는 정보를 자유롭게 기입한다.

• Sequence Data

  • 설명 줄 다음 줄부터 실제 염기 또는 아미노산 서열을 입력한다. 서열은 한 줄로 길게 쓰고나 여러 줄로 나누어 쓸 수 있으며 줄바꿈, 공백, 숫자 등은 프로그램이 알아서 무시한다.

    >gi|568815597|ref|NM_000546.5| Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA
    AGCTTCGACAGCAGCTC ..............

2-2. BLAST (Basic Local Alignment Search Tool)

BLAS는 특정 염기 서열이나 아미노산 서열이 어떤 생물체의 유전자/단백질과 유사한지를 빠르게 검색해준다. 전 세계의 방대한 데이터베이스(ex. NCBI의 GenBank)를 대상으로 내가 가진 서열(Query Sequence)과 비슷한 부분을 찾아준다.

• BLAST 역할

  • 유전자 / 단백질 동정 : 미지의 서열이 어떤 기능을 하는 유전자 혹은 단백질인지 정체를 밝힌다.
  • 종간 유사성 비교 : 비슷한 서열을 가진 다른 종을 찾아내어 진화적인 관계(유연관계)를 추정할 수 있다.

• BLAST 원리와 결과 BLAST는 전체 서열을 하나하나 비교하는 대신, 서열을 짧은 조각(word)으로 나눈 뒤 데이터 베이스에서 이 조각과 일치하는 부분을 먼저 찾아낸다. 그리고 그 주변으로 정렬(alignment)을 확장해 나가며 국소적으로 유사도가 높은 영역(Local Alignment)을 찾아내는 방식이다. 검색이 종료되면 유사도가 높은 순서대로 결과 목록을 보여주는데, 이때 쓰이는 지표가 E-value(Expectation Value) 이다.

• E-value : 현재 검색 결과와 같은 유사도(Score)가, 우연히 발견될 수 있는 기대 횟수를 의미한다. 0에 가까울 수록 우연의 일치일 가능성이 희박하며, 두 서열이 통계적으로 유의미한 관계에 있음 을 시사한다.

  3. BioPython을 활용한 BLAST

  3-1. BioPython

  BioPython 은 생물학 데이터를 위한 거대한 패키지이며, 안에는 다양한 목적을 가진 여러 라이브러리가 있다. 이번에 쓰일 라이브러리인 Entrez와 Blast에 대해서 알아보도록 하자.

  3-2. Bio.Entrez

  Entrez는 미국 국립생물정보센터(NCBI)가 운영하는 방대한 생물학 데이터베이스에 접속하게 해주는 인터페이스 이다.

• Entrez.esearch() : 특정 키워드로 지정된 데이터베이스를 검색해서, 조건에 맞는 데이터들의 고유 ID 목록을 찾아준다. 이때, 데이터를 가져오는것이 아님을 주의하자.

• Entrez.efetch() : 고유 ID를 입력하면 해당 ID의 실제 데이터를 FASTA, GenBank등 원하는 방식으로 가져온다.

3-3. Bio.Blast

Blast는 내가 가진 특정 서열(Query)을 거대한 서열 데이터베이스와 비교하여 유사서열을 찾아주는 프로그램이다.

• NCBIWWW.qblast() : NCBI의 공식 BLAST 서버를 사용하여 검색한다.

3-4. Python에서 해보기

from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import Entrez

# NCBI 데이터베이스에 접속하기 위해 자신의 이메일을 명시합니다. (NCBI 정책)
Entrez.email = "bublman375@knu.ac.kr"

# 1. 데이터 가져오기 (Entrez)
# protein DB에서 인간 TP53(ID: NP_000537.3)의 서열을 FASTA 형식으로 다운로드
handle = Entrez.efetch(db="protein", id="NP_000537.3", rettype="fasta", retmode="text")
fasta_sequence = handle.read()
handle.close()

print("--- 가져온 서열 ---")
print(fasta_sequence)

# 2. BLAST 실행
print("TP53 - NCBI BLAST Starting...")
# blastp (단백질 vs 단백질)로 nr DB에 검색, 상위 200개 결과를 요청
result_handle = NCBIWWW.qblast("blastp", "nr", fasta_sequence, hitlist_size=200)
print("TP53 - NCBI BLAST Ended")

# 3. 결과 파싱 및 필터링
# XML 형식의 BLAST 결과를 파이썬 객체로 변환
blast_record = NCBIXML.read(result_handle)

print("\n--- BLAST 검색 결과 (인간 & 인공 제작물 제외) ---")
count = 0
for alignment in blast_record.alignments:
    # 검색 결과에서 인간(Homo sapiens)과 인공 서열(synthetic construct)은 제외
    if "Homo sapiens" not in alignment.title and "synthetic construct" not in alignment.title:
        for hsp in alignment.hsps:
            count += 1
            print(f"순위 {count}: {alignment.title}")
            print(f"     E-value: {hsp.expect}\n")

    # 상위 10개 결과만 출력 후 종료
    if count >= 10:
        break

result_handle.close()

• Entrez.efetch(db="protein", id="NP_000537.3", rettype="fasta", retmode="text") : protein 데이터베이스에서 고유 ID가 NP_000537.3인 데이터를 요청한다. 이는 인간(Homo sapiens)의 p53이다.
• NCBIWWW.qblast("blastp","nr",fasta_sequence,hitlist_size=200): 일반적으로 qblast는 상위 50개 결과만을 우선적으로 가져오지만hitlist_size=200을 통해 상위 200개 결과를 검색하였다.
• NCBIXML.read를 이용해 BLAST 서버로 받은 XML형식의 데이터를 파이썬이 다루기 쉬운 객체로 변환한다.
• if "Homo sapiens" not in alignment.title and "synthetic construct" not in alignment.title : 인간(Homo sapiens)와 인공제작물(Synthetic construct)를 제외하여 검색한다.

--- 가져온 서열 ---
>NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAA
PPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKT
CPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRN
TFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGR
DRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALEL
KDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD


TP53 - NCBI BLAST Starting...
TP53 - NCBI BLAST Ended

--- BLAST 검색 결과 (인간 & 인공 제작물 제외) --

순위 1: ref|XP_063556659.1| cellular tumor antigen p53 isoform X1 [Gorilla gorilla gorilla]
     E-value: 0.0

순위 2: ref|XP_004058559.3| cellular tumor antigen p53 isoform X2 [Gorilla gorilla gorilla]
     E-value: 0.0

순위 3: ref|XP_002827020.3| cellular tumor antigen p53 [Pongo abelii]
     E-value: 0.0

순위 4: ref|XP_030656345.1| cellular tumor antigen p53 [Nomascus leucogenys]
     E-value: 0.0

순위 5: ref|XP_054313311.2| cellular tumor antigen p53 [Pongo pygmaeus] >gb|PNJ12854.1| TP53 isoform 1 [Pongo abelii] >gb|PNJ12856.1| TP53 isoform 3 [Pongo abelii]
     E-value: 0.0

순위 6: ref|XP_031995228.1| cellular tumor antigen p53 [Hylobates moloch]
     E-value: 0.0

순위 7: ref|XP_055112571.1| cellular tumor antigen p53 isoform X1 [Symphalangus syndactylus]
     E-value: 0.0

순위 8: ref|XP_017706033.1| PREDICTED: cellular tumor antigen p53 isoform X1 [Rhinopithecus bieti]
     E-value: 0.0

순위 9: ref|XP_050621303.1| cellular tumor antigen p53 isoform X1 [Macaca thibetana thibetana]
     E-value: 0.0

순위 10: ref|XP_010360689.1| cellular tumor antigen p53 [Rhinopithecus roxellana] >ref|XP_017706034.1| PREDICTED: cellular tumor antigen p53 isoform X2 [Rhinopithecus bieti]
     E-value: 0.0

4. 결론

결과적으로 이 작성한 코드는 인간의 p53 단백질과 진화적으로 가까운 '사촌' 단백질을 다른 종에서 찾아내라 는 명령을 수행하는 파이프라인이다. 다음 결과에서 1,2위는 Gorilla, 3위는 Orangutan(Pongo abelii) ... 임을 알 수 있다. E-value가 0.0으로 표시되어 있지만 실제로는 0에 무수히 가까운 값이기 때문에 영장류 들이 인간의 유전적 조상임을 알 수 있다.

결정 트리(Decision Tree) 는 '의사결정 규칙'을 'Tree' 자료구조로 표현 한 모델이다. 모델은 데이터에 대해 연속적인 질문을 던져가며 가중치를 보정하고, 끝내 데이터가 어떤 그룹에 속하는지 분류 하거나 어떤 값인지(회귀) 에 대한 결론에 도달한다.

# 필요한 라이브러리 불러오기
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier, plot_tree
import matplotlib.pyplot as plt

# 1. 데이터 준비
iris = load_iris()
X = iris.data
y = iris.target

# 2. 모델 생성 및 학습
dt_clf = DecisionTreeClassifier(max_depth=3, random_state=42)
dt_clf.fit(X, y)

# 3. 학습된 트리 시각화
plt.figure(figsize=(15, 10))
plot_tree(dt_clf, 
          filled=True, # 노드의 클래스에 맞게 색을 채웁니다.
          feature_names=iris.feature_names, # 각 노드가 어떤 특성으로 분기했는지 표시합니다.
          class_names=iris.target_names) # 잎 노드의 클래스 이름을 표시합니다.
plt.title("Decision Tree Visualization")
plt.show()

• 위에 트리를 확인해보면, 모델이 어떤 기준을 통해 붗꽃의 품종을 찾아가는지 한눈에 확인할 수 있다.

2. 과적합과 가지치기

결정 트리의 가장 큰 단점은 훈련 데이터에 과하게 최적화되는 과적합(Overfitting) 이 다. 모델이 훈련 데이터의 사소한 노이즈까지 모두 학습하게 되어 새로운 데이터에 대한 예측 성능이 떨어지게 된다. 이를 해결하기 위해 가지치기(Pruning) 개념이 등장했다. 말 그대로 너무 복잡한 트리의 가지를 쳐내어 모델을 단순화 하는 과정이다. 위 코드의 DecisionTreeClassifier의 파라미터들을 조절하여 가지치기를 할 수 있다.

• max_depth : 트리의 최대 깊이를 제한한다.
• min_samples_split : 노드를 나누기 위해 필요한 최소 데이터 수를 지정한다.
• min_samples_leaf : 마지막 잎 노드가 되기 위해 필요한 최소 데이터 수를 지정한다.

3. 랜덤 포레스트 (Random Forest)

랜덤 포레스트(Random Forest)는 결정 트리의 과적합 문제를 해결하기 위해, 조금씩 다른 여러 개의 결정 트리를 만들어 그 결과를 종합하는 앙상블(Esemble)기법이다.

두가지 핵심적인 '무작위성' 을 통해 다양성을 확보한다.

1. 데이터의 무작위성 (Bagging) : 훈련 데이터에서 무작위로 데이터를 샘플링하여 각. 트리를 학습시킨다.
2. 특성의 무작위성 : 각 트리의 노드를 분할할 때, 전체 특성 중 일부만 무작위로 뽑아 사용한다.

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# 1. 데이터 준비 (훈련/테스트 데이터 분리)
X_train, X_test, y_train, y_test = train_test_split(
    iris.data, iris.target, test_size=0.2, random_state=42
)

# 2. 모델 생성 및 학습
# n_estimators=100 : 숲을 구성할 나무의 개수를 100개로 지정합니다.
rf_clf = RandomForestClassifier(n_estimators=100, random_state=42)
rf_clf.fit(X_train, y_train)

# 3. 예측 및 성능 평가
y_pred = rf_clf.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)

print(f"랜덤 포레스트 모델의 정확도: {accuracy:.4f}")

랜덤 포레스트 모델의 정확도: 1.0000

다음은 100개의 트리중 첫번째 트리만을 시각화 한 결과이다.

다음은 특성 중요도(Feature Importance)를 시각화 한 결과이다.

경사 하강법의 목표는 비용 함수(Cost Function)의 값을 최소로 만드는 파라미터를 찾는 것이다.

1-1. 준비

먼저 최적화 할 대상인 직선의 방정식을 정의하겠다. $$ H(x) = Wx +b $$ 그리고, 비용 함수는 MSE를 사용하여 정의하겠다 . $$ J(W,b) = \frac{1}{m} \sum_{i=1}^{m}(H(x_i)-y_i)^2 $$ 이때 우리가 찾아야할 파라미터는 $W$와$b$이다.

1-2. Gradient

비용함수 $J(W,b)$를 최소화하려면, 현재 파라미터 위치에서 비용이 가장 가파르게 감소하는 방향으로 $W$와$b$를 이동시켜야한다. 수학적으로 가장 가파른 방향은 기울기(Gradient)로 알 수 있으며, 이는 비용 함수를 각 파라미터로 편미분하여 구한다.

• $W$에 대한 편미분 ($W$의 영향력) $$ \frac{∂J}{∂W} = \frac{2}{m}\sum_{i=1}^{m}(H(x_i)-y_i)\times x_i $$

• $b$에 대한 편미분 ($b$의 영향력) $$ \frac{∂J}{∂b} = \frac{2}{m}\sum_{i=1}^{m}(H(x_i)-y_i) $$

1-3. Update

기울기(Gradient)를 구했다면, 이제 기울기의 반대방향으로 파라미터를 업데이트 한다. 이때, 얼마나 큰 간격으로 움직일지를 결정하는 것이 학습률(Learning Rate, $\alpha$) 이다. $$ W : = W-\alpha\frac{∂J}{∂W} \ $$ $$ b := b-\alpha\frac{∂J}{∂b} $$

경사 하강법은 이 과정을 수렴할 때까지 반복하는 것이다.

2. Numpy를 이용한 구현

import numpy as np
import matplotlib.pyplot as plt

# 1. 데이터 준비
num_point = 100

# 모델이 찾아내야 할 실제 정답
true_W = 3.5
true_b = 8.0

# X와 Y 데이터를 랜덤으로 생성
noise = np.random.randn(num_point,1)
X = 2 * np.random.rand(num_point,1)
Y = true_W * X + true_b + noise

# 2. 경사하강법 구현
learning_rate = 0.01 # 학습률
epochs = 1000 # 학습 반복 횟수

# 파라미터를 0으로 초기화
W = 0.0
b = 0.0

m = len(X)
cost_history = [] # 비용의 변화를 기록하기 위한 리스트

for epoch in range(epochs) :
    # 가설 계산
    predictions = W * X + b

    # 비용 계산
    cost = (1/m) * np.sum((predictions - Y) ** 2)
    cost_history.append(cost)

    # 기울기 계산
    W_gradient = (2/m) * np.sum((predictions - Y) * X)
    b_gradient = (2/m) * np.sum(predictions - Y)

    # 업데이트
    W = W - learning_rate * W_gradient
    b = b - learning_rate * b_gradient

# --- 3. 결과 출력 및 시각화 ---
print("--- 학습 완료 ---")
print(f"최종 예측 W: {W:.3f} (실제 값: {true_W})")
print(f"최종 예측 b: {b:.3f} (실제 값: {true_b})")

--- 학습 완료 ---
최종 예측 W: 3.629 (실제 값: 3.5)
최종 예측 b: 7.780 (실제 값: 8.0)

오른쪽 그래프는 학습이 진행됨에 따라 오차가 빠르게 감소하여 특정 값으로 수렴하는 전형적인 경사 하강법의 학습 곡선을 보여준다.

선형 회귀(Linear Regression)이란 머신러닝 중 연속적인 숫자 값을 예측하고 싶을 때 사용하며 주어진 데이터를 가장 잘 설명하는 직선 하나를 찾는 것 이다.

• 독립 변수 (X) : 원인이 되는 데이터 ( ex. 공부 시간 )
• 종속 변수 (Y) : 결과가 되는 데이터 ( ex. 시험 성적 ) 선형 회귀는 이 X가 변할때 Y가 어떻게 변하는지 예측하는 모델을 만드는 과정이다.

2. 핵심 개념

2-1 . 가설 (Hypothesis)

우리가 찾으려는 직선은 다음 방정식으로 나타낼 수 있다. $$ H(x) = Wx + b $$ 이때 $W$는 가중치(Weight)이며 직선의 기울기 이며, $b$는 편향(bias), 즉 y절편이다. 머신 러닝의 목표는 주어진 데이터에 대해 가장 적절한 $W$와 $b$를 찾는 것이다.

2-2. 비용 함수 (Cost Function / Loss Function)

우리가 찾은 직선이 실제 데이터와 얼마나 차이가 나는지 오차를 계산하는 함수이다. 가장 대표적인 비용 함수는 평균 제곱 오차(Mean Squared Error, MSE)이다. $$ Cost = (V_{predict} - V_{real}) ^2의 \ 평균 $$

• 제곱을 하는 이유는 오차의 크기를 극대화하고, 음수 값을 없대기 위함이다.
• 우리의 목표는 이 $Cost$를 최소로 만드는 $W$와$b$를 찾는 것이다.

  해당 내용은 손실 함수 - L1 loss & L2 loss에서 잘 나와있습니다 :)

  2-3. 경사 하강법 (Gradient Descent)

  $Cost$를 최소화하는 $W$와$b$를 어떻게 찾을지에 대한 방법론 중 하나이다.
• 함수의 변화량을 이용하여 계수들 조정*하는 방식으로, 함수가 줄어드는 방향으로 계수를 계속 갱신해 나간다. 우리는 함수가 한번에 얼마나 줄어들게 하는지를 결정하는 파라미터를 정해야하며 이를 학습률이라고 한다. 이때, 파라미터가 너무 작으면 최적화가 오래걸리고, 너무 크면 최솟값을 지나칠 수 있다.

경사하강법에 대한 포스트를 이어서 쓸 것이므로 해당 포스트에서 자세하게 공부해보자.

경사하강법에 대한 포스트 바로가기

3. Scikit-learn 을 이용한 선형 회귀

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# 1. 데이터 준비
num_point = 100

# 모델이 찾아내야 할 실제 정답
true_W = 3.5
true_b = 8.0

# X와 Y 데이터를 랜덤으로 생성
noise = np.random.randn(num_point,1)
X = 2 * np.random.rand(num_point,1)
Y = true_W * X + true_b + noise

# 2. 선형 회귀 모델 생성 및 학습
model = LinearRegression()
model.fit(X,Y)

# 결과 확인
predict_W = model.coef_[0][0]
predict_b = model.intercept_[0]

print(f"true W :{true_W}")
print(f"predict W : {predict_W:.2f}")
print("-" * 30)
print(f"true b : {true_b}")
print(f"predict b {predict_b:.2f}")

true W :3.5
predict W : 3.67
------------------------------
true b : 8.0
predict b 7.81

보이는 것과 같이 열심히 예측하는걸 볼 수 있다! 파이썬의 scikit-learn 라이브러리를 사용하면 Linear Regression() 메서드로 한줄로 선형 회귀를 할 수 있다.

• model.coef_에는 학습된 $W$가 저장되며 2차원 배열 형태로 저장된다.

  • 선형회귀는 여러 개의 특성 ($X$) 를 가질 수 있기때문에(다중 선형 회귀) scikit-learn는 각 특성에 대한 기울기가 별도로 계산되어 배열 형태로 저장된다.

• model.intercept_ 에는 학습된 편향 $b$가 저장된다.

다음 포스트에서 경사 하강법을 씹고 뜯고 맛보도록 하자....!