brainspan_20191017_152_chang.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Analyzing Brain Span Data

Brain Span Data가 만들어지기까지의 과정이 담긴 [PDF](file:///C:/Users/%EC%9E%A5%EC%9D%80%ED%98%95/Downloads/Transcriptome_Profiling%20(8).pdf)이다.
각 단계를요약하자면, 

- developmental stage별로 뇌 조직 샘플을 얻는다.
- Illumina Omni-2.5 million SNP arrays를 이용하여 genotyping 한다.
- gene expression을 알기 위해 total RNA를 추출 후 mRNA를 추출한다.
- mRNA를 이용해 cDNA를 합성하고 sequencing한다. 
- 그 후 reference sequence를 바탕으로 alignment를 진행한다. reference sequence로는 human genome sequence (hg19, GRCh37)를 이용한다.
- alignment output인 BAM을 SAM으로 바꾸고, 최종적으로 MRF(Mapped Read Format)으로 바꾼다.
- MRF를 이용해 RPKM 단위로 Expression을 측정한다.

### 1.Brain Span Data를 다운받고 로드하기

R studio cloud에서 terminal에 wget을 이용하면 다운 받을 수 있다. 혹은 코드 블럭에 {bash}를 써서 다운 받을 수 있다.
```{bash}
wget https://www.dropbox.com/s/u3ii9x4wpnzazm5/data_brainspan_DFC.20190928.Rdata
```

```{r}
load('data_brainspan_DFC.20190928.Rdata')
```


### 2. Data 살펴보기

1. e 데이터의 row에는 52376의 gene, column에는 참가자 35명이 있다. 35명의 gene에 대한 발현량을 나타낸 데이터이다. 발현량의 단위는 RPKM이다. RPKM은 Reads Per Kilobase of transcript per Million mapped reads의 약자이다. 
$\mbox{RPKM = Total fragments / Mapped Reads(millions)*Exon Length(Kb)}$으로 계산한다.
```{r}
head(e)
```


2. g 데이터는 52376개의 유전자에 대한 정보이다. gene id, e데이터와 매칭되는 row number, gene symbol, gene type에 대한 정보를 제공한다.
```{r}
head(g)
```


3. s 데이터에는 유전자를 제공한 donor에 대한 정보가 나와있다. 첫 번째 column에는 e 데이터의 column과 일치하도록 column number가 있다. 두번째 column부터는 donor의 id, name, age, sex(gender라고 나와있지만 sex일 것이다.)가 있고, RNA를 추출한 structure의 id, acronym, name이 있다. 여기에서 structure는 모두 dorsolateral prefrontal cortex이다.
```{r}
s
```


### 3. 분석 주제 정하기(question & aim)

> Expression changes over development by gene type

태아일 때부터 성인이 되기까지 뇌에서 각 gene type의 발현량의 경향성이 어떻게 되는지 알아볼 것이다.

이를 위해 

1. [PDF](file:///C:/Users/%EC%9E%A5%EC%9D%80%ED%98%95/Downloads/Transcriptome_Profiling%20(8).pdf)에 제시된 Developmental Period별로 age를 stage로 구분하여 donor들을 stage별로 묶는다.

2. 각 gene type을 s 데이터의 column에 추가하여 사람마다 각 gene type에 대해 평균 발현량이 얼마인지 계산 후 저장한다.

3. stage 별 평균을 구하고, x축을 stage로 하고 y축을 발현량으로 해서 각 gene type의 발현량을 line plot을 그린다.


### 4. 코드 작성하기

필요한 library를 불러온다.
```{r}
library(dslabs)
library(tidyverse)
library(ggthemes)
```

s 데이터에 stage라는 column을 만들어 각 age에 맞게 stage값이 들어가도록 한다.
```{r}
## stage group 만들기
s<-s %>%
  mutate(stage=case_when(
    age %in% c("8 pcw","9 pcw", "12 pcw") ~ 2,
    age %in% c("13 pcw","16 pcw", "17 pcw") ~ 3,
    age %in% c("19 pcw","21 pcw", "24 pcw") ~ 4,
    age %in% c("26 pcw","37 pcw") ~ 5,
    age %in% c("4 mos") ~ 6,
    age %in% c("10 mos","1 yrs") ~ 7,
    age %in% c("2 yrs","3 yrs", "4 yrs") ~ 8,
    age %in% c("8 yrs","11 yrs") ~ 9,
    age %in% c("13 yrs","18 yrs", "19 yrs") ~ 10,
    TRUE ~ 11))
```

s 데이터에 gene type이 많아서 4가지 큰 분류의 gene type만 분석하기로 했다. Donor마다 4가지 gene type의 평균 발현량을 저장하기 위해 column 4개를 추가한다.
```{r}
## gene type에 따른 발현량 넣을 공간 마련하기
s<-s %>%
  mutate(protein_coding=0 , activated_ig=0, inactivated_ig=0, noncoding_rna=0) 
```

4가지 큰 분류에 해당하는 gene type을 ind에 넣어준다. 이 분류는 [Gencode](https://www.gencodegenes.org/pages/biotypes.html)를 참고하였다.
```{r}
## ind에 g data의 해당 gene type만 추출한 index 넣기

ind_pc<-(g$gene_type == "protein_coding" & !is.na(g$gene_type))

ind_a_ig<-(g$gene_type %in% c("IG_C_gene","IG_D_gene","IG_J_gene","IG_LV_gene","IG_V_gene","TR_C_gene","TR_J_gene","TR_V_gene","TR_D_gene") & !is.na(g$gene_type))

ind_ia_ig<-(g$gene_type %in% c("IG_pseudogene","IG_C_pseudogene","IG_J_pseudogene","IG_V_pseudogene", "TR_V_pseudogene","TR_J_pseudogene") & !is.na(g$gene_type))

ind_nc_rna<-(g$gene_type %in% c("Mt_rRNA","Mt_tRNA","miRNA","misc_RNA","rRNA","scRNA","snRNA","snoRNA","ribozyme","sRNA","scaRNA") & !is.na(g$gene_type))
```

index는 true, false값으로 52376개의 모든 유전자에 대해 값을 가지므로, 각 gene type에 해당하는 gene이 몇 개인지 알기 위해  새로운 변수에 저장해준다. 

row_num값을 저장했지만, 이번 분석에서는 구체적인 row_num을 사용하지 않는다. 단순히 length를 구할 때 사용할 것이므로 굳이 row_num으로 하지 않아도 되지만 다른 분석을 할 때에는 필요할 수 있다.
```{r}
## length 사용해서 해당 gene type의 총 개수를 구하기 위함
pc<-g$row_num[ind_pc]
a_ig<-g$row_num[ind_a_ig]
ia_ig<-g$row_num[ind_ia_ig]
nc_rna<-g$row_num[ind_nc_rna]
```

35명의 donor의 4가지 gene type의 평균 발현량을 구하기 위한 함수와 적용이다.

e 데이터의 column들이 character여서 function을 사용할 때 애를 먹었던 부분이다. 함수를 작성하는 과정에서 column이 character여서 애를 먹을 때 적용할 수 있는 방법을 알게 되었고, paste함수와 column을 부를 때 `$` 말고도 `[[]]`가 있다는 것을 알게 되었다. 

원래의 column이름이 아니라 다른 변수에 저장된 character가 사용될 때는 `$` 가 아닌 `[[]]`을 사용해야 제대로 된 결과가 나왔다. `$`를 쓰면 error가 나지 않으면서 다른 결과값을 주어서 디버깅이 힘들 것 같다. 다음부터 코딩할 때 이점을 유의하면서 코딩을 해야겠다. 둘의 차이점에 대해서 더 구체적으로 알아봐야 겠다.
```{r}
## 35명의 해당 gene type의 평균 발현량 구하기
compute_avg_pc <- function(n){
  a<-as.character(n)
  b<-paste0("X",a)
  s$protein_coding[n]<-sum(e[[b]][ind_pc]) / length(pc)
}

compute_avg_a_ig <- function(n){
  a<-as.character(n)
  b<-paste0("X",a)
  s$activated_ig[n]<-sum(e[[b]][ind_a_ig]) / length(a_ig)
}

compute_avg_ia_ig <- function(n){
  a<-as.character(n)
  b<-paste0("X",a)
  s$inactivated_ig[n]<-sum(e[[b]][ind_ia_ig]) / length(ia_ig)
}

compute_avg_nc_rna <- function(n){
  a<-as.character(n)
  b<-paste0("X",a)
  s$noncoding_rna[n]<-sum(e[[b]][ind_nc_rna]) / length(nc_rna)
}


y<-1:35
s$protein_coding<-sapply(y,compute_avg_pc)
s$activated_ig<-sapply(y,compute_avg_a_ig)
s$inactivated_ig<-sapply(y,compute_avg_ia_ig)
s$noncoding_rna<-sapply(y,compute_avg_nc_rna)
```

드디어 plot을 그리는 단계이다. 

- stage 별로 묶고, 4가지 gene type에 대한 avg를 구한다. 
- plot을 그릴 때 x축을 stage로 하고, y축을 발현량의 avg값으로 하여 총 4개의 line을 그리고, color blind friendly 색을 지정해준다. 각 line에 대해 text가 옆에 오도록 위치시키고 line과 같은 색을 지정해준다.
- log scale을 사용하지 않으면 protein coding을 제외한 gene type의 경향을 파악하기 힘들어서 y축을 log scale로 만들어준다.(log scale을 사용했기 때문에 각 line의 증감률은 서로 비교할 수 있지만 각 line의 증감량은 한눈에 비교할 수 없다.)
- 제목와 축 제목을 설정해준다.
- stage 간의 grid(minor grid)는 의미가 없기 때문에 제거해준다.

```{r}
## stage별 평균 구해서 color blind friendly 색으로 line 그리기
#p<- (ggsave로 저장할 때 사용)
s %>%
  group_by(stage) %>%
  summarise(avg_pc = mean(protein_coding), avg_a_ig = mean(activated_ig), avg_ia_ig=mean(inactivated_ig), avg_nc_rna=mean(noncoding_rna)) %>%
  ggplot() +
  geom_point(aes(stage, avg_pc),col="#CC79A7") +
  geom_point(aes(stage, avg_a_ig),col="#009E73") +
  geom_point(aes(stage, avg_ia_ig),col="#56B4E9") +
  geom_point(aes(stage, avg_nc_rna),col="#E69F00") +
  geom_line(aes(stage, avg_pc),col="#CC79A7") +
  geom_line(aes(stage, avg_a_ig),col="#009E73") +
  geom_line(aes(stage, avg_ia_ig),col="#56B4E9") +
  geom_line(aes(stage, avg_nc_rna),col="#E69F00") +
  geom_text(aes(x=10,y=5,label="Protein coding"),colour="#CC79A7")+
  geom_text(aes(x=10, y=0.4,label="Non-coding RNA"), colour="#E69F00")+
  geom_text(aes(x=10, y=0.07,label="Activated Ig"), colour="#009E73")+
  geom_text(aes(x=10, y=0.003,label="Inactivated Ig"), colour="#56B4E9")+
  scale_y_log10() +
  scale_x_continuous(breaks = seq(0,12,1)) +
  ggtitle("Gene Expressions in Brain by Stage") +
  xlab("Developmental Stages") + 
  ylab("Average Expression 
(RPKM, Log Scale)")+
  theme(panel.grid.minor.x = element_blank())
```


plot만 저장하기 위해서는 아래의 코드를 실행하면 된다.
```{r}
#ggsave("plot.BrainSpan.line_GeneExpressionByStage.20191017.pdf",p, width=9, height=9)
```


### 5. 분석하기(findings)

1. Gene type간의 관계를 파악해보자.
- log scale을 적용하지 않으면 Protein coding만 10 RPKM 정도로 가장 발현량이 많고, 나머지는 0 RPKM 부근으로, 10 RPKM에 가까운 차이가 난다.
- Protein coding을 뺀 나머지 gene type은 모두 0이상 1이하라서 서로 간에 큰 차이가 나지 않지만 그래도 전 stage에 걸쳐 뚜렷하게 상하관계가 유지된다. 전 stage에 걸친 발현량 순서는 Protein coding > Non-coding RNA > Activated Ig > Inactivated Ig 이다. 
- 서로 다른 gene type간에 비슷한 경향성을 보이는 type이 없으므로 뚜렷한 상관관계는 없는 것으로 판단된다. 하지만 stage 8까지 non-coding RNA와 Inactivated Ig가 반비례하는 경향이 보인다.

2. 각 gene type의 stage에 따른 발현량의 양상을 살펴보자.
- Protein coding은 전 stage에 걸쳐 발현량이 10 RPKM 정도로 일정하게 유지된다. 
- Non-coding RNA는 0.1 RPKM 정도의 발현량을 유지하는데, stage 8인 early childhood 시기에 발현량이 증가했다가 다시 감소하는 것을 알 수 있다.
- Activated Ig를 제외한 gene type들은 일정한 증가나 감소 추세를 보이지 않고 일정하게 유지되면서 중간중간에 증감이 나타났다 사라지는 반면, Activated Ig는 상승폭 자체는 작지만 일정하게 상승하는 경향성을 보인다. stage 5와 6 사이에 0.1 RPKM을 넘어선다. 이 시기는 prenatal 단계에서 birth로 넘어가는 시기이다.
- Inactivated Ig는 발현량이 0.01 RPKM 내외로 나타난다. stage 8에 0.001로 급감하고 stage 9에서 0.05로 급증해서, stage 8과 9사이에 약 50배의 차이가 난다. stage 8은 19 months-5yrs의 early childhood 단계이고, stage 9은 6-11yrs의 late childhood 단계이다.