[R] 변이 위치 표시해보기!

아주 오랜만에 포스팅을 하네요  ㅎㅎ

 

샤이니 공부하면서 새로 배운 것들 정리하러 왔습니다.

 

---

 

GWAS를 분석할때나 혹은 변이를 분석하고 논문을 낼 때 참고용 그림으로서

 

변이의 위치가 유전자 어디에 있는지에 대한 그림을 넣는 경우가 있습니다.

 

예전 GWAS분석 후 논문 제출 때에도 비슷한 그림을 그려서 낸적이 있기 때문에 약간의 필요성이 있을 것이라 생각하고  포스팅합니다.

 

물론!! R package에서 더 이쁘게 그려주는 tool이 있을 수 있지만 좀더 커스텀이 가능한게 장점이죠 ㅎㅎ

 

Lolipop plot

 

R 기본 plot으로도 변이 그래프를 그릴 수가 있습니다.

 

먼저 유전자들 전체의 시작과 끝들을 알아야 합니다.

 

이러한 정보들은 UCSC table browser를 통해 얻을 수 있습니다.

 

예시로 다음과 같은 정보를 이용하여 만들어보도록 하겠습니다.

 

# 예시
st=16857 # 시작
ed=19759 # 끝
exon_st<-c(16857, 17232, 17605, 17914, 18267, 18912) # exon 시작
exon_ed<-c(17055, 17368, 17742, 18061, 18379, 19759) # exon 끝
align_id='uc001aai.1' # align ID
symbol='WASH7P' # symbol gene

 

유전자들 몸체를 그려보도록 하겠습니다. 몸체는 rect라는 기능을 이용합니다.

 rect : 사각형을 그려주는 기능을 합니다.

 

plot(NA,xlim = c(st,ed),ylim=c(0,2))
for(i in 1:length(exon_st)){
  rect(xleft = exon_st[i],xright = exon_ed[i],ytop = 1.5,ybottom = .5)
}

 

 

대략적인 몸체가 만들어졌고, 필요없는 axis들 xlab, ylab을 제거해보도록 하겠습니다.

 

plot(NA,xlim = c(st,ed),ylim=c(0,2),axes = F,xlab='',ylab='')
abline(h=1)
for(i in 1:length(exon_st)){
  rect(xleft = exon_st[i],xright = exon_ed[i],ytop = 1.5,ybottom = .5,col='grey')
}

axes = F, xlab,ylab =F 를 이용하여 제거를 해주고, abline을 이용하여 가운데 선을 추가하여 이어짐을 표시해줍니다. rect안에 col을 이용하여 상자 안에 그럴싸한 색을 추가로 넣어 줍니다.

 

고런 다음 원하는 변이를 추가하는데 임의로 두 개의 변이를 추가해줍니다.

 

v=c(19012,19512)

lines(x=c(rep(v[1],2)),y=c(1.5,1.7))
lines(x=c(rep(v[2],2)),y=c(1.5,1.7))
points(x=v,y=c(rep(1.7,2)),pch=16,col='red')

 

 

변이의 크기가 너무 작다 싶으면 points안에 cex를 이용하여 크기를 조절해주면 됩니다.

 

옆에 어떤 유전자인지 써주고 싶다면, par(mar=c(bottom,left,top,right))이용해서 양옆 조절을 해주고 axis를 이용하여 원하는 위치에 label로 글자를 추가해준다.!

 

axis 에서 at, labels를 이용하여 축들의 tick 글자를 변경가능하다!

 

R의 기본적인 기능들로 아주 그럴싸하게 그림을 그려보았습니다.

 

응용으로 TSS 위치 그려주고 상자색 변경해주고 할 수가 있으니 실습해보면서 만들어보세요!

 

---

두번째로! 많이 사용되는 ggplot2를 이용한 그림!

 

d=data.frame(x1=exon_st,x2=exon_ed)
ggplot()+geom_rect(data=d,aes(xmax=x1,xmin=x2,ymin=0.5,ymax=1.5))

 

배경색을 변경하고, 똑같이 xlab, ylab제거 해주겠습니다.

 

theme안에서 axis.text, axis.ticks, acis.line을 설정해주도록 하겠습니다.

 

d=data.frame(x1=exon_st,x2=exon_ed)
ggplot()+geom_rect(data=d,aes(xmax=x1,xmin=x2,ymin=0.5,ymax=1.5))+
  theme_classic()+
  theme(axis.text = element_blank(),axis.ticks = element_blank(),axis.line = element_blank())

 

추가로 가운데 선!과 변이를 추가적으로 넣어보겠습니다.

 

 

d=data.frame(x1=exon_st,x2=exon_ed)
v=c(19012,19512)
v=merge(v,c(1))
ggplot()+geom_rect(data=d,aes(xmax=x1,xmin=x2,ymin=0.5,ymax=1.5))+
  theme_classic()+
  theme(axis.text = element_blank(),axis.ticks = element_blank(),axis.line = element_blank())+
  geom_segment(aes(x=st,xend=ed,y=1,yend=1))+
  geom_point(aes(x=v$x,y=1.7),col='red')+
  geom_segment(aes(x=v$x,xend=v$x,y=v$y,yend=v$y+.7))+xlab("")+ylab("")

 

뭔가 exon들 위에 있는 선이 거슬리고, 변이 위치 표시선이 거슬릴 수 있습니다.

 

해결방법은! ggplot을 만들때 순서!를 변경하면 됩니다! (변경하는 김에 변이 표시 크기도 변경)

 

원하는 변이가 intron일수도 있기 때문에 intron표시로서 점선을 길게 추가해보도록하겠습니다! (추가로 상자 높이도 조절)

ggplot()+
  theme_classic()+
  theme(axis.text = element_blank(),axis.ticks = element_blank(),axis.line = element_blank())+
  geom_hline(yintercept = 1,linetype='dashed')+
  geom_segment(aes(x=st,xend=ed,y=1,yend=1),size=.8)+
  geom_segment(aes(x=v$x,xend=v$x,y=v$y,yend=v$y+.7))+xlab("")+ylab("")+
  geom_point(aes(x=v$x,y=1.7),col='red',size=2)+
  geom_rect(data=d,aes(xmax=x1,xmin=x2,ymin=0.5,ymax=1.5))+
  ylim(-1,2.5)

 

마지막으로 어떤 유전자인지 표시해보도록하겠습니다.

 

ggplot()+
  theme_classic()+
  theme(axis.text.x = element_blank(),axis.ticks = element_blank(),axis.line = element_blank())+
  geom_hline(yintercept = 1,linetype='dashed')+
  geom_segment(aes(x=st,xend=ed,y=1,yend=1),size=.8)+
  geom_segment(aes(x=v$x,xend=v$x,y=v$y,yend=v$y+.7))+xlab("")+ylab("")+
  geom_point(aes(x=v$x,y=1.7),col='red',size=3)+
  geom_rect(data=d,aes(xmax=x1,xmin=x2,ymin=0.5,ymax=1.5))+
  scale_y_continuous(labels = sprintf('%s (%s)',align_id,symbol),breaks = c(1),limits = c(-1,2.5))+
  theme(axis.text.y = element_text(size=15))

 

ylim을 제거를 하고, scale_y_continuous에서 limits를 이용해서 범위조절을 하였습니다.

짠!

 

이렇게 해서 basic 기능과 ggplot2를 이용하여 변이를 표시하는 그림을 그려보았습니다.

 

물론!! 처음에서 이야기했습니다! 더 괜찮은 package가 어딘가에 (github에서도 찾으면 나와요!) 있기 때문에 참고용으로만 사용하세요!

 

 

 오늘 이러한 포스팅을 하는 이유는 shiny이용해서 비슷한 변이 표현해주는 그림을 다운받을 수 있는 사이트를 제작중에 있습니다! (어휴.. 다른건 안끝내고 또 일벌리네.)

 

졸업전에 sql이용도 해볼겸. 사이트도 만들어보고 싶었기 때문에 한번 제작해보도록 하겠습니다. (hosting까지 해보려고 해요!)

 

아무튼! 나중에 만들어지면 놀러와주세요~ ㅎㅎ

 

 

---

유용하셨거나, 잘 보셧다면 주변 광고 한번씩만 클릭 부탁드립니다! 감사합니다!