RNA seq 데이터 처리에 대한 어려움을 느끼는 생물학도들에게

RNA-Seq 데이터 처리는 생물정보학 연구에서 매우 중요한 분야입니다. 저도 RNA-Seq 써먹을 때가 있어서 배우려고 했는데 학부 때 배운 적이 없어서 뒤늦게 배우느라 고생을 많이 했던 기억이 나네요. 교수님도, 선배들도 모르는데 저한테 시키셨던 ㅅㅂ. 저희 학교엔 정보학 수업도 딱히 없어서 맨바닥에 헤딩하는 느낌으로 배웠습니다.

RNA Seq 분석을 위해 저는 R 프로그램을 사용했는데 R이 가장 대중적인 프로그램이기도 하고 RNA-Seq 데이터 처리와 분석에 유용한 도구를 많이 제공하더라고요. R 프로그램을 사용하여 RNA-Seq 데이터를 처리하는 예제를 제공하여, 생물학도 분들이 보다 쉽게 데이터 처리와 분석을 할 수 있도록 도움을 드리고 싶습니다. 저같은 고생 안하시면 좋겠네요.

1. RNA-Seq 데이터 전처리

RNA-Seq 데이터를 분석하기 전에, 먼저 데이터를 전처리하는 과정이 필요합니다. 이 과정은 주로 다음과 같은 단계로 이루어집니다.

1) Quality Control: 데이터의 품질을 검사하고 필요하다면 제거합니다.
2) Read Trimming: Low-quality reads와 adapter contamination 등을 제거합니다.
3) Read Alignment: 참조 시퀀스에 reads를 매핑합니다.
4) Counting: 각 gene 또는 exon에 대한 read 수를 세어줍니다.

이러한 전처리 과정에서 많은 R 패키지들이 사용됩니다. 예를 들어, FastQC 패키지는 Quality Control을 위해, cutadapt 패키지는 Read Trimming을 위해, STAR 패키지는 Read Alignment을 위해, featureCounts 패키지는 Counting을 위해 사용됩니다.

다음은 cutadapt 패키지를 사용하여 RNA-Seq 데이터를 trimming하는 예제 코드입니다.

```
library(cutadapt)

# Read data from input file
input_file <- "path/to/input.fastq"
input_data <- readFastq(input_file)

# Trim adapters from reads
trimmed_data <- adapterTrim(input_data, "AGATCGGAAGAGCACACGTCTGAACTCCAGTCA",
min.length = 50)

# Write output to file
output_file <- "path/to/output.fastq"
writeFastq(trimmed_data, output_file)

```

이 코드에서는 cutadapt 패키지의 **`adapterTrim()`** 함수를 사용하여 adapter를 trimming합니다. 먼저 **`readFastq()`** 함수로 입력 파일을 읽어들인 다음, **`adapterTrim()`** 함수로 adapter sequence를 지정하고 min.length 옵션으로 최소 read 길이를 지정합니다. 마지막으로 **`writeFastq()`** 함수를 사용하여 결과를 출력합니다.

2. RNA-Seq 데이터 분석

RNA-Seq 데이터의 분석에는 주로 다음과 같은 단계가 포함됩니다.

1) Differential Expression Analysis: 여러 조건 또는 그룹 간에 gene expression이 어떻게 다른지 비교합니다.
2) Gene Ontology (GO) Analysis: Differential expression한 gene들을 functional annotation과 함께 분석합니다.
3) Pathway Analysis: Differential expression한 gene들을 기반으로 biological pathway와 interaction network를 분석합니다.

이러한 분석 과정에서 많은 R 패키지들이 사용됩니다. 예를 들어, DESeq2 패키지는 Differential Expression Analysis을 위해, topGO 패키지는 GO Analysis를 위해, pathwayTools 패키지는 Pathway Analysis를 위해 사용됩니다.

다음은 DESeq2 패키지를 사용하여 RNA-Seq 데이터의 Differential Expression Analysis를 수행하는 예제 코드입니다.

```
library(DESeq2)

# Read count data from input file
input_file <- "path/to/input.counts"
count_data <- read.table(input_file, header=TRUE, row.names=1)

# Create DESeqDataSet object
dds <- DESeqDataSetFromMatrix(countData=count_data, colData=col_data,
design=~condition)

# Perform differential expression analysis
dds <- DESeq(dds)

# Extract differentially expressed genes
res <- results(dds, contrast=c("condition", "group1", "group2"))
de_genes <- res[which(res$padj < 0.05 & abs(res$log2FoldChange) > 1),]

# Plot results
plotCounts(dds, gene=de_genes$gene_id, intgroup="condition")

```

이 코드에서는 DESeq2 패키지의 **`DESeq()`** 함수를 사용하여 Differential Expression Analysis를 수행합니다. 먼저 **`read.table()`** 함수로 입력 파일을 읽어들인 다음, **`DESeqDataSetFromMatrix()`** 함수로 count 데이터와 sample metadata를 합쳐 DESeqDataSet object를 생성합니다. **`DESeq()`** 함수로 Differential Expression Analysis를 수행한 다음, **`results()`** 함수를 사용하여 differentially expressed gene을 추출합니다. 마지막으로 **`plotCounts()`** 함수로 결과를 시각화합니다.

또한, GO Analysis와 Pathway Analysis를 위해서는 differential expression한 gene 리스트를 가지고 Enrichment Analysis를 수행해야 합니다. 이 과정에서도 R 프로그램을 활용할 수 있습니다. 예를 들어, topGO 패키지를 사용하여 GO Analysis를 수행하는 예제 코드는 다음과 같습니다.

```
library(topGO)

# Load gene ontology data
obo_file <- "path/to/gene_ontology.obo"
data_file <- "path/to/gene_ontology_annotations.txt"
godata <- read.GOBriefings(oboFile=obo_file, annotFile=data_file)

# Create topGOdata object
geneset <- as.integer(rownames(de_genes))
topGOdata <- new("topGOdata", ontology="BP", allGenes=geneset,
gene2GO=godata, nodeSize=10)

# Run enrichment analysis
result <- runTest(topGOdata, algorithm="weight01", statistic="fisher")

# Plot results
png("GO_enrichment.png", width=800, height=600)
plotEnrichment(result, showCategory=20)
dev.off()

```

이 코드에서는 topGO 패키지의 **`runTest()`** 함수를 사용하여 GO Analysis를 수행합니다. 먼저 **`read.GOBriefings()`** 함수로 gene ontology data를 로딩한 다음, **`new()`** 함수로 topGOdata object를 생성합니다. **`runTest()`** 함수로 Enrichment Analysis를 수행한 다음, **`plotEnrichment()`** 함수로 결과를 시각화합니다.

3. Pathway Analysis를 수행하는 예제 코드는 다음과 같습니다.

```
scssCopy code
library(pathwayTools)

# Load pathway database
db_file <- "path/to/pathway_database.db"
pw_db <- loadDb(db_file)

# Convert gene IDs to Entrez IDs
de_genes_entrez <- convertIDs(de_genes$gene_id, fromType="SYMBOL",
toType="ENTREZID", org="Hs")

# Run pathway analysis
result <- pathwayAnalysis(de_genes_entrez, pw_db)

# Plot results
png("pathway_enrichment.png", width=800, height=600)
plotPathway(result)
dev.off()

```

이 코드에서는 pathwayTools 패키지의 **`pathwayAnalysis()`** 함수를 사용하여 Pathway Analysis를 수행합니다. 먼저 **`loadDb()`** 함수로 pathway database를 로딩한 다음, **`convertIDs()`** 함수로 gene ID를 Entrez ID로 변환합니다. **`pathwayAnalysis()`** 함수로 Enrichment Analysis를 수행한 다음, **`plotPathway()`** 함수로 결과를 시각화합니다.

이처럼 R 프로그램은 RNA-Seq 데이터 처리와 분석에 있어서 매우 강력한 도구입니다. 하지만 처음 시작하는 학생들은 문법적인 어려움과 데이터 처리 방법의 이해가 필요한데 저처럼 배울 곳이 없어서 힘들어하시는 분들도 계시더라고요. 저처럼 헤매지 않으시길 바라며 요새는 온라인 강의도 잘 돼있으니 들어보시면 많은 도움될겁니다.