들어가기 전에

1) R의 꽃, R의 주목적이다보니 분량이 어마무시하게 길다. 네이버에서 글 분량갖고 자르진 않겠지만 그래서 상, 하편으로 나눠서 작성할 예정이다. (cookbook 소챕터 3:3으로 나눔)
2) 그리고 이새기들 자꾸 안 알려주고 뜬금없이 라이브러리 소환하세요 하는데 오늘도 깔아야 할 패키지가 있다. 하편에 나오는 것도 있고 상편에 나오는 것도 있는데 그냥 다 깔고 가자. 

install.packages("car")
install.packages("ggplot2")
install.packages("exact2x2")
install.packages("irr")

3) 처음 보는 분석들이 많아서 나도 이게 뭔지 모른다. 따라서 이론적인 설명은 패스. 

회귀분석과 상관계수

참고로 여기서 만드는 데이터는 진짜 새발의 피수준... 나중 가면 라이브러리 깔아서 불러온다. 

> set.seed(1)
> xvar=1:20+rnorm(20,sd=3)
> zvar=1:20/4+rnorm(20,sd=2)
> yvar=2*xvar+xvar*zvar/5+3+rnorm(20,sd=4)

> dat=data.frame(x=xvar,y=yvar,z=zvar)
> dat
            x          y          z
1  -0.8793614  0.2159694  2.0879547
2   2.5509300  8.1415762  2.0642726
3   0.4931142  6.8627566  0.8991300
4   8.7858424 17.5642539 -2.9787034
5   5.9885233 15.2038911  2.4896515
6   3.5385948  8.2293408  1.3877425
7   8.4622872 23.8173481  1.4384090
8  10.2149741 24.5805906 -0.9415048
9  10.7273441 26.7808960  1.2936999
10  9.0838348 30.7526228  3.3358831
11 15.5353435 52.6505709  5.4673591
12 13.1695297 34.2512053  2.7944245
13 11.1362783 35.6025037  4.0253432
14  7.3559003 18.1851647  3.3923899
15 18.3747928 49.1415013  0.9958809
16 15.8651992 52.7105687  3.1700109
17 16.9514292 47.1691760  3.4614201
18 20.8315086 58.7406015  4.3813732
19 21.4636636 78.0409159  6.9500507
20 21.7817040 74.4542008  6.5263515

회귀분석은 이놈으로 해 보자. 

> head(dat)
           x          y         z
1 -0.8793614  0.2159694  2.087955
2  2.5509300  8.1415762  2.064273
3  0.4931142  6.8627566  0.899130
4  8.7858424 17.5642539 -2.978703
5  5.9885233 15.2038911  2.489651
6  3.5385948  8.2293408  1.387743

> tail(dat)
          x        y         z
15 18.37479 49.14150 0.9958809
16 15.86520 52.71057 3.1700109
17 16.95143 47.16918 3.4614201
18 20.83151 58.74060 4.3813732
19 21.46366 78.04092 6.9500507
20 21.78170 74.45420 6.5263515

참고로 데이터 분량이 씁 아 이건 좀 그렇네 싶으면 head()나 tail()을 써 보자. pandas의 head(), tail()과 똑같이 데이터의 머리꼬리를 뽑아준다. 

> head(dfa,5)
  rater1 rater2
1      3      3
2      3      6
3      3      4
4      4      6
5      5      2

숫자 지정은 이런 식으로 하면 된다. 아무튼 다시 분석하러 가 봅시다. 

Correlation(상관분석)

상관계수는 그냥 이걸로 분석하면 된다. 

> cor(dat$x,dat$y)
[1] 0.9669456

혹시 분석을 뭉텅이로 해야 해서 일일이 못 돌릴 것 같은가? 

> cor(dat)
          x         y         z
x 1.0000000 0.9669456 0.5479301
y 0.9669456 1.0000000 0.6893575
z 0.5479301 0.6893575 1.0000000

아 그럼 행렬 만들면 된다. 

> round(cor(dat),2)
     x    y    z
x 1.00 0.97 0.55
y 0.97 1.00 0.69
z 0.55 0.69 1.00

반올림도 해준다. 

Bradford assay의 경우 무게를 알고 있는 단백질을 Bradford adday를 통해 분석한 다음 그걸로 standard curve를 그리게 되는데, 이 때 r^2를 계산하게 된다. 그리고 r^2가 1에 가까울 수록 잘 된거다. 그 알제곱이 저거다. 아마도? 

Linear regression(선형 회귀분석)

> fit=lm(y~x,data=dat)
> fit

Call:
lm(formula = y ~ x, data = dat)

Coefficients:
(Intercept)            x  
     -1.962        3.172

> fit=lm(dat$y~dat$x)
> fit

Call:
lm(formula = dat$y ~ dat$x)

Coefficients:
(Intercept)        dat$x  
     -1.962        3.172

선형 회귀 분석은 이걸로 하면 된다. 둘 다 같은건데, 위 코드는 데이터프레임 컬럼 불러서 진행한거고 아래는 데이터프레임 컬럼 벡터로 픽해서 했다. Intercept는 절편(여기서는 y 절편). 

> summary(fit)

Call:
lm(formula = dat$y ~ dat$x)

Residuals:
   Min     1Q Median     3Q    Max 
-8.341 -5.304 -1.047  4.505 11.925 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.9619     2.5497  -0.769    0.452    
dat$x         3.1718     0.1971  16.089 3.97e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.894 on 18 degrees of freedom
Multiple R-squared:  0.935,	Adjusted R-squared:  0.9314 
F-statistic: 258.9 on 1 and 18 DF,  p-value: 3.968e-12

자세한 정보가 궁금하다면 summary()를 쓰면 된다. 

네? 여러개요? 

> fit2=lm(y~x+z,data=dat)
> fit2

Call:
lm(formula = y ~ x + z, data = dat)

Coefficients:
(Intercept)            x            z  
     -3.145        2.762        2.190

> fit2=lm(dat$y~dat$x+dat$z)
> fit2

Call:
lm(formula = dat$y ~ dat$x + dat$z)

Coefficients:
(Intercept)        dat$x        dat$z  
     -3.145        2.762        2.190

둘 중 하나 픽해서 하면 된다. (Intercept는 아까도 말했듯 절편)

> summary(fit2)

Call:
lm(formula = dat$y ~ dat$x + dat$z)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4155 -2.8132  0.0799  1.9206  6.6844 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.1445     1.7599  -1.787 0.091817 .  
dat$x         2.7620     0.1610  17.160  3.6e-12 ***
dat$z         2.1896     0.4713   4.646 0.000231 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.025 on 17 degrees of freedom
Multiple R-squared:  0.9714,	Adjusted R-squared:  0.968 
F-statistic: 288.3 on 2 and 17 DF,  p-value: 7.669e-14

근데 summary는 요약 아님? 

Interaction 껴서 분석하기

> fit3=lm(y~x*z,data=dat)
> fit3

Call:
lm(formula = y ~ x * z, data = dat)

Coefficients:
(Intercept)            x            z          x:z  
     2.8057       2.1931      -0.4598       0.1950

> fit3=lm(y~x+z+x:z,data=dat)
> fit3

Call:
lm(formula = y ~ x + z + x:z, data = dat)

Coefficients:
(Intercept)            x            z          x:z  
     2.8057       2.1931      -0.4598       0.1950

Interaction 껴서 분석하는 것도 된다. 근데 이게 뭐 하는건지는 나도 모른다. 

> summary(fit3)

Call:
lm(formula = y ~ x + z + x:z, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5313 -2.5372  0.1603  1.8645  6.7627 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.80570    2.41253   1.163  0.26189    
x            2.19312    0.22715   9.655 4.47e-08 ***
z           -0.45978    0.94445  -0.487  0.63299    
x:z          0.19497    0.06346   3.072  0.00729 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.29 on 16 degrees of freedom
Multiple R-squared:  0.982,	Adjusted R-squared:  0.9786 
F-statistic: 290.7 on 3 and 16 DF,  p-value: 3.673e-14

물론 얘도 summary()로 자세히 볼 수 있다. 

T-test

> x=floor(runif(25)*100)
> y=floor(runif(25)*105)
> x
 [1] 99 49 48 17 75 45 51 20 22 59 57  7  3 64 92 59 56 52 98 50 68 60 23 25 72
> y
 [1] 47 18 78 11 90 64 58 34 47 52 18 55  7 29 22 29 93 46 81 92 43  6 35 75 35

데이터 가져오기 귀찮으면 만들면 된다. 

> data=data.frame(x,y)
> head(data)
   x  y
1 99 47
2 49 18
3 48 78
4 17 11
5 75 90
6 45 64

오케이 가즈아! 

평범하게 T-test 웰치스? 

> t.test(x,y)

	Welch Two Sample t-test

data:  x and y
t = 0.56293, df = 47.977, p-value = 0.5761
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.90445  19.38445
sample estimates:
mean of x mean of y 
    50.84     46.60

이거 치면 끝. 

> t.test(data$x,data$y)

이것도 먹힌다. 

> t.test(x,y,var.equal=TRUE)

	Two Sample t-test

data:  x and y
t = 0.56293, df = 48, p-value = 0.5761
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.90427  19.38427
sample estimates:
mean of x mean of y 
    50.84     46.60

var.equal=TRUE 옵션을 주면 스튜던트 T-test를 진행할 수 있다. default는 웰치스 T-test. 잠깐만 웰치스라고 했음? 아니 그 음료수 아니... 스펠이 왜 똑같지? 

Paired T-test

> t.test(x,y,paired=TRUE)

	Paired t-test

data:  x and y
t = 0.6259, df = 24, p-value = 0.5373
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.741436 18.221436
sample estimates:
mean of the differences 
                   4.24

paired=TRUE 옵션을 주면 paired T-test도 가능하다. 

One sample T-test

엥? 근데 저 데이터 하난데 어캄? 

> t.test(x,mu=0)

	One Sample t-test

data:  x
t = 9.6518, df = 24, p-value = 9.776e-10
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 39.9686 61.7114
sample estimates:
mean of x 
    50.84

> t.test(y,mu=0)

	One Sample t-test

data:  y
t = 8.6554, df = 24, p-value = 7.617e-09
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 35.48812 57.71188
sample estimates:
mean of x 
     46.6

단식도 된다. mu=0은 모평균을 0으로 놓고 본다는 얘기. 모평균은 뮤, 모표준편차는 시그마 제곱이다. 

빈도 분석

> data <- read.table(header=TRUE, text='
+  condition result
+    control      0
+    control      0
+    control      0
+    control      0
+  treatment      1
+    control      0
+    control      0
+  treatment      0
+  treatment      1
+    control      1
+  treatment      1
+  treatment      1
+  treatment      1
+  treatment      1
+  treatment      0
+    control      0
+    control      1
+    control      0
+    control      1
+  treatment      0
+  treatment      1
+  treatment      0
+  treatment      0
+    control      0
+  treatment      1
+    control      0
+    control      0
+  treatment      1
+  treatment      0
+  treatment      1
+ ')
> head(data)
  condition result
1   control      0
2   control      0
3   control      0
4   control      0
5 treatment      1
6   control      0

데이터는 고마운데 분량 이거 ㄷㄷ

이럴 땐 이런 분석! 이라고 올려뒀더라고... 

X^2(카이스퀘어, Chi-square)

> ct=table(data$result)
> ct

 0  1 
17 13

참고로 빈도 분석은 어떤 형태로든 contingency table을 먼저 만들고 시작한다. 여기서는 글 분량때문에 같은 종류의 테이블에 대해서는 코드 생략하고 뭘로 만들었는지만 기재한다. 

> chisq.test(ct)

	Chi-squared test for given probabilities

data:  ct
X-squared = 0.53333, df = 1, p-value = 0.4652

평범한 카이스퀘어

> pt=c(.75,.25) # 확률값
> chisq.test(ct,p=pt)

	Chi-squared test for given probabilities

data:  ct
X-squared = 5.3778, df = 1, p-value = 0.02039

기대 빈도수가 다르다고 전제하고 계산하는 카이스퀘어

> chi_res=chisq.test(ct,p=pt)
> str(chi_res)
List of 9
 $ statistic: Named num 5.38
  ..- attr(*, "names")= chr "X-squared"
 $ parameter: Named num 1
  ..- attr(*, "names")= chr "df"
 $ p.value  : num 0.0204
 $ method   : chr "Chi-squared test for given probabilities"
 $ data.name: chr "ct"
 $ observed : 'table' int [1:2(1d)] 17 13
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr [1:2] "0" "1"
 $ expected : Named num [1:2] 22.5 7.5
  ..- attr(*, "names")= chr [1:2] "0" "1"
 $ residuals: 'table' num [1:2(1d)] -1.16 2.01
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr [1:2] "0" "1"
 $ stdres   : 'table' num [1:2(1d)] -2.32 2.32
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr [1:2] "0" "1"
 - attr(*, "class")= chr "htest"

얘는 카이스퀘어 돌린 거 변수 지정하고 str()로 봐야 자세히 보인다. 

Exact binomial test

> binom.test(ct,p=0.5)

	Exact binomial test

data:  ct
number of successes = 17, number of trials = 30, p-value = 0.5847
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.3742735 0.7453925
sample estimates:
probability of success 
             0.5666667

요거는 위에서 만든 contingency table을 바탕으로 Exact binomial test를 한 결과. 저거 근데 확률 안 써도 값 같더만... (심지어 확률을 바꿔도 값이 같음) 아, p-value는 좀 다르긴 함. 

> binom.test(ct,p=0.25)

	Exact binomial test

data:  ct
number of successes = 17, number of trials = 30, p-value = 0.0002157
alternative hypothesis: true probability of success is not equal to 0.25
95 percent confidence interval:
 0.3742735 0.7453925
sample estimates:
probability of success 
             0.5666667

카이스퀘어와 달리 두 번째 컬럼값의 확률만 기재한다. 그러면 알아서 앞에껏도 계산해주는 듯. 

> bin_res=binom.test(ct,p=0.5)
> str(bin_res)
List of 9
 $ statistic  : Named num 17
  ..- attr(*, "names")= chr "number of successes"
 $ parameter  : Named num 30
  ..- attr(*, "names")= chr "number of trials"
 $ p.value    : Named num 0.585
  ..- attr(*, "names")= chr "0"
 $ conf.int   : num [1:2] 0.374 0.745
  ..- attr(*, "conf.level")= num 0.95
 $ estimate   : Named num 0.567
  ..- attr(*, "names")= chr "probability of success"
 $ null.value : Named num 0.5
  ..- attr(*, "names")= chr "probability of success"
 $ alternative: chr "two.sided"
 $ method     : chr "Exact binomial test"
 $ data.name  : chr "ct"
 - attr(*, "class")= chr "htest"

얘도 자세히 보려면 카이스퀘어처럼 변수 지정해주고 str()로 봐야 한다. 

X^2(카이스퀘어)

> ct=table(data$condition,data$result)
> ct
           
             0  1
  control   11  3
  treatment  6 10

이번꺼는 그룹간의 독립성을 보는 거라 테이블이 조금 다르다. 위 테이블은 결과값만 가지고 만는거지만, 이 테이블은 결과값과 컨디션으로 들어간다. 

> chisq.test(ct)

	Pearson's Chi-squared test with Yates' continuity correction

data:  ct
X-squared = 3.593, df = 1, p-value = 0.05802

아이고 어슨이형 또뵙네 

> chisq.test(ct,correct=FALSE)

	Pearson's Chi-squared test

data:  ct
X-squared = 5.1293, df = 1, p-value = 0.02353

correct=FALSE 옵션을 주면 Yates's correction for continuity가 빠진다. ...근데 그거 막 그렇게 빼도 돼? 

Fisher’s exact test

> fisher.test(ct)

	Fisher's Exact Test for Count Data

data:  ct
p-value = 0.03293
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.966861 45.555016
sample estimates:
odds ratio 
  5.714369

이건 피셔의 정확도 검정? 그건데 샘플 수가 작을 때 사용하는 분석법이라고 한다. 사실 통계가 데이터가 많을수록 좋은 것도 사실이지만, 그게 항상 그렇게 되리라는 법이 없다. 예를 들자면 희귀병 환자라던가... 

Cochran-Mantel-Haenszel test

fish <- read.table(header=TRUE, text='
  Location Allele   Habitat Count
 tillamook     94    marine    56
 tillamook     94 estuarine    69
 tillamook non-94    marine    40
 tillamook non-94 estuarine    77
   yaquina     94    marine    61
   yaquina     94 estuarine   257
   yaquina non-94    marine    57
   yaquina non-94 estuarine   301
     alsea     94    marine    73
     alsea     94 estuarine    65
     alsea non-94    marine    71
     alsea non-94 estuarine    79
    umpqua     94    marine    71
    umpqua     94 estuarine    48
    umpqua non-94    marine    55
    umpqua non-94 estuarine    48
')
> head(fish)
   Location Allele   Habitat Count
1 tillamook     94    marine    56
2 tillamook     94 estuarine    69
3 tillamook non-94    marine    40
4 tillamook non-94 estuarine    77
5   yaquina     94    marine    61
6   yaquina     94 estuarine   257

코크란-만텔-헨젤 검정은 또 표가 다른 게 들어간다. 살려줘. 

> ct=xtabs(Count~Allele+Habitat+Location,data=fish)
> ct
, , Location = alsea

        Habitat
Allele   estuarine marine
  94            65     73
  non-94        79     71

, , Location = tillamook

        Habitat
Allele   estuarine marine
  94            69     56
  non-94        77     40

, , Location = umpqua

        Habitat
Allele   estuarine marine
  94            48     71
  non-94        48     55

, , Location = yaquina

        Habitat
Allele   estuarine marine
  94           257     61
  non-94       301     57

얘는 그냥 contengency table로 만들면 3차원이 된다. 지역별로 테이블이 다 나오기때문에 

> ftable(ct)
                 Location alsea tillamook umpqua yaquina
Allele Habitat                                          
94     estuarine             65        69     48     257
       marine                73        56     71      61
non-94 estuarine             79        77     48     301
       marine                71        40     55      57

이걸 2차원으로 다시 합칠거다. 

> mantelhaen.test(ct)

	Mantel-Haenszel chi-squared test with continuity correction

data:  ct
Mantel-Haenszel X-squared = 5.0497, df = 1, p-value = 0.02463
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.6005522 0.9593077
sample estimates:
common odds ratio 
         0.759022

그리고 분석하면 된다... 분석 개빡시네 이거. 

McNemar’s test

> data <- read.table(header=TRUE, text='
+  subject time result
+        1  pre      0
+        1 post      1
+        2  pre      1
+        2 post      1
+        3  pre      0
+        3 post      1
+        4  pre      1
+        4 post      0
+        5  pre      1
+        5 post      1
+        6  pre      0
+        6 post      1
+        7  pre      0
+        7 post      1
+        8  pre      0
+        8 post      1
+        9  pre      0
+        9 post      1
+       10  pre      1
+       10 post      1
+       11  pre      0
+       11 post      0
+       12  pre      1
+       12 post      1
+       13  pre      0
+       13 post      1
+       14  pre      0
+       14 post      0
+       15  pre      0
+       15 post      1
+ ')
> head(data)
  subject time result
1       1  pre      0
2       1 post      1
3       2  pre      1
4       2 post      1
5       3  pre      0
6       3 post      1

Mcnemar's test는 또 다른 데이터로 진행한다. 살려줘요. 

> library(tidyr)
> data_wide <- spread(data, time, result)
> head(data_wide)
  subject post pre
1       1    1   0
2       2    1   1
3       3    1   0
4       4    0   1
5       5    1   1
6       6    1   0

어? 이거 이렇게 처리하는 거 어디서 봤는데? 

> ct=table(data_wide[,c("pre","post")])
> ct
   post
pre 0 1
  0 2 8
  1 1 4

아무튼 테이블을 만들어주고... 

> mcnemar.test(ct)

	McNemar's Chi-squared test with continuity correction

data:  ct
McNemar's chi-squared = 4, df = 1, p-value = 0.0455

분⭐석

> library(exact2x2)
필요한 패키지를 로딩중입니다: exactci
필요한 패키지를 로딩중입니다: ssanv
필요한 패키지를 로딩중입니다: testthat

다음의 패키지를 부착합니다: ‘testthat’

The following object is masked from ‘package:tidyr’:

    matches

> mcnemar.exact(ct)

	Exact McNemar test (with central confidence intervals)

data:  ct
b = 8, c = 1, p-value = 0.03906
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   1.072554 354.981246
sample estimates:
odds ratio 
         8

excat버전을 쓰려면 exact2x2라는 패키지를 깔아야 한다. (위에 썼던 그거 맞다)

ANOVA

들어보긴 했는데 이거 뭐임? SW 이름 아니었어? 이사람

> data <- read.table(header=TRUE, text='
+  subject sex   age before after
+        1   F   old    9.5   7.1
+        2   M   old   10.3  11.0
+        3   M   old    7.5   5.8
+        4   F   old   12.4   8.8
+        5   M   old   10.2   8.6
+        6   M   old   11.0   8.0
+        7   M young    9.1   3.0
+        8   F young    7.9   5.2
+        9   F   old    6.6   3.4
+       10   M young    7.7   4.0
+       11   M young    9.4   5.3
+       12   M   old   11.6  11.3
+       13   M young    9.9   4.6
+       14   F young    8.6   6.4
+       15   F young   14.3  13.5
+       16   F   old    9.2   4.7
+       17   M young    9.8   5.1
+       18   F   old    9.9   7.3
+       19   F young   13.0   9.5
+       20   M young   10.2   5.4
+       21   M young    9.0   3.7
+       22   F young    7.9   6.2
+       23   M   old   10.1  10.0
+       24   M young    9.0   1.7
+       25   M young    8.6   2.9
+       26   M young    9.4   3.2
+       27   M young    9.7   4.7
+       28   M young    9.3   4.9
+       29   F young   10.7   9.8
+       30   M   old    9.3   9.4
+ ')
> head(data)
  subject sex age before after
1       1   F old    9.5   7.1
2       2   M old   10.3  11.0
3       3   M old    7.5   5.8
4       4   F old   12.4   8.8
5       5   M old   10.2   8.6
6       6   M old   11.0   8.0

Age에 old랑 young은 뭔 기준인지 모르겠지만 일단 패스. 

> data$subject=factor(data$subject)

이 데이터는 ANOVA 들어가기 전에 저기 subject를 팩터화해줘서 연속적인 데이터 취급 받는 사태를 방지해야 한다. 참고로 ANOVA는 one-way랑 two-way가 있다. 

one-way ANOVA

> aov1=aov(before~sex,data=data)
> aov1
Call:
   aov(formula = before ~ sex, data = data)

Terms:
                     sex Residuals
Sum of Squares   1.52861  74.70105
Deg. of Freedom        1        28

Residual standard error: 1.633369
Estimated effects may be unbalanced
> summary(aov1)
            Df Sum Sq Mean Sq F value Pr(>F)
sex          1   1.53   1.529   0.573  0.455
Residuals   28  74.70   2.668
> model.tables(aov1,"means")
Tables of means
Grand mean
         
9.703333 

 sex 
     F      M
    10  9.532
rep 11 19.000

two-way ANOVA

> aov3=aov(after~sex*age,data=data)
> aov3
Call:
   aov(formula = after ~ sex * age, data = data)

Terms:
                      sex       age   sex:age Residuals
Sum of Squares   16.07755  38.96110  89.61138 103.51164
Deg. of Freedom         1         1         1        26

Residual standard error: 1.995299
Estimated effects may be unbalanced
> summary(aov3)
            Df Sum Sq Mean Sq F value  Pr(>F)    
sex          1  16.08   16.08   4.038  0.0550 .  
age          1  38.96   38.96   9.786  0.0043 ** 
sex:age      1  89.61   89.61  22.509 6.6e-05 ***
Residuals   26 103.51    3.98                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> model.tables(aov3,"means")
Tables of means
Grand mean
         
6.483333 

 sex 
         F      M
     7.445  5.926
rep 11.000 19.000

 age 
       old  young
     7.874  5.556
rep 12.000 18.000

 sex:age 
     age
sex   old    young 
  F    6.260  8.433
  rep  5.000  6.000
  M    9.157  4.042
  rep  7.000 12.000

> aov3=aov(after~sex+age+sex:age,data=data)
> summary(aov3)
            Df Sum Sq Mean Sq F value  Pr(>F)    
sex          1  16.08   16.08   4.038  0.0550 .  
age          1  38.96   38.96   9.786  0.0043 ** 
sex:age      1  89.61   89.61  22.509 6.6e-05 ***
Residuals   26 103.51    3.98                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

둘이 같은 코드. 

Tukey HSD post-hoc test

> TukeyHSD(aov1)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = before ~ sex, data = data)

$sex
          diff       lwr       upr     p adj
M-F -0.4684211 -1.736038 0.7991961 0.4554053

One-way

> TukeyHSD(aov3)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = after ~ sex + age + sex:age, data = data)

$sex
         diff       lwr        upr     p adj
M-F -1.519139 -3.073025 0.03474709 0.0549625

$age
              diff       lwr        upr     p adj
young-old -2.31785 -3.846349 -0.7893498 0.0044215

$`sex:age`
                      diff        lwr        upr     p adj
M:old-F:old      2.8971429 -0.3079526  6.1022384 0.0869856
F:young-F:old    2.1733333 -1.1411824  5.4878491 0.2966111
M:young-F:old   -2.2183333 -5.1319553  0.6952887 0.1832890
F:young-M:old   -0.7238095 -3.7691188  2.3214997 0.9138789
M:young-M:old   -5.1154762 -7.7187601 -2.5121923 0.0000676
M:young-F:young -4.3916667 -7.1285380 -1.6547953 0.0008841

Two-way

ANOVA 데이터를 바탕으로 tukey HSD 분석을 할 수도 있다. 

ANOVAs with within-subjects variables

> data_long <- gather(data, time, value, before:after)
> head(data_long)
  subject sex age   time value
1       1   F old before   9.5
2       2   M old before  10.3
3       3   M old before   7.5
4       4   F old before  12.4
5       5   M old before  10.2
6       6   M old before  11.0
> data_long$subject <- factor(data_long$subject)

subject variable를 끼고 하려고 표를 길게 만들었단다. (아까 팩터화한 그거 맞음)

> aov_time=aov(value~time+Error(subject/time),data=data_long)
> aov_time

Call:
aov(formula = value ~ time + Error(subject/time), data = data_long)

Grand Mean: 8.093333

Stratum 1: subject

Terms:
                Residuals
Sum of Squares   261.2473
Deg. of Freedom        29

Residual standard error: 3.001421

Stratum 2: subject:time

Terms:
                   time Residuals
Sum of Squares  155.526    63.144
Deg. of Freedom       1        29

Residual standard error: 1.475595
Estimated effects are balanced

one-way

> aov_age_time=aov(value~age*time+Error(subject/time),data=data_long)
> aov_age_time

Call:
aov(formula = value ~ age * time + Error(subject/time), data = data_long)

Grand Mean: 8.093333

Stratum 1: subject

Terms:
                      age Residuals
Sum of Squares   24.44011 236.80722
Deg. of Freedom         1        28

Residual standard error: 2.908161
1 out of 2 effects not estimable
Estimated effects are balanced

Stratum 2: subject:time

Terms:
                   time age:time Residuals
Sum of Squares  155.526   18.769    44.375
Deg. of Freedom       1        1        28

Residual standard error: 1.258897
Estimated effects may be unbalanced

mixed design

일반적으로 우리가 쓰는 숫자는 10진수가 맞는데, 컴퓨터는 손가락이 두 개라 이진법을 쓴다. 예, 그겁니다. 일반적으로 10진법을 2진법으로 변환할때는 2로 나누는데, 13을 예시로 들자면 

1. 13/2=6...1
2. 6/2=3...0
3. 3/2=1...1

이렇게 구한 다음 3번 결과의 몫 1부터 시작해 3번 결과의 나머지-2번 결과의 나머지-1번 결과의 나머지 순으로 올라가서 1101이 된다. 1의 보수는 여기서 0을 1로, 1을 0으로 바꾼 0010. 2의 보수는 1의 보수에 1을 더하면 된다. (0011)

참고로 이진수로 변환된 숫자는 8, 16진수와 상호변환이 가능한데 

1. 110010010을 8진수로 바꾸려면 110 010 010으로 세자리씩 끊어서 10진법으로 바꾼다. (622)
2. 110010010을 16진수로 바꾸려면 1 1001 0010으로 네자리씩 끊어서 10진법으로 바꾼다. (192)
3. FF(16)를 8진수로 바꾸려면 2진수로 바꾼 다음 세자리씩 끊어서 10진법으로 바꾼다. (11 111 111->377)
4. 77(8)를 16진수로 바꾸려면 2진수로 바꾼 다음 네자리씩 끊어서 10진법으로 바꾼다. (11 1111->3F)

a = int(input())
bin_list = []
while a >= 1:
    bin_list.append(a%2)
    a = int(a/2)
bin = bin_list[::-1]
bin = ''.join(map(str,bin))
print(bin)
# 입력한 숫자를 이진수로 바꿔준다. (물론 수동 방식이다)

수동변환

for i in range(len(bin_list)):
    if bin_list[i] == 1:
        bin_list[i] = 0
    else:
        bin_list[i] = 1
bin = bin_list[::-1]
bin = ''.join(map(str,bin))
print(bin)
# 1의 보수(2의 보수는 이 방식으로 할 경우 처리가 되게 애매해지는 문제가 있다)

1의 보수

b=int(comp_1,2)
b=b+1
print(b)
# 그래서 2의 보수는 일단 10진수로 바꾼 다음, 1을 더하고 다시 바꿀 예정. 
comp2_list=[]
while b >= 1:
    comp2_list.append(b%2)
    b = int(b/2)
comp_2 = comp2_list[::-1]
comp_2 = ''.join(map(str,comp_2))
print(comp_2)
# 아, 0 빠진건 알아서 붙이세요.

2의 보수

a = int(input())
hex_list = []
while a >= 1:
    if a % 16 == 10: # 16진수에서 10~15까지는 각각 A~F로 표기한다. 즉, 여기에 대한 처리를 따로 진행해야 한다. 
        hex_list.append("A")
    elif a % 16 == 11:
        hex_list.append("B")
    elif a % 16 == 12:
        hex_list.append("C")
    elif a % 16 == 13:
        hex_list.append("D")
    elif a % 16 == 14: 
        hex_list.append("E")
    elif a % 16 == 15:
        hex_list.append("F")
    else: 
        hex_list.append(a % 16)
    a = int(a / 16)
hex = hex_list[::-1]
hex=''.join(map(str,hex))
print(hex)

자매품(16진수)

16진수의 경우 10~15까지를 A~F로 표기하기때문에 그걸 변환해줘야 한다. 

풀긴 풀었는데 이거 로직이 뭐가 문제인건지 제시한거랑 출력값이 다름. 

일단 문제가 뭐냐... 예를 들어서 아래와 같은 텍스트가 있다면 

츄라이츄라이

출력값은 

000123

이 된다. 네 번째 글자 츄를 기준으로 했을 때 

츄라이에서 츄/이츄/라이츄를 찾는건데 저기서 일치하는 게 츄 하나고 그게 한글자짜리거든. 

다섯번째 글자 라를 기준으로 하면 츄라이츄에서 라/츄라/이츄라/라이츄라 이렇게 찾게 되는거고 그렇게 되면 일치하는 문자열 중 제일 긴 게 츄라라서 2. 

이해하는 데 하루 걸렸다더니 실화냐고요? 네. 잠깐 뇌에 블루스크린 버프 와서요. 

text=input("Insert text: ")
# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다.
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다.
find_list=[]
length=len(text)
for i in range(1,length+1):
    if len(text[0:i]) == 1:
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위
        max_list=[]
        for j in range(1,len(find)+1):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                max_list.append(len(subset))
            else:
                max_list.append(0)
            max_list=set(max_list)
            max_list=list(max_list)
            max_values=max(max_list) # 리스트에서 최대값을 추출
        find_list.append(max_values) # 했으면 넣어주세요
# 리스트 출력이 뭔가 이상하다.
# set을 적용하게 되면 각 iteration별로 고유값이 나오는 게 아니라 엉뚱한 값이 나오게 된다.
# append의 대상이 되는 리스트는 밖에 있지만, append는 안에 있어서 또 애매하고... append를 밖으로 뺄 수도 없다.
find_text=''.join(map(str,find_list))
# 리스트 형태로 처리한 것을 문자열로 붙여서 출력
print(find_text)

일단 전체 코드는 이거.

​

사실 문제가 정말 궁서체로 개빡세다. 어제 추워서 쉬긴 했지만 그래도 이틀 걸렸다. (주말에는 공부 안 함) 코딩은 컴퓨터 시켜먹으려면 일단 사람이 이해를 해야 하는데 이해하고 착수하는데 꽤 걸림+이게 내 맘대로 안돼서 걸림...

​

위에 설명한 게 문제인데, 그럼 필요한 기능이 뭐냐...

1. 리스트의 인덱싱/슬라이싱(각각 필요한 부분을 단식/범위로 찾는 기능)

2. For(반복문)

3. if(제어문)

4. 특정 문자를 찾는 기능(대충 R의 grep같은 거)

5. 입력(이건 맨 나중이라 우선순위 미뤄도 된다)

크게 이렇게 필요하다.

​

그럼 차근차근 게임코딩빨로 해봅시다. 이사람 대체 근데 솔직히 두뇌풀가동 이후로 엑스트라 밀고 겜코딩 켠 적 없음

Indexing/slicing 일반화 및 찾을 영역 정의하기 

일단 본인은 별찍기때도 그러했듯 small scale(고정값 혹은 실행하는 데 오래 안 걸리는 작은 범위)에서 해보고 일반화 과정을 거친다. 대충 수열로 치자면 어떤 규칙으로 변하는지 몇 개 찍어보고 일반항 공식 도출하는 방식이다. (별찍기 그렇게 한 다음 while로 찍었더니 while이 더 편하더라...)

어떤 n글자짜리 문자열이 있을 때 

-이 문자열의 인덱싱 주소는 0부터 (n-1)이고
-이 문자열의 슬라이싱 주소는 0부터 n까지이다. 이게 왜 그러냐... 아래 텍스트를 예시로 들어보자. 

>>> text="HAPPY DAYS!"
>>> len(text)
11
>>> print(text[0:11])
HAPPY DAYS!

H,A,P,P,Y, ,D,A,Y,S,!, 해서 11자(중간에 오타가 아니라 공백이 두 개 이고 이 때 인덱싱 주소는 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10이다. 근데 왜 11까지냐... 0에서 10까지니까 슬라이싱을 0:10으로 해 보면 

>>> print(text[0:10])
HAPPY DAYS

뭐야 느낌표 돌려줘요!!! 

이게 왜 이렇게 되는거냐면 파이썬 슬라이싱은 [부터:미만]으로 들어가기 때문이다. 그래서 0:10으로 하면 0부터 10 미만, 즉 9까지 뽑아준다. 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        print(find) # 이게 찾을 범위

이게 찾는 범위를 지정하는 코드. 왜 1부터냐... 저게 [부터:미만]이라고 했는데 일단 첫 빠따는 뭐가 들어가던 0이다. 왜냐, 찾을 영역은 지정할 수 있는데 찾을 글자가 없다. 뭔 소리냐... 

[0,1,2,3,4,5]

위와 같은 리스트가 있을 때, 찾을 범위는 [0], [0,1],...,[0,4]까지이고 찾아야 할 글자는 범바범인데 범위가 [0]일 때 [1], [0,1]일 때 [2],[1,2] 이런 식으로 들어간다. 즉 찾는 영역은 앞에서부터 순차적으로 들어가고 찾아야 하는 영역은 뒤에서부터 순차적으로 들어간다. 

일차적인 for문에서 range가 1,len(text)+1인 이유도 그것때문이다. 저거 없이 그냥 박아버리면 0부터 시작하기때문에 나중에 음수 인덱싱을 하게 될 경우 문제가 생긴다. (-0이나 0이나 그게 그거라...) 그리고 length+1 해야 슬라이싱하면 length값까지 잡아준다. 

찾아야 할 부분집합 슬라이싱하기

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        subset_init=text_sub[i-1]
        print(subset_init) # 첫 서브셋이자 시발점

이 부분이 1차 난관이었다. 찾아야 하는 부분집합의 개수는 찾아야 할 영역의 크기에 따라 순차적으로 줄어들고, 찾아야 할 영역의 크기는 문자열의 길이에 따라 달라진다. 위에서 깜빡하고 안 썼는데, 이것은 마치 김밥 한 줄을 썰어서 왼쪽 꼬다리에서부터 오른쪽 꼬다리를 찾는 것과 비슷한 이치. 그러면 범위는 일단 김밥 한 줄로 잡아야 한다. 

이 때 전체 텍스트는 김밥 한 줄, 찾을 부분은 왼쪽 꼬다리라고 보면 된다. 찾는 부분집합은 찾는 영역에 따라 다르지만 김밥 한 줄 기준으로 오른쪽 꼬다리부터 하나씩 시작이다. (찾을 범위에서 오른쪽 꼬다리가 빠지고, 부분집합을 만드는 범위에서 왼쪽 꼬다리가 빠진다) 

그러니까 김밥 한 줄에서 오른쪽 꼬다리부터 순차적으로 지정해야 해서 원래 저 안에 for문이 또 들어간다. 그래서 전체 코드를 보면 for문이 두 개 들어간 것. ...저 시발점을 잘못 잡은 게 오른쪽 꼬다리면 무조건 -1번 인덱싱(맨 뒤)인데... 

제어문

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_subset=[]
max_subset=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_subset.append(0)
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)+1):
            subset=text_sub[1:j+1] # 찾을 텍스트
            # 텍스트 유무에 따른 처리는 했는데, 문제는 텍스트가 있을 때 길이가 아니라 리스트 자체가 출력된다... 이거 어쩔겨... 
            if subset in find: 
                print(subset)
                find_subset.append(1)
            else: 
                find_subset.append(0)
print(find_subset)

[0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0]

(마른세수) 

출력이 개판인데 그거는 이따가 또 얘기합시다... 진짜 출력 잡는것도 대환장파티였어요... 

아무튼 저기 안쪽에 for문 하나 더 들어간 게 찾아야 할 부분집합이라는 건 알겠는데, if가 왜 있느냐... 일단 저 코드는 일치하는 게 있으면 1, 없으면 0을 끼워 넣으라는 의미긴 하다. 그러니까 쉽게 말하자면 일치하는 게 있냐 없냐를 먼저 보면서 범위부터 손볼 요량이었다. 

if subset in find: 
    print(subset)
    max_subset.append(1)
    find_subset.append(len(max_subset))
else: 
    find_subset.append(0)

제어문은 일단 이런 식. 해당하는 텍스트가 있을 때 일치하는 텍스트의 '길이'를 넣는다. 

if subset in find:
    max_list.append(len(subset))
else:
    max_list.append(0)
max_list=set(max_list)
max_list=list(max_list)
max_values=max(max_list) # 리스트에서 최대값을 추출

전체 코드에서는 이런 식으로 처리한다. 

출력 처리

여기서 정말 마른세수 여러번 했다... (마른세수)

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_subset=[] # 서브셋들의 최대값을 담아서 최종적으로 출력하는 리스트
max_subset=[] # 서브셋들의 최대값을 매기기 위한 리스트
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_subset.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)):
            print(j,"th iteration")
            subset=text_sub[j:] # 찾을 텍스트
            print(find,len(find))
            print(subset,"*")
            if subset in find:
                max_subset.append(len(subset))
                max_values=max(max_subset)
                print(max_subset,max_values)
                find_subset.append(max(max_subset))
                print(find_subset)
            else:
                find_subset.append(0)
                print(find_subset)
# 리스트 출력이 뭔가 이상하다. 이거 중복값 처리를 어떻게 해야 하는거지?

이놈은 돌려봤더니 중복값이 고대로 나오고... 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_list=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                find_list.append(len(subset))
            else: 
                find_list.append(0)
            find_list=set(find_list)
            find_list=list(find_list)
# 리스트 출력이 뭔가 이상하다.
print(find_list)

이놈은 돌렸더니 중복값이 다 사라져서 글자수랑 안 맞고... 

일단 iteration의 기준이 뭐냐면, 찾을 범위를 하나 지정하고 거기서 문자열 차즌ㄴ 게 하나의 iteration이다. 즉, 첫번째 for문을 돌 때마다 값이 하나씩 저장되어야 한다. 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_list=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        max_list=[]
        for j in range(1,len(find)+1):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                max_list.append(len(subset))
            else: 
                max_list.append(0)
            max_list=set(max_list)
            max_list=list(max_list)
            max_values=max(max_list)
        find_list.append(max_values)
# 리스트 출력이 뭔가 이상하다.
# set을 적용하게 되면 각 iteration별로 고유값이 나오는 게 아니라 엉뚱한 값이 나오게 된다. 
# append의 대상이 되는 리스트는 밖에 있지만, append는 안에 있어서 또 애매하고... append를 밖으로 뺄 수도 없다. 
find_text=''.join(map(str,find_list))
print(find_text)

근데 이게 이거 뺀다고 해결되는거 실화냐... 아, 저기 join은 출력을 리스트 형태로 하기 위해 넣은거다. 

AGTC AGTC CGTG ATCT AGCT AGCT AGTC GCTG CATC AGTA C
0000 1234 1121 1121 1212 3456 7834 2232 2223 3452 1 # 돌려서 나온 결과
0000 1234 0000 1000 1200 1200 1234 0000 0100 1231 0 # 문제에서 제시한 결과

근데 이거 왜 결과가 이렇게 나오지...? 

역시나 회사명은 밝히지 않음. 

어제 면접보면서 코테 얘기가 나와서 블로그에 올려도 되겠느냐고 여쭤봤고, 그 쪽이 커리어에 도움이 될 테니 올려도 된다+올리면 안 되는 거면 따로 언질을 주겠다는 답변을 받았음. 

문제가 이거 말고 하나 더 있는데 그거는 바이오파이썬을 쓸 일이 없고(이것도 안썼지만...), 뇌에 블루스크린이 떠서 이해하느라 시간이 좀 걸려서 아마 구현은 빨라야 내일... 잘못하면 다음주까지도 갈 것 같음. 그래도 두번째껀 실행하는 시간은 빨라서 망정이지 얘는... 저 모듈 Jupyter에서 불러오지도 못함 OTL 

전체 코드

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
Gene=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
     INFO_get=record.INFO.get('CLNSIG')
     Gene_get=record.INFO.get('GENEINFO')
     if INFO_get:
          Chromosome.append(record.CHROM)
          CLNSIG.append(str(INFO_get))
     else:
          Chromosome.append(record.CHROM)
          CLNSIG.append("No data")
     # CLNSIG에 키가 없으면 처리하는거
     if Gene_get:
          Gene.append(Gene_get)
     else:
          Gene.append("No data")
     # Gene도 몇개 없더라...
# 이거 자체는 간단하다. 키가 없으면 NaN으로 넣는다는 얘기.
Chr_series=pd.Series(Chromosome)
CLN_series=pd.Series(CLNSIG)
Gene_series=pd.Series(Gene)
# 제발 시리즈 한번에 나오게 해주소서
Gene_df=pd.DataFrame({"Chromosome":Chr_series, "Gene":Gene_series,"CLNSIG":CLN_series})
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서
Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
# Groupby로 시원하게 묶어보자
# 근데 이거 object라 정렬이 1, 2, 3, 4로 안된다...
Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자

문제

Q1. VCF파일 내 데이터 중 CLNSIC 데이터를 염색체별로 분류할 것
Q2. Pathogenic한 Gene Top 30 꼽기

공통 절차-분석을 위한 Dataframe 만들기

import vcf
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
     print(record)
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...

일단 저걸 하려면 VCF 파일을 열어야 한다. VCF file의 경우 ClinVar에서 받을 수 있고, ftp에 일주일마다 업그레이드가 된다. 아무튼... 저거는 바이오파이썬으로 불러오는 게 아니라 PyVCF라는 모듈이 또 있어서 설치를 해야 하는데 여기서부터 난관이었다. 

1. 설치를 했는데 Jupyter에서 못 불러온다 
2. 파이참에서도 못 불러온다 
3. IDLE은 되네? 
4. 어? 갑자기 파이참에서 되네? (Jupyter는 안됨)

참고로 코드 실행 속도가 정말 어마무시하게 느리다. 데이터가 방대하니 어쩔 수 없지... 

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
     print(record.INFO['CLNSIG'])
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...
     CLNSIG_frame=pd.DataFrame(record.INFO['CLNSIG'])

그리고 그냥 이렇게 Dataframe 만들려고 하면 Key error가 난다. 찾아보니 해당하는 키가 없으면 나는 에러라고 한다. 

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
     INFO_get=record.INFO.get('CLNSIG')
     if INFO_get:
          Chromosome.append(record.CHROM)
          CLNSIG.append(INFO_get)
          print(record.CHROM,INFO_get)
     else:
          Chromosome.append(record.CHROM)
          CLNSIG.append("NaN")
          print(record.CHROM,"No keys")
# 이거 자체는 간단하다. CLNSIG에서 Key error가 나는데, 없으면 NaN으로 쓴다는 얘기.

그래서 이렇게 if함수를 넣어줍니다. Key가 있으면 해당하는 값을 넣고, 없으면 Nan(현재는 No data)을 넣으라는 얘기. 

Chr_series=pd.Series(Chromosome)
CLN_Series=pd.Series(CLNSIG)
# 심보러시여 제발 시리즈 한번에 나오게 해주소서

에이 솔직히 시리즈는 잘 돼요... 면접관님 보시는데 무슨 심보러 드립을 치고 있냐 

CLN_df=pd.DataFrame({"Chromosome":Chr_series, "CLNSIG":CLN_Series})
print(CLN_df)
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서

데이터프레임...... (마른세수) 

중간에 개고생하긴 했으나 해결. 

Chromosome                    CLNSIG
0                     1  [Uncertain_significance]
1                     1           [Likely_benign]
2                     1                  [Benign]
3                     1           [Likely_benign]
4                     1  [Uncertain_significance]
...                 ...                       ...
1105836              MT  [Uncertain_significance]
1105837              MT  [Uncertain_significance]
1105838              MT                  [Benign]
1105839              MT  [Uncertain_significance]
1105840  NW_009646201.1                  [Benign]

[1105841 rows x 2 columns]

Process finished with exit code 0

이렇게 데이터프레임이 나왔다. (참고로 최종 Dataframe은 2번 문제도 하나의 데이터프레임으로 풀기 위해 Gene 컬럼도 추가되어있음)

CLNSIG 집계

CLN_df2=CLN_df.groupby('Chromosome').count()
print(CLN_df2)

CLNSIG
Chromosome            
1                89545
10               41010
11               65665
12               49168
13               29113
14               33996
15               38488
16               60126
17               76248
18               18665
19               51709
2               112429
20               20230
21               12514
22               21575
3                56286
4                34468
5                58692
6                48407
7                51289
8                36992
9                50025
MT                2838
NW_009646201.1       1
X                46316
Y                   46

참고로 저게 Object라 정렬이 저렇게 된다... 예전에 scatter plot 그릴때도 저랬는데, 그때는 결측값 지우고 object에서 numeric으로 바꿔서 해결 봤다... 그건 축 값이 순수하게 숫자라 그런거고 저거는 염색체라 지우는 게 안돼요... (주륵) 

CLN_df2=CLN_df.groupby('Chromosome').aggregate(['CLNSIG'])
print(CLN_df2)
# 염색체별로 묶어보자.

이거 뭐 계속 바꿔서 해봤더니 아 리스트라 안된다고!!! 하더라... 그 답이 전체 코드에 있다. 

CLNSIG.append(str(INFO_get))

정확히는 이거. 저 INFO에 CLNSIG 말고 데이터가 많은데 그걸 다 리스트로 불러와서 저런 사태가 일어난 것이라고 한다. 딕셔너리 키값은 튜플이나 문자열같이 절대불변인 것들로 해야 하기 때문. 그래서 if문 돌리고 리스트에 넣을 때 문자열로 변환했다. 

Gene_df2=Gene_df.groupby(['CLNSIG','Chromosome']).count()
print(Gene_df2)
# Groupby(CLNSIG, Chromosome 순)

Gene
CLNSIG          Chromosome      
No data         1             56
                10            29
                11            57
                12            18
                13            10
...                          ...
['risk_factor'] 6             30
                7             16
                8             21
                9             20
                X             10

[527 rows x 1 columns]

니네 근데 진짜 이걸로 해결되는거 실화냐고... (No data: CLNSIG 데이터가 없음)

Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
print(Gene_df2)
# Groupby

Gene
Chromosome CLNSIG                                   
1          No data                                56
           ['Affects']                            19
           ['Benign', '_other']                    1
           ['Benign']                          15397
           ['Benign/Likely_benign', '_other']      2
...                                              ...
Y          ['Benign']                              4
           ['Likely_benign']                       5
           ['Likely_pathogenic']                   2
           ['Pathogenic']                         28
           ['Uncertain_significance']              7

[527 rows x 1 columns]

난 저거 말고 염색체별로 볼 건데? 하면 순서를 바꾸시면 됩니다. 

Top 30 pathogenic gene

이거 Top 30 뽑는거는 간단하다. ascending에 0 주고 head(30)으로 뽑으면 된다. 

Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자

Chromosome  CLNSIG
Gene                                          
BRCA2:675                        13480   13480
BRCA1:672                        11565   11565
TTN:7273|TTN-AS1:100506866        9999    9999
APC:324                           8901    8901
NF1:4763                          7658    7658
TTN:7273                          7601    7601
TSC2:7249                         6392    6392
ATM:472                           6260    6260
MSH6:2956                         5591    5591
POLE:5426                         4939    4939
FBN1:2200                         4796    4796
RYR2:6262                         4358    4358
MSH2:4436                         4262    4262
RYR1:6261                         4063    4063
DMD:1756                          4063    4063
ATM:472|C11orf65:160140           3829    3829
PALB2:79728                       3788    3788
NEB:4703                          3706    3706
SYNE1:23345                       3486    3486
DICER1:23405                      3403    3403
MLH1:4292                         3354    3354
BRIP1:83990                       3350    3350
USH2A:7399                        3300    3300
PLEC:5339                         3267    3267
SMARCA4:6597                      3024    3024
PMS2:5395                         2949    2949
LDLR:3949                         2842    2842
TSC1:7248                         2776    2776
POLD1:5424                        2736    2736
CDH1:999                          2715    2715

이거는 유전자 전체 통틀어서 출력해주는거고 

Gene_df3=Gene_df.groupby(['Gene','CLNSIG']).count()
Gene_df3=Gene_df3.sort_values('Chromosome',ascending=0)
print(Gene_df3.head(30))
# 그냥 Gene으로 정렬하는 코드는 이거

Chromosome
Gene                       CLNSIG                                
BRCA2:675                  ['Uncertain_significance']        5467
APC:324                    ['Uncertain_significance']        4835
TTN:7273|TTN-AS1:100506866 ['Uncertain_significance']        3933
BRCA2:675                  ['Pathogenic']                    3590
TTN:7273                   ['Uncertain_significance']        3238
ATM:472                    ['Uncertain_significance']        3207
NF1:4763                   ['Uncertain_significance']        3101
POLE:5426                  ['Uncertain_significance']        3093
BRCA1:672                  ['Uncertain_significance']        3012
MSH6:2956                  ['Uncertain_significance']        2944
BRCA1:672                  ['Pathogenic']                    2892
TTN:7273|TTN-AS1:100506866 ['Likely_benign']                 2680
BRCA1:672                  ['not_provided']                  2654
BRCA2:675                  ['Likely_benign']                 2445
TTN:7273                   ['Likely_benign']                 2441
TSC2:7249                  ['Uncertain_significance']        2212
RYR2:6262                  ['Uncertain_significance']        2123
MSH2:4436                  ['Uncertain_significance']        2002
ATM:472|C11orf65:160140    ['Uncertain_significance']        1995
PALB2:79728                ['Uncertain_significance']        1956
BRIP1:83990                ['Uncertain_significance']        1917
NF1:4763                   ['Pathogenic']                    1895
RYR1:6261                  ['Uncertain_significance']        1757
DICER1:23405               ['Uncertain_significance']        1745
SYNE1:23345                ['Uncertain_significance']        1712
PLEC:5339                  ['Uncertain_significance']        1670
BRCA1:672                  ['Likely_benign']                 1617
NEB:4703                   ['Likely_benign']                 1617
NF1:4763                   ['Likely_benign']                 1596
RECQL4:9401                ['Uncertain_significance']        1577

Process finished with exit code 0

이렇게 하면 유전자별로 가장 높은 CLNSIG만 나온다. 

이놈들아 이것도 되면 좀 된다고 말좀 해줘... 

참고로 이거 어떻게 알았냐면 면접보는 회사에서 발표주제 중 하나가 저 두놈이었는데 찾다보니 NCBI에서 만든거네? -> Entrez에 있네? -> 비켜봐 시켜볼 게 있어(주섬주섬 파이참을 켠다) 가 된 거임. 

dbSNP

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com" # 내가 누구인지 말해주는 과정이 필요하다고...
# 이메일은 자기꺼 그냥 쓰세요
handle = Entrez.esearch(db="snp", term="EGFR", retmax="40" )
record = Entrez.read(handle)
print(record)

참고로 db에는 snp라고 써야지 dbsnp라고 쓰면 안된다. 

Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="snp", term="rs1050171", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])
IdList=list(record['IdList'])
# dbSNP+esearch
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="snp", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records)
# 일단 esummary 써서 다 가져왔다.

일단 워드클라우드 코드를 좀 응용해서 있는걸 다 털었는데 분량이 장난없어요... (마른세수) 이와중에 5000자 넘는다고 네이버 코드블록이 짤랐어 또... 

{'DocumentSummarySet': DictElement({'DocumentSummary': [DictElement({'SNP_ID': '1050171', 'ALLELE_ORIGIN': '', 'GLOBAL_MAFS': [{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}, {'STUDY': 'ALSPAC', 'FREQ': 'G=0.424235/1635'}, {'STUDY': 'Estonian', 'FREQ': 'G=0.435714/1952'}, {'STUDY': 'ExAC', 'FREQ': 'G=0.477223/57869'}, {'STUDY': 'FINRISK', 'FREQ': 'A=0.486842/148'}, {'STUDY': 'GnomAD', 'FREQ': 'G=0.483832/67753'}, {'STUDY': 'GnomAD_exomes', 'FREQ': 'G=0.476248/119743'}, {'STUDY': 'GoESP', 'FREQ': 'G=0.457558/5951'}, {'STUDY': 'GoNL', 'FREQ': 'G=0.380762/380'}, {'STUDY': 'HapMap', 'FREQ': 'A=0.405488/133'}, {'STUDY': 'KOREAN', 'FREQ': 'A=0.133219/389'}, {'STUDY': 'Korea1K', 'FREQ': 'A=0.126092/231'}, {'STUDY': 'MGP', 'FREQ': 'G=0.440075/235'}, {'STUDY': 'NorthernSweden', 'FREQ': 'G=0.483333/290'}, {'STUDY': 'PRJEB36033', 'FREQ': 'G=0.319149/30'}, {'STUDY': 'PRJEB37584', 'FREQ': 'A=0.153944/121'}, {'STUDY': 'PharmGKB', 'FREQ': 'G=0.344828/40'}, {'STUDY': 'Qatari', 'FREQ': 'G=0.439815/95'}, {'STUDY': 'SGDP_PRJ', 'FREQ': 'G=0.302083/116'}, {'STUDY': 'Siberian', 'FREQ': 'G=0.368421/14'}, {'STUDY': 'TOMMO', 'FREQ': 'A=0.159327/2670'}, {'STUDY': 'TOPMED', 'FREQ': 'G=0.482935/127828'}, {'STUDY': 'TWINSUK', 'FREQ': 'G=0.416667/1545'}, {'STUDY': 'Vietnamese', 'FREQ': 'A=0.285016/175'}, {'STUDY': 'ALFA', 'FREQ': 'G=0.427251/38873'}], 'GLOBAL_POPULATION': '', 'GLOBAL_SAMPLESIZE': '0', 'SUSPECTED': '', 'CLINICAL_SIGNIFICANCE': 'likely-benign,benign', 'GENES': [{'NAME': 'EGFR', 'GENE_ID': '1956'}, {'NAME': 'EGFR-AS1', 'GENE_ID': '100507500'}], 'ACC': 'NC_000007.14', 'CHR': '7', 'HANDLE': 'GENOMED,GRF,BUSHMAN,PJP,EVA_SAMSUNG_MC,SEATTLESEQ,CSHL,ILLUMINA,WUGSC_SSAHASNP,COMPLETE_GENOMICS,CSHL-HAPMAP,EGP_SNPS,KHV_HUMAN_GENOMES,EVA,EVA_FINRISK,BCMHGSC_JDW,EGCUT_WGS,JMKIDD_LAB,JJLAB,CGAP-GAI,CANCER-GENOME,EVA_UK10K_TWINSUK,ACPOP,TOMMO_GENOMICS,BGI,TISHKOFF,SGDP_PRJ,GNOMAD,SSMP,NHLBI-ESP,1000GENOMES,WEILL_CORNELL_DGM,ILLUMINA-UK,HUMANGENOME_JCVI,EVA-GONL,GMI,SI_EXO,EVA_EXAC,PACBIO,ENSEMBL,HAMMER_LAB,EVA_UK10K_ALSPAC,DDI,TOPMED,USC_VALOUEV,ABI,EVA_MGP,PERLEGEN,BIOINF_KMB_FNS_UNIBA,CLINSEQ_SNP,URBANLAB,KRGDB,PHARMGKB_PAAR-UCHI,HUMAN_LONGEVITY,LMM-PCPGM,LEE,OMUKHERJEE_ADBS,HGSV,CORNELL,SWEGEN,KOGIC,FSA-LAB,CPQ_GEN_INCA,EVA_DECODE', 'SPDI': 'NC_000007.14:55181369:G:A,NC_000007.14:55181369:G:C', 'FXN_CLASS': 'non_coding_transcript_variant,coding_sequence_variant,missense_variant,synonymous_variant,genic_downstream_transcript_variant', 'VALIDATED': 'by-frequency,by-alfa,by-cluster', 'DOCSUM': 'HGVS=NC_000007.14:g.55181370G>A,NC_000007.14:g.55181370G>C,NC_000007.13:g.55249063G>A,NC_000007.13:g.55249063G>C,NG_007726.3:g.167339G>A,NG_007726.3:g.167339G>C,NM_005228.5:c.2361G>A,NM_005228.5:c.2361G>C,NM_005228.4:c.2361G>A,NM_005228.4:c.2361G>C,NM_005228.3:c.2361G>A,NM_005228.3:c.2361G>C,NM_001346899.2:c.2226G>A,NM_001346899.2:c.2226G>C,NM_001346899.1:c.2226G>A,NM_001346899.1:c.2226G>C,NM_001346900.2:c.2202G>A,NM_001346900.2:c.2202G>C,NM_001346900.1:c.2202G>A,NM_001346900.1:c.2202G>C,NM_001346941.2:c.1560G>A,NM_001346941.2:c.1560G>C,NM_001346941.1:c.1560G>A,NM_001346941.1:c.1560G>C,NM_001346898.2:c.2361G>A,NM_001346898.2:c.2361G>C,NM_001346898.1:c.2361G>A,NM_001346898.1:c.2361G>C,NM_001346897.2:c.2226G>A,NM_001346897.2:c.2226G>C,NM_001346897.1:c.2226G>A,NM_001346897.1:c.2226G>C,NR_047551.1:n.1201C>T,NR_047551.1:n.1201C>G,NP_005219.2:p.Gln787His,NP_001333828.1:p.Gln742His,NP_001333829.1:p.Gln734His,NP_001333870.1:p.Gln520His,NP_001333827.1:p.Gln787His,NP_001333826.1:p.Gln742His|SEQ=[G/A/C]|LEN=1|GENE=EGFR:1956,EGFR-AS1:100507500', 'TAX_ID': '9606', 'ORIG_BUILD': '86', 'UPD_BUILD': '155', 'CREATEDATE': '2000/10/05 19:10', 'UPDATEDATE': '2021/04/26 11:38', 'SS': '1524562,4415417,14121776,16263880,17933589,23461029,24778961,43094183,52067157,65740167,66861655,74802379,77320037,85017386,86269798,93684223,98280473,104430036,105110075,112044216,116097662,142664329,143002862,159714800,159903775,162355178,166625159,203360709,212046487,223092668,233988791,240940108,279318753,285632487,293873711,342235809,479680970,482283259,485456160,490945938,491906494,534578333,560020602,654383955,779491194,781710098,834961318,836317040,974464379,984303175,1067488412,1074630855,1325217930,1431130856,1584052633,1593883812,1618262516,1661256549,1688745702,1711163743,1805015052,1927546725,1966658513,2024460166,2152655954,2294217561,2463237960,2634609883,2708323052,2736449362,2747825596,2853377601,3001152883,3023062922,3026026429,3347597637,3530770690,3530770691,3629823416,3632517137,3636855046,3646352982,3648637018,3669078603,3719737966,3734653562,3766593686,3785824290,3791125584,3796005564,3809753344,3824277165,3825524129,3825539990,3825720179,3830585782,3838783084,3844235500,3867317995,3914396600,3961507699,3984367631,3984367632,3984588810,3985298625,3986039617,3986382885,4746722701,5183240605,5236855147,5236861046,5237033464,5237196188', 'ALLELE': 'V', 'SNP_CLASS': 'snv', 'CHRPOS': '7:55181370', 'CHRPOS_PREV_ASSM': '7:55249063', 'TEXT': 'MergedRs=1050171', 'SNP_ID_SORT': '0001050171', 'CLINICAL_SORT': '1', 'CITED_SORT': '', 'CHRPOS_SORT': '0055181370', 'MERGED_SORT': '1'}, attributes={'uid': '58115520'})], 'DbBuild': 'Build211004-0955.1'}, attributes={'status': 'OK'})}

(대충 하나 뽑은게 이정도)

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="snp", term="rs1050171", retmax="40" )
record = Entrez.read(handle)
IdList=list(record['IdList'])
# dbSNP+esearch
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="snp", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records['DocumentSummarySet']['DocumentSummary'][0]['GLOBAL_MAFS'][0])
# 일단 esummary 써서 다 가져왔다. 근데 딕셔너리인데 왜 픽이 안될까...
# 그 전에 저렇게 꼭 4단픽까지 가야겠냐... 아니 진짜 너무한 거 아닙니까...

사실 5단이다. 딕셔너리가 아주 이중삼중이여... 

/home/koreanraichu/PycharmProjects/pythonProject/venv/bin/python "/home/koreanraichu/PycharmProjects/pythonProject/dbSNP in Entrez.py"
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}

다 똑같아보이는데 제대로 가져온 거 맞냐고요? 네, 맞습니다. 

ClinVar

dbSNP는 그 뭐지... single nucleotide polymoorphism, 그러니까 point mutation으로 인한 변이들을 기록해 둔 데이터베이스고 ClinVar는 거기에 대한 보고서에 접근할 수 있게 해 주는 public archive다. 

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="clinvar", term="EGFR", retmax="40" )
record = Entrez.read(handle)
print(record)
# clinVar+esearch

기본적으로는 이렇게 가져오고 

from Bio import Entrez
handle = Entrez.esearch(db="clinvar", term="L858R", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])
# clinVar+esearch
IdList=list(record['IdList'])
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="clinvar", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records)
# 일단 esummary 써서 다 가져왔다.

일단 다 쓸어오자... (저 for문은 wordcloud 만들 때 썼던 걸 응용한거다)

{'DocumentSummarySet': DictElement({'DocumentSummary': [DictElement({'obj_type': 'single nucleotide variant', 'accession': 'VCV001314051', 'accession_version': 'VCV001314051.', 'title': 'NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)', 'variation_set': [{'measure_id': '1304312', 'variation_xrefs': [], 'variation_name': 'NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)', 'cdna_change': 'c.2744T>G', 'aliases': [], 'variation_loc': [{'status': 'current', 'assembly_name': 'GRCh38', 'chr': '1', 'band': '1q44', 'start': '247444052', 'stop': '247444052', 'inner_start': '', 'inner_stop': '', 'outer_start': '', 'outer_stop': '', 'display_start': '247444052', 'display_stop': '247444052', 'assembly_acc_ver': 'GCF_000001405.38', 'annotation_release': '', 'alt': '', 'ref': ''}, {'status': 'previous', 'assembly_name': 'GRCh37', 'chr': '1', 'band': '1q44', 'start': '247607354', 'stop': '247607354', 'inner_start': '', 'inner_stop': '', 'outer_start': '', 'outer_stop': '', 'display_start': '247607354', 'display_stop': '247607354', 'assembly_acc_ver': 'GCF_000001405.25', 'annotation_release': '', 'alt': '', 'ref': ''}], 'allele_freq_set': [], 'variant_type': 'single nucleotide variant', 'canonical_spdi': 'NC_000001.11:247444051:T:G'}], 'trait_set': [{'trait_xrefs': [{'db_source': 'MedGen', 'db_id': 'CN517202'}], 'trait_name': 'not provided'}], 'supporting_submissions': {'scv': ['SCV002002472'], 'rcv': ['RCV001771282']}, 'clinical_significance': {'description': 'Uncertain significance', 'last_evaluated': '2021/01/14 00:00', 'review_status': 'criteria provided, single submitter'}, 'record_status': '', 'gene_sort': 'NLRP3', 'chr_sort': '01', 'location_sort': '00000000000247444052', 'variation_set_name': '', 'variation_set_id': '', 'genes': [{'symbol': 'NLRP3', 'GeneID': '114548', 'strand': '+', 'source': 'submitted'}], 'protein_change': 'L801R, L858R, L915R, L917R', 'fda_recognized_database': ''}, attributes={'uid': '1314051'})], 'DbBuild': 'Build211201-1855.1'}, attributes={'status': 'OK'})}

(마른세수)

그래도 얘는 짤리진 않았음... 하지만 딕셔너리가 대체 몇중이지 하면서 쌍욕박는 건 똑같습니다... 

결국 차근차근 겜코딩 하던 머리로 차근차근 딕셔너리 마킹해가면서 찾음... (딕셔너리 다중일때는 ['상위 키']['하위 키'] 이런 식으로 찾으면 된다)

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="clinvar", term="L858R", retmax="40" )
record = Entrez.read(handle)
# clinVar+esearch
IdList=list(record['IdList'])
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="clinvar", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records['DocumentSummarySet']['DocumentSummary'][0]['title'])
# 일단 esummary 써서 다 가져왔다.
# 야이 근데 무슨 이건 진짜 할 말을 잃었음

NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)
NM_000245.4(MET):c.2573T>G (p.Leu858Arg)
NM_000222.3(KIT):c.2585T>G (p.Leu862Arg)
NM_005228.5(EGFR):c.2573_2574delinsGT (p.Leu858Arg)
NM_005228.5(EGFR):c.2572_2573inv (p.Leu858Arg)
NM_005228.5(EGFR):c.2369C>T (p.Thr790Met)
NM_005228.5(EGFR):c.2573T>G (p.Leu858Arg)

EGFR쪽에 보면 Leu858Arg 있는데 이게 L858R, Thr790Met이 T790M이다. (A는 알라닌)

Notes: 파이썬 모듈이나 버전은 따라해도 되는데 OS는 본인 편한거 쓰는게 좋음. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 한데 맥이나 리눅스는 윈도우와 달라서 무턱대고 쓰기 좀 불편하다. 나 고딩때 애들 커맨드 컨트롤 헷갈리는것만 여러번 봤음 진짜. (맥에서는 복사가 Ctrl+C가 아니라 Command+C) 그리고 어지간한 개발툴은 윈도 맥 리눅스 다 지원한다. 카톡 빼고 그건 개발툴 아닌데 아 그래서 지원 안하나

​

참고로 카톡이 리눅스에서 안돼서(와인 깔면 된다는데 실패함...) 본인 본의아니게 공부하는동안 연락 안됐음... 리눅스에서는 디코로만 연락 가능합니다.

Laptop

Lenovo thinkpad E470 (현재 단종) +4GB RAM

​

내가 리눅스 쓸거다 근데 특히 우분투 쓸거다 그러면 우분투에서 여기다 우리꺼 깔면 잘됨 보장함 땅땅땅 있는 노트북 중에서 사는게 좋음.

​

일단 본인은 노트북 1호/3호와 달리 얘는 들고 다닐 거 감안해서 사양보다 휴대성에 중점을 뒀는데, 그러면 보통 그램을 삼. 1킬로도 안되는게 되게 가볍다더라고... 얘네 진짜 외계인 고문하나 근데 우분투 측에서 여기다 우리꺼 깔면 잘됨 땅땅땅 한 기종이 아니라 걸렀음... 예산 범위도 오버였고... 100만원 안으로 잡아놨는데 그램은 가격부터 오버여... 그램이 개발할 때 안 좋다는 얘기가 아니라 OS와 예산때문에 거른겁니다. 오해 ㄴㄴ요. 맥북도 예산때문에 걸렀음. 맥이 솔직히 편하긴 한게 고등학생때 써봤거든... 본인 놋북빼고 패드 폰 다 애플기기고. 근데 비싸요...

​

리눅스는 공짜기때문에 free-DOS 사서 직접 설치하면 됨. 부팅디스크 만드는 법도 찾아보면 나오고 설치법도 나오고... 뭐 우분투의 경우 쓰다보면 터미널에서 돌리는게 좀 많긴 한데 설치까지 터미널로 하는 건 아니고 윈도처럼 GUI로 돌아감. (리눅스도 종류가 많다... 리눅스 민트 센토스 뭐... 우분투는 그 중 하나.)

​

이 노트북은 원래 내장된 램이 4기가인데, 개발하다보면 생각보다 램이 많이 필요하다고 한다. 당장 노트북 1호가 게이밍이었던 이유도 MATLAB으로 영상 처리를 해야 하기 때문이고(영상처리 생각보다 고사양 필요함)... 그래서 4기가 엑스트라로 박았는데도 파이참 좀 버벅거림... (spyder나 jupyter는 제깍제깍 열린다) 저거 살 때는 파이썬 안했고 자바 했는데 이클립스로 돌아가야 하다 보니 저사양에서 버벅거리는 걸 몇 번 보기도 했고. 놋북 CPU는 i5고... 본인은 아무래도 잘알까지는 아니라서 CPU 사양을 크게 신경쓰는 건 아닌데 일단 i 뒤에 붙는 숫자는 세대가 같으면 큰게 좋고, 본인은 일단 셀러론 거름. 그래서 마지노선이 i3~i5단.

​

내가 셀러론에 안좋은 추억이 있어... 게임하다 랙걸려서 사리적립만 몇년 했는지 몰라...OTL 무슨 메이플 하는데 랙이 걸리냐고...

OS

Ubuntu 20.04 LTS(2022년 8월 21일 현재는 22.04 LTS)

​

커널은 올렸는데 LTS가 새로운게 안 풀린건지 내가 저장소를 따로 추가해야 하는 건지... 리누스형 이러지 말자 진짜

​

리눅스는 맥과 함께 유닉스 계열 OS고 악마의 명령어 sudo rm -rf/가 먹힌다 우분투의 경우 GUI가 있지만 터미널을 통해 CUI로 해야 하는 게 많다. 대표적인 예가 글꼴 뭉텅이로 설치할 때 복사할 글꼴이 있는 폴더로 가서 sudo cp *.ttf /usr/share/fonts 같은 거(cp=copy). 윈도우의 관리자 권한으로 실행같은 게 sudo인데 생각보다 이거때문에 터미널에서 해야 하는 게 좀 있다. 윈도보다 편한 부분도 있는데 윈도보다 불편한 부분도 존재함. 맥도 그렇잖아요? 대신 우분투 ~설치 이런걸로 찾아보면 정보는 많이 나오는 편.

​

그리고 얘 종특인지는 모르겠는데 한영키를 한영키가 아니라 오른쪽 알트로 인식해서 그거 셋업하고 한글 언어팩 깔아야 한글 입력 됨. 나는 우분투 깔면 제일 먼저 하는게 한글 깔고-mozc 설정하고(일어 자판)-무선랜 드라이버 추가하는거임. 그 뒤에 부차적으로 필요한 거 깔고 구축함.

​

대신 터미널에서 하는거에 적응하면 터미널 하나 켜놓고 거기서 다운로드 설치 실행 이동 복사 삭제가 다 가능하다. 본인 깃헙 세팅 터미널에서 하고 깃크라켄 까니까 알아서 셋업 인식하더만. 파이썬도 터미널에서 깔고 필요한 모듈도 일단은 터미널에서 깔았고 아나콘다 베이스도 다 터미널에서 돌리고(conda update 이런거) 파이참도 바로가기 없어서 터미널에서 켰음. (파이썬 IDLE은 바로가기가 있는데 왜 파이참만...)

​

아 그리고 바탕화면 그림은 바꿀 수 있는데 윈도처럼 바탕화면에 따라 작업표시줄 색 바꾸는건 안됨. 리누스씨 이 기능 빨리 넣어줘요

언어

Python 3.8

언어는 가급적 최신꺼 쓰는 편. (R도 깔려는 있음)

Module

-numpy

-sympy(+mpmath)

-scipy

-pandas

-matplotlib(부를때는 matplotlib.pyplot)

-Biopython

-sklearn

-tensorflow

-anaconda

(2022년 8월 21일 현재 request, Flask, pymongo 추가됨)

​

아 NetworkX도 깔았음. 그래프(이산수학 그래프... 막대그래프는 matplotlib이나 seaborn 쓰세요) 그려준대서...

​

아나콘다는 파이썬에서 소환하는 모듈은 아니고 뭐라고 해야되나... 저거 베이스로 깔아야 하는 모듈이나 IDE가 또 있음. sympy도 아나콘다 끼고 깔았던 것 같고... 심파이는 방정식 해주는건가 그렇고 scipy는... 뭐지 저거? 이봐요

​

넘파이는 배열 만드는거고(가끔 바이오파이썬 하다보면 마려운 그거 맞음) 판다스는 데이터프레임, 그러니까 표 만드는겁니다. csv나 엑셀파일 불러와서 표 만드는것도 가능하고 불러와서 계산하고 피벗테이블마냥 그룹별로 묶는것도 됨. 쓰다보면 진짜 편합니다. Biopython은 바이오파이썬이고 matplotlib은 그래프 그려주는거고... (seaborn이 좋다고 하는데 안써봐서 모름)

​

sklearn, tensorflow는 각각 머신러닝, 딥러닝 관련인데 이쪽은 뭐... 쓸 일이 있을 지 모르겠음. 사용법은 둘째치고 머신러닝이나 딥러닝이 어려워서 내장 데이터도 제대로 못 불렀음. ChEMBL에서 가져온 데이터 처리해서 써보려고 했더니 ChEMBL ID가 문자라서 인식을 못한다는데 ID도 빼줘야되냐고... OTL

IDE

-pycharm

-Jupyter(lab, notebook)

-spyder

-R studio(교육때문에 깔았지만 딱히 쓰지는 않음)

(2022년 8월 21일 현재 VScode도 설치해서 사용중)

​

spyder는 에리본씨 추천인데 한영키 안먹어서 뭐야 이거 왜이래요 했던거 드디어 해결봤고... ㄹㅇ 디자인 내 취향이었는데 거기서 막혀서 개고생했고... 바이오파이썬 할 때 고정적으로 쓰는건 파이참, 간단한건 IDLE(파이썬 깔면 IDLE도 따라옴), 간단한 코딩인데 IDLE에서 안될 것 같은건(예: 판다스 쿼리검색) jupyter로 돌린다. 저거 랩이랑 놋북 다 깔아야됩니다. (주파이터 랩 놋북 세개 깔아야됨)

​

리눅스는 바로가기고 뭐고 공통적으로 터미널에서 소환하면 알아서 열리는데(jupyter notebook 치면 실행됨) 파이참은 pycharm 치면 안되고 파이참 설치 경로 찾아가서 쉘 스크립트 실행해줘야 한다(설치도 다른 IDE랑 좀 다름). 즉 다른 IDE는 터미널 키고

​

spyder

jupyter lab

jupyter notebook

​

이렇게 치면 바로 열리는데 파이참은

​

cd /파이참 설치 경로/bin

sh pycharm.sh

​

해야 열린다. 그나마 탭키 누르면 자동완성 돼서 망정이지... (IDLE은 바로가기도 있고 진짜 간단한거 아니면 IDLE은 안씀) 그래서 편한거는 스파이더랑 주파이터가 편함. 쳐주면 알아서 실행도 되겠다... 근데 spyder는 설치를 pip가 아니라 anaconda 경유해서 해야 하고 한글 입력하려면 셋업을 또 따로 해줘야 한다. 코드는 영어 아닌가 주석도 영어로 다시게요? 파이참은 바로가기 아이콘이 없는 문제가 있지만 일단 설치하고 쓰는 데 크게 지장이 있는 건 아니다. 뭐야 아이콘 돌려줘요

​

셋업 면에서는 파이참 주파이터가 편하지만 실행 면에서는 스파이더 파이참이 편하다. 그리고 똑같은 다크테마여도 스파이더 다크테마가 내 취향임. (Jupyter는 다크모드 없는거같던디...) 주파이터는 랩보다 노트북을 더 많이 쓰는 편이다.

​

참고로 파이참은 터미널에서 인스톨하는 게 아니라 파이참 페이지에서 tar파일 받아서 압축풀고 거기 경로 들어가서 쉘 스크립트 실행하는 게 끝이다. 그래서 버전업 할거면 새 파이참 tar파일 받아서 압축 풀고 거기 가서 쉘 스크립트 돌리면 된다.

​

대부분 Jupyter 주피터로 읽는데 본인은 주파이터로 읽는다. python을 파이썬이라고 하잖음. 개발자는 뭐라고 읽나 궁금하네

​

여담이지만 처음에 Jupyter notebook이라고 하길래 인공지능 돌릴 용도로 노트북을 따로 준비해야 하나 했음... 나중에 교육 실습에 주파이터 나와서 어? 저거? 하다가 찾아봤더니 Jupyter notebook이래서 충격받았다. 그게 그거였냐고 근데 인공지능은 서버에서 돌아가는데 노트북을 별도로 준비해야됨? 아니 귀하신 분이라 별도로 모셔뒀나 했지

기타 툴

-ClustalW

-MUSCLE

-Gitkraken

-ATOM

(2022년 8월 21일 현재 MongoDB가 추가되었음)

​

노션은 바이오파이썬 Q&A에서도 언급했듯 리눅스에서는 설치가 아니라 웹노션을 쓴다. (윈도우와 맥 버전은 설치 가능) 그래서 쓰고 있지만 설치한 건 아니므로 생략. 코딩한 거 정리할 때는 노션이 짱이다. 코드블럭이 있어...

​

ClustalW와 MUSCLE은 MSA 툴인데 Biopython에서 돌리려면 깔아야 한대서 깔았더니 안돼서 걍 터미널로 돌리는 중. 역시 spyder나 jupyter 시리즈처럼 터미널에 clustalw, muscle 치면 실행된다. (muscle은 명령어 어떻게 치라고 나온다)

​

Gitkraken은 본진에도 설치한건데, 깃헙 관련된 건 다 이쪽에서 처리한다. 저장소에는 파일 이동만 하고 풀 커밋 푸시는 깃크라켄 통해서 하는 중. 이것도 okky에서 추천받은 툴이다. 리눅스에서 풀 커밋 푸시는 터미널로 해도 되는데 그거보단 얘가 편함. 이사람 깃헙 있음

​

ATOM에서 코드 쳐도 된다는데 본인은 오타 교정+자동완성이 되는 IDE가 편하다보니 그냥 글 쓸 때만 쓴다.

Q&A지만 자문자답이다. 어쨌든 질답은 맞음 

Q1. MSA의 그 clustalW랑 MUSCLE은 어찌됐나요? 

A1. 그거 둘다 깔아야 됩니다. OS 박고 경로 박아서 돌리는 거 나오긴 했는데 트라이 해보려고 했더니 윈도 기준이네... 머슬은 경로 박아서 해봤는데 커맨드만 나와서 MSA는 터미널로 돌리고 있습니다. 리눅스는 clustalw와 MUSCLE 둘 다 일단 설치해두면 터미널에서 돌릴 수 있습니다. (py파일로 아웃풋 설정도 가능)

Q2. 일부 건너뛴 챕터들이 있던데...? 

A2. 실습용 자료 구하기가 빡세거나(얘네도 다 안올려줌...) MSA처럼 할 수 없는 여건인 경우 건너뜁니다. 그래서 16 17(이건 하려고 했는데 wordcloud 하느라 시간 다 잡아먹음) 건너뛰고 케그 했지... 

Q3. 하면서 제일 재미없었을 때+제일 재밌었을 때는? 

A3. 실습할 거리가 많으면 많을수록 꿀잼입니다. clustering은 실습할 거리도 없는데다가 초장에 이론 나와서 노잼... 

Q4. 쿡북 보면서 할 때 제일 빡쳤던 부분은? 

A4. 처음에는 설명도 막 디테일하게 하고 이거이거이거 다돼! 이러더니 뒤로 갈수록 코드에 오타나고 모듈 빼먹고... 예제대로 따라하면 에러가 반기는... 심지어 그 에러 내가 찾아서 다 수정했어...... 쓰다가 귀찮았는지는 모르겠는데 귀차니즘이 거기서 터지면 안되지 이놈들아... 

Q5. 전체적인 개발 환경은? 

A5. 이건 따로 글 올려드리겠음. 여기다 풀면 분량이 길어져서;; 

Q6. 가끔 보다보면 응용편이 있던데 그건 대체 어떻게 떠올리신거죠 

A6. 화장실에서 X싸다가요. 이게 측간신 버프인가 그건가 그냥 코드 보면 떠오를 때도 있습니다. for문으로 돌리는거 while로 돌려본다던가... 

Q7. 저도 할수 있을까요? 

A7. 저도 파이썬 잘 못하는데 하잖음. 

Q8. IDE 배경색 뭐쓰세요?

A8. 원래 밝은 색 썼었는데 지금은 걍 깜장배경 씁니다. 개발자들은 깜장배경이 국룰인가... Jupyter는 흰 바탕입니다. 

Q9. 코드블럭에 입력하는 거 귀찮지 않으세요? 

A9. 노션은 코드블럭이나 있으니까 그나마 복붙하면 장땡이지 미디움은 코드블럭도 없어서 스크린샷 찍어서 붙입니다. 네이버도 초창기에 코드블럭 찾기 전에 쓴 건 스크린샷이었는데... 근데 결과 길면 복붙도 귀찮음 진짜... 

Q10. 코딩은 어디서 하세요?

A10. 방해받기 싫어서 스터디카페 가서 합니다. 군자역 열공다방 짱짱맨. 근데 밥 먹을데가 없음... 그리고 집중하다보면 끼니도 잊어버려서 음료수로 배 채우다가 집 가서 밥 먹습니다... 그래서 아침먹고 저녁먹고 야식먹고 근데 왜 살이 빠졌지...

Q11. 코딩할 때 중요한건 뭐라고 생각하십니까?

A11. 막 한번에 완성하려는 생각은 버리시고, 코드 입력하기 전에 대충 어떻게 돌아가면 좋을지 한번 그려보세요. 기능 구현하는 방법은 모르면 찾으면 되지만 그것도 무슨 기능을 어떤 식으로 구현할 지 알아야 찾아요. 아, 한가지 더. 큰 그림을 그릴 때는 일단 차근차근 해봅시다. 작은 단위(의 고정된 수치)로 먼저 해보고 점차 규모를 키워나가던가 일반화해야지 첫빠따에 일반화하려고 하면 에러가 여러분을 반깁니다. 내가 그래서 에러랑 많이 만났음 

Q12. 개발자한테 한마디?

A12. 문의메일 왜 안받냐. 심지어 내가 썬더버드로 보내서 메일 지메일인데. (리눅스 썬더버드에 지메일 연결함) 

Q13. 운영체제는 꼭 리눅스여야 하나요? 

A13. 윈도 편하면 그거 쓰세요. 리눅스는 tkinter인가 안됨. 저도 처음엔 그랬는데 이거 터미널 적응하는거랑 키세팅도 빡셉니다. 깔자마자 한글 입력 안돼서 그거 셋업해야되고 저는 놋북 무선랜 드라이버도 따로 넣어줘야되고 개귀찮음 진짜. 그리고 어지간한 개발툴은 윈도 맥 리눅스 지원하니까 편한거 쓰면 됩니다. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 하지만, 그냥 많이 쓴다고 무작정 그거 쓰는것보단 자기 편한거 쓰는 게 짱입니다. 디코도 지원하는 리눅스를 지원 안 하는 카톡은 뭘까 사랑해요 디코 

Q14. Notion 좋아요?

A14. 전직장 개발자님이 전파한건데 이거 ㄹㅇ 와 진짜 개편합니다. 달력 생성도 돼서 스케줄러에 메모장으로도 쓰고 있어요. 실행도 빠릿빠릿하고 (인터넷 연결은 필요하지만) 싱크도 빠르고... 패드에서 그림 올리면 PC에 금방 뜨고, 그림을 옮기거나 글을 쓰면 그것도 빠릿빠릿하게 나타나고... 리눅스에서는 설치형은 없고 웹노션이 있습니다. 노션도 지원하는 웹을 지원 안하는 카톡은 뭘까 

Q15. 워드프레스에는 따로 연재 안하시나요?

A15. 거기는 이미지도 HTML로 올려야돼서 귀차니즘이 배가됩니다. 간단한거(별찍기)는 이미지 많이 안 올려도 돼서 워드프레스에 올리지만... 

Q16. 책 추천해주세요!

A16. Biopython 책 산 건 있는데 한 번 보고 안읽은듯... 

Q17. 이거 재밌나요?

A17. 저는 재밌습니다. Wordcloud도 재밌었고... 근데 균 학명에 미역김치는 진짜 어디서 나온걸까... 정답: 미역김치 

Q18. 코드는 따로 올려두는 공간이 있나요? 여기에는 py파일이 없는 거 같은데...

A18. 깃헙에 올려두고 있습니다. 그거 덕분에 파이썬이 매트랩 분량 이겼음. 

쿡북 분량 개짧음 진짜 이거보다 짧을수가 없음. 

KEGG? 

https://www.genome.jp/kegg/

Kyoto Encyclopedia of Genes and Genomes의 준말. 그렇다, 이름에 교토가 들어간 걸 보면 아시겠지만 일제 DB다. 

여기가 메인페이지

여기가 KEGG brite. 본인이 자주 가는 곳이다. 생각보다 쏠쏠한 정보가 많고 KEGG brite의 경우 golden standard dataset이라고 해서 야먀니시가 인공지능 학습시킬 때 쓴 데이터셋(GPCR, 효소, 핵 리셉터, 트랜스포터) 분류별로 약물 타겟을 알려준다. 

Parsing

파싱할거면 일단 파일이 당신 컴퓨터에 있어야 한다. 

급한대로 이거라도 쓰자. (아님 밑에 쿼링에서 갖고오든가...)

from Bio.KEGG import Enzyme
records = Enzyme.parse(open("/home/koreanraichu/ec_5.4.2.2.txt"))
record = list(records)[0]
print(record)

전체를 뽑고 싶으면 이렇게 하고 

from Bio.KEGG import Enzyme
records = Enzyme.parse(open("/home/koreanraichu/ec_5.4.2.2.txt"))
record = list(records)[0]
print(record.classname)

['Isomerases;', 'Intramolecular transferases;', 'Phosphotransferases (phosphomutases)']

특정 부분만 보고 싶으면 이렇게 한다. 근데 쟤 뭐 하는 효소냐... 

아, parse 친구 read도 있다. 데이터가 단식이면 리드, 뭉텅이면 파스. 

Querying

쿼리에 ing면 쿼링인가... 

Enzyme

from Bio.KEGG import REST
from Bio.KEGG import Enzyme
request = REST.kegg_get("ec:7.1.2.2")
# 참고: ATP synthase의 EC번호이다
open("ec_7.1.2.2.txt", "w").write(request.read())
records = Enzyme.parse(open("ec_7.1.2.2.txt"))
record = list(records)[0]
print(record)

참고: EC 7.1.2.2는 ATP synthase다 

파일이 생긴 것을 볼 수 있다. (본인은 Pycharm으로 해서 파이참 경로에 생겼다)

Pathway

사실상 KEGG의 근본은 pathway data다. 정말 세상 온갖가지 pathway가 다 있음. 진짜로. 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
print(human_pathways)

뭔지 모르겠을 땐 일단 다 불러오는 게 국룰이다. 밑에 보여주겠지만 pathway만 써도 불러오는 건 된다. 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(entry)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드

['path:hsa03410', 'path:hsa03420', 'path:hsa03430']

다 좋은데 그걸 꼭 hsa로 뽑아야겠냐... 알아서 찾아라 이거지 이럴거면 KEGG가지 뭐하러 여기서 찾냐 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(description)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드

['Base excision repair - Homo sapiens (human)', 'Nucleotide excision repair - Homo sapiens (human)', 'Mismatch repair - Homo sapiens (human)']

아니 그럼 출력을 description으로 바꾸면 되잖아. 와 천잰데? 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "*sis" in description:
        repair_pathways.append(description)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드

[]

와일드카드는 안먹거나 *이 아닌 다른 기호인 듯 하다. 

Gene

Pathway에서 이어진다. 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(entry)
# 조건부로 Repair pathway만 보는 코드
# 여기까지는 파일에 있었으니까 다들 아시리라 믿고...
repair_genes = []
for pathway in repair_pathways:
    pathway_file = REST.kegg_get(pathway).read()  # query and read each pathway

    # iterate through each KEGG pathway file, keeping track of which section
    # of the file we're in, only read the gene in each pathway
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()  # section names are within 12 columns
        if not section == "":
            current_section = section

        if current_section == "GENE":
            gene_identifiers, gene_description = line[12:].split("; ")
            gene_id, gene_symbol = gene_identifiers.split()

            if not gene_symbol in repair_genes:
                repair_genes.append(gene_symbol)

print(
    "There are %d repair pathways and %d repair genes. The genes are:"
    % (len(repair_pathways), len(repair_genes))
)
print(", ".join(repair_genes))

일단 코드가 긴데 

1. Pathway 찾아주는 코드(그래서 entry로 쳐야 먹힌다)
2. 에서 Gene 섹션 찾아서 출력해준다. (hsa 번호가 맞는 pathway에서 gene을 찾아준다)

There are 3 repair pathways and 78 repair genes. The genes are:
OGG1, NTHL1, NEIL1, NEIL2, NEIL3, UNG, SMUG1, MUTYH, MPG, MBD4, TDG, APEX1, APEX2, POLB, POLL, HMGB1, XRCC1, PCNA, POLD1, POLD2, POLD3, POLD4, POLE, POLE2, POLE3, POLE4, LIG1, LIG3, PARP1, PARP2, PARP3, PARP4, FEN1, RBX1, CUL4B, CUL4A, DDB1, DDB2, XPC, RAD23B, RAD23A, CETN2, ERCC8, ERCC6, CDK7, MNAT1, CCNH, ERCC3, ERCC2, GTF2H5, GTF2H1, GTF2H2, GTF2H2C_2, GTF2H2C, GTF2H3, GTF2H4, ERCC5, BIVM-ERCC5, XPA, RPA1, RPA2, RPA3, RPA4, ERCC4, ERCC1, RFC1, RFC4, RFC2, RFC5, RFC3, SSBP1, PMS2, MLH1, MSH6, MSH2, MSH3, MLH3, EXO1

그래서 이렇게 나온다. 

근데 솔직히 인간적으로 저것만 찾으면 개노잼 아니냐. 그래서 준비했다. 

이건 산화적 인산화(Oxidative phosphorylation)인데, ATP를 얻기 위한 과정 중 하나이다. 여러분들이 일반적으로 알고 있는 ATP 생산 절차는

​

1. 우리가 먹은 포도당을 해당과정(Glycolysis)을 통해 피루브산 두 개를 만들고

2. TCA cycle(크렙 회로나 시트르산 사이클이라고도 한다) 들어가고(여기서부터 미토콘드리아 영역)

3. Oxidative phosphorylation을 통해 미토콘드리아가 ATP를 만든다(여기서 겁나 나온다)

​

이건데... Oxidative phosphorylation은 마지막 절차이다. 보통 TCA cycle부터 미토콘드리아 영역이고 피루브산이 TCA cycle로 가면 돌아올 수 없는 강을 건너게 된다. (피루브산 상태에서 포도당신생합성(Gluconeogenesis)을 돌리기도 한다. 번역명 왜이래요. 번역 간지나네)

​

여담이지만 과당도 해당과정 중간에 낄 수 있다. (Galactose도 낀다...젖당이었나 갈락이었나 아무튼 낌)

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
Oxidative_pathways = []
Oxidative_pathways_name =[]
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "Oxidative" in description:
        Oxidative_pathways.append(entry)
        Oxidative_pathways_name.append(description)
print("Entry no: ",Oxidative_pathways,"description is: ",Oxidative_pathways_name)
# 조건부로 Repair pathway만 보는 코드
# 여기까지는 파일에 있었으니까 다들 아시리라 믿고...
Oxidative_genes = []
for pathway in Oxidative_pathways:
    pathway_file = REST.kegg_get(pathway).read()  # query and read each pathway

    # iterate through each KEGG pathway file, keeping track of which section
    # of the file we're in, only read the gene in each pathway
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()  # section names are within 12 columns
        if not section == "":
            current_section = section

        if current_section == "GENE":
            gene_identifiers, gene_description = line[12:].split("; ")
            gene_id, gene_symbol = gene_identifiers.split()

            if not gene_symbol in Oxidative_genes:
                Oxidative_genes.append(gene_symbol)

print(
    "There are %d pathways and %d genes. The genes are:"
    % (len(Oxidative_pathways), len(Oxidative_genes))
)
print(", ".join(Oxidative_genes))

Entry no:  ['path:hsa00190'] description is:  ['Oxidative phosphorylation - Homo sapiens (human)']
There are 1 pathways and 134 genes. The genes are:
ND1, ND2, ND3, ND4, ND4L, ND5, ND6, NDUFS1, NDUFS2, NDUFS3, NDUFS4, NDUFS5, NDUFS6, NDUFS7, NDUFS8, NDUFV1, NDUFV2, NDUFV3, NDUFA1, NDUFA2, NDUFA3, NDUFA4, NDUFA4L2, NDUFA5, NDUFA6, NDUFA7, NDUFA8, NDUFA9, NDUFA10, NDUFAB1, NDUFA11, NDUFA12, NDUFA13, NDUFB1, NDUFB2, NDUFB3, NDUFB4, NDUFB5, NDUFB6, NDUFB7, NDUFB8, NDUFB9, NDUFB10, NDUFB11, NDUFC1, NDUFC2, NDUFC2-KCTD14, SDHA, SDHB, SDHC, SDHD, UQCRFS1, CYTB, CYC1, UQCRC1, UQCRC2, UQCRH, UQCRHL, UQCRB, UQCRQ, UQCR10, UQCR11, COX10, COX3, COX1, COX2, COX4I2, COX4I1, COX5A, COX5B, COX6A1, COX6A2, COX6B1, COX6B2, COX6C, COX7A1, COX7A2, COX7A2L, COX7B, COX7B2, COX7C, COX8C, COX8A, COX11, COX15, COX17, CYCS, ATP5F1A, ATP5F1B, ATP5F1C, ATP5F1D, ATP5F1E, ATP5PO, ATP6, ATP5PB, ATP5MC1, ATP5MC2, ATP5MC3, ATP5PD, ATP5ME, ATP5MF, ATP5MG, ATP5PF, ATP8, ATP6V1A, ATP6V1B1, ATP6V1B2, ATP6V1C2, ATP6V1C1, ATP6V1D, ATP6V1E2, ATP6V1E1, ATP6V1F, ATP6V1G1, ATP6V1G3, ATP6V1G2, ATP6V1H, TCIRG1, ATP6V0A2, ATP6V0A4, ATP6V0A1, ATP6V0C, ATP6V0B, ATP6V0D1, ATP6V0D2, ATP6V0E1, ATP6V0E2, ATP6AP1, ATP4A, ATP4B, ATP12A, PPA2, PPA1, LHPP

물량 쩌네. 

Appendix. 뭐 딴거 없어요? 메인페이지에 개많던데? 

나도 궁금해서 KEGG 메인에 있는거 다 돌려봤음. 

from Bio.KEGG import REST
human_pathways = REST.kegg_list('brite').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('pathway').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('module').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('enzyme').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('genome').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('glycan').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('compound').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('reaction').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('disease').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('drug').read()
print(human_pathways)

from Bio.KEGG import REST
human_pathways = REST.kegg_list('network').read()
print(human_pathways)

저기서 되는거 Genes빼고 어지간한건 다 된다. 

1. KEGG Brite
2. KEGG Pathway
3. KEGG Module
4. KEGG Enzyme
5. KEGG Compound
6. KEGG Genome Gene은 안되는데 Genome은 되네? 
7. KEGG Glycan 
8. KEGG Reaction
9. KEGG Network
10. KEGG Drug
11. KEGG Disease

Appendix 2. 응용편

from Bio.KEGG import REST
Disease = REST.kegg_list('disease').read()
# 가져올 수 있는 것: brite, pathway, genome(gene은 안됨), module, enzyme, glycan, compound, reaction, network, drug, disease
disease_leukemia = []
for line in Disease.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "leukemia" in description:
        disease_leukemia.append(description)
print(disease_leukemia)

다 불러올 수 있으면 아까 그 쿼링 코드 써서 검색도 되나 궁금해서 찾아봤음. 

from Bio.KEGG import REST
drug = REST.kegg_list('drug').read()
# 가져올 수 있는 것: brite, pathway, genome(gene은 안됨), module, enzyme, glycan, compound, reaction, network, drug, disease
drug_list = []
drug_name=input("찾고 싶은 약과 관련된 것을 입력해주세요: ")
for line in drug.rstrip().split("\n"):
    entry, description = line.split("\t")
    if drug_name in description:
        drug_list.append(description)
print(drug_list)

이건 나름 업그레이드 버전인데 Drug에서 사용자가 입력한 항목을 찾아준다. (와일드카드는 안 먹는다)

이것도 Project Wordcloud의 초기 코드이다. Text with wordcloud로, 이 당시에는 영문만 됐었다. (한글도 되긴 한데 접속사가 안떼짐)

from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
from PIL import Image
import numpy as np  
# Summon module
text = []
input_text= input("wordcloud로 만들 텍스트를 입력해주세요. ")
text.append(input_text)  
yorn = input("더 추가할 텍스트가 있나요?")
# 일단 입력을 받는다. (없으면 n, 있으면 n 말고 다른거)
while yorn != "n":
    yorn = input("더 추가할 텍스트가 있나요? ")
    if yorn == "n":
        text = ''.join(text)
    else:
        input_text= input("wordcloud로 만들 텍스트를 입력해주세요 ")
        text.append(input_text)
# 추가로 입력이 있을 경우 입력이 '없다'고 할 떄까지 입력을 받고 입력이 더 이상 들어오지 않으면 join한다. (입력 받을때마다 join하면 에러남)
image = np.array(Image.open("/home/koreanraichu/600px-793Nihilego.png"))
# 이미지를 부르고
font_path = '/usr/share/fonts/Maplestory Light.ttf'
wordcloud = WordCloud(font_path = font_path,background_color="#ffffff",colormap="PuBu_r",width = 800, height=800, mask=image)
wordcloud = wordcloud.generate_from_text(text)
# 그리고
plot.figure(figsize=(15,15))
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# 출력

전체 코드(n or N으로 줘봤는데 안먹힘...)

결과(릴리요 도감 설명) 

Project wordcloud의 극초반 코드. 

from Bio import Entrez
# 논문 긁어올 때 필요한거
from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
# Wordcloud 그릴 때 필요한거
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="pubmed", term="Arabidopsis[title]", retmax="15")
record = Entrez.read(handle)
IdList=record['IdList']
# 일차적으로 Pubmed에서 논문을 찾을 수단인 PMID를 입수한다.
# 날짜검색이 안돼서 키워드에 '애기장대'가 있는 논문 10개를 찾음.
for i in range(len(IdList)):
    handle = Entrez.esummary(db="pubmed", id=IdList[i], retmode="xml")
    records = Entrez.parse(handle)
    title=[]
    for record in records:
        title.append(record['Title'])
# 타이틀 뽑았다.
# wordcloud에서 뺄 단어들(접속사 이런거)
font_path = '/usr/share/fonts/truetype/nanum/BinggraeMelona.ttf'
wordcloud = WordCloud(stopwords=STOPWORDS,   font_path = font_path, width = 800,
    height = 800,background_color="#ffffff")
wordcloud = wordcloud.generate_from_text(record['Title'])
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# wordcloud generate

코드 구성은 크게 

1. 모듈(wordcloud용 모듈과 Entrez 접속용 모듈)
2. Entrez에서 긁어오기(PMID-title)
3. Wordcloud 그리기

로 나누어진다. 

주석에는 10개로 되어 있었는데, 이게 40개를 가져오다가 서버 에러가 나서 규모를 좀 줄였다. (현재 15개) 코드를 일반화함에 따라 중간에 for문도 바꼈다. 

워드클라우드의 경우 Stopword라고 해서 접속사나 감탄사같은 것들을 다 빼야 하는데(in, at, the같은 거) 영문의 경우 STOPWORDS에 저장되어 있다. 

1차 결과물(10개)

2차 결과물(15개, 폰트 변경)

근데 내가 원하는 결과는 대충 위에꺼인 거 같다... 

from Bio import Entrez
# 논문 긁어올 때 필요한거
from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
# Wordcloud 그릴 때 필요한거
import numpy as np
from PIL import Image
# 이미지 마스킹할 때 필요한거
image=np.array(Image.open("/home/koreanraichu/alice.png"))
font_path = '/usr/share/fonts/dovemayo.otf'
# 마스킹할 이미지
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="pubmed", term="EGFR[title]", retmax="20")
record = Entrez.read(handle)
IdList=record['IdList']
# 일차적으로 Pubmed에서 논문을 찾을 수단인 PMID를 입수한다.
# 숫자를 20개로 늘려서 시간이 좀 걸린다.
title=[]
for i in range(len(IdList)):
    handle = Entrez.esummary(db="pubmed", id=IdList[i], retmode="xml")
    records = Entrez.parse(handle)
    for record in records:
        title.append(record['Title'])
title=''.join(title)
# 타이틀 뽑았다.
# 아 이거 근데 묶어줘야되더라...
wordcloud = WordCloud(stopwords=STOPWORDS,font_path = font_path, width = 800,height = 800,
                      background_color="#ffffff",colormap='spring',mask=image)
wordcloud = wordcloud.generate_from_text(record['Title'])
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# wordcloud generate

PS. 마스크 기능 추가하는 김에 확인해봤는데 하나밖에 안되길래 리스트업한거 join해버림. 

PS2. +코드의 wordcloud = wordcloud.generate_from_text(record['Title'])부분을 wordcloud = wordcloud.generate_from_text(Title)로 수정해줘야 제대로 나온다. 

뭐야 왜이렇게 작아요 

하는김에 김치로 찾음. 

Arabidopsis+2021[Years]로 준 것