R

[R] Neuroimaging Analysis 2021.02.01
[R] 국민건강영양조사 분석 2021.02.01 1
[R] Function (사용자 정의 함수) 2020.02.03
[R] Logistic regression (로지스틱 회귀분석) 2020.01.21
[R] Multiple linear regression (다중회귀분석) 2020.01.14 1
[R] ChIP-seq 분석 2020.01.05 1
[R] Circos plot 2019.12.30
[R] ggplot2 2019.12.16
[R] 상자그림(Box plot) 2019.12.10
[R] 파이차트(Pie plot) 2019.12.03

[R] Neuroimaging Analysis

2021. 2. 1. 19:39

[R] Neuroimaging Analysis

Start

BioinformaticsAndMe

Neuroconductor

: Neuroconductor는 재현 가능한 이미징 소프트웨어의 오픈 소스 플랫폼
: 이미지 분석 전용 R 소프트웨어의 중앙 저장소를 제공
: 프로그래밍 언어 R을 기반

: 시각화, 데이터 처리 및 저장, 통계적 추론을 포함한 이미징 영역 관련 수많은 패키지들을 보유

: Neuroconductor는 공식적인 테스트를 거쳐 새로운 R 패키지 submit을 허용함

1. Packages download

oro.nifti - reading/writing NIfTI images

*NIfTI(Neuroimaging Informatics Technology Initiative): 자기 공명 영상 기법을 사용하여 얻은 뇌 영상 데이터 형식
neurobase - extends oro.nifti and provides helpful imaging functions

install.packages("oro.nifti")
install.packages("neurobase")
library(oro.nifti)
library(neurobase)

2. Display neuroimaging

연습용 데이터 다운로드

ortho2 - displays nifti objects in 3 different planes (sagittal, coronal, transverse)

t1 = neurobase::readnii("./Raw_data/training01_01_t1.nii.gz")
neurobase::ortho2(t1)

이미지를 밝게 설정

ortho2(robust_window(t1, probs = c(0, 0.975)))

3. Display multiple neuroimaging

double_ortho - 동일한 크기의 이미지 2 개를 나란히 표시

double_ortho(t1, y = t1 > 300, col.y = "white")

4. Display specific neuroimaging

slice - 개별 슬라이스를 플롯팅 (아래는 sagittal까지 지정)

oro.nifti::slice(t1, z = 125, plane = "sagittal")

5. Function 정리

ortho2 - show orthographic images (and with overlays)
image - shows multiple slices of an image
slice - shows only specified slices
slice_overlay - similar to image but with an overlay
double_ortho - similar to ortho2 but side-by-side
robust_window - good for setting high values to not so high

#Reference

1) http://www.kjg.or.kr/journal/view.html?uid=3528&vmd=Full neuroconductor.org/

2) johnmuschelli.com/imaging_in_r/

[R] Neuroimaging Analysis

End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] 국민건강영양조사 분석 (1)	2021.02.01
[R] Function (사용자 정의 함수) (0)	2020.02.03
[R] Logistic regression (로지스틱 회귀분석) (0)	2020.01.21
[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] ChIP-seq 분석 (1)	2020.01.05

[R] 국민건강영양조사 분석

2021. 2. 1. 12:31

[R] 국민건강영양조사 분석 Start

BioinformaticsAndMe

국민건강영양조사

: 음주, 영양, 만성질환 등 500여 개 보건지표를 산출하는 국가 건강통계조사로 1998년에 도입하여 매년 1만여 명을 대상으로 실시

: 국민건강영양조사 18년도 원시자료 (v. 2020-09-01)

참여대상(row): 7,992명

변수(column): 785개

: 아래 홈페이지에서 기본DB의 SAS 다운로드

https://knhanes.cdc.go.kr/knhanes/sub03/sub03_02_02.do

1. Data loading

# 국민건강영양조사 18년도 원시자료 (v. 2020-09-01) 로딩
install.packages('data.table') ; install.packages('sas7bdat')
library(data.table)
library(sas7bdat)
sasr <- read.sas7bdat("./sas_class/hn18_all.sas7bdat")
data_sasr <- data.table(sasr)
dim(data_sasr)
[1] 7992  785

2. Data overview

# Gender (Pie chart)
mytable <- table(data_sasr$sex)
names(mytable) <- c('Man', 'Woman')
lbls <- paste(names(mytable), " : ", mytable, sep="")
pie(mytable, labels = lbls, col=c("#56B4E9", "#E69F00"), main="Pie Chart of Gender\n (with sample sizes)")

# Age (Histogram)
hist(data_sasr$age, freq=T, col="gray", xlab="Age")

3. Data filtering

: 흡연량과 폐기능검사 사이의 간단한 분석을 위해 데이터 필터링

: 6개변수 selection

가. BS3_2 (하루평균 흡연량)

#888 : 비해당(문항3-③⑧)

#999 : 모름, 무응답

나. HE_fev1fvc (1초간 노력성 호기량 / 노력성 폐활량)

다. sex (성별)

#1 : 남자

#2 : 여자

라. age (나이)

#1~79 : 1~79세

#80 : 80세이상

마. HE_BMI (체질량지수)

바. HE_COPD (폐기능 판정결과)

#1 : 폐기능정상

#2 : 제한성환기장애

#3 : 폐쇄성환기장애

#9 : 판정불능

# Data filtering
New_data_sasr <- data_sasr[, c('BS3_2','HE_fev1fvc', 'sex', 'age', 'HE_BMI', 'HE_COPD')]
Parsed_data_sasr <- New_data_sasr[ BS3_2 != 999 & BS3_2 != 888 & is.na(HE_fev1fvc) != T,] 
dim(Parsed_data_sasr)
[1] 644   6

#지표설명

FVC (Forced Vital Capacity; 노력성 폐활량)

-최대로 숨을 들이쉰 다음, 최대 노력으로 끝까지 내쉬었을 때 공기량

FEV1 (Forced Expiratory Volume in One second; 1초간 노력성 호기량)

-첫 1초간 얼마나 빨리 숨을 내쉴 수 있는지 보는 지표

FEV1/FVC ratio

-기관지폐쇄 유무를 확인하는 지표

-0.7 이하인 경우 숨을 내쉬는데 장애, 즉 기도폐쇄가 있음을 의미

4. Correlation test

# Correlation test (Raw data)

cor.test(Parsed_data_sasr$BS3_2, Parsed_data_sasr$HE_fev1fvc, method='pearson')

Pearson's product-moment correlation

data:  Parsed_data_sasr$BS3_2 and Parsed_data_sasr$HE_fev1fvc
t = -2.8655, df = 642, p-value = 0.004299
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.18800432 -0.03542489
sample estimates:
      cor 
-0.112377

5. 계층화 추출법 후, correlation test

# Plotting
Temp_bin_frame <- Parsed_data_sasr[,Bin_group:= cut(BS3_2, c(-Inf, 5, 10, 15, 20, 25, 30, Inf), paste0(1:7))]
boxplot(Temp_bin_frame$HE_fev1fvc ~ Temp_bin_frame$Bin_group, outline=F, las=1, xlab='', ylab='')#, varwidth=T )
Bin_frame <- aggregate(Temp_bin_frame, by=list(Temp_bin_frame$Bin_group), mean)
plot(Bin_frame$BS3_2, Bin_frame$HE_fev1fvc, pch=19, las=1, xlab='', ylab='')

# Correlation test (After filtering)

cor.test(Bin_frame$BS3_2, Bin_frame$HE_fev1fvc, method='pearson')

Pearson's product-moment correlation

data:  Bin_frame$BS3_2 and Bin_frame$HE_fev1fvc
t = -3.3085, df = 5, p-value = 0.02127
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9739242 -0.2006617
sample estimates:
       cor 
-0.8285219

#Reference

1) https://www.khidi.or.kr/nutristat

2) https://github.com/seung-00/classification_obesity_risk_groups

3) https://m.blog.naver.com/PostView.nhn?blogId=i-doctor&logNo=220605517140&proxyReferer=https:%2F%2Fwww.google.com%2F

4) https://rsas.tistory.com/250

[R] 국민건강영양조사 분석 End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Neuroimaging Analysis (0)	2021.02.01
[R] Function (사용자 정의 함수) (0)	2020.02.03
[R] Logistic regression (로지스틱 회귀분석) (0)	2020.01.21
[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] ChIP-seq 분석 (1)	2020.01.05

[R] Function (사용자 정의 함수)

2020. 2. 3. 15:00

[R] Function (사용자 정의 함수) Start

BioinformaticsAndMe

R function

: 프로그래밍 언어에서 함수(Function)는 반복적으로 사용될 수 있는 문장 블록의 형태

: R은 여러 내장 함수(Built-In function)을 제공하며, 사용자가 직접 자신의 함수(User defined function)를 정의할 수 있음

1. 함수 생성 및 실행

# function 커맨드를 사용하여, 새로운 함수(myfunc)를 생성
myfunc <- function() {
print('Hello, World!')
}

# 정의된 함수(myfunc)를 실행
myfunc()
[1] "Hello, World!"

2. 함수에 특정 인수값 전달하기

: 사용자는 인수(Argument)를 통해, 특정 값을 함수에 전달할 수 있음

: 인수는 function 커맨드 다음에 있는 괄호 안에서 선언됨

*쉼표로 구분하여 지정된 개수만큼 인수 전달 가능

# x, y 인수를 갖는 함수 생성
sum <- function(x, y) {
    x + y
}

# sum 함수에 특정 인수값을 전달하여 실행
sum(2, 3)
[1] 5

3. 기본값(Default value) 지정하기

: 함수에서 인수는 기본값을 지정할 수 있음

: 사용자가 인수 입력없이 함수를 호출하면 기본값이 사용됨

# 인수 y에 기본값 3을 지정
pow <- function(x, y=3) {
    x ^ y
}

# y 인수값이 입력되지 않았으므로, 기본값 'y=3'이 사용됨
pow(2)
[1] 8

# y 인수값이 입력됐으므로, 'y=4'이 사용됨
pow(2, 4)[1] 16

4. 함수값 반환하기

: 함수에서 특정값을 반환하기 위해서, return 커맨드를 사용

# return 커맨드를 사용하여, 함수값 반환
sum <- function(x, y) {
    return(x + y)
}

sum(2, 3)
[1] 5

# 결과를 벡터(or 리스트)에 저장하여, 여러 값을 동시에 반환할 수 있음
math <- function(x, y) {
    add <- x + y
    sub <- x - y
    mul <- x * y
    div <- x / y
    c(addition = add, subtraction = sub, multiplication = mul, division = div)
}

math(6, 3)
      addition    subtraction multiplication       division 
             9              3             18              2

5. 인수의 개수가 가변적인 상황

: 함수 생성시 생략부호(...)를 사용하여, 개수가 가변적인 인수를 수용함

: 사용자가 함수에 전달할 인수의 개수를 미리 알지 못하는 상황에서 유용함

# 아래 myfunc 함수는 생략부호(...)를 사용하여 , 모든 인수를 수용할 수 있음
myfunc <- function(x, ...) {
print(x)
summary(...)
}
v <- 1:10

myfunc("Summary of v:", v)
[1] "Summary of v:"
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    3.25    5.50    5.50    7.75   10.00

#Reference

1) https://www.learnbyexample.org/r-functions/

[R] Function (사용자 정의 함수) End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Neuroimaging Analysis (0)	2021.02.01
[R] 국민건강영양조사 분석 (1)	2021.02.01
[R] Logistic regression (로지스틱 회귀분석) (0)	2020.01.21
[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] ChIP-seq 분석 (1)	2020.01.05

[R] Logistic regression (로지스틱 회귀분석)

2020. 1. 21. 14:42

[R] Logistic regression (로지스틱 회귀분석) Start

BioinformaticsAndMe

Logistic regression (로지스틱 회귀분석)

: 로지스턱 회귀분석은 종속변수(반응변수)가 범주형 데이터인 경우에 사용되는 회귀 분석법

: 종속변수 y는 '성공(1) 및 실패(0)'의 두 가지 값(이항변수)을 갖음

*환자사망여부/전염병발병여부/교통사고발생여부 등

: 로지스티 회귀분석은 지도 학습으로 분류되며, 특정 결과의 분류 및 예측을 위해 사용됨

일반화선형모형 (Generalized linear model)

: 일반화선형모형은 정규분포를 따르지 않는 종속변수의 선형 모형 확장으로, 로지스틱회귀 또는 포아송회귀에 사용됨

: 일반화선형모형은 R의 내장함수인 glm()함수를 사용

: 로지스틱 회귀분석에서는 glm()함수에 'family=binomial' 인수를 지정해야함

1. 실습 대장암 데이터 로딩

# survival 패키지의 1858명 colon 데이터

install.packages(“survival”) library(survival) str(colon)

'data.frame':   1858 obs. of  16 variables:
 $ id      : num  1 1 2 2 3 3 4 4 5 5 ...
 $ study   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ rx      : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
 $ sex     : num  1 1 1 1 0 0 0 0 1 1 ...
 $ age     : num  43 43 63 63 71 71 66 66 69 69 ...
 $ obstruct: num  0 0 0 0 0 0 1 1 0 0 ...
 $ perfor  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ adhere  : num  0 0 0 0 1 1 0 0 0 0 ...
 $ nodes   : num  5 5 1 1 7 7 6 6 22 22 ...
 $ status  : num  1 1 0 0 1 1 1 1 1 1 ...
 $ differ  : num  2 2 2 2 2 2 2 2 2 2 ...
 $ extent  : num  3 3 3 3 2 2 3 3 3 3 ...
 $ surg    : num  0 0 0 0 0 0 1 1 1 1 ...
 $ node4   : num  1 1 0 0 1 1 1 1 1 1 ...
 $ time    : num  1521 968 3087 3087 963 ...
 $ etype   : num  2 1 2 1 2 1 2 1 2 1 ...

2. 로지스틱 회귀분석 수행

■반응변수 - status(대장암 재발 또는 사망인 경우 1)

■예측변수

- obstruct : 종양에 의한 장의 폐쇄 (obstruction)

- perfor : 장의 천공 (perforation)

- adhere : 인접장기와의 유착 (adherence)

- nodes : 암세포가 확인된 림프절 수

- differ : 암세포의 조직학적 분화 정도 (1=well, 2=moderate, 3=poor)

- extent : 암세포가 침습한 깊이 (1=submucosa, 2=muscle, 3=serosa, 4=인접장기)

- surg : 수술 후 등록까지의 시간 (0=short, 1=long)

# 로지스틱 회귀분석에서 'family=binomial'로 지정

colon1<-na.omit(colon) result<-glm(status ~ sex+age+obstruct+perfor+adhere+nodes+differ+extent+surg, family=binomial, data=colon1) summary(result)

Call:
glm(formula = status ~ rx + sex + age + obstruct + perfor + adhere + 
    nodes + differ + extent + surg, family = binomial, data = colon1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.575  -1.046  -0.584   1.119   2.070  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.430926   0.478301  -5.082 3.73e-07 ***
rxLev       -0.069553   0.122490  -0.568 0.570156    
rxLev+5FU   -0.585606   0.124579  -4.701 2.59e-06 ***
sex         -0.086161   0.101614  -0.848 0.396481    
age          0.001896   0.004322   0.439 0.660933    
obstruct     0.219995   0.128234   1.716 0.086240 .  
perfor       0.085831   0.298339   0.288 0.773578    
adhere       0.373527   0.147164   2.538 0.011144 *  
nodes        0.185245   0.018873   9.815  < 2e-16 ***
differ       0.031839   0.100757   0.316 0.752003    
extent       0.563617   0.116837   4.824 1.41e-06 ***
surg         0.388068   0.113840   3.409 0.000652 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2461.7  on 1775  degrees of freedom
Residual deviance: 2240.4  on 1764  degrees of freedom
AIC: 2264.4

Number of Fisher Scoring iterations: 4

3. 유의한 변수 선택

: backward elimination방법으로 stepwise logistic regression 수행

*backward elimination 참고 - https://bioinformaticsandme.tistory.com/290

# 유의하지 않은 변수를 누락하고 로지스틱 회귀모형을 새롭게 정의

reduced.model=step(result, direction = "backward") summary(reduced.model)

Call:
glm(formula = status ~ rx + obstruct + adhere + nodes + extent + 
    surg, family = binomial, data = colon1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5583  -1.0490  -0.5884   1.1213   2.0393  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.30406    0.35138  -6.557 5.49e-11 ***
rxLev       -0.07214    0.12221  -0.590 0.554978    
rxLev+5FU   -0.57807    0.12428  -4.651 3.30e-06 ***
obstruct     0.22148    0.12700   1.744 0.081179 .  
adhere       0.38929    0.14498   2.685 0.007251 ** 
nodes        0.18556    0.01850  10.030  < 2e-16 ***
extent       0.56510    0.11643   4.854 1.21e-06 ***
surg         0.38989    0.11371   3.429 0.000606 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2461.7  on 1775  degrees of freedom
Residual deviance: 2241.5  on 1768  degrees of freedom
AIC: 2257.5

Number of Fisher Scoring iterations: 4

4. 예측 인자들의 Odds ratio 구하기

: 예측 변수들의 오즈비 계산

# 오즈비 출력 함수 정의
ORtable=function(x,digits=2){
    suppressMessages(a<-confint(x))
    result=data.frame(exp(coef(x)),exp(a))
    result=round(result,digits)
    result=cbind(result,round(summary(x)$coefficient[,4],3))
    colnames(result)=c("OR","2.5%","97.5%","p")
    result
}

ORtable(reduced.model)
              OR 2.5% 97.5%     p
(Intercept) 0.10 0.05  0.20 0.000
rxLev       0.93 0.73  1.18 0.555
rxLev+5FU   0.56 0.44  0.72 0.000
obstruct    1.25 0.97  1.60 0.081
adhere      1.48 1.11  1.96 0.007
nodes       1.20 1.16  1.25 0.000
extent      1.76 1.41  2.22 0.000
surg        1.48 1.18  1.85 0.001

# Odds ratio 시각화
install.packages(“moonBook”)
library(moonBook)
odds_ratio = ORtable(reduced.model)
odds_ratio = odds_ratio[2:nrow(odds_ratio),]
HRplot(odds_ratio, type=2, show.CI=TRUE, cex=2)

#Reference

1) https://www.tech-quantum.com/classification-logistic-regression/

2) https://rstudio-pubs-static.s3.amazonaws.com/41074_62aa52bdc9ff48a2ba3fb0f468e19118.html

3) http://www.dodomira.com/2016/02/12/logistic-regression-in-r/

4) https://link.springer.com/chapter/10.1007/978-1-4842-4470-8_20

[R] Logistic regression (로지스틱 회귀분석) End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] 국민건강영양조사 분석 (1)	2021.02.01
[R] Function (사용자 정의 함수) (0)	2020.02.03
[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] ChIP-seq 분석 (1)	2020.01.05
[R] Circos plot (0)	2019.12.30

[R] Multiple linear regression (다중회귀분석)

2020. 1. 14. 16:32

[R] Multiple linear regression (다중회귀분석) Start

BioinformaticsAndMe

Multiple linear regression (다중회귀분석)

: 다중회귀분석은 예측변수(독립변수;설명변수)가 2개 이상인 회귀분석

: 다중회귀분석에서 예측변수 개수가 많을 경우, 적절한 회귀모형 선택이 필요함

: 회귀모형에 포함되는 예측변수의 선정 기준

ㄱ) 반응변수(종속변수)와 높은 상관관계

ㄴ) 선택된 예측변수들은 서로 낮은 상관관계를 보임 (다중공선성 문제를 회피)

ㄷ) 예측변수의 개수는 적을수록 유리

Feature selection (변수 선택법)

1) All possible regressions

- 변수들의 가능한 모든 조합들로부터 최적의 모형을 찾아냄

- 유의한 변수가 누락되지 않는 안전한 방법

- 변수가 많을수록 탐색 시간이 급증함

2) Forward stepwise selection (Forward selection)

- 기여도가 높은 유의한 변수부터 하나씩 추가하는 기법

- 빠른 계산이 장점

- 이미 선택된 변수는 다시 제거되지 않음

3) Backward stepwise selection (Backward elimination)

- 모든 변수를 포함한 상태에서 불필요한 변수를 제거해나가는 방법

- 중요한 변수가 제외될 가능성이 매우 적음

- 이미 제외된 변수는 다시 선택되지 않음

1. 실습데이터 로딩

# MASS 패키지의 birthwt 데이터

library(MASS) str(birthwt)

'data.frame':   189 obs. of  10 variables:
 $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
 $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
 $ race : int  2 3 1 1 1 3 1 3 1 1 ...
 $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
 $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
 $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
 $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

2. 다중회귀분석 수행

■반응변수 - bwt(출생시 체중)

■예측변수 - age(엄마나이), lwt(엄마몸무게), race(인종), smoke(임신중흡연상태), ptl(과거조산횟수), ht(고혈압기왕력), ui(uterine irritability 존재여부)

# race는 연속형 변수로 인식하므로, factor(race)로 사용해야함
model_1=lm(bwt ~ age+lwt+factor(race)+smoke+ptl+ht+ui, data=birthwt)
summary(model_1)
Call:
lm(formula = bwt ~ age + lwt + factor(race) + smoke + ptl + ht + 
    ui, data = birthwt)

Residuals:
    Min      1Q  Median      3Q     Max 
-1838.7  -454.5    57.6   465.1  1711.0 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2934.132    311.450   9.421  < 2e-16 ***
age             -4.093      9.440  -0.434 0.665091    
lwt              4.301      1.722   2.497 0.013416 *  
factor(race)2 -488.196    149.604  -3.263 0.001318 ** 
factor(race)3 -353.334    114.319  -3.091 0.002314 ** 
smoke         -351.314    106.180  -3.309 0.001132 ** 
ptl            -47.423    101.663  -0.466 0.641443    
ht            -586.836    200.841  -2.922 0.003925 ** 
ui            -514.937    138.483  -3.718 0.000268 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 648.7 on 180 degrees of freedom
Multiple R-squared:  0.2424,	Adjusted R-squared:  0.2087 
F-statistic: 7.197 on 8 and 180 DF,  p-value: 2.908e-08

3. 유의하지 않은 변수 제거

: 위 다중회귀분석 결과에서 age와 ptl 변수가 유의하지 않은 것으로 확인됨

: 아래의 3가지 방법 중 하나로, 유의하지 않은 변수를 제거할 수 있음

1) 유의하지 않은 변수를 누락하고 회귀모형을 새롭게 정의
model_2 = lm(bwt~ lwt+factor(race)+smoke+ht+ui, data=birthwt)

2) update 함수를 사용하여, 기존 회귀모형에서 유의하지 않은 변수 제거
model_2 = update(model_1, .~. -age-ptl)

3) step 함수를 사용하여, 기존 회귀모형에서 유의하지 않은 변수를 제거해나감
model_2 = step(model_1, direction = "backward")

summary(model_2)
Call:
lm(formula = bwt ~ lwt + factor(race) + smoke + ht + ui, data = birthwt)

Residuals:
     Min       1Q   Median       3Q      Max 
-1842.14  -433.19    67.09   459.21  1631.03 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2837.264    243.676  11.644  < 2e-16 ***
lwt              4.242      1.675   2.532 0.012198 *  
factor(race)2 -475.058    145.603  -3.263 0.001318 ** 
factor(race)3 -348.150    112.361  -3.099 0.002254 ** 
smoke         -356.321    103.444  -3.445 0.000710 ***
ht            -585.193    199.644  -2.931 0.003810 ** 
ui            -525.524    134.675  -3.902 0.000134 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 645.9 on 182 degrees of freedom
Multiple R-squared:  0.2404,	Adjusted R-squared:  0.2154 
F-statistic:   9.6 on 6 and 182 DF,  p-value: 3.601e-09

4. 회귀모델 비교 및 평가

: F-test 분산분석으로 두 회귀모형의 설명력을 비교하여, 첫 번째 회귀모형에서 제거된 변수들의 기여도를 평가함

# anova( 제거 전 모형, 제거 후 모형 )
anova(model_1, model_2)Analysis of Variance Table

Model 1: bwt ~ age + lwt + factor(race) + smoke + ptl + ht + ui
Model 2: bwt ~ lwt + factor(race) + smoke + ht + ui
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    180 75741025                           
2    182 75937505 -2   -196480 0.2335  0.792

: F-test 결과 p-value가 0.792로 매우 크므로, 앞서 제거된 두 변수는 회귀모형에 대한 기여도가 적음을 알 수 있음

# summary 함수로 최종 회귀모형을 평가

summary(model_2)

Call:
lm(formula = bwt ~ lwt + factor(race) + smoke + ht + ui, data = birthwt)

Residuals:
     Min       1Q   Median       3Q      Max 
-1842.14  -433.19    67.09   459.21  1631.03 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2837.264    243.676  11.644  < 2e-16 ***
lwt              4.242      1.675   2.532 0.012198 *  
factor(race)2 -475.058    145.603  -3.263 0.001318 ** 
factor(race)3 -348.150    112.361  -3.099 0.002254 ** 
smoke         -356.321    103.444  -3.445 0.000710 ***
ht            -585.193    199.644  -2.931 0.003810 ** 
ui            -525.524    134.675  -3.902 0.000134 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 645.9 on 182 degrees of freedom
Multiple R-squared:  0.2404,	Adjusted R-squared:  0.2154 
F-statistic:   9.6 on 6 and 182 DF,  p-value: 3.601e-09

: R² 값이 0.2404 → 최종 회귀모형이 예측변수들의 24.04% 설명함

5. 다중공선성(Multicollinearity) 확인

: 다중공선성은 분산팽창지수(Variation Inflation Factor;VIF)라는 통계량을 사용하여 계산 가능

: 한 예측변수에 대해 VIF의 제곱근은 다중공선성의 정도를 나타냄

: VIF가 10을 넘으면(= GVIF^(1/(2*Df)) 값이 2를 넘으면), 다중공선성 문제를 보임

# VIF값은 car 패키지의 vif() 함수로 계산
install.packages('car')
library(car)
vif(model_2)
                 GVIF Df GVIF^(1/(2*Df))
lwt          1.182676  1        1.087509
factor(race) 1.257687  2        1.058993
smoke        1.154755  1        1.074595
ht           1.073548  1        1.036122
ui           1.036844  1        1.018255

: vif 함수로 완성된 회귀모형을 검사했을 때, 다중공선성 문제는 없는 것으로 확인됨

6. 변수의 상대적 중요성

: 현재까지 사용한 변수들 중에서, 결과 예측에 가장 중요한 변수를 서열화할 수 있음

# 아래 코드는 Dr. Johnson 논문에 소개된 스크립트로, 변수의 상대적 중요도를 출력함
relweights <- function(fit,...){
         R <- cor(fit$model)
         nvar <- ncol(R)
         rxx <- R[2:nvar, 2:nvar]
         rxy <- R[2:nvar, 1]
         svd <- eigen(rxx)
         evec <- svd$vectors
         ev <- svd$values
         delta <- diag(sqrt(ev))
         lambda <- evec %*% delta %*% t(evec)
         lambdasq <- lambda ^ 2
         beta <- solve(lambda) %*% rxy
         rsquare <- colSums(beta ^ 2)
         rawwgt <- lambdasq %*% beta ^ 2
         import <- (rawwgt / rsquare) * 100
         import <- as.data.frame(import)
         row.names(import) <- names(fit$model[2:nvar])
         names(import) <- "Weights"
         import <- import[order(import),1, drop=FALSE]
         dotchart(import$Weights, labels=row.names(import),
            xlab="% of R-Square", pch=19,
            main="Relative Importance of Predictor Variables",
            sub=paste("Total R-Square=", round(rsquare, digits=3)),
            ...)
         return(import)
}

# 변수의 상대적 중요도를 시각화
model_3 = lm(bwt~ lwt+race+smoke+ht+ui, data=birthwt)
result = relweights(model_3, col="blue")
result

# ggplot2을 사용하여, 변수의 상대적 중요도를 시각화
library(ggplot2)
plotRelWeights=function(fit){
  data<-relweights(fit)
  data$Predictors<-rownames(data)
  p<-ggplot(data=data,aes(x=reorder(Predictors,Weights),y=Weights,fill=Predictors))+ 
    geom_bar(stat="identity",width=0.5)+
    ggtitle("Relative Importance of Predictor Variables")+
    ylab(paste0("% of R-square \n(Total R-Square=",attr(data,"R-square"),")"))+
    geom_text(aes(y=Weights-0.1,label=paste(round(Weights,1),"%")),hjust=1)+
    guides(fill=FALSE)+
    coord_flip()
  p
}

plotRelWeights(model_3)

#Reference

1) http://rstudio-pubs-static.s3.amazonaws.com/189354_277dfb3a83a34a2abaae855b90fcf269.html

2) http://brokerstir.com/multiple-linear-regression-intuition/

3) https://quantifyinghealth.com/stepwise-selection/

4) https://slidesplayer.org/slide/11708770/

5) https://rpubs.com/cardiomoon/190997

6) http://leoslife.com/archives/3874

7) Johnson(2000, Mutilvariate BehavioralResearch, 35,1-19)

[R] Multiple linear regression (다중회귀분석) End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Function (사용자 정의 함수) (0)	2020.02.03
[R] Logistic regression (로지스틱 회귀분석) (0)	2020.01.21
[R] ChIP-seq 분석 (1)	2020.01.05
[R] Circos plot (0)	2019.12.30
[R] ggplot2 (0)	2019.12.16

[R] ChIP-seq 분석

2020. 1. 5. 15:51

[R] ChIP-seq 분석 Start

BioinformaticsAndMe

ChIP-seq

: ChIP-seq(ChIP-sequencing)은 Chromatin ImmunoPrecipitation sequencing의 약자로, DNA에 결합하는 단백질 분석에 사용되는 방법

: 특정 단백질과 결합된 DNA를 면역침강방법으로 분리하여 해당 서열을 확인

*Microarray로 서열 확인 → ChIP-chip

*NGS로 서열 확인 → ChIP-seq

: ChIP-seq은 전사인자 또는 염색질관련단백질들이 Phenotype(주로 발현)에 미치는 기작을 연구하기 위해 사용

: 'ChIPseeker'는 ChIP-seq 분석을 간단하게 수행할 수 있는 R 패키지로 높은 인용수를 보임

1. 패키지 설치 및 로딩

#ChIPseeker 설치
source("https://bioconductor.org/biocLite.R")
biocLite("ChIPseeker")
biocLite("clusterProfiler")

library(ChIPseeker)
library(clusterProfiler)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

2. ChIP profiling

1) 분석에 사용될 예제 샘플들 확인
files <- getSampleFiles()
peak <- readPeakFile(files[[4]])
peak

2) ChIP peaks coverage plot
covplot(peak, weightCol="V5")
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))

3) Transcription Start Site(TSS)의 ChIP peak 결과를 heatmap으로 보여주기
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

3. Peak annotation

1) ChIP-seq peak에 annotation 하기 (파이 차트)
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)

2) TSS에 결합하는 Transcription Factor(TF)의 상대적 위치 그리기
plotDistToTSS(peakAnno)

4. Functional enrichment analysis

1) Peak 부분에 대한 생물학적 패스웨이 분석 시각화
source("https://bioconductor.org/biocLite.R")
biocLite("ReactomePA")
library(ReactomePA)

pathway1 <- enrichPathway(as.data.frame(peakAnno)$geneId)
dotplot(pathway1)

5. ChIP peak data set comparison

1) 여러 데이터 셋의 ChIP peak 결과를 동시 비교 시각화
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)

2) 여러 데이터 셋의 ChIP peak annotation 결과를 동시 비교 시각화
peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb, tssRegion=c(-3000, 3000), verbose=FALSE)
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes) = sub("_", "\n", names(genes))
compKEGG <- compareCluster(geneCluster=genes, fun="enrichKEGG", pvalueCutoff=0.05, pAdjustMethod="BH")
plot(compKEGG, showCategory=15, title="KEGG Pathway Enrichment Analysis")

6. Data Mining with ChIP seq data

1) GEO(Gene Expression Omnibus)의 ChIP 데이터 수집
getGEOspecies()
getGEOgenomeVersion()
hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)

2) GEO ChIP 데이터 다운로드
downloadGEObedFiles(genome="hg19", destDir="hg19")
gsm <- hg19$gsm[sample(nrow(hg19), 10)]
downloadGSMbedFiles(gsm, destDir="hg19")

#Reference

1) https://en.wikipedia.org/wiki/ChIP-sequencing

2) https://www.insilicogen.com/blog/tag/CHIP-Seq

3) https://www.ncbi.nlm.nih.gov/pubmed/25765347

4) http://bioconductor.org/packages/devel/bioc/vignettes/ChIPseeker/inst/doc/ChIPseeker.html

[R] ChIP-seq 분석 End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Logistic regression (로지스틱 회귀분석) (0)	2020.01.21
[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] Circos plot (0)	2019.12.30
[R] ggplot2 (0)	2019.12.16
[R] 상자그림(Box plot) (0)	2019.12.10

[R] Circos plot

2019. 12. 30. 14:34

[R] Circos plot Start

BioinformaticsAndMe

Circos plot

: Circos plot은 보유한 데이터를 종합적으로 보여주기 위한 시각화 방법

: 데이터를 원형 레이아웃으로 시각화하여, 특정 정보 사이의 위치 관계를 탐색하기에 유리함

: R 패키지 'OmicCircos'를 통해, 오믹스 데이터 간의 관계를 Circos plot으로 나타냄

: http://bioconductor.org/packages/release/bioc/html/OmicCircos.html

1. 데이터 준비

1) 패키지 설치 및 데이터 로딩
options(stringsAsFactors = F)
source("https://bioconductor.org/biocLite.R")
biocLite("OmicCircos")
library("OmicCircos")

data("UCSC.hg19.chr")
data("TCGA.BC.gene.exp.2k.60")
dat <- UCSC.hg19.chr

dat$chrom <- gsub("chr", "",dat$chrom)
dat[dat[,1]==1,]

2) 시뮬레이션 데이터 생성
colors <- rainbow(10, alpha = 0.8)    #alpha : 투명도
lab.n <- 50
cnv.n <- 200
arc.n <- 30
fus.n <- 10

arc.d <- c()
for(i in 1:arc.n){
  chr <- sample(1:19, 1)
  chr.i <- which(dat$chrom == chr)
  chr.arc <- dat[chr.i,]
  arc.i <- sample(1:nrow(chr.arc), 2)
  arc.d <- rbind(arc.d, 
                 c(chr.arc[arc.i[1], c(1,2)], 
                   chr.arc[arc.i[2], c(2,4)]))
}
colnames(arc.d) <- c("chr", "start", "end", "value")

2-1) Gene funsion 시뮬레이션 데이터
fus.d <- c()
for(i in 1:fus.n){
  chr1 <- sample(1:19, 1)
  chr2 <- sample(1:19, 1)
  chr1.i <- which(dat$chrom == chr1)
  chr2.i <- which(dat$chrom == chr2)
  chr1.f <- dat[chr1.i,]
  chr2.f <- dat[chr2.i,]
  fus1.i <- sample(1:nrow(chr1.f), 1)
  fus2.i <- sample(1:nrow(chr2.f), 1)
  n1 <- paste0("geneA", i)
  n2 <- paste0("geneB", i)
  fus.d <- rbind(fus.d, c(
    chr1.f[fus1.i, c(1,2)], n1,
    chr2.f[fus2.i, c(1,2)], n2
  ))
}
colnames(fus.d) <- c("chr1","po1","gene1","chr2","po2","gene2")

2-2) Copy number 시뮬레이션 데이터
cnv.i <- sample(1:nrow(dat), cnv.n)
vale <- rnorm(cnv.n)
cnv.d <- data.frame(dat[cnv.i,c(1,2)], value=vale)

2-3) Genomic position 시뮬레이션 데이터
gene.pos <- TCGA.BC.gene.exp.2k.60[,1:3]

2-4) Gene expression 시뮬레이션 데이터
gene.exp <- TCGA.BC.gene.exp.2k.60

2-5) p-value 시뮬레이션 데이터
gene.pos$p <- rnorm(250,0.01,0.001)*
  sample(c(1,0.5,0.01,0.001,0.0001),250,replace=TRUE)

2. Circos plot 그리기

1) 염색체 그리기
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=T, print.chr.lab = T)    #R : 반지름,   W : 원형 너비

2) 같은 높이의 bar chart 그리기
par(mar=c(1,1,1,1))
plot(c(1,800),c(1,800),type="n",axes=FALSE,xlab="",ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=T, print.chr.lab = T)
circos(R=355, cir="hg19", type="b3", W=40, mapping=cnv.d, B=T, col=colors[7])    #b3 : 같은 높이의 형태로 bar chart 생성

3) 고정된 반지름을 점으로 표시
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=TRUE, print.chr.lab=TRUE)
circos(R=355, cir="hg19", type="b3", W=40, mapping=cnv.d, B=TRUE, col=colors[7])
circos(R=355, cir="hg19", type="s2", W=40, mapping=cnv.d, B=FALSE, col=colors[1], cex=0.5)    #s2 : 고정된 반지름 점 표시

4) 고정된 반지름으로 원호 그리기
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=TRUE, print.chr.lab=TRUE, cex=14)
circos(R=355, cir="hg19", type="b3", W=40, mapping=cnv.d, B=TRUE, col=colors[7])
circos(R=355, cir="hg19", type="s2", W=40, mapping=cnv.d, B=FALSE, col=colors[1], cex=0.5)
circos(R=320, cir="hg19", type="arc2", W=35, mapping=arc.d, B=TRUE, col=colors, lwd=10, cutoff=0)
#color는 matrix의 row 순서대로 들어감

5) cutoff value 반대편에 bar chart 그리기
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=TRUE, print.chr.lab=TRUE)
circos(R=355, cir="hg19", type="b3", W=40, mapping=cnv.d, B=TRUE, col=colors[7])
circos(R=355, cir="hg19", type="s2", W=40, mapping=cnv.d, B=FALSE, col=colors[1], cex=0.5)
circos(R=320, cir="hg19", type="arc2", W=35, mapping=arc.d, B=TRUE, col=colors, lwd=10) 
circos(R=280, cir="hg19", type="b2", W=40, mapping=cnv.d, B=TRUE, col=colors[c(5,9)], lwd=2, cutoff=-0.2, col.v=3)
#col.v : 값이 있는 칼럼 위치,    cutoff : +이상값제외/-이하값제외

6) type 파라미터를 계속해서 추가해주면서, Circos plot을 그려나감
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=TRUE, print.chr.lab=TRUE)
circos(R=355, cir="hg19", type="b3", W=40, mapping=cnv.d, B=TRUE, col=colors[7])
circos(R=355, cir="hg19", type="s2", W=40, mapping=cnv.d, B=FALSE, col=colors[1], cex=0.5)
circos(R=320, cir="hg19", type="arc2", W=35, mapping=arc.d, B=TRUE, col=colors, lwd=10, cutoff=0)
circos(R=280, cir="hg19", type="b2", W=40, mapping=cnv.d, B=TRUE, col=colors[c(7,9)], lwd=2, cutoff=-0.2, col.v=3)
circos(R=240, cir="hg19", type="arc", W=40, mapping=arc.d, B=TRUE, col=colors[c(1,7)], lwd=4, scale=TRUE, col.v=4)
circos(R=200, cir="hg19", type="box", W=40, mapping=cnv.d, B=TRUE, col=colors[1], lwd=0.1, scale=TRUE, col.v=3)
circos(R=160, cir="hg19", type="h", W=40, mapping=cnv.d, B=FALSE, col=colors[3], lwd=0.1, col.v=3)
circos(R=120, cir="hg19", type="link", W=10, mapping=fus.d, col=colors[c(1,7,9)], lwd=2)

7) 바깥쪽으로 Label을 두고 그리기
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=300, cir="hg19", type="chr", W=10, scale=FALSE, print.chr.lab=FALSE)
circos(R=310, cir="hg19", type="label", W=40, mapping=gene.pos, col=c("black","blue","red"), cex=0.4, side="out")
circos(R=250, cir="hg19", type="b3", W=40, mapping=cnv.d, B=TRUE, col=colors[7])
circos(R=250, cir="hg19", type="s2", W=40, mapping=cnv.d, B=FALSE, col=colors[1], cex=0.5)
circos(R=220, cir="hg19", type="arc2", W=30, mapping=arc.d, B=TRUE, col=colors, lwd=10, cutoff=0)
circos(R=190, cir="hg19", type="b2", W=30, mapping=cnv.d, B=TRUE, col=colors[c(7,9)], lwd=2, cutoff=-0.2, col.v=3)
circos(R=160, cir="hg19", type="arc", W=30, mapping=arc.d, B=TRUE, col=colors[c(1,7)], lwd=4, scale=TRUE, col.v=4)
circos(R=130, cir="hg19", type="box", W=30, mapping=cnv.d, B=TRUE, col=colors[1], lwd=0.1, scale=TRUE, col.v=3)
circos(R=100, cir="hg19", type="h", W=30, mapping=cnv.d,B=FALSE, col=colors[3], lwd=0.1, col.v=3)
circos(R=90, cir="hg19", type="link", mapping=fus.d, col=colors[c(1,7,9)], lwd=2)

6) Heatmap circos plot 그리기
par(mar=c(1,1,1,1))
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="")
circos(R=400, cir="hg19", type="chr", W=10, scale=FALSE, print.chr.lab=T)
circos(R=300, cir="hg19", type="heatmap2", W=100, mapping=gene.exp, col.v=4, cluster=F, col.bar=T, lwd=0.1, col="green")
circos(R=200, cir="hg19", type="s", W=100, mapping=gene.pos, col.v=4, col=colors, scale=TRUE, B=TRUE)
sig.gene <- gene.pos[gene.pos$p<0.000001,]
circos(R=190, cir="hg19", type="label", W=40, mapping=sig.gene, col=c("black","blue","red"), cex=0.4, side="in")

#Reference

1) http://circos.ca/

2) Hu, Ying, et al. "OmicCircos: a simple-to-use R package for the circular visualization of multidimensional omics data." Cancer informatics 13 (2014): CIN-S13495.

3) http://felixfan.github.io/circos/

4) http://bioconductor.org/packages/release/bioc/html/OmicCircos.html

[R] Circos plot End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Multiple linear regression (다중회귀분석) (1)	2020.01.14
[R] ChIP-seq 분석 (1)	2020.01.05
[R] ggplot2 (0)	2019.12.16
[R] 상자그림(Box plot) (0)	2019.12.10
[R] 파이차트(Pie plot) (0)	2019.12.03

[R] ggplot2

2019. 12. 16. 15:52

[R] ggplot2 Start

BioinformaticsAndMe

ggplot2

: ggplot2은 그래픽 문법에 기반한 R plotting 패키지로, 복잡한 그래픽을 쉽게 표현할 수 있는 강력한 툴

: 뉴질랜드 통계학자 Hadley Wickham이 개발

*Hadley는 강연의 참석자로부터 'R basic graphic에서 ggplot2으로 바꿔야 하나?'라는 질문을 받고,

'R basic graphic은 그림 그리기에 유리한 툴이지만, ggplot2은 그려진 데이터를 쉽게 이해할 수 있는 훌륭한 시각화 툴이다'라고 답함

: ggplot2 패키지는 R뿐만 아니라, Python에서도 plotnine 패키지를 통해 ggplot2 사용 가능

: ggplot2의 시각화 작업은 손으로 직접 그래프를 그리는 과정과 흡사

1) ggplot() 메소드로 축을 그림

2) geom_boxplot()/geom_bar() 등의 메소드로 그래프를 그림

3) xlabs()/ylabs() 등의 메소드로 기타 옵션을 조절

R ggplot2 예제

1. 산점도(Scatterplot)
# install.packages("ggplot2")                #ggplot2가 설치되지 않았다면, '#'을 지우고 설치
options(scipen=999)                           # 1e+48'과 같은 과학적 기수법을 사용하지 않음
library(ggplot2)                             
theme_set(theme_bw())                       #bw theme를 사전 설정
data("midwest", package = "ggplot2")    #'midwest' 데이터 로딩

#ggplot2 시작 (Scatterplot)
gg <- ggplot(midwest, aes(x=area, y=poptotal)) + 
  geom_point(aes(col=state, size=popdensity)) + 
  geom_smooth(method="loess", se=F) + 
  xlim(c(0, 0.1)) + 
  ylim(c(0, 500000)) + 
  labs(subtitle="Area Vs Population", 
       y="Population", 
       x="Area", 
       title="Scatterplot", 
       caption = "Source: midwest")

plot(gg)

2. 분기막대(Diverging bars)
library(ggplot2)
theme_set(theme_bw())  

data("mtcars")                                        #'mtcars' 데이터 로딩
mtcars$`car name` <- rownames(mtcars)      #'car names'의 새로운컬럼 생성 
mtcars$mpg_z <- round((mtcars$mpg - mean(mtcars$mpg))/sd(mtcars$mpg), 2)        #표준화 작업
mtcars$mpg_type <- ifelse(mtcars$mpg_z < 0, "below", "above")
mtcars <- mtcars[order(mtcars$mpg_z), ]  # sort
mtcars$`car name` <- factor(mtcars$`car name`, levels = mtcars$`car name`)               #플로팅에 나타날 정렬화 작업

#ggplot2 시작 (Diverging bars)
ggplot(mtcars, aes(x=`car name`, y=mpg_z, label=mpg_z)) + 
  geom_bar(stat='identity', aes(fill=mpg_type), width=.5)  +
  scale_fill_manual(name="Mileage", 
                    labels = c("Above Average", "Below Average"), 
                    values = c("above"="#00ba38", "below"="#f8766d")) + 
  labs(subtitle="Normalised mileage from 'mtcars'", 
       title= "Diverging Bars") + 
  coord_flip()

#아래 링크에서 ggplot2으로 그리는 50가지 종류의 R 소스코드 제공

http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

#Reference

1) http://freesearch.pe.kr/archives/3134

2) https://wikidocs.net/33913

3) http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

[R] ggplot2 End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] ChIP-seq 분석 (1)	2020.01.05
[R] Circos plot (0)	2019.12.30
[R] 상자그림(Box plot) (0)	2019.12.10
[R] 파이차트(Pie plot) (0)	2019.12.03
[R] 생존분석(Survival analysis) (0)	2019.11.25

[R] 상자그림(Box plot)

2019. 12. 10. 16:46

[R] 상자그림(Box plot) Start

BioinformaticsAndMe

1. 상자그림(Box plot)

: 상자그림은 특정한 수치 값을 기반으로 그려진, 자료 특성이 요약된 그래프

: 사분위수범위(Inter-Quartile Range;IQR) = Q1~Q3 = 상자그림몸통

: 중앙값(Median) = 상자그림몸통가운데선 = Q2(두번째사분위수)

: 최댓값(Maximum) = 위쪽 수염의 끝부분 (수염은 상자그림몸통 끝에 수직을 이룬 선)

: 최솟값(Minimum) = 아래쪽 수염의 끝부분

: 이상점(Outlier) = 최대최소를 벗어난 값 = 울타리(Fence) 바깥의 값 (Fence는 IQR x 1.5로 지정)

: 상자그림은 정규분포와 흡사

: 전체 분포 범위의 50%가 상자그림몸통에 포함됨

: 상자그림을 통해, 데이터의 분포를 확인하고 이상점을 처리하기 쉬움

2. R 상자그림 예제

# 뉴욕의 공기 퀄리티 예제 로딩

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 ...

# 오존값으로 박스플롯 그리기
boxplot(airquality$Ozone)

# 상자그림 파라미터 지정

boxplot(airquality$Ozone, main = "Mean ozone in parts per billion at Roosevelt Island", xlab = "Parts Per Billion", ylab = "Ozone", col = "orange", border = "brown", horizontal = TRUE, #수평 상자그림 notch = TRUE #Notch(파인 형태) 상자그림 )

# 다중 상자그림
boxplot(Temp~Month,
            data=airquality,
            main="Different boxplots for each month",
            xlab="Month Number",
            ylab="Degree Fahrenheit",
            col="orange",
            border="brown"
)

#아래 링크에서 'R 기본상자그림' 및 'ggplot2 상자그림'의 다양한 스크립트를 제공

https://www.r-graph-gallery.com/boxplot.html

#Reference

1) https://www.datamentor.io/r-programming/box-plot/

2) https://www.r-graph-gallery.com/boxplot.html

3) https://namu.wiki/w/%EC%83%81%EC%9E%90%20%EC%88%98%EC%97%BC%20%EA%B7%B8%EB%A6%BC

[R] 상자그림(Box plot) End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] Circos plot (0)	2019.12.30
[R] ggplot2 (0)	2019.12.16
[R] 파이차트(Pie plot) (0)	2019.12.03
[R] 생존분석(Survival analysis) (0)	2019.11.25
[R] 히스토그램(Histogram) (0)	2019.11.18

[R] 파이차트(Pie plot)

2019. 12. 3. 16:13

[R] 파이차트(Pie plot) Start

BioinformaticsAndMe

1. R pie chart

: R에는 파이 차트를 그릴 수 있는 '내장 함수 pie'가 존재

# 기본 pip 차트
apple <- c(10,20,30,40)
pct <- round(apple/sum(apple)*100, 1)        # 비율 값 지정
lab <- paste(pct,"%")
pie(apple, init.angle=90, col=rainbow(length(apple)), label=lab)        # init.angle 옵션은 파이 차트의 시작 각도 지정
legend("topright", ,c("1주","2주","3주","4주"), cex=0.8, fill=rainbow(length(apple)) )

2. ggplot2 파이 차트

: R에서 제공하는 기본 pie 함수는 시각적/기능적 한계를 가지므로, ggplot2 함수로 성능 확장

# 예제 데이터 생성
df <- data.frame(
  group = c("Male", "Female", "Child"),
  value = c(25, 25, 50)
  )
head(df)
   group value
1   Male    25
3  Child    50
2 Female    25

# ggplot2 로딩 후, barplot 그리기
library(ggplot2)
bp<- ggplot(df, aes(x="", y=value, fill=group))+
geom_bar(width = 1, stat = "identity")
bp

# Pie chart 그리기
pie <- bp + coord_polar("y", start=0)
pie

# Pie chart 색상 지정
pie + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))

# Pie chart 팔레트 색상 사용
pie + scale_fill_brewer(palette="Blues") + theme_minimal()

#Reference

1) http://www.sthda.com/english/wiki/ggplot2-pie-chart-quick-start-guide-r-software-and-data-visualization

[R] 파이차트(Pie plot) End

BioinformaticsAndMe

저작자표시

'R' 카테고리의 다른 글

[R] ggplot2 (0)	2019.12.16
[R] 상자그림(Box plot) (0)	2019.12.10
[R] 생존분석(Survival analysis) (0)	2019.11.25
[R] 히스토그램(Histogram) (0)	2019.11.18
[R] Excel 읽기/쓰기 (0)	2019.11.05

PREV 1 2 3 4 NEXT

BioinformaticsAndMe

R

[R] Neuroimaging Analysis

'R' 카테고리의 다른 글

[R] 국민건강영양조사 분석

'R' 카테고리의 다른 글

[R] Function (사용자 정의 함수)

'R' 카테고리의 다른 글

[R] Logistic regression (로지스틱 회귀분석)

'R' 카테고리의 다른 글

[R] Multiple linear regression (다중회귀분석)

'R' 카테고리의 다른 글

[R] ChIP-seq 분석

'R' 카테고리의 다른 글

[R] Circos plot

'R' 카테고리의 다른 글

[R] ggplot2

'R' 카테고리의 다른 글

[R] 상자그림(Box plot)

'R' 카테고리의 다른 글

[R] 파이차트(Pie plot)

'R' 카테고리의 다른 글

+ Recent posts

티스토리툴바