TCGA Data structure & survival analysis

/*

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

TCGA 사이트에서 데이터를 받아봅시다.

http://cancergenome.nih.gov/

오른쪽 위의 Launch Data Portal 버튼을 누르셔서 Data download 탭에서 원하는 데이터를 받으시면 됩니다.

BRCA1_Methyl.txt

Clinical.csv

Expression.txt

Mutation.txt

R 프로그래밍과 예제로 위의 데이터를 이용하여 데이터 구조를 분석해보겠습니다.

> setwd("D:/Analysis_data") # 데이터가 있는 경로로 R의 현재풀더 바꾸기

> dir()

[1] "BRCA1_Methyl.txt" "Clinical.csv"     "Expression.txt"  

[4] "Mutation.txt"     "lecture03_code.txt"   

> cli <- read.csv("Clinical.csv", header=T, na.strings="null")

> mut <- read.delim("Mutation.txt", header=T)

> met <- read.delim("BRCA1_Methyl.txt",header=T)

> exp <- read.delim("Expression.txt", header=T)

# mutation raw data 구조

> head(mut)

    SampleName GeneName  rsID              type

1 TCGA-36-1569    AARS2 novel Missense_Mutation

2 TCGA-23-1024    AARS2 novel Missense_Mutation

3 TCGA-25-1325    ABCA1 novel       Splice_Site

(다른 데이터는 row의 데이터 갯수가 많기 때문에 생략)

#BRCA1/2 mutation을 가진 데이터 샘플 추출

unique(x, incomparables = FALSE, ...)

unique returns a vector, data frame or array like x but with duplicate elements/rows removed.

>brca1.mut <- unique(as.character(mut[which(mut[,2] =="BRCA1" ),1]))  #mutation data에서 gene 이름이 brca1인 샘플 추출

>brca2.mut <- unique(as.character(mut[which(mut[,2] =="BRCA2" ),1])) #mutation data에서 gene 이름이 brca2인 샘플 추출

> brca2.mut

 [1] "TCGA-13-0792" "TCGA-24-1555" "TCGA-23-1120" "TCGA-04-1331"

 [5] "TCGA-23-1030" "TCGA-09-2050" "TCGA-13-0890" "TCGA-13-0885"

 [9] "TCGA-24-1103" "TCGA-13-1481"

# 마찬가지로 brca1.mut 데이터도 확인 가능

# BRCA1 expression level 데이터 추출하기 (microarray 데이터)

brca1.exp <- as.numeric(exp[which(exp[,1]=="BRCA1"),-1]) # BRCA1 gene에 대한 expression level 추출

head(brca1.exp)

# BRCA1 expression과 상관관계가 높은 methylation probe 추출

exp.cor <- apply(met[,-1], 1, function(x)cor(x, brca1.exp))

#apply 함수를 이용하여 BRCA1 expression level과 각각의 methylation probe의 beta value간의 pearson correlation coefficient를 구한다.

met[,-1]은 met[,1]을 제외한 모든 데이터에 대해서 다음의 function을 호출하라는 뜻

위의 apply 함수는 아래의 for loop로 구현할 수 있습니다.

cor.result<-list()

 for (i in c(1:9)){

   cor.result[i] <- cor(as.numeric(met[i,-1]), brca1.exp)

 }

# BRCA1 expression이 달라진 methylation prove를 찾는다.

TestMatrix <- met[which(abs(exp.cor) > 0.3),-1] 

# methylation pattern이 다른 샘플을 찾기 위해heatmap을 그린다.

heatmap(as.matrix(TestMatrix))

# cluster dendrogram을 이용해 BRCA1 methylation 이 달라진 샘플 찾기###

AA<- hclust(dist(t(TestMatrix)))

plot(AA)

#R graphic 창에 나타난 dendrogram에서 가장 왼쪽에 차이나는 샘플들을 선택(왼쪽클릭)한 후 stop(오른쪽클릭해서 나타나는 stop선택)

> M<-identify(AA)  

> brca1.met<- names(met[,as.numeric(M[[1]])])  #BRCA1 methylation 이 다른 sample을 변수에 할당

> brca1.met  #확인

 [1] "TCGA.04.1542" "TCGA.09.0367" "TCGA.09.1665" "TCGA.09.1666"

 [5] "TCGA.10.0926" "TCGA.10.0931" "TCGA.13.0714" "TCGA.13.0717"

 [9] "TCGA.13.0723" "TCGA.13.0730" "TCGA.13.0760" "TCGA.13.0765"

[13] "TCGA.13.0792" "TCGA.13.0795" "TCGA.13.0801" "TCGA.13.0802"

[17] "TCGA.13.0885" "TCGA.13.1412" "TCGA.13.1483" "TCGA.13.1489"

[21] "TCGA.13.1494" "TCGA.13.1497" "TCGA.13.1498" "TCGA.13.1505"

[25] "TCGA.13.1506" "TCGA.24.1105"

# survival analysis 하기

type <- c(rep('wild', 317)) # type이 wild 인 변수 생성

names(type)<-cli[,1] #column을 환자 id로 바꿔줌

brca1.met<-gsub(".","-",brca1.met, fixed=T) #샘플의 이름 표기가 서로 다르므로 이를 바꿔준다.

brca1.mut #분석할 테이블 확인

brca2.mut

brca1.met #교재 오타

type[brca1.mut]   #샘플별 type을 바꿔줌

type[brca1.met] = 'brca1.met'

type[brca1.mut] = 'brca1.mut'

type[brca2.mut] = 'brca2.mut'

library(survival) #survival 분석을 위한 함수를 불러온다.

out <- survfit(Surv(cli$days_to_last_followup) ~ type)

color <- c("green","blue", "red","black") #4개 군을 나타낼 색깔을 정해둔다.

plot(out, col=color , main="Association of BRCA1/2 Mutations with Survival", xlab="Time,days", ylab="Proportion") #plot

legend(4200, 0.9, c("BRCA1 mutation","BRCA1 methylation","BRCA2 mutation","wild"), col=color, lty=1, lwd=3)