Gene Expression data로부터 PCA 분석하기

 

/*

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

1. Gene Expression Data 구조

2. ALL 공개데이터를 이용한 PCA 분석

3. TCGA BRCA RNAseq V2를 이용한 PCA 분석

 

 

 

 

1. Gene Expression Data 구조

 

Gene expression 데이터는 대표적으로 Microarray와 RNA-seq이 있습니다.

일반적으로 이 데이터는 matrix 구조로 되어있으며, 구조는 다음과 같습니다.

 

column: sample

row: gene

value: expression

 

 

 - TCGA Gene Expression 데이터 예

 

특정 유전자가 각 환자별로 얼마나 많이 발현되었는지 알 수 있습니다.

 

TCGA 등의 공개데이터나 혹은 가지고 계신 expression data가 있다면

read.table 등의 함수를 이용하여 불러오셔서 분석하시면 되겠습니다.

 

데이터가 없는데 연습하고자 하는 경우 아래에서 ALL 데이터를 가져오시면 됩니다.

 

 

 

2. ALL 공개데이터를 이용한 PCA 분석

 

> source("http://www.bioconductor.org/biocLite.R")

> biocLite("ALL")
> library("ALL")
> data("ALL")

> exp <- exprs(ALL)

 

> dim(exp)
[1] 12625   128

 

# 12625개의 gene과 128개의 sample로 이루어진 데이터입니다.

# 데이터에 대한 자세한 이해는 ?ALL 를 이용해 확인하세요.

 

 

> heatmap(exp[1:100,])

 

# 기본적으로 microarray data의 전체적인 양상을 보기위해서는 heatmap을 자주 이용합니다.

# 각 샘플과 gene의 expression 값을 바로 볼 수 있기 떄문에 직관적입니다.

 

 

 

 

# expression data로부터 Principle component를 구해보겠습니다.

# PCA 분석은 prcomp라는 함수가 있습니다.

 

> exp.obj <- prcomp(exp,center=TRUE,scale=TRUE)

> str(exp.obj)
List of 5
 $ sdev    : num [1:128] 10.945 1.101 0.932 0.753 0.629 ...
 $ rotation: num [1:128, 1:128] -0.0897 -0.0881 -0.0889 -0.0872 -0.0884 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:128] "01005" "01010" "03002" "04006" ...
  .. ..$ : chr [1:128] "PC1" "PC2" "PC3" "PC4" ...
 $ center  : Named num [1:128] 5.63 5.65 5.63 5.63 5.63 ...
  ..- attr(*, "names")= chr [1:128] "01005" "01010" "03002" "04006" ...
 $ scale   : Named num [1:128] 1.88 1.85 1.89 1.86 1.97 ...
  ..- attr(*, "names")= chr [1:128] "01005" "01010" "03002" "04006" ...
 $ x       : num [1:12625, 1:128] -11.77 3.77 10.55 -2.19 -1.17 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...
  .. ..$ : chr [1:128] "PC1" "PC2" "PC3" "PC4" ...
 - attr(*, "class")= chr "prcomp"
>

> plot(exp.obj$x[,1], exp.obj$x[,2])
# PC1, PC2로 plot 그리기

 

# 위 그림에서 결국 exp.obj의 x 변수들이 PC에 해당됩니다.

 

# plot을 그릴 때 환자군, 일반군이나 다른종류의 group별로 색을 이용하여 점을 찍고,

# 그것들을 유의하게 나눠지는 Principle component가 찾아지면 정말 이상적이겠지요.

 

 

 

 

 

 

3. TCGA BRCA RNAseq V2를 이용한 PCA 분석

 

BRCA_RNASeqV2.z01

BRCA_RNASeqV2.z02

BRCA_RNASeqV2.zip

 

 

#데이터 읽기

> BRCA <- read.table("BRCA_RNASeqV2.txt", header=TRUE, row.names=1)

 

 

#용량이 큰 matrix 데이터의 경우 col 수가 워낙 많기 때문에, head()로 데이터를 보더라도 페이지가 넘어가는 경우가 많아서

#저는 직접 범위를 입력해서 출력합니다.

> BRCA[1:10,1:10]

 

 

첫번쨰 row는 데이터가 아니라 group이네요. column의 수와 level은 다음과 같습니다.

1-100: N

101-200: S1

201-300: S2

301-400: S3

 

# 데이터 전체에서 그룹정보가 담긴 row가 있으면 계산이 불가능하기 때문에 다른 변수에 저장해줍니다.

# factor 데이터 구조로 바꿔줍니다. 일단 list형태에서 matrix로 구조를 바꿔주어야 factorize가 가능합니다.

> group <- BRCA[1,]

 

 

> group <- factor(as.matrix(group))

> levels(group)
[1] "N"  "S1" "S2" "S3"

 

# BRCA에서 group을 나타내는 row를 제거합니다.

 

> BRCA <- BRCA[-1,]

> BRCA[1:5,1:10]

 

 

# PCA 분석을 합니다.

> prcomp(BRCA,scale=TRUE)
다음에 오류가 있습니다colMeans(x, na.rm = TRUE) : 'x'는 반드시 수치형이어야 합니다

 

# 기본적으로 list 형태를 바로 prcomp 함수에 넣을 수 없기 때문에,

# list를 vector로 return하는 sapply함수를 이용하여 numeric으로 형변환을 해준 후에 matrix로 바꿔줍니다.

> BRCAnm <- as.matrix(sapply(BRCA, as.numeric))

 

# PCA 분석 후에 PC1, PC2로 plot 그려보기

> BRCAnm.obj <- prcomp(BRCAnm,scale=TRUE)

> plot(BRCAnm.obj$x[,1], BRCAnm.obj$x[,2])

 

 

 

 

# log transform 후에 PCA analysis 하기

# 0인 경우 문제가 발생하기 때문에 1을 더해준 후에 log transform

> BRCAnm_log <- log2( BRCAnm+1 )
> BRCAnm_log.obj <- prcomp(BRCAnm_log,scale=TRUE)
> plot(BRCAnm_log.obj$x[,1], BRCAnm_log.obj$x[,2])