BioinformaticsAndMe

[R] 생존분석(Survival analysis) Start

BioinformaticsAndMe

#생존분석에 대한 개념적 설명은 아래 포스팅 참조

https://bioinformaticsandme.tistory.com/223

1. 생존분석 라이브러리 로딩

: 생존분석에 필요한 패키지(survival, survminer) 및 데이터프레임 관리 패키지(dplyr) 로딩

: 참고로, 패키지 설치 에러에는 여러 이유가 있겠으나, 일반적으로 R 최신 버전 상태에서 패키지를 설치하면 해결됨

*아래 튜토리얼은 'R version 3.6.1 (2019-07-05)'에서 수행됨

# 패키지 인스톨 후, 로딩
install.packages(c('survival', 'survminer', 'dplyr'))
library(survival)
library(survminer)
library(dplyr)

2. 예제 데이터셋 로딩

: 난소암(Ovarian cancer) 데이터를 생존분석의 예제로 사용

# glimpse 함수로 난소암 데이터의 내용/구조/타입을 파악
data(ovarian)
glimpse(ovarian)
Observations: 26
Variables: 6
$ futime   <dbl> 59, 115, 156, 421, 431, 448, 464, 475, 477, 563, 638,...  #생존 시간
$ fustat   <dbl> 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,...  #Censored data 여부
$ resid.ds <dbl> 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1,...  #종양 퇴행
$ age      <dbl> 72.3315, 74.4932, 66.4658, 53.3644, 50.3397, 56.4301,...  #환자 나이
$ rx       <dbl> 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, 1,...  #치료 요법 
$ ecog.ps  <dbl> 1, 1, 2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1,...  #환자 전신 활성도(ECOG)

# 데이터 Labeling
ovarian$rx <- factor(ovarian$rx, 
                     levels = c("1", "2"), 
                     labels = c("A", "B"))
ovarian$resid.ds <- factor(ovarian$resid.ds, 
                           levels = c("1", "2"), 
                           labels = c("no", "yes"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps, 
                          levels = c("1", "2"), 
                          labels = c("good", "bad"))

# 히스토그램으로 나이 데이터 분포 확인
hist(ovarian$age) 

# 위 이항분포 히스토그램 결과에 근거하여, '50 years' 기준의 두 그룹으로 나눔
ovarian <- ovarian %>% mutate(age_group = ifelse(age >=50, "old", "young"))
ovarian$age_group <- factor(ovarian$age_group)

3. Kaplan-Meier(카플란-마이어) 생존분석

# Surv 함수로 Survival object 생성 (Survival object는 survfit 함수에 인식되도록 특정칼럼을 컴파일한 형태)
# surv_object 결과의 '+' 표시는 Censored data를 의미
surv_object <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv_object 
[1]   59   115   156   421+  431   448+  464   475   477+  563   638 
[12]  744+  769+  770+  803+  855+ 1040+ 1106+ 1129+ 1206+ 1227+  268 
[23]  329   353   365   377+

# survfit 함수를 통해, Kaplan-Meier 곡선 Fitting
fit1 <- survfit(surv_object ~ rx, data = ovarian)
summary(fit1)
Call: survfit(formula = surv_object ~ rx, data = ovarian)

               rx=A 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  59     13       1    0.923  0.0739        0.789        1.000
  115     12       1    0.846  0.1001        0.671        1.000
  156     11       1    0.769  0.1169        0.571        1.000
  268     10       1    0.692  0.1280        0.482        0.995
  329      9       1    0.615  0.1349        0.400        0.946
  431      8       1    0.538  0.1383        0.326        0.891
  638      5       1    0.431  0.1467        0.221        0.840

               rx=B 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  353     13       1    0.923  0.0739        0.789        1.000
  365     12       1    0.846  0.1001        0.671        1.000
  464      9       1    0.752  0.1256        0.542        1.000
  475      8       1    0.658  0.1407        0.433        1.000
  563      7       1    0.564  0.1488        0.336        0.946

# ggsurvplot 함수로 생존분석 결과 시각화 (추가로 p value 표시)
# 그래프 내에 있는 수직바는 Censored data를 의미함
ggsurvplot(fit1, data = ovarian, pval = TRUE)

# resid.ds 기준으로 나뉜 두 그룹의 생존분석은 p-value가 약 0.05로 유의한 결과를 냄
fit2 <- survfit(surv_object ~ resid.ds, data = ovarian)
ggsurvplot(fit2, data = ovarian, pval = TRUE)

4. Cox proportional hazards model (콕스 비례위험모델)

: Hazard Ratio(HR)는 상대적 위험도이므로, 그 위험비는 항상 reference에 대한 특정값임

: HR>1(사망 위험 증가), HR<1(사망 위험 감소)

: 아래 그림에서 0.25는 Treatment B를 받은 환자들이 Treatment A를 받은 환자들에 비해 사망률이 0.25배 감소했다는 의미

*Treatment A는 reference로 사용됨

# coxph 및 ggforest 함수를 통해, Hazard ratio 결과 출력
fit.coxph <- coxph(surv_object ~ rx + resid.ds + age_group + ecog.ps, 
                   data = ovarian)
ggforest(fit.coxph, data = ovarian)

#생존분석에 대한 개념적 설명은 아래 포스팅 참조

https://bioinformaticsandme.tistory.com/223

#Reference

1) http://www.gums.ac.ir/Upload/Modules/Contents/asset68/Medical%20Statistic%20Made%20Easy.pdf

2) https://www.datacamp.com/community/tutorials/survival-analysis-R#fourth

3) https://ko.wikipedia.org/wiki/%EC%83%9D%EC%A1%B4%EB%B6%84%EC%84%9D

4) https://namu.wiki/w/%EC%83%9D%EC%A1%B4%20%EB%B6%84%EC%84%9D

5) https://www.partek.com/webinar/survival-analysis-with-partek-genomics-suite-software/

6) https://ko.wikipedia.org/wiki/%EC%A4%91%EB%8F%84%EC%A0%88%EB%8B%A8

7) https://ko.wikipedia.org/wiki/%EB%A1%9C%EA%B7%B8%EC%88%9C%EC%9C%84%EB%B2%95

8) https://www.youtube.com/watch?v=czQ3l0QXxnA

9) https://rexsoft.org/?page_id=485

10) http://www.e-urol-sci.com/viewimage.asp?img=UrolSci_2018_29_5_223_240363_t6.jpg

[R] 생존분석(Survival analysis) End

BioinformaticsAndMe

[R] 히스토그램(Histogram) Start

BioinformaticsAndMe

1. 히스토그램(Histogram)?

: 히스토그램(Histogram)은 표로 되어 있는 도수분포표를 그래프로 나타낸 것

2. 히스토그램 예제

# R에 내장된 airquality 데이터셋
str(airquality)
'data.frame':	153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

# 히스토그램 생성을 위한 hist 함수 사용
Temperature <- airquality$Temp
hist(Temperature)

# hist 함수 파라미터 정보 추가
hist(Temperature,
     main="Maximum daily temperature at La Guardia Airport",
     xlab="Temperature in degrees Fahrenheit",
     xlim=c(50,100),
     col="darkmagenta",
     freq=FALSE
)

# text 함수를 사용하여, 특정값 Labeling
h <- hist(Temperature, ylim=c(0,40))
text(h$mids, h$counts, labels=h$counts, adj=c(0.5, -0.5))

# break 옵션으로 히스토그램 구간을 나눔
hist(Temperature, breaks=4, main="With breaks=4")
hist(Temperature, breaks=20, main="With breaks=20")

#Reference

1) https://terms.naver.com/entry.nhn?docId=707162&cid=42318&categoryId=42318

2) https://www.datamentor.io/r-programming/histogram/

[R] 히스토그램(Histogram) End

BioinformaticsAndMe

[R] Excel 읽기/쓰기 Start

BioinformaticsAndMe

1. R xlsx 패키지 설치

: R에서 엑셀 데이터의 읽기/처리/쓰기 작업을 위해, 'xlsx' 패키지 설치

*참고로 엑셀 데이터를 처리하는 여러 R 패키지들이 존재 (xlsx, readxl, XLConnect 등)

# xlsx 패키지
install.packages("xlsx")
library("xlsx")

2. Read an Excel file

: read.xlsx 함수를 사용하여, 엑셀 데이터 로딩

Expression_example.xlsx

# 예제 파일 경로 지정 후, 데이터 로딩
file <- 'Expression_example.xlsx'
res <- read.xlsx(file, 1)       # 첫번째 시트 읽기
res

3. Write data to an Excel file

: write.xlsx 함수를 사용하여, 엑셀 데이터 생성

# R 예제 데이터를 엑셀 형태로 파일 생성
write.xlsx(USArrests, file="myworkbook.xlsx",sheetName="USA Arrests")

#Reference

1) http://www.sthda.com/english/wiki/r-xlsx-package-a-quick-start-guide-to-manipulate-excel-files-in-r

[R] Excel 읽기/쓰기 End

BioinformaticsAndMe

[R] Merge Start

BioinformaticsAndMe

R Merge?

: R에서는 두 개의 data.frame을 병합하기 위해, Merge 함수를 사용

*데이터프레임들은 최소 1개 이상의 동일한 column(열)을 가져야 함

: merge(x, y, by, all)

* x - 데이터프레임1

* y - 데이터프레임2

* by - 공통된 column (동일 변수)

* all - Merging의 형태 (default값=FALSE)

# 데이터 준비 (첫번째 data.frame)
df1 = data.frame(CustomerId = c(1:6), Product = c(rep("Oven", 3), rep("Television", 3)))
df1
  CustomerId    Product
1          1       Oven
2          2       Oven
3          3       Oven
4          4 Television
5          5 Television
6          6 Television

# 데이터 준비 (두번째 data.frame)
df2 = data.frame(CustomerId = c(2, 4, 6), State = c(rep("California", 2), rep("Texas", 1)))
df2
  CustomerId      State
1          2 California
2          4 California
3          6      Texas

Merge 형태

1) Natural join: 두 데이터프레임에서 공통된 행만 유지 (all=FALSE)

2) Full outer join: 두 데이터프레임에서 모든 행을 유지 (all=TRUE)

3) Left outer join: 두 데이터프레임 중 좌측 x 프레임의 모든 행을 유지 (x=TRUE)

4) Right outer join: 두 데이터프레임 중 우측 y 프레임의 모든 행을 유지 (y=TRUE)

1) Natural join (Inner join)
df<-merge(x=df1,y=df2,by="CustomerId")
df
  CustomerId    Product      State
1          2       Oven California
2          4 Television California
3          6 Television      Texas

2) Full outer join 
df<-merge(x=df1,y=df2,by="CustomerId",all=TRUE)
df
  CustomerId    Product      State
1          1       Oven       <NA>
2          2       Oven California
3          3       Oven       <NA>
4          4 Television California
5          5 Television       <NA>
6          6 Television      Texas

3) Left outer join
df<-merge(x=df1,y=df2,by="CustomerId",all.x=TRUE)
df
  CustomerId    Product      State
1          1       Oven       <NA>
2          2       Oven California
3          3       Oven       <NA>
4          4 Television California
5          5 Television       <NA>
6          6 Television      Texas

4) Right outer join
df<-merge(x=df1,y=df2,by="CustomerId",all.y=TRUE)
df
  CustomerId    Product      State
1          2       Oven California
2          4 Television California
3          6 Television      Texas

#Reference

1) http://www.datasciencemadesimple.com/join-in-r-merge-in-r/

[R] Merge End

BioinformaticsAndMe

[R] Heatmap 시각화 Start

BioinformaticsAndMe

Heatmap?

: 히트맵(heatmap)은 2차원의 데이터를 계층적 클러스터링으로 시각화하는 방법

: 행과 열에 존재하는 데이터 값의 칼라 스케일을 조절하면서, 눈에 띄는 데이터 패턴을 찾고자 함

1) heatmap:  R base function, stats package

2) heatmap.2:  gplots R package

3) pheatmap:  pheatmap R package

4) d3heatmap:  d3heatmap R package

# 데이터 준비
df <- scale(mtcars)     #scale: 데이터 표준화
df
                         mpg        cyl        disp         hp       drat           wt       qsec
Mazda RX4          0.1508848 -0.1049878 -0.57061982 -0.5350928  0.5675137 -0.610399567 -0.7771651
Mazda RX4 Wag      0.1508848 -0.1049878 -0.57061982 -0.5350928  0.5675137 -0.349785269 -0.4637808
Datsun 710         0.4495434 -1.2248578 -0.99018209 -0.7830405  0.4739996 -0.917004624  0.4260068
Hornet 4 Drive     0.2172534 -0.1049878  0.22009369 -0.5350928 -0.9661175 -0.002299538  0.8904872
Hornet Sportabout -0.2307345  1.0148821  1.04308123  0.4129422 -0.8351978  0.227654255 -0.4637808
Valiant           -0.3302874 -0.1049878 -0.04616698 -0.6080186 -1.5646078  0.248094592  1.3269868

1) heatmap() function

: stat 패키지에 내장된 R의 기본 히트맵함수

: heatmap(x, scale = "row")

* x → 숫자 matrix

* scale → 숫자 값들을 row 혹은 column 방향으로 표준화 (Default 값은 row)

# 기본 히트맵 그리기
heatmap(df, scale = "none")     # none은 scale 과정 수행 x
#높은 값은 Red 계열, 낮은 값은 Yellow 계열

2) heatmap.2() function

: gplots 패키지로 이전 heatmap() 보다 확장성이 늘어난 히트맵함수

# heatmap.2 히트맵 그리기
install.packages("gplots")
library("gplots")
heatmap.2(df, scale = "none", col = bluered(100), 
          trace = "none", density.info = "none")

3) pheatmap() function

: pheatmap 패키지로 시각적으로 더 선명하고 깔끔해진 히트맵함수

# pheatmap 히트맵 그리기
install.packages('pheatmap')
library("pheatmap")
pheatmap(df, cutree_rows = 4)

4) d3heatmap() function

: d3heatmap 패키지로 Interactive하게 결과를 확인할 수 있는 히트맵함수

# d3heatmap 히트맵 그리기
install.packages('d3heatmap')
library("d3heatmap")
d3heatmap(scale(mtcars), colors = "RdYlBu",
          k_row = 4,     # row 그룹의 숫자
          k_col = 2      # column 그룹의 숫자
)

#Reference

1) https://www.datanovia.com/en/lessons/heatmap-in-r-static-and-interactive-visualization/

[R] Heatmap 시각화 End

BioinformaticsAndMe

[R] transform (데이터열생성 함수) Start

BioinformaticsAndMe

transform (데이터열생성 함수)

: transform 함수는 데이터프레임에 새로운 열(변수)을 생성하는 함수

: 아래 예제는 iris 데이터의 Sepal.Length x 4 값을 갖는 새로운 변수 Sepal_multiply 열을 생성

# transform 함수 실행
result_1 <- transform( iris, Sepal_multiply = 4*Sepal.Length )
head(result_1)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species Sepal_multiply
1          5.1         3.5          1.4         0.2  setosa           20.4
2          4.9         3.0          1.4         0.2  setosa           19.6
3          4.7         3.2          1.3         0.2  setosa           18.8
4          4.6         3.1          1.5         0.2  setosa           18.4
5          5.0         3.6          1.4         0.2  setosa           20.0
6          5.4         3.9          1.7         0.4  setosa           21.6

# cbind 함수를 사용해 동일한 결과를 얻을 수 있음
result_2 <- cbind(iris, data.frame(Sepal_multiply = 4*iris$Sepal.Length))
head(result_2)

[R] transform (데이터열생성 함수) End

BioinformaticsAndMe

[R] subset (데이터추출 함수) Start

BioinformaticsAndMe

subset (데이터추출 함수)

: subset 함수는 조건에 맞는 matrix 혹은 data.frame 결과를 추출

: 추출하고자 하는 조건이 복잡해질수록 구문이 길어지는 문제점을 보완하기 위한 함수

: 결측치가 포함된 데이터에서 에러를 쉽게 피할 수 있음

# iris 데이터에서 'setosa' 종 추출
result_1 <- subset(iris, Species=='setosa')
head(result_1)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

# Sepal.Length가 5를 넘고, 'setosa' 종 추출
result_2 <- subset(iris, Sepal.Length>5 & Species=='setosa')
head(result_2)
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1           5.1         3.5          1.4         0.2  setosa
6           5.4         3.9          1.7         0.4  setosa
11          5.4         3.7          1.5         0.2  setosa
15          5.8         4.0          1.2         0.2  setosa
16          5.7         4.4          1.5         0.4  setosa
17          5.4         3.9          1.3         0.4  setosa

# Sepal.Length가 5를 넘고, 'setosa' 종 추출한 데이터에서 select된 정보만 출력
result_3 <-subset(iris, Sepal.Length>5 & Species=='setosa', select=c('Petal.Length', 'Petal.Width', 'Species') )
head(result_3)
  

[R] subset (데이터추출 함수) End

BioinformaticsAndMe

[R] p value 보정 Start

BioinformaticsAndMe

P value 보정 (본페로니, FDR)

# Raw p value
pvals = c(.001, .002, .003, .01, .02, .05, .22, .59, .87, .88)

# Bonferroni correction
BONF = p.adjust(pvals, "bonferroni")

# FDR (False discovery rate)
FDR = p.adjust(pvals, "fdr")

# 결과 비교를 위한 데이터프레임 생성
res = data.frame(pvals=pvals, BONF=round(BONF, 3), FDR=round(FDR, 3))
res

#FDR 보정 과정 참고

https://bioinformaticsandme.tistory.com/129

[R] p value 보정 End

BioinformaticsAndMe

[R Shiny] 샤이니 소개 Start

BioinformaticsAndMe

What is R Shiny?

: R Shiny(샤이니)는 R 사용자들에게 Interactive web app 제작을 가능케하는 R 패키지

: Shiny code로 HTML, CSS로 제작되는 웹앱을 동등하게 구현

: 웹앱 제작에 대한 지식이 부족해도 R script로 간단하게 웹앱 개발

: R 사용자들은 분석 과정에서 보유한 데이터를 쉽게 어플리케이션에 연동 가능

Installing R Shiny

# R Shiny 패키지 설치
install.packages('shiny')

# Shiny 패키지 로딩
library(shiny)

Structure of a Shiny app

1. User Interface (UI)

: 사용자 인터페이스(UI)는 앱의 레이아웃과 모양을 결정

: 앱에서 사용되는 Input(입력)/Ouput(출력)이 포함

: 탭, 메뉴 등 앱 내부의 각 요소들을 정의

2. Server

: 샤이니 앱의 서버 로직을 정의

: 입력에서 출력까지 생성하는 어플리케이션의 전반적인 과정 포함

: 샤이니 앱을 로딩할 때 서버 기능이 호출됨

#Reference

1) https://www.edureka.co/blog/r-shiny-tutorial/

2) https://deanattali.com/blog/building-shiny-apps-tutorial/

[R Shiny] 샤이니 소개 End

BioinformaticsAndMe

Hierarchical clustering (R 계층적 군집화) Start

BioinformaticsAndMe

Hierarchical clustering (R 계층적 군집화)

: 주어진 데이터에서 모든 두 군집간의 거리를 계산하는 알고리즘 (데이터의 이진트리 구성)

ㄱ) 가장 가까운 거리에 있는 데이터를 서로 묶음 (반복 수행)

ㄴ) 최종적으로 하나의 클러스터로 합쳐질 때까지 진행

ㄷ) Dendrogram 형태의 자연적인 계층 구조로 정렬

# 1~150 사이에서 랜덤하게 40개 추출
idx <- sample(1:nrow(iris), 40)   
     
# 랜덤하게 뽑혀진 40개 숫자 인덱스에 대응하는 iris 행 추출
irisSample <- iris[idx,]

# Species 칼럼 제거
irisSample$Species <- NULL

# dist 함수는 데이터 간 거리 계산 (거리가 멀수록 차이가 큼)
hc <- hclust(dist(irisSample), method="average")
Call:
hclust(d = dist(irisSample), method = "average")

Cluster method   : average 
Distance         : euclidean 
Number of objects: 40 

# 'hang=-1' 은 모든 데이터의 y축(height)을 0부터 시작하게 함 
plot(hc, hang=-1, labels=iris$Species[idx])

# 'k=3' 은 3개의 클러스터 생성을 의미
rect.hclust(hc, k=3)

groups <- cutree(hc, k=3) 
groups
 98  73  34  68 133  79  95 150 136  71 123  36 118  42  38  49 117  97  19 129 120  53 149  75 115  89 
  1   1   2   1   1   1   1   1   3   1   3   2   3   2   2   2   1   1   2   1   1   1   1   1   1   1 
143  47  18 147  55  57  10 111 134  54  17  48 100 102 
  1   2   2   1   1   1   2   1   1   1   2   2   1   1

#Reference

1) https://rstudio-pubs-static.s3.amazonaws.com/249084_09c0daf4ceb24212a81ceddca97ba1ea.html

2) http://ropatics.com/data-mining/r-and-data-mining/RDM-Clustering.html

3) https://www.brandidea.com/hierarchicalclustering.html

4) http://blog.naver.com/PostView.nhn?blogId=wjddudwo209&logNo=220046493579&parentCategoryNo=&categoryNo=42&viewDate=&isShowPopularPosts=true&from=search

5) http://girke.bioinformatics.ucr.edu/GEN242/pages/mydoc/Rclustering.html

6) https://bcho.tistory.com/1204

Hierarchical clustering (R 계층적 군집화) End

BioinformaticsAndMe