2017년 12월 18일 월요일

독서 기록 - 통계의 힘(니시우치 히로무 지음/신현호 옮김)

2014년 12월에 구입해서 거의 다 읽었다가 마지막에 집중력이 흐려져서 독서를 마무리하지 못했던 책을 다시 붙들고 처음부터 꼼꼼하게 읽었다. 책의 정식 제목은 '통계의 힘'이지만 부제를 붙이면 '빅테이터를 지배하는 통계의 힘'이다. 요즘은 AI(인공지능)·4차 산업혁명에 약간 가리워진 느낌이 없지는 않으나 방대한 데이터로부터 의미를 추출하는 과정에서 통계의 중요성은 말할 나위가 없다.


이 책은 통계학 입문서나 교재가 절대로 아니다. 일반 독자에게 통계한 이런 것임을 보여주기 위한 교양 서적이다. 그럼에도 불구하고 회귀분석, 일반화 선형모델(generalized linear model), 빈도론자와 베이즈론자의 대립 등 그동안 피상적으로만 이해했던 중요한 개념들을 알기 쉽게 설명하였다. 특히 현대 통계학의 기반을 마련한 로널드 피셔와 그의 유명한 밀크티 일화가 소개되어 흥미를 돋우고 있다. 밀크티를 만들 때 홍차에 우유를 넣었는지, 혹은 그 반대로 우유에 홍차를 넣었는지를 정확하게 알아맞힌다는 어떤 부인의 주장을 검증하기 위하여 피셔는 세계 최초의 임의화 비교실험을 수립하였다. 그 과정은 피셔가 1935년에 쓴 책 <실험계획법(The Design of Experiments)>에 잘 소개되었다고 한다. 피셔는 비료가 작물의 생산성이 효과를 주는지를 알아보기 위하여 경작지를 어떻게 구획하고 각각 어떻게 처리를 해야 하는지를 상세하게 연구하였다고 한다. 말하자면 현대 통계학은 생물통계학의 역사와 그 뿌리가 같은 셈이다. 통계학이 먼저 확립되고 난 다음에 생물학 분야에 적용된 것이 아니라는 뜻이다.

55쪽의 도표 5에서는 빅데이터 관련 전문용어를 설명하고 있는데, 맨 마지막 칸에 보인 R언어에 대한 설명이 재미있다. '오픈소스의 통계해석용 언어'. 유료 소프트웨어를 살 수 없는 비교적 가난한 학자들이 사용하는 언어인데, 최근 갑자기 주목을 받고 있다. 엑사데이터나 그린플럼, 나아가 SPSS로부터도 직접 R 라이이브러리를 호출할 수 있게 되었다.'

이 책의 결론은 다음의 도표 한 장으로 나타낼 수 있다. 설명변수가 몇개인지 여부에 상관없이 일반화 선형모델을 사용하여 (거의) 모든 통계적 분석을 감당할 수 있다는 것이다.


문맹(文盲)보다 더 무서운 것이 통계맹이라는 말도 있다. 요즘은 새로운 주장을 하기 위한 근거로서 여러 수치 자료를 제시하는 경우가 많은데, 통계학적 관점에서 볼 때 오류가 있거나 심지어는 의도적으로 자료의 의미를 왜곡하는 일이 벌어지고는 한다. 오죽하면 '세상에는 두 가지 거짓말이 있으니 하나는 새빨간 거짓말이요, 나머지는 통계'라는 말까지 있지 않은가. 통계는 정확히 쓰여야 하고 어떤 주장을 일방적으로 지지하기 위한 목적으로 악용되어서도 안된다. 또한 단순한 집계와 통계를 구별할 줄 알아야만 한다. 수학적 리터러시(literacy - 읽고 쓰는 능력 또는 특정 분야의 지식), 특히 통계 리터러시가 있어야 이러한 근거를 이용한 잘못된 주장을 제대로 파악하고 대처할 수 있는 것이다.

데이터가 말하려는 것은 무엇인가? 이것을 정확하게 파악하는 것이 통계학의 중요한 목표 중 하나라고 본다. 데이터가 말하는 바를 귀담아 듣고 판단을 내리는 것은 결코 쉬운 노릇은 아니다. 그래서인지 지금은 이 과정 자체를 AI가 하도록 내버려 두는 것 같다. 가까운 예로 각종 인터넷 구매 사이트에서 내 구매 성향을 파악하여 쇼핑할 거리를 제안하는 서비스가 그렇다. 이렇게 판단 자체를 인공지능에 맡기고 남는 시간에는 더 가치있는 일을 하겠다는 멋진 전제와 함께 말이다. 그런데 과연 미래는 커즈와일이 예측한 대로 그렇게 장미빛으로만 돌아갈까? 최근에 읽은 소설 <밈:언어가 사라진 세상>(독후감은 곧 블로그에 쓸 예정)과 니콜라스 G. 카의 주장(IT doesn't matter (2003) - 하바드 비즈니스 리뷰, 2016년 블로터의 논평; 2014년 <유리 감옥, 자동화와 우리> 출간 직후의 인터뷰)을 보면 걱정이 밀려온다.
이 책(유리 감옥)은 소프트웨어 개발자를 포함 많은 직종이 '디스킬(Deskill, 일을 처리하는데 필요한 숙련도 수준이 낮아짐)' 될 수 있다고 경고한다. 업무 프로세스의 상당 부분이 기계로 넘어간다는 이유에서다.
그렇다고 러다이트(Luddites, 신기술 반대자)를 옹호하는 책은 아니다. 이보다는 우리를 대체할 시스템이 초래할 결과와 파장을 훌륭히 탐구하고 있는 책이다. 연방 항공청(FAA)의 항공 자동화에 대한 경고, 전자 의료 기록이 초래할 수 있는 비용 상승과 건강 악영향 등에서 카가 우려하는 부분을 확인할 수 있다.

통계학 자체는 가치 중립적인 도구이다. 이 점은 명심해야 한다. 통계적 결론을 악용하는 것, 그리고 최종 판단을 내리는 주체 역할을 하지 않으려는 것, 이 두가지를 경계해야 할 것이다.

Statistics for Biologists

과학잡지 Nature에서는 생물학자를 위한 통계학 정보 콜렉션을 운영한다.


Statistics in biology, practical guides, points of significants 및 other resources로 구성되어 있으니 통계학 입문서 단행본을 읽다가 지루해지면 종종 방문해 보도록 하자.

2017년 12월 15일 금요일

DataCamp의 온라인 데이터 과학 강좌

해들리 위컴의 강의를 온라인으로 들으라고? R 관련 강좌를 웹에서 둘러보다가 DataCamp라는 곳을 발견했다. R, python, ggplot2, machine learning 등 총 98개의 강좌가 서비스되는 중이다. 한 달에 29 달러를 내면 모든 강좌를 무제한으로 수강할 수 있으며, 연 단위로 결제하면 한 달에 25달러로 할인이 적용된다. 영어 수강에 자신이 있다면 충분히 들을 만한 강좌가 아닐까 한다.


늘어나는 데이터 과학자의 수요에 맞춘 유용한 서비스라고 생각된다.

대표적인 MOOC 서비스로 잘 알려진 Coursera에서도 데이터 과학 분야 아래에 data analysis, machine learning 및 probability & statistics 소분야로 나누어서 261개 강좌를 서비스한다. 존스홉킨스대에서는 Coursera 내에 총 10개 강좌로 구성된 Data Science Specialization를 제공한다.

MOOC에 대해서는 2015년 내 블로그에 작성한 글이 있다.

UST 교수법 워크숍을 다녀와서 - MOOC와 Flipped learning

한국형 MOOC를 표방하는 K-MOOC에서도 비슷한 성격의 강좌를 다수 서비스한다. 2015년 10개 대학이 참여하여 27개 강좌로 시작하여 2017년 12월 현재 개설된 강좌는 총 570개로 꽤 큰 규모로 성장하였다. Coursera는 2,000개가 넘는 강좌를 제공하고 있다. K-MOOC의 강좌 분류는 Coursera와 약간 달라서 '데이터 과학'과 같이 트렌디한 분류명은 없다.


데이터 과학에 해당하는 K-MOOC의 강좌는 어떤 대분야에 속해 있을까? 몇 가지를 조사해 보았다. 생각보다 다양한 대분야로 분산된 모습을 볼 수 있다.

  • [A.I. Series] R을 이용한 통계학개론 -> 자연
  • Machine learning for data science -> 공학
  • 인공지능 및 기계학습 -> 공학
  • 빅데이터의 세계, 원리와 응용 -> 사회과학
  • 빅데이터 첫걸음 -> 공학
  • 통계학의 이해 -> 자연
  • 데이터 과학을 위한 Python 입문 -> 공학
  • 바이오헬스 빅데이터 마이닝 -> 인문과학
  • 텍스트 마이닝 실전 및 분석 -> 공학
이제는 학교나 학원에 다니지 않고서도 이렇게 학습을 할 수 있는 기회가 널리 제공되고 있으니 기초 통계학이나 스크립트 언어 사용법을 배울 기회가 없었다고 변명을 하기는 점점 곤란해질 것이다. 앞으로 KOBIC의 차세대생명정보강좌에서도 이러한 점을 감안해야 할 것이다. 그리고 기존에 IT나 통계 실무 분야에서 교육 서비스를 제공하던 회사와의 상생 방안도 마련되어야 할 것으로 보인다.

2017년 12월 14일 목요일

[하루에 한 R] 매트릭스에서 특정 값을 갖는 row 추출하기

[작성 후 추가한 글] 제목이 잘못되었다. 엄밀히 말해서 매트릭스가 아니라 데이터프레임이다.

해들리 위컴(Hardley Wickham)은 뉴질랜드의 통계학자 및 데이터 사이언티스트이자 저명한 R 개발자이다. RStudio의 Chief Scientist이기도 한 그는 ggplot2(), dplyr 등 데이터 처리와 시각화에 탁월한 R 패키지를 다수 개발하였고 활발한 저술과 강연활동도 펼치고 있다. 그가 개발한 패키지들은 GitHub 사이트에서 볼 수 있다. Quora에는 이런 질문도 있었다. '어떻게 해들리 위컴은 패키지 개발을 통해 그렇게 R에 많이 기여할 수 있었나요?'

이 질문에 대해서 위컴 자신이 답을 달았다. 다른 사람이 쓴 답이 무슨 소용이 있으랴.

  • 나는 많이 씁니다. 침대에서 일어나자마자 매일 아침 60-90분 정도를 씁니다(코드? 아티클?)
  • 나는 많이 읽습니다. 300개 정도의 블로그를 팔로우하고, 트위터와 스택 오버플로우에서 R이라는 태그가 달린 글을 계속 봅니다. 물론 모든 글을 깊이있게 읽지는 못하고 훑어보는 수준이지만 나에게 큰 도움이 됩니다.
  • Chunking(덩이 짓기?). '컨텍스트 스위칭'은 비용이 많이 듭니다(여기에서 context switching이란 소프트웨어 개발에 관한 용어가 아니라 어떤 일을 하다가 다른 일로 전환하는 것을 의미하는 것으로 보인다). 동시에 여러 패키지를 작업하면 아무것도 되지 않습니다. 새로운 기능에 대한 문제점과 아이디어를 꾸준히 축적하고 있다가 임계 질량에 다다르면 이 패키지에 며칠을 보내는 식으로 일을 합니다.
한마디로 말해서 참 멋진 사람이다. dplyr 패키지에 대한 국문 소개는 양우성님의 웹사이트를 참조하면 좋다. 내가 여기까지 흘러들어오게 된 것은 인하대 유동현 교수님 덕분이다.

오늘은 데이터프레임에서 특정 조건에 맞는 행을 뽑아오는 방법에 대해 공부해 보려고 한다. 숫자를 가지고 조건에 맞는 행을 뽑는 것이 아니라, 캐릭터 값을 가지고 일을 하려는 것이다. 몇 천 유전자에 대한 blast 결과를 테이블로 출력하여 R로 읽어들여 데이터프레임으로 만든 다음 내가 관심을 갖고 있는 몇십, 혹은 몇백개의 유전자에 해당하는 hit row를 쏙 빼내려는 것이다. Genes of interest는 별도의 벡터에 저장된 상태라고 가정하자. Perl에서는 hash를 사용하여 if (exists $seen{$item}) {...}으로 늘 하던 일이다. 여담이지만 내가 Perl애소 이 기법을 처음 접한 것은 phred/phrap/consed 패키지의 phredPhrap 스크립트를 뜯어보면서였다.

dplyr 패키지를 쓰면 오늘 논하는 것을 포함하여 더욱 정교한 작업을 할 수 있을 것이다. 그러나 이 글을 작성하는 목적은 R base package의 기본 기능을 최대한 활용해 보자는 것이다. 핵심이 되는 것은 which() 함수와 %in% 연산자이다. 이 특수 연산자에 대한 도움말을 띄워 보자. 함수에 대한 도움말을 찾을 때만 help() 또는 ? 명령어를 쓰라는 법은 없다.
> ?"%in%"   # 또는 help("%in%")
a %in% b라고 입력하면 a 벡터의 모든 원소에 대해서 이 값이 b 벡터에도 존재하는지를 TRUE or FALSE로 출력하는 것이다. 실제 사례를 보자.
> a = c("A","B","C","D","E")
> b = c("C","D","E","F","G")
> a %in% b
[1] FALSE FALSE  TRUE  TRUE  TRUE
which() 함수를 이용하면 a 벡터의 원소 중 b에도 있는 것, 그리고 a에만 있는 것을 다음과 같이 출력할 수 있다.
> a[which(a %in% b)]
[1] "C" "D" "E"
> a[-which(a %in% b)]
[1] "A" "B"
벡터에 대하여 방법을 알아보았으니 어제 글에서 활용한 샘플, 즉 blast 결과 파일을 파싱한 데이터프레임에 대해서 이 일을 해 보자. 'd'라고 명명한 데이터프레임의 구조는 다음과 같다. head() 함수를 이용하여 앞부분만 발췌하였다.


추출할 유전자는 벡터 e(for extract)에 넣어두었다. 몇 개 되지 않으면 손으로 입력해도 되지만, 그 수가 많다면 파일로부터 읽어들이는 것이 현명할 것이다.
> e
[1] "BfmR" "CsuB" "PgaB" "EntA" "PlcD"
d$V1 컬럼의 값 중 e 벡터에 존재하는 것을 찾아서 그 row를 추출하면 된다. 코드는 다음의 한 줄이다.
> d[which(d$V1 %in% e), ]
실제 실행 화면을 보자.


너무나 간단해서 내가 다 미안할 지경이다. 다음에는 데이터프레임의 row 혹은 column을 그룹 단위로 조작하는 방법에 대해서 공부해 보련다.

2017년 12월 13일 수요일

독서 기록 - 1코노미 외 4 권


1코노미(부제: 1인 가구가 만드는 비즈니스 트렌드)

이준영 지음. 1인 가구가 보편화하면서 나타날 새로운 비즈니스 기회에 대하여 쓴 책이다. 나는 이에 대해서는 다소 부정적이다. 1인 가구 구매력의 국가적 총합은 과연 과거 다인 가구의 그것을 월등히 초과할 수 있을까? 주택 구매나 육아 및 교육에 돈을 쓰지 않으니 '나'를 위해서 더 많은 지출을 할 수 있을까? 인구가 줄면 경제 전반에 활력이 떨어지고 소득 자체가 늘어나기 힘들다. 더군다나 사회 구성원이 재생산되지 않음으로 인해 획기적인 경제 체질 개선이 일어나지 않으면 더욱 밝은 전망을 하기 어렵다. 1인 가구가 늘어나는 것은 자발적인 선택도 있지만 가정을 꾸릴 경제적 여건이 되지 못하여 불가피하게 이를 선택하는 면도 있다. 이러한 현상을 골고로 고려하지 않고 지나치게 낙관적으로 쓴 책이라는 생각이 들었다. 그리고 지나치게 많은 신조어를 소개함으로써 책을 읽는 것이 아니라 인터넷 기사를 넘기는 느낌이었다.

고독력, 관태기(관계권태기), 자뻑, 포미족, 탕진잼, 0.5가구, 혼행, DD(Do not Disturb족), 혼놀족...

커피가 죄가 되지 않는 101가지 이유

로잔느 산토스, 다르시 리마 지음. 김정운 옮김. 커피(카페인)만큼 건강에 미치는 긍정적 효과와 부정적 효과에 대한 상반되는 연구 결과가 꾸준히 발표되는 기호식품도 없을 것이다. 커피 수십 잔 분량에 해당하는 카페인을 갑자기 섭취하면 당연히 좋을 리가 없다. 카페인은 인체에 가벼운 흥분을 일으키는 물질임은 분명하지만, 커피에는 폴리페놀의 일종인 클로로겐산(chlorogenic acid) 등 유익한 성분은 더욱 많다. 커피 원산지의 농부들이 과연 합당한 노동의 댓가를 받고 있는지는 좀 더 고민해 봐야 되겠지만, 커피를 마시는 것이 그들에게 도움이 됨은 분명하다. 즐겨 마시자! 단, 탄산음료나 에너지 드링크 등에 합성 카페인을 지나치게 많이 넣는 것은 법으로 규제를 하면 좋겠다.

어떤 경제를 만들 것인가(지금의 시대정신은 '행복한 경제 만들기'다)

김동열 지음. 34개 OECD 회원국 중 우리나라 국민의 행복감은 꼴찌 수준이다. 경제 규모는 세계 11위 수준이라고 하지만 고용 불안, 노후 불안, 소득 불평등이 우리를 짓누르고 있다. 한국 경제가 추구해야할 새로운 모델을 알아본 책이다. 이에 대한 방안으로서 통일 대비, 한일 FTA, 연금의 확충, 주택연금,세금 마일리지 등을 예로 들었다.

정서적 협박에서 벗어나라(부제: 내 마음을 옭아매는 영혼의 감옥)

저우무쯔 지음, 하은지 옮김. '네가 어떻게 나한테 이럴 수 있어!' 가장 가까운 사이에서 오히려 흔한 정서적 협박은 개인을 파멸에 이르게 한다. 일단 그 상황을 벗어나고, 거절할 줄 알아야 한다. 호의를 마땅한 권리라고 여기는 사람들에게 정서적 경계선을 세우고 단호하게 대해야 한다.

리얼리스트를 위한 유토피아 플랜(부제: 우리가 바라는 세상을 현실에서 만드는 법)

뤼트허르 브레흐만 지음, 안기순 옮김. 단연코 이번 독서의 하이라이트, 먼저 지난 10월 한겨레 신문 인터뷰 기사를 보자("주 15시간 일한다는 게 왜 비현실적인 거죠?"). 1988년 생으로 이제 서른살도 안된 젊은이가 그동안 우리가 굳건하게 믿어왔던 자본주의의 토대 - 가난은 나태·노력 부족에서 기인한다 - 를 허무는 급진적인 주장을 한다면 우리 사회에서는 이를 과연 받아들일까? 인류가 경험한 많은 실험 사례에서 기본 소득은 충분히 긍정적 효과가 있음이 입증되었다. '돈을 이리저리 옮기면서' 돈을 버는 금융 부문은 과감히 개혁하고, 노동 시간을 줄이고, 가치 있는 직업을 갖도록 노력하고, 기본 소득 지급을 통해서 현실 세계에서의 유토피아에 다가갈 수 있을 것으로 저자는 내다보고 있다. 다시 한번 더 읽고 싶은 책이다.

[하루에 한 R] BLAST tabular output에서 best hit 뽑아내기

BLAST+(최신 버전)가 완전히 정착된 것이 언제인데 나는 아직도 'legacy' BLAST(v2.2.26)를 즐겨 사용한다. blastall 커맨드 라인 사용법에 너무나 길들여져 있기 때문이다. BLAST+의 옵션 설명서를 인쇄하여 책상 앞에다 붙여놓아야 겨우 익숙해질까?

예전, 하지만 그렇게 오래 전도 아니던 시절, 미생물 유전체의 functional annotation을 하려면 Swiss-Prot이나 NCBI NR 데이터베이스를 로컬 컴퓨터에 받은 다음 수천 개의 query를 blastp로 날리는 일을 수시로 했었다. 요즘은 RAST server에 올려 버리거나, 로컬 머신을 쓴다 하여도 Prokka를 쓰면 단 몇 분에 해결이 되이 큰 용량의 공개 단백질 서열 DB에 대해서 BLAST를 직접 돌릴 일이 많지 않다. 또 그럴 일이 있으면 유료 서비스에 가입한 Blast2GO에서 검색을 수행하면 된다.

그래도 간혹 연구 목적을 위해 특별한 목적으로 만든 서열 데이터베이스에 대해서 BLAST를 실행할 일이 있다. 그 결과물을 파싱하는 프로그램을 직접 구현하는 것이 예전에는 매우 보편적인 생명정보학 실습 과제이기도 했다. 누가 그랬더라? 이제 제발 그런 거 하지 말라고...

사람의 눈으로만 보기 편한 BLAST 결과물을 컴퓨터에게 읽히려고 하니 프로그래밍으로 처리하기에 얼마나 불편하겠는가? 사실 내가 tabular output format(blastall -m 8 또는 -m 9; BLAST+에서는 -outfmt 6)에 관심을 가진 것도 얼마 되지 않았다. 이를 엑셀로 불러들이면 완벽하지는 않더라도 대량의 결과를 처리하기에 좋으니까 말이다. 엑셀의 귀재라면 각 query에 대한 hit 중에서 best를 골라내는 방법을 함수로 구현할 수 있을지도 모른다.

R을 사용하면 BLAST의 tabular format output으로부터 각 query에 대한 best hit을 쉽게 뽑을 수 있을 것이라는 생각을 오랫동안 하고 있었다. 사실 R 환경 자체에서 BLAST를 실행할 수 있는 rBLAST라는 패키지가 있으니 이 환경 안에서 best hit을 찾을 수 있을 것이고 또 더 검색을 해 보면 다른 도구가 있을 것이다. 범위를 좀 더 넓힌다면 reciprocal best hit을 찾아주는 유틸리티도 있을 것이다. 하지만 나는 좀 더 '하등한' 수준의 해결 방법을 알고 싶었다.

[SEQanswers] How to output only best-hit result using standalone blast <= 심지어 sort 명령어를 사용한 방법도 있다. one-liner를 즐기는 사람에게는 좋은 소일거리가 될 것이다.

그러면 본론으로 들어가서, blastall -m 9(tabular with comment lines; #으로 시작) 명령으로 작성한 결과물을 R에서 읽어들여서 best hit을 찾는 방법을 공부해 보자. 예전에 마이크로어레이 자료 분석을 공부하면서 control spot 때문에 한 유전자에 대한 형광 측정값이 여럿 존재하는 경우 이것의 median 값을 취하도록 R 코드를 짰던 것 같은데 도저히 기억이 나질 않는다. 아마 apply() 혹은 by() 계열의 함수를 쓰지 않았었나 생각한다.

그런데 구글링을 계속 해 보니 base R만을 이용한 아주 원초적인 방법을 제시한 글이 있어서 소개하고자 한다. 이 안에는 외부 패키지와 함수를 이용한 복잡한 방법도 포함되어 있다. 공부를 위해서 나머지 방법도 따라서 해 보아야 한다.

[Stack Overflow] Extract the maximum value within each group in a dataframe

blastall -m 9 옵션을 사용하여 출력 파일(tblastn.m9,.txt)을 만든 뒤 각 query에 대한 best hit을 뽑아보자.
> d = read.table("tblastn.m9.txt",sep="\t",comment.char="#",header=F)
> d
             V1   V2     V3  V4 V5 V6  V7  V8   V9  V10      V11    V12
1          BfmR bfmR 100.00 238  0  0   1 238    1  714 8.0e-156  423.0
2          BfmR bfmS  30.56  36 25  0  21  56 1437 1330  2.3e-01   22.3
3          BfmS bfmS  99.82 549  1  0   1 549    1 1647  0.0e+00 1083.0
4          BfmS bfmS  27.78  54 39  1 133 186  658  533  7.1e+00   18.9
5          BfmS  nfu  38.46  13  8  0 117 129  442  480  8.1e+00   18.5
6          CsuA csuA 100.00 192  0  0   1 192    1  576 2.0e-129  353.0
7          CsuA csuE  43.75  16  9  0 177 192  970 1017  4.5e-01   20.8
8          CsuB csuB 100.00 172  0  0   1 172    1  516 4.0e-114  312.0
...중간 생략
> do.call(rbind, lapply(split(d,d$V1), function(x) {return(x[which.max(x$V12),])}))
                       V1   V2     V3  V4 V5 V6 V7  V8 V9  V10    V11  V12
BfmR                 BfmR bfmR 100.00 238  0  0  1 238  1  714 8e-156  423
BfmS                 BfmS bfmS  99.82 549  1  0  1 549  1 1647  0e+00 1083
CsuA                 CsuA csuA 100.00 192  0  0  1 192  1  576 2e-129  353
CsuB                 CsuB csuB 100.00 172  0  0  1 172  1  516 4e-114  312
CsuC                 CsuC csuC 100.00 277  0  0  1 277  1  831  0e+00  568
CsuD                 CsuD csuD 100.00 832  0  0  1 832  1 2496  0e+00 1650
CsuE                 CsuE csuE 100.00 339  0  0  1 339  1 1017  0e+00  666
EntA                 EntA entA 100.00 256  0  0  1 256  1  768  0e+00  502
EpsA                 EpsA epsA  99.73 366  1  0  1 366  1 1098  0e+00  727
NfuA                 NfuA  nfu 100.00 212  0  0  1 212  1  636 2e-162  438
OmpA_partial OmpA_partial ompA 100.00 315  0  0  1 315  1  945  0e+00  592
pgaA                 pgaA pgaA 100.00 747  0  0  1 747  1 2241  0e+00 1544
PgaB                 PgaB pgaB  99.80 510  1  0  1 510  1 1530  0e+00  980
PgaC                 PgaC pgaC 100.00 392  0  0  1 392  1 1176  0e+00  805
PgaD                 PgaD pgaD 100.00 154  0  0  1 154  1  462 6e-101  278
PlcD                 PlcD plcD 100.00 487  0  0  1 487  1 1461  0e+00 1018
Ptk                   Ptk  ptk 100.00 727  0  0  1 727  1 2181  0e+00 1372
V1(query)의 값이 동일한 hit group에 대해서 V12(bit score)가 최대인 것을 출력하는 것이다. 어떻게 이보다 더 간단명료할 수 있겠는가?

[하루에 한 R] order() 함수의 이해

어떤 벡터가 있을 때 sort() 함수를 적용하면 그 원소들을 다음과 같이 오름차순으로 정리하여 출력한다. 내림차순으로 정동하고 싶으면 decreasing=T를 설정하면 된다.
> x = c(1,5,3,2,4)
> sort(x)
[1] 1 2 3 4 5
order() 함수는 벡터의 원소를 이 인덱스 번호순으로 놓으면 값이 정렬된다는 것을 의미한다. Character 타입의 벡터를 사용하여 설명하는 것이 이해하기 쉬울 것이다. 이를 이해한다면 아래에서 보인 마지막 명령이 무엇을 뜻하는지도 납득할 수 있다.
> y = c("a","c","e","b","c")
> sort(y)
[1] "a" "b" "c" "c" "e"
> order(y)
[1] 1 4 2 5 3ㅁㄴㅇ
> y[order(y)]
[1] "a" "b" "c" "c" "e"
order(y)를 사용하여 알 수 있는 것은 벡터 y의 원소를 값에 따라서 순서대로 나열하려면 첫번째, 네번째, 두번째, 다섯번째, 그리고 세번째 원소의 순서대로 놓으면 된다는 것이다.

order() 함수를 이용하면 행이나 열을 name에 따라서 정렬하는 것도 가능하다. 연습용 4x3 매트릭스를 만들고 colnames()로 각 열에 이름을 붙인 뒤 이것에 따라서 컬럼을 정렬하는 연습을 해 본다. 컬럼의 위치만 바뀔 뿐 각 컬럼 내의 값은 원래의 행 번호를 그대로 지킨다.
> z =  matrix(c(4:1,12:9,5:8),4,3)
> z
     [,1] [,2] [,3]
[1,]    4   12    5
[2,]    3   11    6
[3,]    2   10    7
[4,]    1    9    8
> colnames(z) = c("A","C","B")
> z
     A  C B
[1,] 4 12 5
[2,] 3 11 6
[3,] 2 10 7
[4,] 1  9 8
> z[, order(colnames(z))]
     A B  C
[1,] 4 5 12
[2,] 3 6 11
[3,] 2 7 10
[4,] 1 8  9
apply()와 join() 계열의 함수를 잘 써야 하는데... 아직도 멀었다!

2017년 12월 12일 화요일

Pseudomolecule을 만들 때 사용할 space sequence

Pseudomolecule 혹은 pseudo-molecule이란 contig sequence의 사이에 적당한 spacer sequence를 넣어서 이어 붙인(concatenation) 염기서열을 말한다. 말하자면 가상의 염색체에 해당하는 염기서열인 셈이다. Contig들의 연결 순서는 reference sequence와 비교했거나 혹은 pair end sequencing read의 alignment를 이용하여 사실에 가깝게 결정해 놓은 다음 연결하는 것이 가장 바람직할 것이다. Spacer의 길이는 실제 sequence gap을 반영하도록 설계하면 더 좋을 것이다. 그러나 항상 이러한 행운이 따르는 것은 아니라서, 연결 정보가 마련되지 않은 contig 서열을 그냥 이어 붙이는 경우도 종종 있다. Contig의 순서는 어떻게 할 것인가? 특별히 따라야 할 규칙은 없지만 가장 긴 서열을 앞에 오도록 정렬한 다음 붙이기도 한다.

가장 바람직한 spacer sequence는 무엇일까? 이것에 대해서 특별히 심각하게 생각해 본 일은 없다. Pseudomolecule이 필요할 때 그저 50~100개 정도의 N을 삽입하여 연결하는 것이 다반사였다. 하지만 컴퓨터를 이용하여 pseudomolecule로부터 gene prediction을 하게 되면 spacer를 내부에 갖는 무의미한 유전자가 예측될 수도 있다.

논문을 읽다가 6개 프레임 전부에서 stop codon을 만나게 되는 36-bp 서열을 spacer sequence로 쓴 사례를 보게 되었다. 대단히 현명하다. 이렇게 pseudomolecule을 만들면 최소한 spacer를 사이에 두고 양 contig를 가로지르는 무의미한 유전자는 나타나지 않을 것이다. 이 방법을 인용한 다른 논문도 최소한 두 편을 확인하였다(아르헨티나 국적의 동일 저자 논문임).
The contig sequences of each strain were concatenated with the sequences NNNNNCACACACTTAATTAATTAAGTGTGTGNNNNN, which puts stop codons in all six reading frames...
[출처] Pyrosequencing-based comparative genome analysis of the nosocomial pathogen Enterococcus faecium and identification of a large transferable pathogenicity island. BMC Genomics 2010 11:239 https://doi.org/10.1186/1471-2164-11-239. Pan-genome 개념을 창시한 H. Tettlin이 이 논문의 공저자이다.

하지만 실제로 해당 genome에 존재하지 않는 서열을 spacer sequence로 사용한다는 점이 약간 불편하다. N을 제외한 중간 염기서열(palindromic 26 bp)을 NCBI BLAST site에서 검색해 보았다. DB는 RefSeq Genome으로 택하였다. 이 서열과 100% 동일한 것을 보유한 genome이 54개였고, Cronobacter의 어느 스트레인을 시퀀싱한 것에는 여러 차례 출현한다. 이 서열을 spacer로 사용하여 이어 붙인 pseudomolecule을 NCBI에 등록했을 가능성이 크다. 실제로 하나를 클릭하여 GenBank record에서 match가 일어나 서열을 보니 좌우에 N이 5개씩 존재한다.

대규모의 comparative genomic analysis를 하는 사람은 이것을 기억해 두는 것이 좋겠다.