RNA-seq 데이터로 DEG 구하기 - R

/*

 * http://sosal.kr/
 * made by so_Sal
 */

DEG (Differential Expression Gene). 차별 유전자 발현

Definitions:

1. gene expression that responds to signals or triggers; a means of gene regulation,

    effects of certain hormones on protein biosynthesis.

실험군이 대조군에 비해 특별히 expression이 많이 일어나거나 적게 일어나는 gene들을 DEG라고 합니다.

예를들어, 10명에 해당하는 환자들을 모아, 그들의 tumor cell과 normal cell을 분석하여,

DEG 분석을 통해 tumor cell에서 expression value가 상당히 높거나 낮은 gene들의 목록을 찾는다면

그 gene들이 혹시 환자들의 질병에 영향을 미치지 않았을까, 연구의 폭을 줄여볼 수 있을것입니다.

http://www.ncbi.nlm.nih.gov/books/NBK10061/

- 예제데이터

cancer

normal

100개의 gene에 대한 10명의 cancer cell과 10명의 normal cell의 sample을 TCGA에서 긁어왔습니다.

Cancer는 LUAD 데이터이며, RNAseq V2 Expression data 입니다.

대충 긁어오느라 각 파일의 columns (환자목록)은 다르지만, DEG를 구하는데는 큰 문제 없습니다.

Figure 1.0 - 데이터 구조

rows: Gene list

columns: Patient list.

 

 

# bioconductor 및 DEGseq library

source("http://bioconductor.org/biocLite.R")
biocLite("DEGseq")
library(DEGseq)

 

-R Programming

cancer <- readGeneExp(file="cancer", geneCol=1, valCol=c(2:11))

normal <- readGeneExp(file="normal", geneCol=1, valCol=c(2:11))

# normal, cancer 파일 불러오기

DEGexp(geneExpMatrix1 = normal, geneCol1 = 1, expCol1 = c(2:11), groupLabel1 = "Normal", geneExpMatrix2 = cancer, geneCol2=1, expCol2=c(2:11), groupLabel2 = "Cancer", method="MARS", output='output')

DEGexp 함수는 아래와 같이 cancer, normal 데이터의 비교와 DEG 목록들을 파일로 만들어준다.

output 파일 안에는 아래 5개의 figure가 들어있고, output.html로 보여준다.

output_score.txt는 각 gene의 DEG score를 보여준다.

 

boxplot, barplot 등 다양한 figure들을 통해 두 군의 비교 결과를 보기 쉽게 전달해준다.

fold change log값과, normalization의 fold_change log값, 그리고 z-score 값을 줍니다.

z-score를 정렬하여 프로그램을 출력해주며, 어디까지 DEG인지 자르는것은 데이터를 분석하는분의 몫이 되겠습니다.

유의값을 0.01로 한다면 이경우 gene이 100개이기 때문에 ABCA3|21 하나의 gene만이 DEG로 선택될것입니다.

- 간단 해석

LUAD tumor cell 10개와, normal cell 10개의 샘플 속에, 100개의 gene에 대한 RNA-seq 데이터를 분석하였고, 그 데이터에서 찾은 DEG  ABCA3|21 gene, ABCC3 등이 되겠습니다.

ABCA3|21 gene은 tumor cell이 normal cell에 비해 유의하게 RNA-seq 발현값이 증가하였고,

ABCC3|8714 gene은 normal cell이 tumor cell에 비해 유의하게 RNA-seq 발현값이 증가한 gene이라고 해석할 수 있겠습니다.

이것은 간단한 예제이므로 실제로 발현이 증가하였다고 보긴 어렵고,

TCGA에서 제공하는 모든 LUAD data의 비교를 통해 실제 실험을 하시기 바랍니다.