BioinformaticsAndMe

R apply 함수 Start.

BioinformaticsAndMe

# R apply() Function

- apply() : 데이터프레임 또는 매트릭스에 함수를 적용하여 행 또는 열 단위의 계산 수행.

apply( X, MARGIN, FUN)

- X : matrix

- MARGIN : 1(row) , 2(col)

- FUN : row/column 단위로 적용할 함수. 사용자가 만든 함수도 가능함.

EX1) apply(X, 1, sum)

# 각 특성(Column)별로 평균을 구한다

> apply(patients.sub, 2, mean, na.rm=TRUE)

weight height

124.72 165.00

# 각 셀에 2를 곱해준다

MulTwo <- function(x){ return(2*x) }

> apply(patients.sub, c(1,2), MulTwo)

     weight height

[1,]  130.8    340

[2,]  110.0    310

[3,]  760.0     NA

[4,]  144.4    346

#다양한 apply 함수군이 존재한다.

lapply: 결과를 리스트 형태로 반환

sapply: 벡터, 또는 행렬의 형태로 반환 (s: simplify)

tapply: 입력값을 index에 지정한 factor 값으로 분류(그룹화)하여 매개변수로 넘어온 function을 적용하는 함수다.

# lapply

R apply 함수 End.

BioinformaticsAndMe

R 상관분석 (R correlation test) Start.

BioinformaticsAndMe

# R을 이용한 상관분석(correlation test)을 시작해보자.

(상관분석에 대한 개념정리는 아래의 'Statistic' 카테고리에 있으니, 먼저 선행하고 오면 좋을 듯하다..!)

http://bioinformaticsandme.tistory.com/58?category=808983

# Pearson 상관분석은 변수들이 얼마나 직선적인 관계를 가지는지 분석하는 기법으로 상관계수를 이용하여 측정한다.

# 상관계수: Correlation coefficient

# 아이리스 예제데이터 불러오기

> attach(iris)

# cor() 함수 사용하여 상관계수 확인

> cor(Sepal.Length, Petal.Width)

[1] 0.8179411

# Pearson 상관계수: 0.8179

> cor.test(Sepal.Length, Petal.Width)

        Pearson's product-moment correlation

data:  Sepal.Length and Petal.Width

t = 17.2965, df = 148, p-value < 2.2e-16

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

 0.7568971 0.8648361

sample estimates:

      cor 

0.8179411 

# iris 데이터의 4가지 변수에 대해서 (Species를 제외한) 상관계수를 구해보자.

R 상관분석 (R correlation test) End.

BioinformaticsAndMe

막대그래프 (Barplot) Start.

BioinformaticsAndMe

# barplot( ) 함수를 이용해 막대 그래프를 그려보자.

.

- 위 파라미터들은 barplot 에 사용되는 옵션이니 하나씩 테스트해보면 이해가 빠르겠다.

# 기본 bar 그래프

> x <- c(1,2,3,4,5)

> barplot(x)

# 그래프 가로로 출력 (horiz=T)

> x <- c(1,2,3,4,5)

> barplot(x, horiz=T)

# 그룹으로 묶어서 가로로 출력(beside, horiz=T)

> x <- matrix(c(1,2,3,4),2,2)

> barplot(x, beside=T, horiz=T, names=1:2)

> barplot(x, horiz=T, names=1:2)

# 여러 가지 옵션을 추가하여 그래프 그리기

> apple <- c(100,120,140,130,150)

> barplot(apple, main= "Apple", xlab="요일", ylab="금액", names.arg=c("월","화","수","목","금"),  border="green", density=c(10,20,30,40,50))

# 여러 막대그래프를 그룹으로 묶어서 한꺼번에 출력

> apple <- c(100,120,140,130,150)

> peach <- c(180,200,210,190,170)

> berry <- c(110,150,160,90,120)

> fruits <- data.frame(APPLE=apple, PEACH=peach, BERRY=berry)

> fruits

  APPLE PEACH BERRY

1   100   180   110

2   120   200   150

3   140   210   160

4   130   190    90

5   150   170   120

> barplot(as.matrix(fruits), main="FRUITS", ylab="수량", beside=T, col=rainbow(5), ylim=c(0,400))

> legend("topright", c("월","화","수","목","금"), cex=0.8, fill=rainbow(5))

# 하나의 막대그래프에 여러가지 내용 출력하기

> barplot(t(fruits), main="FRUITS", ylab="판매량", ylim=c(0,900), col=rainbow(3), space=0.1, cex.axis=0.8, las=1, 

names.arg=c("월","화","수","목","금"), cex=0.8)

> legend(3.5, 900, names(fruits), cex=0.8, fill=rainbow(3))

# 조건을 주고 그래프 그리기

> peach <- c(180,200,210,190,170)

> colors <- c( )

> for( i in 1:length(peach)){

if (peach[i] >= 200) {colors <- c(colors, "red")}

else { colors <- c(colors, "blue")}

   }

> barplot(peach, main="PEACH", xlab="요일", ylab="수량", names.arg=c("월","화","수","목","금"), col=colors)

막대그래프 (Barplot) End.

BioinformaticsAndMe

R, 결측치 처리 (Missing value, NA) Start.

BioinformaticsAndMe

#R 작업시 발생하는 결측치 (missing value)를 다뤄보자.

-R 프로그래밍에서 결측지(missing value)는 NA (Not Available) 라는 문자로 처리해야 한다. NaN (Not a Number)는 분모를 0으로 나누는 것과 같이 계산이 불가능 할 경우 출력되는 문자다.

> y <- c(1,2,3, NA)

> y

[1]  1  2  3 NA

#is.na()는 벡터의 결측지가 존재할 경우 true

> is.na(y)

[1] FALSE FALSE FALSE  TRUE

> summary(y)

#특정 값을 NA로 바꾸기 (-999 -> NA )

> ages <- c(48, 78, 56, 88, -999, 13, 26, -999)

> ages[ ages == -999] <- NA

> ages

[1] 48 78 56 88 NA 13 26 NA

#결측지(missing value)가 하나라도 포함된 데이터가 존재할 경우 연산의 결과 역시 NA가 된다. 따라서 함수 역시 아래와 같이 NA가 결과로 나온다.

> sum(ages)

[1] NA

> mean(ages)

[1] NA

#NA 데이터를 제외하고 연산하고 싶을 경우  na.rm = TRUE 매개변수를 넣어주면 된다.

> sum(ages, na.rm = TRUE)

[1] 309

> mean(ages, na.rm=TRUE)

[1] 51.5

#결측치 (missing value) 실습 예제

> weight <- c(65.4, 55, 380, 72.2, 51, NA)

> height <- c(170, 155, NA, 173, 161, 166)

> gender <- c("M", "F","M","M","F","F")

> testDate <- c("2013/09/01", "2013/09/01", "2013/09/05", "2013/09/14", "2013/10/11", "2013/10/26")

> patients <- data.frame( weight = weight, height=height, gender=gender, testDate=testDate)

# complete.cases(patients) # NA가 존재하는 경우 FALSE를 리턴

[1]  TRUE  TRUE FALSE  TRUE  TRUE FALSE

#patients[complete.cases(patients),] # complete.cases 함수 리턴값이 참인 경우를 출력하는 방법

  weight height gender   testDate

1   65.4    170      M 2013/09/01

2   55.0    155      F 2013/09/01

4   72.2    173      M 2013/09/14

5   51.0    161      F 2013/10/11

R, 결측치 처리 (Missing value, NA) End.

BioinformaticsAndMe

R, Command line interface Ⅱ Start.

BioinformaticsAndMe

파트 1 에 이어서,

R 의 기본 명령어와 Component 를 다뤄보자.

8. 매트릭스에 row/column 추가하기

#Row 추가하기

> mat = matrix(1:20, ncol=4, nrow=5)

> rbind(mat, c(26:29) )

     [,1] [,2] [,3] [,4]

[1,]    1    6   11   16

[2,]    2    7   12   17

[3,]    3    8   13   18

[4,]    4    9   14   19

[5,]    5   10   15   20

[6,]   26   27   28   29

> colnames(CountTable) <- c("No heart attack", "Heart attack")

> CountTable

        No heart attack Heart attack

Placebo           189        10845

Aspirin             104        10933

> CountTable["Placebo",]

No heart attack    Heart attack 

            189           10845 

> colnames(CountTable)

[1] "No heart attack" "Heart attack"

13. 데이터프레임 예제

#변속기가 자동(am == 0)이고 & 실린더가 4개, 6개 (cyl == c(4, 6)) 인 자동차들의 연비(mpg) 평균(mean())는?

R, Command line interface Ⅱ End.

BioinformaticsAndMe

R, Command line interface Ⅰ Start.

BioinformaticsAndMe

R 의 기본 명령어와 Component 를 다뤄보자.

1. 변수 할당

-알파벳, 숫자, _(언더바), .(마침표)로 구성된다.

- '-' (하이픈)은 사용불가.

-첫글자는 알파벳 또는 .(마침표)로만 시작해야 한다.

2. Data Type & Structure

6. 벡터의 조건문

EX1) x <- c(11, 12, 13, 14, 19, 20)

R, Command line interface Ⅰ End.

BioinformaticsAndMe

R, RStudio 설치 Start.

BioinformaticsAndMe

오픈소스 통계 프로그래밍 언어 R 과

R 통합개발 분석 환경 (IDE, Integrated Development Environment)인 RStudio 설치를 해보자.

1. R 설치

R은 Windows, MAC OS X, Linux 환경 모두에서 설치 가능하다. 여기서는 Window를 기준으로 설치를 진행한다.

R을 설치하려면 아래 보이는 공식사이트(https://www.r-project.org/)에서 R소프트웨어를 다운로드 한다. 
현재 최신버전은 3.5.1 (2018년07월02일) 이다.

아래의 download R을 클릭하면 R을 설치할 수 있는 Mirror 사이트로 연결된다. 

RStudio는 코드 직접실행, 구문강조, 괄호 자동입력지원, 명령어 완성, 다양한 단축키, 데이터 보기 및 가져오기, 그래픽 조작, 프로젝트 관리, 버전 관리 등의 다양한 기능을 제공한다. RStudio를 설치하기 위해서 RStudio 홈페이지로 이동한다. 

https://www.rstudio.com/

공식홈페이지에서 Download RStudio를 누르자.

클릭해서 들어가면 처음에 무슨 달러 써있고 이래서 놀랄 수 있는데, 가볍게 지나치고

스크롤을 내리면 아래 같은 화면이 나오고 해당 OS에 맞게 설치해주면 끝난다.

모든 설치가 끝나면 4개의 영역으로 구성된 RStudio를 볼 수 있겠다.

1. 소스 편집기 및 데이터view

2. R 콘솔

3. 작업환경과 히스토리

4. 파일, 플롯, 패키지, 도움말

R, RStudio 설치 End.

BioinformaticsAndMe

Permutation test (순열검정법) Start.

BioinformaticsAndMe

Permutation test 는 t-test 등의 일반적인 통계 검정을 수행할 만큼 샘플의 수가 크지 않은 경우에 사용할 수 있는 검정 방법. 이 경우 주어진 샘플을 무작위로 추출하여 인공적으로 샘플 숫자를 늘림으로써 전체 모수를 통계 검정이 가능한 크기만큼 키운 다음, 원래 주어진 샘플의 통계 값(ex. 평균, 분산 등)이 전체 모수와 비교하여 얼마나 유의하게 차이 나는지를 검정하는 방법이다.

#‘저체중아의 산모’ vs ‘저체중아x의 산모’ 의 체중 차이를 Permutation test 해보자

2. 정상군과 실험군 분류

normal = birthwt[birthwt[,"low"]==0, "lwt"]

normal

case = birthwt[birthwt[,"low"]==1, "lwt"]

case

t.test(normal, case)

3. 두 그룹의 산모 체중에 대한 t 검정 값

real_test = t.test(normal, case)$statistic 

real_test

4. 두 그룹간의 permutation test

permfunc.R

source("permfunc.R") #첨부파일 다운하여 실행

tperm = perm.test(normal, case, n.perm=1000) #1000번 Permutation을 통해 1000개의 t value 생성

hist(tperm)

abline(v=abs(real_test), lty=2, col=2) #실제 t value 가 분포의 극단치에서 보임 (우연이 아닐 것이라고 예상)

# 위 내용은 BITEC (Biomedical Informatics Training and Education Center) 교육내용을 참고하였다.

6. 실습 Example (위와 같은 맥락이지만, 연습삼아 해보셔도 좋을듯하다)

1) make up some ‘true’ data

carrier <- rep(c(0,1), c(100,200))

null.y <- rnorm(300)

alt.y <- rnorm(300, mean=carrier/2)

2) t-test

t.test(null.y~carrier, var.equal=TRUE)

t.test(alt.y~carrier, var.equal=TRUE)

3) permutation test

Permutation test (순열검정법) End.

BioinformaticsAndMe

Cogena-2 (CoExpression 분석) Start.

BioinformaticsAndMe

Cogena-1 (CoExpression 분석) 에 이어지는 내용이다.

8. Drug enrichment analysis

Cogena 에서는 expression에 따른 클러스터에 관련된 Drug을 매칭해주는 작업을 말한다

# Drug data set

cmapDn100_cogena_result <- clEnrich_one(genecl_result, "pam", "10",

sampleLabel=sampleLabel, annofile=system.file("extdata", "CmapDn100.gmt.xz", package="cogena") )

summary(cmapDn100_cogena_result)

CmapDn100.gmt 을 보면 아래와 같은 테이블 형태이다

9. Drug repositioning

Drug repositioning : 임상에서 실패한 약물 or 시판되는 의약품을 재평가하여 새로운 약효를 찾는 방법

(Drug repositioning 은 깊게 들어가면 끝도 없으므로, 나중에 새로운 칼럼으로 다루겠다..)

Cogena-1 (CoExpression 분석) Start.

BioinformaticsAndMe

array 분석을 하다보면, 때로 수많은 DEG(Deferentially Expressed Gene) 형님들을 만나게 된다.

몇 천..?개는 다른 보정이 필요하겠고, 사실 몇 백개라도 골치 아프다.

아.. 수많은 DEG가지고 DAVID 돌려버리면 뭔가 두리뭉실한 기작만 나온다. GG

이럴때 나는

'Co-Expression 분석' 을 한다.

뭐 여러 정의가 있겠지만

간단히 말해서 유전자 발현(up, down) 양상이 비슷한 것 끼리 분석합시다.. 라고 나는 정의하겠다.

다시말해, 수많은 DEG들을 비슷한 발현양상을 갖는 그룹으로 클러스터링하여, 클러스터링된 DEG들의 기작을 살펴보는 것이다.

이것을 쉽게하기 위한 아래 R package를 추천한다.

'Cogena'

Co-expressed gene-set enrichment analysis 

논문에 나와있는 위 그림이 Cogena 를 단번에 설명해준다.

1. DEG 나옴.

2. 발현에따라 클러스터링 해줌 (클러스터링방법은많다).

3. 클러스터링 된 그룹에 무슨 Drug, Pathway가 Enrichment(빵빵)한지 분석을 해줌.

#############################################################################

BioinformaticsAndMe script

#############################################################################

1. cogena 패키지 설치

source("https://bioconductor.org/biocLite.R")

biocLite("cogena")

library(cogena)

# 예제(건선) 데이터 로딩

data(Psoriasis)

# 데이터 확인

ls()

3. Pathway 데이터 로딩

# KEGG Pathway gene set

annoGMT <- "c2.cp.kegg.v5.0.symbols.gmt.xz“

# annotation 경로 지정

annofile <- system.file("extdata", annoGMT, package="cogena")

# GO 분석을 하고 싶다면,

# GO biological process gene set

# annoGMT <- "c5.bp.v5.0.symbols.gmt.xz"

아래표는 Cogena 가 가진 Annotation set이다. 잘 정리했군..

4. 파라미터 지정 (몇개 클러스터로 할건지, 무슨 clustering method를 쓸건지 등등)

# the number of clusters. It can be a vector.

# nClust <- 2:20

nClust <- 10

# Making factor "ct" , "Psoriasis"

sampleLabel <- factor(sampleLabel, levels=c("ct", "Psoriasis"))

# the clustering methods

# clMethods <- ("hierarchical","kmeans","diana","fanny","som","model",

# "sota","pam","clara","agnes") # All the methods can be used together.

clMethods <- c("hierarchical","pam")

# the distance metric

metric <- "correlation“

# the agglomeration method used for hierarchical clustering

method <- "complete"

5. Co-expression Analysis 시작

genecl_result <- coExp(DEexprs, nClust=nClust, clMethods=clMethods,

                       metric=metric, method=method)

summary(genecl_result)

# Enrichment (Pathway) analysis for the co-expressed genes

clen_res <- clEnrich(genecl_result, annofile=annofile, sampleLabel=sampleLabel)

summary(clen_res)

6. 클러스터링 결과 Heatmap 확인

# pam 클러스터링을 통해 10개의 클러스터 제작

# Always make the number as character, please!

# heatmap 그리기

heatmapCluster(clen_res, "pam", "10", maintitle="Psoriasis")

크... 일단 이쁘다.. 10개의 클러스티렁, up|down DEG, two샘플 그룹으로 heatmap이 만들어진다.

아! Pathway 봐야지!

7. 클러스터 Pathway 확인

# The enrichment score for 10 clusters, together with Down-regulated,Up-regulated and All DE genes.

heatmapPEI(clen_res, "pam", "10", maintitle="Pathway analysis for Psoriasis")

유전자 클러스터마다 어떤 pathway가 enrichment한지 보인다 (물론, 전체 only up, only down 도 볼 수있다).

옵션 조정해주면 특정 클러스터의 score 내림차순 기준으로 그림을 다시그려준다. 클러스터9 를 기준으로 다시그려보았다.

그림 왜이렇게 커...;;

Cogena-2 (CoExpression 분석) 에서 이어서 살펴보자.

Cogena-1 (CoExpression 분석) End.

BioinformaticsAndMe