Differential Expression Gene(DEG) with R

---

DEG(차등 발현 유전자 분석)은 microarray나 RNA-seq을 이용한 분석방법

 

같은 유전자에 대해 sample별로 발현량을 비교하는 방법이다.

 

DEG 분석 후 GO, KEGG pathway 분석 등을 통해 차등 발현한 유전자의 fuctional analysis도 함께 진행한다.

---

목차

 

1. DEG 분석 (차등 발현된 유전자 골라내기)

 

2. 간단한 GO, KEGG-pathway

 

3. web을 이용한 방법

---

1. DEG 분석 (차등 발현된 유전자 골라내기)

GEO에서 받은 데이터와 R limma package를 이용하여 분석한다.

# limma 설치 및 로드

BiocManager::install('limma')

# RMA 전처리한 data

GSEdata

# limma package를 이용하기 위해 dataframe으로 변형

eset <- data.frame(exprs(GSEdata))

# featureData

fset <- fData(GSEdata)

featuredata를 이용하여 Gene_symbol을 확인한다. (여러 gene 중 같은 것을 묶기 위해서)

 

eset data에 Gene_symbol 칼럼을 만들고 같은 유전자끼리 값을 평균 낸다.

# featureData에서 gene_symbol 가져오기 (있는 경우에)

gene_id = fset[,'Gene_Symbol']

# eset에 gene_id 컬럼을 만들고 평균

eset$gene_id = gene_id

eset <- aggregate(.~gene,data=eset,mean)

rownames(eset) <- eset$gene_id

eset <- eset[,-1]

phenoData를 이용하여 sample을 구분한다.

# phenoData에서 sample을 구분할 수 있는 값으로 이용한다.

grd <- pset$tissue

# limma에서 활용할 수 있도록 행렬로 만든다.

design <- model.matrix(~ 0 + grd)

colnames(design) <- c('ctrl','test') # 간편하게 사용하기 위해서

DEG 분석, 대비 행렬(test - ctrl 값)을 갖는 행렬을 만들고 통계분석을 진행한다.

# sample group화 한 design과 eset 연결

fit <- lmFit(eset, design)

# 대비행렬 만들기 test - ctrl

cont <- makeContrasts(test-ctrl, levels=design)

# fit 데이터와 cont data 연결

fit_cont <- contrasts.fit(fit, cont)

# 통계분석

fit_cont <- eBayes(fit_cont)

# 유전자 별 통계분석 결과를 알려주는 Table

result <- topTable(fit_cont, number=Inf, lfc=1, adjust='BH')

p-value와 FC값 기준으로 원하는 유전자만 골라낸다.

# p-value 0.05, logFC 2를 기준으로 필터링

target <- result[result$P.Value <0.05 & result$logFC >2,]

두 sample 간 차등 발현하는 유전자에 대한 정보를 얻을 수 있다.

 

▶ Volcanoplot을 이용하여 유전자들의 p-value, FC 분포를 확인한다.

library(EnhanceVolcano)

EnhancedVolcano(result,lab=rownames(result)
                      ,x = 'logFC'
                      ,y = 'P.Value'
                      ,xlab = bquote(~Log[2]~ 'fold change')
                      ,pCutoff = 0.05
                      ,FCcutoff = 2.0
                      ,pointSize = 2.0
                      ,labSize = 5.0
                      ,colAlpha = 1
                      ,legendPsition = 'right'
                      ,legendLabSize = 14
                      ,legendIconSize = 5.0)

---

2. 간단한 GO, KEGG-pathway 분석

- clusterprofiler, enrichplot, AnnotationDBi, org.Hs.eg.db를 이용한다.

BiocManager::install('clusterProfiler')

library(clusterProfiler)

# 위 4개 모두 import

DEG 분석에서 사용했던 gene_symbol을 기준으로 org.Hs.eg.db에서 ENTREZID를 찾는다. (여러 가지 ID를 이용할 수 있다.)

# org.Hs.eg.db : 사람 유전자 db

eg <- bitr(rownames(target),fromType='Symbol',toType='ENTREZID',orgDb='org.Hs.eg.db')

# 저장(다양하게 활용 ex. 웹에서)

write.csv(eg,file='./target.csv')

GO 분석 후 dotplot 그리기

go <- enrichGO(eg$ENTREZID,OrgDb = 'org.Hs.eg.db',ont='all')

# dotplot

dotplot(go, split='ONTOLOGY',showCategory=10) + facet_grid(ONTOLOGY~.,scale='free')

KEGG 분석 후 dotplot 그리기

kk <- enrichKEGG(gene=eg$ENTREZID, organism='hsa',pvalueCutoff = 0.05)

dotplot(kk)

 

---

3. web을 이용한 방법

1. DAVID

 

DAVID 페이지에 접속한다.

Tools에서 fuction annotation을 클릭한다.

upload 란에 1번 칸에 위에서 저장한 csv파일에서 ENTREZID를 복사 붙여 넣기 한다.(모두 입력)

2번 칸에서 ID 종류를 입력한다. 이번 예시에서는 ENTREZID를 사용했다.

gene list를 클릭하고 submit 한다.

결과 페이지에서 원하는 분석을 클릭하면 완료

2. string

 

string 페이지에 접속한다.

search, multiple protein을 클릭한다.

1번 칸에 ENTREZID를 모두 붙여넣기 한 후

2번 칸에서 organism을 선택한다.