RNA-seq 분석 환경 구축 및 데이터 가져오기

학습 목표

  • RNA-seq와 차등 발현 유전자 분석의 전체 흐름 파악
  • 실험 설계의 핵심 고려사항 이해
  • R 언어를 활용한 데이터 분석 방법 습득

1. RNA-seq 개요

최근 10년간 RNA-seq은 전사체 차등 발현 및 mRNA 선택적 스플라이싱 분석의 핵심 기술로 자리 잡았다. 특정 조건 하에서 어떤 유전자나 전사체의 발현량이 변하는지를 정확히 식별하는 것은 생물학적 반응 메커니즘을 이해하는 데 필수적이다.

이 토리얼에서는 여러 R 패키지를 활용하여 RNA-seq 분석의 전 과정을 다룬다. 데이터 불러오기부터 가상 카운트의 정규화, 품질 평가와 샘플 간 관계 탐색, 차등 발현 분석, 결과 시각화를 거쳐 다운스트림 기능 분석까지 수행한다.

2. 실험 데이터

Kenny PJ et al, Cell Rep 2014의 일부 데이터를 사용한다. HEK293F 세포에서 MOV10 유전자를 과발현시키거나 녹다운하여 발현 변화를 유도한 실험이다.

처리 조건123
MOV10 유전자과발현녹다운무관 siRNA
설명오버익스프레션유전자 침묵대조군

해당 데이터셋의 원시 서열은 Sequence Read Archive(SRA)에서 다운로드받아 Linux/Unix 환경에서 파이프라인 도구로 전처리하였다.

MOV10: microRNA 경로와 관련된 RNA 헬리케이스로, 성별 관련 발달에 연관되어 있다.

3. 연구 질문

  1. MOV10 발현 변화가 세포에 미치는 영향은 무엇인가?
  2. 각 변화 간 공통적인 특성이 존재하는가?

4. 작업 환경 설정

RStudio에서 새 프로젝트를 생성한다.

  1. File 메뉴 → New Project 선택
  2. New Directory → DEanalysis 디렉토리 생성
  3. RStudio가 해당 프로젝트를 자동으로 연다

getwd()로 작업 디렉토리를 확인하면 .../DEanalysis 형태로 출력되어야 한다. 이 디렉토리 내에 meta와 results 두 개의 하위 폴더를 추가로 만든다.

분석용 파일(MOV10 데이터)을 다운로드하여 압축 해제하면 data 폴더가 생성되며, 각 샘플별 하위 디렉토리가 포함된다. 또한 전사체 식별자를 유전자명으로 변환할注해석注 파일도 다운로드한다. 이 파일은 R의 AnnotationHub 패키지에서 추출한 것이다.

프로젝트 폴더에서 de_script.R 파일을 새로 만들어 다음 주석을 작성하고 저장한다.

## DESeq2를 이용한 유전자 수준 차등 발현 분석

5. 패키지 로드

분석에 필요한 패키지를 불러온다. 일부는 CRAN, 나머지는 Bioconductor에서 설치한다.

library(DESeq2)
library(tidyverse)
library(RColorBrewer)
library(pheatmap)
library(DEGreport)
library(tximport)
library(ggplot2)
library(ggrepel)

6. 데이터 불러오기

Salmon의 주요 산출물인 quant.sf 파일은 각 샘마다 하나씩 존재하며, 다음 정보를 담고 있다.

  • 전사체 식별자
  • 전사체 길이
  • 유효 길이(effective length)
  • TPM(transcripts per million, 유효 길이 기반 산출)
  • 추정 카운트

유효 길이(effective length): 서열 구성에 따라 동일한 실제 길이의 전사체라도 샘플링 확률이 달라진다. 샘플링 가능성이 높은 전사체는 더 큰 유효 길이를 갖게 되며, 이는 서열 특이적 편향과 GC 편향을 보정한 길이이다.

tximport 패키지로 quant.sf 파일을 DESeq2용 형태로 준비한다. 먼저 각 파일 경로를 변수에 저장하고, 출력 행렬에서 샘플을 구분할 수 있도록 이름을 부여한다.

## 파일 목록 작성
sample_dirs <- list.files(path = "./data", full.names = TRUE, pattern = "salmon$")

## 경로 벡터 생성
quant_files <- file.path(sample_dirs, "quant.sf")

## 샘플별 이름 지정
names(quant_files) <- str_replace(sample_dirs, "./data/", "") %>% 
                     str_replace("\\.salmon$", "")

Salmon 색인은 Ensembl ID 기반으로 구축되었으므로, tximport가 어느 유전자에 속하는지 알 수 있도록注해석注 정보가 필요하다.

# 注해석注 파일 로드
gene_map <- read.delim("tx2gene_grch38_ens94.txt")

# 구조 확인
gene_map %>% View()

gene_map은 전사체 ID, 유전체 ID, 유전체 symbol 세 열로 구성된 데이터프레임이다.

tximport()는 Salmon, Kallisto 등 외부 도구의 전사체 수준 카운트를 불러와 유전체 수준으로 집계하거나 전사체 수준 행렬을 출력한다. quant.sf의 풍부도 추정치나 대안값 계산 여부를 옵션으로 지정할 수 있다.

DESeq2 분석에는 유전체 수준의 비정규화된 "원시" 카운트가 필요하다. 기본값이 유전체 수준 카운트 행렬이므로, countsFromAbundance 매개변수만 조정하여 "원시" 카운트를 획득한다.

옵션설명
noTPM(스케일링 값)과 NumReads("원시" 카운트)를 유전체 수준으로 축소
scaledTPMTPM을 라이브러리 크기만큼 확장하여 "원시" 카운트로 사용
lengthScaledTPMTPM × featureLength × 라이브러리 크기로 "원시" 카운트換算, 원시 카운트와 동일한 척도

TPM 계산 과정:

  1. RPK(reads per kilobase): 리드 카운트를 킬로베이스 단위 유전체 길이로 나눔
  2. "per million" 스케일링 팩터: 모든 RPK 합계를 1,000,000으로 나눔
  3. 최종 TPM: 각 RPK를 스케일링 팩터로 나눔
# tximport 실행
txi <- tximport(
  quant_files, 
  type = "salmon", 
  tx2gene = gene_map[, c("tx_id", "ensgene")], 
  countsFromAbundance = "lengthScaledTPM"
)

7. 데이터 탐색

txi 객체는 풍부도(abundance), 카운트(counts), 길이(length) 행렬을 포함한 리스트이다. countsFromAbundance 요소에는 사용된 문자 인자가 저장된다. 길이 행렬은 각 유전체의 평균 전사체 길이를 담아 유전체 수준 분석의 오프셋으로 활용할 수 있다.

attributes(txi)

$names
[1] "abundance"           "counts"              "length"             
[4] "countsFromAbundance"

txi 객체를 DESeq2 입력으로 그대로 사용하되, 다음 단계까지 저장해둔다. 카운트 행렬을 살펴보면 소수점 값이 포함되어 있으므로, 반올림하여 정수로 변환하고 데이터프레임으로 저장한다.

# 카운트 확인
txi$counts %>% View()

# 데이터프레임으로 변환
expr_data <- txi$counts %>% 
  round() %>% 
  as.data.frame()

태그: RNA-seq DESeq2 tximport Salmon 차등발현

7월 26일 01:58에 게시됨

인기 태그

©2026 괴물 클럽