SAMtools BCFtools 1.8버전 설치

다운로드 링크

http://www.htslib.org/download/

설치 절차

1. 위 링크에서 우분투 운영체제 wget을 통해 samtools, bcftools tar.gz2 파일을 다운받는다.

2. tar -xvf [파일명] 을 통해 압축을 해제한다.

3. 이 때, 압축 해제는 아무데나 해도된다.

cd samtools-1.x    # and similarly for bcftools and htslib
./configure --prefix=/where/to/install
make
make install

4. 이후 /where/to/install/bin에 samtools binary 파일이 생성된다.

5. /usr/bin 에 심볼릭 링크를 만들어준다.

sudo ln -s /where/to/install/bin /usr/bin/samtools

6. samtools --version  명령어를 통해 잘 설치 되었는지 확인한다.

BCFTools도 이와 마찬가지 방법이다. 폴더 명만 바꿔주면된다.

부모-자식 데이터를 통한 유전율 추정

본 포스팅은 부모-자식 데이터를 통해 heritability를 추정하는 것을 포스팅한다. 아래 데이터는 parent와 child의 height를 나타내는 수치이다. 옥수수의 키라고 생각하자. 알고 싶은 것은 child의 키가 유전적 요소로 설명되는 정도이다. 우리는 이 수치를 heritability(h^2)라는 값으로 표현한다. 데이터는 이곳에서 다운로드 할 수 있다.

\[ h^2 = \frac{V(G)}{V(Y)}\]

library(data.table)

## Warning: package 'data.table' was built under R version 3.4.3

load("C:\\workspace\\r\\genetics\\data\\mod7_data.R")
  head(Galton)

##     parent child
  ## 1 71.40788  61.7
  ## 2 68.57945  61.7
  ## 3 64.33681  61.7
  ## 4 62.92260  61.7
  ## 5 62.21549  61.7
  ## 6 67.16524  62.2

데이터에 약간의 random error를 주어 아래의 그래프를 통해 표현해보자.

plot(Galton$parent+rnorm(928,0,.1),Galton$child+rnorm(928,0,.1),xlab="Father's Height",ylab="Child's Height")
  abline(v=mean(Galton[,1]),col=2,lty=3,lwd=3)
  abline(h=mean(Galton[,2]),col=2,lty=3,lwd=3)
  abline(a=0,b=1,col=2,lty=3,lwd=3)

summary(lm(Galton$child~Galton$parent))

## 
  ## Call:
  ## lm(formula = Galton$child ~ Galton$parent)
  ## 
  ## Residuals:
  ##     Min      1Q  Median      3Q     Max 
  ## -7.8050 -1.3661  0.0487  1.6339  5.9264 
  ## 
  ## Coefficients:
  ##               Estimate Std. Error t value Pr(>|t|)    
  ## (Intercept)   36.87187    1.98827   18.55   <2e-16 ***
  ## Galton$parent  0.45700    0.02909   15.71   <2e-16 ***
  ## ---
  ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## 
  ## Residual standard error: 2.239 on 926 degrees of freedom
  ## Multiple R-squared:  0.2105, Adjusted R-squared:  0.2096 
  ## F-statistic: 246.8 on 1 and 926 DF,  p-value: < 2.2e-16

회귀계수가 0.457이다. 이 때, 유전율은 0.457*2 = 약 92%로 추정된다. 왜 그런지 살펴보기 위해 우선 beta의 추정치를 아래와 같이 회귀식의 양변에 Cov(, X)를 취해 구해보자.

\[ Y = \alpha + \beta X \] \[ Cov(Y, X) = \beta Cov(X,X) \] \[ \hat{\beta} = \frac{Cov(Y,X)}{Var(X)} \]

Phenotype Y에 대해 아래와 같은 일반적인 Phenotype에 대한 모델을 세우자. 이 때, A,D,C,E 각각의 컴포넌트들에 대한 것은 나중에 살펴보고 i, j 두 개체의 Covariance만 살펴본다. A : additive effect D : dominance effect C : shared environment effect E : residuals \[ Y = A + D + C + E\] r은 relatedness, k2는 genome상의 IBD2의 비율, a는 shared environment를 나타내는 수치이다. 만약 shared environment이면 a=1, 아니면 a=0이다. 

이 Covariance 모델은 어디까지나 가정이다. 이 가정 하에서 각각의 Variance를 추정하고 이를 통해 heritability를 추정하게 된다.

\[ E[Cov(Y_i, Y_j)] = rVar(A) + k_2Var(D) + aVar(C)\]

이제 현재 데이터에서 위 식을 적용하자. 이 때, r=0.5, k2=0, a=0을 가정하자. (부모-자식 간의 IBD2=0이므로 k2=0이다. ) 그러면 랜덤으로 i번째 부모-자식 pair를 골랐을 때, Covariance 모델은 아래와 같다. (P : Parent, C : Child의 약자)

이러면 Var(D)가 사라지기 때문에 우리는 쉽게 Var(A)를 데이터를 통해 추정할 수 있다.

\[ E[Cov(P_i, C_i)] = 0.5Var(A)\] \[ Var(A) = 2E[Cov(P_i, C_i)] \]

print(var(Galton$parent))

## [1] 6.389121

print(var(Galton$parent, Galton$child)) # Sample Covariance

## [1] 2.919806

따라서 위 데이터에는 D(Dominance effect가 없으므로 아래와 같이 유전율을 추정할 수 있다.) \[ h^2 = \frac{V(G)}{V(Y)} = \frac{V(A)}{V(P)} \]

따라서 2.91*2 = 5.82 = Var(A) 이다. Var(P) = 6.39 (Phenotype의 분산을 V(P)로 추정) 이므로 이 모델하에서 추정한 유전율은 아래와 같다.

\[ \frac{V(A)}{V(P)} = 5.82/6.39 = 0.91 \]

print(paste0("Heritability : ", 2*var(Galton$parent, Galton$child) / var(Galton$parent)))

## [1] "Heritability : 0.913992906289406"

따라서 부모-자식 데이터에 이러한 모델을 통해 유전율을 추정할 때 유전율은 회귀계수의 2배를 해주면 된다.

참고

https://jvanderw.une.edu.au/AGSCcourse.htm

인공지능, 머신러닝은 더 이상 현실 세계와 동떨어진 이야기가 아니다. 헬스케어 분야를 리드하는 기업들은 인공지능, 머신러닝 기술들을 자신의 분야에 이미 적용시키고 있다. 이는 AI 기술들이 헬스케어 분야에 실용적으로 적용되고 있다는 것을 뜻한다. 회사들은 AI를 활용하여 자신들의 서비스를 개선하고, 수익을 증대하며 발생할 수 있는 리스크를 줄인다. 본 포스팅에서는 AI가 헬스케어 분야에 적용되어 불러올 수 있는 사회적 가치와 응용 사례를 살펴보려고 한다.

AI FOR GOOD

기술은 그 자체로도 중요하지만 사회적으로 그 기술을 활용하여 무언가의 이득을 얻을 때 그 가치가 가장 크다고 생각한다. 예를 들어, 기존에 해결할 수 없던 문제를 해결한다거나, 어려웠던 문제를 더 쉽게 만든다거나, 돈을 더 벌어준다더가 하는 것이다. AI FOR GOOD 이라는 말은 이러한 여러 사회적 가치 중 "공익"에 관한 것이다. 인공지능이 헬스케어분야에 적용하는 것의 가장 큰 가치는 바로 공익의 증진이다. AI의 헬스케어 적용은 기존의 헬스케어 분야에서 기존의 해결하기 어렵거나, 해결하는데 시간이 오래걸렸던 문제를 해결하고, 가속화하여 공익의 증가에 기여한다. 

AI-based Medicine을 통한 의학 분야의 돌파구

Accenture에 따르면, 현재 가장 큰 가치를 갖고 있는 기술들은 로봇 보조 수술, 가상 환자 돌봄 서비스(virtual nursing assistant), 행정 보조 기술 등이다. 이는 빅데이터, 인공지능 등의 키워드가 뜨기전부터 의학 분야에 적용되어온 기술들이다. 하지만, 의영상 분야도 큰 가치를 갖고 있다. 왜냐하면 취약 계층의 경우, 질병을 진단 받을 때, 의영상 진단 기술의 도움을 받기 힘들다. 전문가의 진단은 비싸며, 일정 수준 이상의 수요를 필요로하기 때문이다. 이러한 취약 계층에 AI 기반의 의영상 진단 기술이 적용될 경우, 공익에 큰 기여를 할 수 있다. 본 포스팅에서는 AI 기반 의료의 실제 사례를 몇 가지 살펴보려고 한다. 

AI 기반 의료의 실제 사례

- 개발 도상국 국가에서 X-RAY를 결핵 환자 발견

의영상 이미지에서 무언가의 패턴을 찾는 것은 가장 유망하며, 실제로 AI가 많이 적용되고 있는 분야이다.  많은 연구자들이 이러한 AI 기반의 의영상 진단 기술을 연구하고 있다. (원논문) 이러한 기술들이 개발 도상국 국가에서 활용된다면, 영상의학 전문의의 도움 없이도 결핵에 걸린 환자를 발견해 낼 수 있다. 또, 한 국가의 여러 지역에서 이런 기술을 적용한다면, 결핵 취약 지역을 발견해 낼 수 있다. 만약 결핵 예방, 치료 사업을 한다고 하면, 이러한 방식의 접근 법을 통해 정책의 우선순위를 결정할 수도 있다.

- 응급실에서의 환자 뇌출혈 발견

이스라엘의 회사인 MedyMatch와 IBM Watson은 응급실에서 환자의 뇌출혈을 확인하는데 인공지능을 활용하고 있다. 

- 피부 사진을 통한 암 진단

암진단에도 종종 AI가 응용된다. CT, MRI, X-ray 등이 암진단을 위한 의영상 기술로 적용되고 있지만, 이 기술들을 통해 암 진단을 하지 못할 경우, 치명적이다. 암환자가 아닌데 암환자로 진단하는 경우 (false positive, 1종 오류) 보다 암환자인데 암환자로 진단하지 못한 경우(false negative, 2종 오류), 환자의 생명에 영향을 주는 치명적인 판단 오류라고 볼 수 있다. 또한 false negative를 위해 false positive를 늘리는 소위 과잉진료로 인한 문제도 심각하다. 이는 환자에 생명에 영향을 주진 않지만, 잘못된 진단 결과로 인해 의료 비용을 높이고 환자의 심적 부담을 증가시킨다. 따라서 정확한 암 진단 기술이 중요하며, 이 진단 과정에 AI가 적용되어 진단의 속도와 정확성을 높이려는 시도가 이루어지고 있다. 

Stanford의 AI 팀은 피부 이미지를 통해 피부암을 진단하는 AI 기술을 내놓았다. 이 AI 기술은 21명의 피부과 의사들과 비교해 진단 능력이 떨어지지 않았다. 

- CT 사진을 통한 폐암 진단

Enlitic은 회사는 폐 CT 사진을 통해 폐암을 진단하는 AI 기술을 개발하였다. 딥러닝을 활용해 병변의 존재와 그 위치를 파악하는 것이다. 회사는 정확도가 영상의학전문의의 50% 이상이라고 주장한다. 

역사상 많은 기술들이 과대 선전되어 왔다. 누군가는 인공지능 기술이 과대 포장되고 실체는 별로 없는 기술이라고 생각할 수 있다. 또한 누군가는 인공지능이 실제로 산업에 적용되어 가치를 불러올 수 있을까하는 의문을 제기할 수 있다. 인공지능이 인간을 완전히 대체하는 것은 어려울지도 모른다. 오로지 데이터 기반 접근법이라는 문제로 인해, 적절한 수의 훈련 데이터 없이는 동작할 수 없으며, 과적합의 가능성이 있다. 또한 AI의 경우 기기별 차이를 감안하여 모든 기계에 대해 robust한 모델을 만드는 것이 아직까지는 완벽하지 않은 경우가 많다. 하지만, 확실한 것은 특정한 제한된 범위의 태스크에서는 AI가 인간보다 나을 수 있다는 것이다.  현재 의료 산업 분야에 활발히 적용되고 있는 인공지능 기술과, 인공지능 스타트업에 투자되고 있는 금액의 규모가 인공지능의 의료분야에서의 가능성을 보여준다.

참고

https://sigmoidal.io/artificial-intelligence-and-machine-learning-for-healthcare/

GATK4 설치

이곳에서 원하는 버전 zip파일 다운로드

https://github.com/broadinstitute/gatk/releases

예를 들어, 4.0.1.1 버전이면 다운로드를 받은 후,

unzip gatk-4.0.1.1 

명령어를 통해 압축을 풀면 설치가 끝난다.

만약에 binary 파일을 다른 곳에 위치시키고 싶다면(예를 들어 /usr/bin), gatk 파일을 카피해서 원하는 폴더 (/usr/bin)으로 옮긴다.

cp gatk /usr/bin

하지만 binary 파일이 다른 파일들과 떨어져있을 경우, 기존 폴더의 경로를 명시적으로 지정해주어야한다. 이는 환경변수를 지정함으로써 할 수 있는데 환경변수 $ GATK_LOCAL_JAR 을 gatk-package-4.0.1.1-local.jar 파일의 위치로 지정하면 된다.

예를 들어, export GATK_LOCAL_JAR=/home/gatk-package-4.0.1.1-local.jar 을 통해 환경변수를 지정하면, binary 파일이 떨어져있어도 실행 가능하다.

Solar 설치

홈페이지

http://www.solar-eclipse-genetics.org/

설명

Heritability 분석을 위한 프로그램이다.

1. 이 링크(https://www.nitrc.org/frs/?group_id=558)에 들어가 아래 압축 파일을 다운받는다.

SOLAR Eclipse General Version 8.1.1

2. 다운 받은 압축 파일의 압축을 푼다.

tar -zxvf 파일명.tar.gz -C /원하는위치

3. 아래 명령어로 인스톨 한다. 

./install_solar <solar-base> <solar-script-dir>

파라미터에 대한 설명 (readme 파일 열어보면 나옴)

<solar-base> is the directory name where the bulk of this release of SOLAR

is to be installed.  To keep it separate from other releases, the release

number should be included in this directory name.  For example, we use:

/opt/appl/solar/8.1.1

<solar-script-dir> is the existing directory to which the 'solar' script

that starts SOLAR is to be put.  /usr/local/bin is a typical choice.

This should be in the PATH of everyone who will use solar.

You can move this script later if desired.  If you already have an older

script named 'solar' there, you must rename or delete it first.

이 때 solar-base에는 실제 실행 파일이 위치한 곳, solar-script-dir에는 실행파일을 연결할 파일의 경로를 써준다. 일반적으로 solar-script-dir은 /usr/local/bin 파일에 위치하도록 하고, solar-base는 아무 곳에나 놔두어도 되지만 /usr/local/bin에 위치해도 크게 상관은 없다. 

***********************************************************

Ignore error messages (if any) from the tests below

Ignore error messages (if any) from the tests above

***********************************************************

    *** Successful Installation ***

SOLAR has been installed with the command name solar

The new documentation directory is /mnt/c/Users/JuYoungAhn/Downloads/solar/8.1.1/doc

For command line editing, give command ./install_rlwrap

root@DESKTOP-F14K9HF:/mnt/c/Users/JuYoungAhn/Downloads/solar_linux# solar

/mnt/c/Users/JuYoungAhn/Downloads/solar/8.1.1/bin/rlwrap: 1: /mnt/c/Users/JuYoungAhn/Downloads/solar/8.1.1/bin/rlwrap: Syntax error: "(" unexpected

SOLAR Eclipse version 8.1.1 (General), last updated on January 29, 2018

Copyright (c) 1995-2018 Texas Biomedical Research Institute

Enter help for help, exit to exit, doc to browse documentation.

위와 같은 메시지가 뜨면 설치 완료된 것이고, /usr/local/bin에 script 파일을 넣은 경우 환경변수를 설정하지 않아도 커맨드 창에 solar를 치면, solar가 실행된다.

어떻게 세포는 자신들이 필요한 유전자를 발현 시킬까?

위에서는 유전자가 단백질로 변화하는 과정 속에서 유전자 조절이 언제 이루어지는지 알아보았습니다. 그렇다면 구체적으로 어떻게 세포는 원하는 유전자만 발현시켜 원하는 기능을 수행할 수 있을까요? 이는 각각의 세포들이 제각기 다른 전사 조절 인자들을 갖고 있기 때문입니다. 이러한 조절 인자들은 전사를 촉진(activate)시키기도 하고 방해(repress)하기도 합니다. 전사 조절 인자는 단백질이며, 이 전사조절인자도 유전자들이 코딩하고 있습니다.

일반적으로 전사는 RNA 중합효소가 소위말하는 프로모터 시퀀스에 결합함으로써 시작됩니다. 이 시퀀스는 전사 시작 지점에 근접한 upstream 방향에 존재합니다. (5' 쪽 방향) 하지만 downstream 방향에도 존재할 수 있습니다. (3' 쪽 방향) 비교적 최근에 연구자들은 인핸서 시퀀스라는 것을 발견하였습니다. 인핸서는 전사 조절 단백질들이 결합할 수 있는 결합 위치(binding site)를 제공함으로써 전사에서 중요한 역할을 하는 시퀀스입니다. 이 인핸서의 전사 조절 단백질이 붙게되면 염색질 구조가 변화되어 RNA 중합효소나 조절 단백질의 결합을 촉진하거나 억제하는 역할을 하게됩니다. 이 때 염색질 구조를 open chromatin structure라고 하는데, 이는 유전자 전사가 활성화 된다는 것과 연관이 있습니다. 반면 염색질의 구조가 빽빽한 경우, 전사가 억제되어 있는 것과 연관이 있습니다.

몇몇 조절 단백질은 다양한 유전자의 전사에 영향을 미칩니다. 이는 전사 조절 단백질 결합 위치 (regulatory protein binding site 혹은 transcription factor binding site라고도 함)가 다양한 곳에 존재하기 때문입니다. 결과적으로 조절 단백질은 다양한 유전자의 걸쳐 다양한 역할을 합니다. 조절 단백질의 역할은 단지 어떤 한 유전자의 발현을 조절하는 것이 아니라 다양한 유전자의 발현을 조절합니다. 이것이 각각의 세포가 한 번에 많은 수의 유전자를 조절할 수 있는 하나의 메커니즘입니다.

Figure 2: 전사의 조절

활성 단백질이 인핸서 시퀀스에 결합하게되면 RNA 중합효소를 활성화 시키는 프로모터 부분에 단백질을 끌어들일 수 있고, 이로 인해 전사가 촉진된다. DNA는 위와 같이 굽어져 활성자 단백질은 RNA 중합효소의 활동을 중재하는 다른 단백질들과의 상호작용을 하게 된다.

우선 원핵 세포에 대해서 알아보면, 원핵생물에서 조절 단백질은 종종 영양소의 이용가능성에 의해 조절됩니다. 이는 박테리아와 같은 생물이 환경 조건에 반응하여 전사 패턴을 빠르게 조절할 수 있게 합니다. 덧붙여, 원핵생물의 조절 부위는 프로모터와 가깝게 위치합니다.

Figure 3: 프로모터 주위에서의 전사 조절

특정 단백질은 RNA 중합효소에 간섭을 함으로써 전사를 조절한다. 불활성화 상태로 존재하는 억제 단백질(repressor)은 다른 분자에의해 활성화 될 수 있으며, 활성화된 상탱서 operator라 불리는 부위에 결합하여 RNA 중합효소가 프로모터에 결합하는 것을 방해한다. RNA 중합효소가 프로모터 부위에 결합하여야 전사가 개시됨으로 이는 전사를 효과적으로 억제한다. 

활성자(activator)는 프로모터 주위에 조절 부위에 결합하여 RNA 중합효소의 활동을 촉진합니다. 억제자(repressor)는 조절부위에 결합하여 RNA 중합효소의 결합을 방해합니다. 

진핵생물의 유전자 발현 조절은 원핵생물에 비해 복잡합니다. 위에서 기술한 것과 마찬가지로 기본적으로는 프로모터 주위 조절 부위에 결합하는 활성자, 억제자에 의해 RNA 중합효소의 활동량이 조절되어, 전사 과정에서의 유전자 조절이 일어납니다. 하지만 진핵생물에서는 그 이상의 많은 수의 조절 단백질이 존재하며, 조절 단백질의 결합 부위는 프로머터와 멀리 떨어져 있는 경우도 많습니다. 이로 인해 유전자 발현의 조절이 더욱 유연하게 됩니다.

Figure 4: 많은 수의 전사 조절 인자

전사 조절 인자들은 각각 다른 역할을 갖고 있다. 위의 3가지 전사 조절인자는 Mediator 복합 단백질과 각각 다르게 상호작용하여 전사를 조절한다. 각각의 조절단백질이 있고 없음, 또한 이들이 어떻게 조합되느냐에 따라 유전자의 발현이 달라진다. 진핵 생물의 유전자 발현의 특징은 이러한 복합한 조절 과정을 통해 같은 유전자라도 다르게 번역될 수 있다는 것이다.

다른 세포 유형은 특징적인 전사 조절인자를 갖고 있습니다. 다세포 생물에서 다른 세포는 각각 다른 조절 인자들의 조합을 갖고 있습니다. 이로인해 각기 다른 기능을 하는 다양한 세포가 생성되고, 기능할 수 있는 것입니다.

Figure 5: 전사 조절인자가 세포 유형을 결정한다.

세포 발달 과정에 따라 전사 조절 인자가 달라지는 것을 나타내는 그림

세포 유형의 다양성은 다른 전사조절인자들의 활동에 의한 것이다.

참고 : Nature Education Gene expression

Chip-seq 데이터를 통한 binding motif의 분석

Chip-seq 을 통해 대략적인 protein binding site (혹은 histon modification)의 시퀀스를 알고, 이 시퀀스들 중에서 비슷한 시퀀스를 찾아내어 해당 protein에 specific한 시퀀스가 무엇인지를 알아내는 분석이다. 이를 수행하는 bioinformatics Tool은 매우 많으며, [Evaluating tools for transcription factor binding site prediction] 논문에서 좋은 평가를 받았던 rGADEM을 이용하여 binding motif를 찾는 분석을 해보도록 한다. rGADEM은 GADEM 방법을 r로 구현한 것이며, GADEM은 motif를 찾는데 Genetic 알고리즘과 EM 알고리즘을 결합하여 사용한다. 자세한 방법론은 본 포스팅에서 다루지 않는다. 본 포스팅은 이 튜토리얼을 참고하였다.

# rGADEM Tutorial
# reference : https://bioconductor.org/packages/devel/bioc/vignettes/rGADEM/inst/doc/rGADEM.pdf

source("https://bioconductor.org/biocLite.R")
# biocLite("rGADEM", lib="C:/Users/Documents/R/win-library/3.4")
# biocLite("BSgenome.Hsapiens.UCSC.hg19", lib="C:/Users/Documents/R/win-library/3.4")
# biocLite("DNAStringSet", lib="C:/Users/R/win-library/3.4")

# install.packages("RCurl")

library(RCurl)
library(rGADEM)
library(BSgenome.Hsapiens.UCSC.hg19)

rGADEM은 bioLite 패키지에 있으며, 위와 같이 설치할 수 있다. rGADEM, BSgenome..., DNAStringSet 3개의 bioLite 패키지를 설치한 후, 이와 의존성이 있는 RCurl 패키지까지 설치하면 라이브러리 로드가 완료된다.

데이터

Test_100.bed

Test_100.fasta

rGADEM을 설치하면 위 예제 데이터가 패키지 안에 기본으로 들어있다. 이 예제 데이터를 통해 실습을 해보도록 한다.

chipseq 데이터로 binding motif를 분석하기까지의 pipeline은 QC, peak calling 등의 과정을 거쳐야한다. 하지만 이번 분석에서는 peak calling까지 끝난 데이터로 분석을 한다고 가정한다. peak calling이 완료되면 protein binding이 많이 일어난 부분의 시퀀스를 얻어낼 수 있다. 이 데이터는 주로 BED 혹은 FASTA 파일로 주어진다.

1. 데이터가 BED 포맷인 경우

BED <- read.table("C:\\Test_100.bed",header=FALSE,sep="\t")
BED <- data.frame(chr=as.factor(BED[,1]),start=as.numeric(BED[,2]),end=as.numeric(BED[,3]))

BED

rgBED <-IRanges(start=BED[,2],end=BED[,3])
Sequences <- RangedData(rgBED,space=BED[,1])

rgBED
Sequences

> BED

chr     start       end
1  chr1 145550911 145551112
2  chr1 112321064 112321265
3  chr1 120124183 120124384
4  chr1 114560780 114560981
5  chr1   8044002   8044203
6  chr1 154708057 154708258
7  chr1  21537833  21538034
8  chr1 117197889 117198090
9  chr1  22547577  22547778
10 chr1 170982215 170982416
11 chr1 203938925 203939126
12 chr1 198607893 198608094
13 chr1 244838696 244838897
14 chr1  22613354  22613555
15 chr1  36725231  36725432

50 chr1  35024860  35025061

> rgBED
IRanges object with 50 ranges and 0 metadata columns:
           start       end     width
       <integer> <integer> <integer>
   [1] 145550911 145551112       202
   [2] 112321064 112321265       202
   [3] 120124183 120124384       202
   [4] 114560780 114560981       202
   [5]   8044002   8044203       202
   ...       ...       ...       ...
  [46]  24574361  24574562       202
  [47]  46110458  46110659       202
  [48] 114377432 114377633       202
  [49] 152657429 152657630       202
  [50]  35024860  35025061       202

> Sequences
RangedData with 50 rows and 0 value columns across 1 space
       space                 ranges   |
    <factor>              <IRanges>   |
1       chr1 [145550911, 145551112]   |
2       chr1 [112321064, 112321265]   |
3       chr1 [120124183, 120124384]   |
4       chr1 [114560780, 114560981]   |
5       chr1 [  8044002,   8044203]   |
6       chr1 [154708057, 154708258]   |
7       chr1 [ 21537833,  21538034]   |
8       chr1 [117197889, 117198090]   |
9       chr1 [ 22547577,  22547778]   |
...      ...                    ... ...
42      chr1 [221975293, 221975494]   |
43      chr1 [ 59123669,  59123870]   |
44      chr1 [ 76297843,  76298044]   |
45      chr1 [188010599, 188010800]   |
46      chr1 [ 24574361,  24574562]   |
47      chr1 [ 46110458,  46110659]   |
48      chr1 [114377432, 114377633]   |
49      chr1 [152657429, 152657630]   |
50      chr1 [ 35024860,  35025061]   |

2. GADEM 실행

gadem <- rGADEM::GADEM(Sequences,verbose=1,genome=Hsapiens)

interation을 돌며 GADEM 알고리즘이 실행된다.

총 2개의 motif 가 생성되었다.

> length(gadem@motifList)
[1] 2

이 중에서 첫 번째 motif는 아래와 같이 생겼다.

plot(gadem@motifList[[1]])

> gadem@motifList[[1]]@alignList
[[1]]
An object of class "align"
Slot "seq":
[1] "GCCCCGACCTCTTATCTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 100

Slot "pval":
[1] 2.209965e-09

Slot "fastaHeader":
[1] 38


[[2]]
An object of class "align"
Slot "seq":
[1] "ACCCCAACCTCTTATCTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 58

Slot "pval":
[1] 2.365196e-09

Slot "fastaHeader":
[1] 38


[[3]]
An object of class "align"
Slot "seq":
[1] "ACCCCAACCCCTTATTTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 121

Slot "pval":
[1] 5.472386e-09

Slot "fastaHeader":
[1] 38

GADEM 알고리즘은 이 motif를 설명할 수 있는 총 26개의 시퀀스를 찾았다. 이 중 3개의 시퀀스가

GCCCCGACCTCTTATCTCTG
ACCCCAACCTCTTATCTCTG
ACCCCAACCCCTTATTTCTG

이다. 이러한 비슷한 시퀀스를 tfbs 라고 하며, tf protein이 이 곳에 달라붙어 regulatory role을 수행하게 된다.

[[26]]
An object of class "align"
Slot "seq":
[1] "CCTTCCTTCCTTTCTTCCTT"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 200254533

Slot "end":
[1] 200254734

Slot "strand":
[1] "-"

Slot "seqID":
[1] 36

Slot "pos":
[1] 193

Slot "pval":
[1] 0.0001878598

Slot "fastaHeader":
[1] 36

위 결과처럼 - strand의 결과도 볼 수 있는데, 이 경우에는 BED파일에서 reverse strand를 보면 된다. 하지만 위 결과는 자동으로 reverse strand를 보여준다. http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:200254533,200254734를 확인하면 위 시퀀스의 reverse인 aaggaagaaaggaaggaagg 를 확인할 수 있다.

3. INPUT이 FASTA 파일인 경우

# 2. Input이 FASTA 파일일 경우
path <- system.file("extdata/Test_100.fasta", package="rGADEM")
FastaFile <- paste(path,sep="")
Sequences <- readDNAStringSet(FastaFile, "fasta")
Sequences

Fasta파일의 경우 위와 같이 전처리 하면 되며, 이후의 절차는 BED 파일과 동일하다.

> Sequences
  A DNAStringSet instance of length 49
     width seq                                                                                                      names               
 [1]   202 CTACAGCTGTTCCTTGTCATCAGCCTGGGGGTGGGTAGTATTTTGATCTTA...GTCCCCATGTGACTGTAGGGTTCCTGAAACCTGGCAGGCCACTCTGCTTG FOXA1 _ 1
 [2]   202 TTATTCTGATGTGGTTTTGCGGTTATACAGTAAGCAGCACTGCTTATGTGG...ACACCAGAGGCCACCAGGACCAGAATGTTTACCAATGTAGGCAGTCACTA FOXA1 _ 10
 [3]   202 AAAGGAGAAACACAGCCAAATAATAAAACAATATCTTCTGTAAGTAAAGAG...CAAAGGTGCAAAGCCCTCTTTCAGATCCATCTCCACCATTTCCCTTCAGG FOXA1 _ 11
 [4]   202 TGTACCCCCCCAATATTTCATGATATTTACATGTTTGCATAGTACTTTCCT...CACAGGCAGAGTCAACATTGGAACTCGGAACACTGAGCCTGGCATTCCAA FOXA1 _ 12
 [5]   202 TTTAAGACTGCCACCTGAAATCAAGTCCAGTGGCTGTTCCTGTCTTCCCGT...TTATTTAACCTTTTGCCATTAGTTTACAAGAAGACGGGTTGAGCGAGTCA FOXA1 _ 13
 ...   ... ...
[45]   202 ATTAGTTCATGCCAGGCAGGGTTTGACCAAAACCTTGGACTCACTCCTGTC...TGCTGCTGCTAAGCCAGTGAGATAGATGGTGGTTTGGAAACACCCTCATG FOXA1 _ 5
[46]   202 CCAGAGCCACCACAGCCAGGCCTCTGCGGCCAGTGTCAACAAGGGGCCAGG...TGGGGTTTTGTTTACACACCTAGGGTCACCTGAAAACACCTGCTGTTTAA FOXA1 _ 6
[47]   202 CAGATCAGAGCCTGGGAGCGGGCCATGTGCAGCTGGATGGGCAGCTGGAGG...GCCCACCCTCCTGCCTTCCCAGCCAAGGAGGGAAGGGAGCTGTGGGAGGA FOXA1 _ 7
[48]   202 GGAATGTGATTTACCCAGATAAATCATCAGCTCAAGGGACTGCTTGGAGAG...TCCAGAGCTAGACGCAGAGGAAGGGGTCAAACATGTCCACATGGAACCTG FOXA1 _ 8
[49]   202 GGCTTTTCCAAGAACAATAGTGTTCTCCTAACAAACATTCGTTCCATCAGC...GGTGTTAGGGACAGAGTCCCAGGTGGCATAAAGTCGGGTTGGTCCTTGGC FOXA1 _ 9

예제 Fasta 파일의 경우 아래와 같은 motif를 찾았고, 총 74개의 sequence가 이 consensus에 부합했다. 즉, FOXA1라는 transcription factor의 binding motif (혹은 transcription factor binding site)의 consensus를 위의 Fasta 파일의 데이터를 통해 예측해볼 수 있다.

plot(gadem@motifList[[1]])

[[72]]
An object of class "align"
Slot "seq":
[1] "CTATTGACTTT"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "-"

Slot "seqID":
[1] 4

Slot "pos":
[1] 93

Slot "pval":
[1] 0.0001517677

Slot "fastaHeader":
[1] 4


[[73]]
An object of class "align"
Slot "seq":
[1] "ATATTTACTTG"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "+"

Slot "seqID":
[1] 42

Slot "pos":
[1] 135

Slot "pval":
[1] 0.0001719873

Slot "fastaHeader":
[1] 42


[[74]]
An object of class "align"
Slot "seq":
[1] "CTATTTGCTGG"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "+"

Slot "seqID":
[1] 41

Slot "pos":
[1] 53

Slot "pval":
[1] 0.0001754087

Slot "fastaHeader":
[1] 41

Gene 구성요소의 5가지 분류

1. Exon

- 폴리펩타이드를 코딩하는 DNA 서열이다.

2. Intron

- 단백질을 코딩하지 않는 DNA 서열이다. * 하지만 miRNA 등의 non-coding RNA를 코딩할 수 있는 것으로 보인다.

3. Transcription Start Site

- RNA 중합효소2(polymerase)가 붙어서 전사가 시작되는 위치이다. (RNA 중합효소2는 12개의 단백질의 중합체로 mRNA를 합성해내는 물질이다.)

4. Promoter

- Core Promoter : TATA 박스(타타 박스라고 읽음)라는 것이 있는데, 이 곳은 실제로 Transcription Factor와 기타 여러 Regulatory Protein이 실제로 붙는 곳이다. 이를 통해 50여개의 단백질 complex가 만들어지며, 이것이 실제로 trascription이 일어나는데 필요한 물질들이다.

- Upstream Promoter : transcription factor 와 기타 조절 단백질(regulatory protein)이 달라붙는 부분이다. Upstream promoter는 gene마다 갯수와 타입이 다르다. transcription factor는 붙어서 해당 유전자를 활성화하거나 비활성화한다.

Promoter는 gene expression을 "총괄" 한다고 보면 되며, 이를 Enhancer가 돕는다고 볼 수 있다. 참고로 잘 알려지지 않은 것 중에 Silencer 라는 것도 있다. 이는 Histon modification과 같은 방식으로 Promoter의 반대 기능 gene expression을 억제하는기능을 수행한다.

5. Enhancer

- 보통 upstream 또는 downstream 방향으로 수천 bp 떨어진 곳에 위치하며, (또는 gene 안에 존재할 수도 있다.) transcription factor가 달라붙게 되면 promoter 쪽으로 굽어서 유전자 발현을 조절한다. 예를 들어, promoter에게 "여기는 brain이니까 이 gene을 더 발현해라!" 라고 말하는 것으로 보면 된다.

유전체에서 유전자 발현 조절 과정 (Regulation)을 수행하는 요소를 찾는 것은 중요한데, 이러한 유전체상의 조절부위(Regulatory region)을 찾는 기본적인 원리는 이 부위에 transcription factor가 달라붙는다는 성질을 이용하여 antibody를 이용한 assay를 하는 것이다.

참고그림 - 출처 위키피디아

한국어 설명 자료

https://blog.naver.com/david6703/220973718228

Gene Regulation 관련 읽어보면 좋을 글

http://genetics.thetech.org/original_news/news14

★ https://www.nature.com/scitable/topicpage/gene-expression-14121669

https://www.khanacademy.org/science/biology/gene-regulation/gene-regulation-in-eukaryotes/a/overview-of-eukaryotic-gene-regulation

GWAS, QTL, Linkage study 셋의 차이점

Linkage Study : 가족데이터 - 멘델 유전 이용, 알려진 마커로 high-penetrance, single gene disease 멘델리안 질병에 대해서 연구

QTL : 가족데이터 - 멘델 유전 이용, polygenic & quantitative trait (ex. height)

GWAS : 인구집단 데이터 - LD를 이용, polygenic, categorical trait 둘다

출처

https://www.biostars.org/p/47725/

Variant Calling

정의

variant call 은 어떤 개인의 genome혹은 transcriptome에서 reference와 뉴클레오타이드 서열의 차이가 있는지를 결정하는 것이다. 이 과정은 variant frequency와 confidence measure가 고려된다.

유즈케이스

DNA-seq : variants

• 질병과의 유전적 연관성
• 암에서의 돌연변이 알아내기
• 이형접한 인구 빈도 알아내기

RNA-seq : allele-speicifc expression

• Allelic imbalance
• Association with isoform usage (splicing QTLs)
• RNA editing (allele absent from genome)

Chip-seq : allele-specific binding

• protein binding이 allele-specific하게 일어나는 현상

핵심은 variant call을 통해 어떤 locus에 'variant가 있다' 라는 정보를 알 수 있기 때문에 이를 통해 위와 같은 연구를 할 수 있음

Variant Call은 Genotype 보다 더 일반적인 접근이다.

genotype은 어떤 locus의 allele의 집합을 나타낸다. 근데 배수성을 보통 2n으로 가정하고 한다.  또한 genotype은 주로 SNP에 대해서만 고려된다. 하지만 RNA-seq 같은 경우 2n이 아니며 genotype이 없다. 또 암유전체학 등에서는 copy number variation이 중요하기 때문에 variantl call이 genotype보다 일반적인 접근법이라 할 수 있다.

참고

https://www.bioconductor.org/help/course-materials/2014/CSAMA2014/3_Wednesday/lectures/VariantCallingLecture.pdf

파일첨부

VariantCallingLecture.pdf

참고할 사이트들

Best Practices workflow

https://software.broadinstitute.org/gatk/best-practices/bp_3step.php?case=GermShortWGS