- 분노가 사람을 만든다 - 왜 배워야 하는가
- 존잘이 너무 많다 - 시퀀싱을 이용한 실험들
- 까만 바탕에 하얀 글씨-리눅스에서 작업 하기
- 1주차 숙제에 대한 간략하면서도 쓸데없이 긴 해설
- 염기서열을 다뤄보기-포맷 익히기
- 2주차 숙제에 대한 간략하면서도 쓸데없이 긴 해설
- 코딩을 하지 않는 방법
- 공공 자료에 빨대 꽂기
- 어떤 유전자가 열심히 일하나
이 자료는 다양한 분노가 쌓여 작성됐습니다.
제 지인 중 하나는 컴퓨터공학을 공부했고 생명쪽 자료를 잘 다룬다는 것을 늘상 자랑했습니다. 멍청하고 경험 없던 저는 그 말에 홀랑 넘어갔고, 위대하신 박사 지도교수 빛☆준☆호 선생님께 그 인간을 소개시켜드리는 큰 죄를 저질러 버리고 말았습니다. 그 중생은 일을 전혀 하지 않으면서 미루기만 했고, 참다참다 이제는 못 참겠다 싶어 제가 직접 분석하기 시작했습니다. 이제 막 공부하는 과정이었는데도 시작하고 보니 그 새기가 1-2년을 끌어온 일을 한두 달만에 끝낼 수 있게 됐고, 그 새기가 싼 똥 치우고 새로 똥을 쌀 때마다 쌍욕 오지게 하면서 경험치를 쌓게 됐습니다. 그리고 시퀀싱 자료 분석을 통해 만들어낸 유전학 논문으로 박사 졸업을 하게 됐습니다.
당시 저는 정신이 피폐해 주변 생물학자들을 만날 때마다 이런 유형을 신나게 욕하고 다녔는데, 저와 비슷하게 고통을 받았다는 이들이 한둘이 아니었습니다. "데이터 분석 잘 한다고 해서 공동1저자/공동교신저자 걸고 일을 시작했는데 정작 분석을 너무 늦게 해서 진행이 안 됐다." "분석 열심히 했다고 하면서 한참 뒤에 결과 보내줬는데, 그만큼 시간 걸릴 일인 건지 모르겠다." "그림 하나 그려놓고 공동1저자 달라고 하는데, 그만한 분석인지 모르겠다." 세상에 날로 먹으려는 훌륭한 분들이 이렇게나 많았다니! 역시 분석을 할 줄 아느냐가 중요한 게 아니었습니다.
그 연구가 어떤 목적으로 계획됐고, 어떤 실험을 통해 논리를 만들고 있으며, 해당 대용량 자료 분석 결과에 대해 어떤 맥락에서, 어느 정도의 증거력을 기대하면서, 무슨 결론을 내야 하는지, 제대로 이해할 수 있는 사람만이 자료를 정확하게 분석할 수 있다는 지극히 당연한 사실을 깨닫게 됐습니다.
바로 연구자 본인 말입니다.
반대로 제가 공동연구를 하면서 제대로 생산되지 않은 시퀀싱 자료가 차고 넘친다는 것도 알게 됐습니다. 시퀀싱 또한 하나의 실험이고 가설을 검증하는 과정이기 때문에 제대로 된 대조군을 시험해야 하고, 시퀀싱으로 할 수 있는 어떤 실험이 있는지 명확하게 알아야 하며, 각 실험을 대체할 수 있는 것은 또 무엇이 있는지 알아야 한다는 걸 이해하게 된 것입니다. 각 시퀀싱 실험을 할 때 어떤 대조군이 있어야 더 명확하고 선명한 결과를 얻을 수 있는지 알려주는 교육이 없으니 당연한 결과 아니겠습니까? 저만 해도 학부 때나 대학원 때 분석 관련 수업을 배울 수 없었고, 비슷한 수업을 들어본 친구들 얘기로는 자기가 원하는 자료 분석과는 거리가 먼 수업이었다고도 합니다. 결국 독학할 수밖에 없는 형편이었습니다.
그래서 연구실과 공동연구하는 동료들을 위해 자료를 만들기 시작했고, 그렇게 만들어낸 임시방편 교육 자료가 이 교육 자료입니다. 앞으로 여러분은 다음 내용을 공부하게 될 것입니다. 각 자료는 읽는 데 30분에서 1시간, 실습하는 데 2-6시간 정도 들이면 충분합니다. 난이도는 들쭉날쭉합니다만, 제 어머니나 할머니는 못하더라도 지도교수님은 따라할 수 있는 수준(Even my PI!)으로 잡았습니다. 또 이 세상에 존재하는 모든 분석법을 익히는 것을 목표로 하진 않습니다. 그보다는 굉장히 제한적인 DNA와 RNA 시퀀싱 자료 분석법을 익히게 될 것입니다. 대신, 다 끝낸 뒤에는 본인이 원하는 실험에 적절한 새로운 방법론을 스스로 찾고 스스로 분석할 수 있도록 경험을 쌓는 것을 목표로 하고 있습니다. 그때가 되면 적어도 "이게 하면 되는 거였구나!"라는 자신감이 붙게 될 테고, 몇몇은 지금까지 공동연구자들이나 회사가 똑바로 하지 않았다는 것을, 또는 기대 이상으로 정말 잘해주고 있었다는 것을 깨닫게 될지 모르겠습니다.
세상에는 위대한 분들이 넘쳐나고, 그들이 만들어낸 온갖 시퀀싱 실험 기법과 분석 기법이 넘쳐 납니다. 시퀀싱은 그저 DNA나 RNA 등이 어떤 서열을 지니고 있는지 빠르게 읽어내는 방법에 지나지 않지만, 이 두 분자가 워낙에 중요해야죠. DNA 서열 정보, RNA 서열 정보만으로도 알아낼 수 있는 것이 정말 많습니다. 이번 글에서는 널리 쓰이는 시퀀싱 실험 기법과 각 기법으로 살펴볼 수 있는 생명현상이 무엇인지 가볍게 알아보려고 합니다.
현존하는 시퀀싱 기법을 나누는 방법이야 여러 가지가 있겠지만, 저는 크게 짧은 리드(read)와 긴 리드 시퀀싱으로 나눠서 설명하는 편입니다. 안타깝게도 현존하는 기법으로는 염색체 하나를 처음부터 끝까지 통째로 읽어내는 것이 불가능하기 때문에, 특정 길이로 잘린 DNA 등을 읽어내는 수밖에 없습니다. 이렇게 한번에 읽어낸 DNA 또는 RNA를 리드라고 부릅니다. 짧은 리드 시퀀싱은 값이 훨씬 싸고 정확한 대신 100-200 bp 정도로 리드 길이가 짧습니다. 긴 리드 시퀀싱은 값도 비싸고 상대적으로 부정확하지만, 리드 길이가 평균 10,000-100,000 bp에 이를 만큼 훨씬 더 긴 길이를 한번에 읽어낼 수 있습니다. 이 때문에 몇 가지 장점을 지니는데, 복잡한 염색체 지역의 정보를 확인하거나 유전자의 서로 다른 아이소폼(isoform)을 확인할 때 특히 더 좋습니다.
BRIC View, 길이가 긴 리드를 이용한 염기서열분석의 응용, 아래 논문 요약본
Piercing the dark matter: bioinformatics of long-range sequencing and mapping
Long-read human genome sequencing and its applications
한번에 염색체 처음부터 끝까지 이어붙일 수 없다면 우리가 쓰고 있는 표준 유전체는 어떻게 만들어낼 수 있었던 걸까요? 같은 자리를 여러 번 읽고, 퍼즐 맞추듯 겹치는 부분을 이어붙여서 만들어냈습니다. 최근 긴 리드 시퀀싱 기법이 발전하면서 이 과정이 훨씬 더 쉬워졌는데, 이 때문에 약 5조 원 정도를 들여 지구에 살고 있는 온갖 진핵생물의 표준 유전체를 만들겠다는 야심찬 프로젝트가 시작되기도 했습니다. 제 밥벌이 방법론이니 공동연구 필요하신 분은 연락주시면 되겠습니다.
Earth BioGenome Project 홈페이지
Earth BioGenome Project: Sequencing life for the future of life
Darwin Tree of Life 홈페이지
정방향 유전학(forward genetics)은 20세기 중후반부터 온갖 모델생물에 적용되며, 각 형질에 영향을 주는 수많은 유전자를 알 수 있도록 해준 대애단한 방법론입니다. 이는 돌연변이 유발 물질을 처리해 온갖 유전자에 돌연변이를 만들어내고, 그 중 원하는 표현형 변이를 지니는 개체를 골라내, 그 개체가 지니고 있는 돌연변이 중 어떤 유전자에 생긴 돌연변이가 그 표현형 변이를 만들어냈을지 찾아내는 과정을 포함합니다. 이 원인 돌연변이를 찾는 일은 몇 년 전까지만 해도 정말 고통스러운 일이었지만, 예쁜꼬마선충 기준, 작정하고 실험한다고 치면, 이제는 1-2주 정도면 그 원인 유전자를 매우 빠르게 좁혀낼 수 있습니다. 이해하기 가장 쉬운 예는 recessive이고 penetrance가 100%이며 표현형이 non-lethal or non-sterile일 때입니다. 돌연변이를 만들어낸 계통(A)과 유전적으로 거리가 먼 야생 계통(B)을 짝짓기 시켜 F1을 얻어내고, F2를 엄청나게 얻어 그 중 표현형이 나타난 개체만 골라내 하나로 섞어 한번에 DNA를 뽑습니다. 이러면 대부분의 locus는 A 유형과 B 유형이 마구잡이로 뒤섞여 시퀀싱 후 각각이 50%에 가깝게 나타날 테지만, 그 표현형에 관련된 돌연변이가 있는 locus는 모두 똑같이 A 유형으로 선택됐을 테니 특정 지역은 A가 100%로 나타날 것입니다. 만약 recombination이 자주 일어나는 종이라면 F2를 50-100마리만 골라내도 순식간에 원인 돌연변이를 찾아낼 수 있게 됩니다.
관련 분석을 자동화해 그림까지 다 그려주는 MiModD
Caenorhabditis elegans mutant allele identification by whole-genome sequencing (original paper)
정방향 유전학처럼 효과 크기(effect size)가 크고 단일 유전자 변이로 인해 생기는 표현형 변이를 연구할 수도 있지만, 이미 자연에 넘쳐나는 자연 변이(natural variation)를 이용해 해당 표현형 관련 유전자와 그 변이를 찾아낼 수도 있습니다. 가장 대표적인 예는 유전형과 표현형이 조금씩 다른 다양한 야생 개체를 채집하여 표현형을 측정하고 시퀀싱으로 유전형을 확인한 뒤 둘 사이의 상관관계를 확인하는 genome-wide association mapping 방법입니다. 연구실에서 쉽게 키울 수 없는 생물에서도 쉽게 쓸 수 있다 보니 다양한 작물이나 사람에게 널리 적용되고 있는 방법론입니다.
샘플 수가 적은 실험을 계획할 때 RNA 정량화만큼 예쁜 그림 만들 수 있는 실험도 드물 겁니다. 각 유전형에 따라 다르게 나타나는 RNA 양은 그 자체로 좋은 표현형이고, 특정 유전형에 따라 유전자들의 발현이 어떻게 바뀌는지 살핌으로써 해당 유전형이 끼치는 영향까지도 알아낼 수 있습니다. 게다가 나이, 약품 처리 여부, 스트레스 여부 등 온갖 조건을 바꿔가며 손쉽게 살펴볼 수 있다는 장점도 있습니다. 물론 제대로 된 RNA 정량 비교 실험은 그만큼 비싸기도 해서 보통 마지막 그림을 아름답게 장식해서 논문 값을 좀 더 높여주거나 새로운 가설을 제안할 수 있는 용도로 많이 쓰이곤 합니다. 요새는 개별세포 하나하나에서 각 유전자의 발현 정도를 비교할 수 있는 기법도 활발하게 개발되고 있죠. 워낙 많이 쓰는 방법이라 RNA 정량 및 통계처리에서 더 자세하게 다룰 예정입니다.
쉽고도 강력한 기법인 ATAC-seq은 히스톤(histone) 단백질로 감겨 있지 않은 열린 크로마틴 지역을 탐색할 수 있게 해줍니다. 이 실험은 Tn5라는 독특한 transposase를 이용하는데, Tn5는 크로마틴으로 감기지 않은 DNA 서열로 끼어들어가면서 원래 트랜스포존이 그렇듯 양끝에 특정한 DNA 서열을 붙여줍니다(Tagmentation). 히스톤으로 감기지 않은 지역에는 동시에 여러 Tn5 효소가 끼어들어갈 수 있게 되고, 이들이 양끝에 남긴 특정 서열을 이용해 곧바로 시퀀싱할 수 있게 됩니다. 워낙 실험이 편한 데다 가성비가 좋아서 개별세포 수준에서도 다양한 기법이 개발돼 있습니다.
Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position
Tn5 transposase and tagmentation procedures for massively scaled sequencing projects
염색체는 핵 속에서 복잡한 구조(conformation)로 접히고 묶여 있습니다. 이런 구조는 특정 유전자가 더 잘 작동하게 하는 데에 도움을 주기도 하며 핵 속에서 염색체가 작동하는 방식에 대한 이해를 높여줄 것이기 때문에 다양한 실험기법이 개발되고 있습니다. 가장 대표적인 것은 Hi-C라 부르는 실험 기법으로, (1) 핵에 화학 물질을 처리해 가까운 DNA 분자끼리 교차 결합(crosslink)을 형성하도록 하고, (2) 제한효소(restriction enzyme)를 처리해 DNA를 조각 낸 뒤, (3) 라이게이션(ligation)으로 교차결합을 형성한 DNA 분자끼리 이어지도록 만들어 시퀀싱하는 기법입니다. 이런 실험 기법으로 인해 1차원 DNA 서열에서는 분명히 멀리 떨어져 있던 두 지역이 함께 시퀀싱되는 결과가 나타나며, 이를 통해 어떤 DNA 서열끼리 핵 안에서 더 가까이 있는지 추정할 수 있습니다. 노이즈가 매우 심하고 여전히 값이 꽤 비싼 실험이라서, 지금도 더 발전시키려는 노력이 계속 되고 있습니다.
Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome
앞서 언급한 실험 기법들 외에도 세상에는 엄청나게 많은 시퀀싱 실험 기법이 존재하며, 이를 통해 다양한 생명 현상을 high throughput으로 연구할 수 있게 됐습니다. 또 이런 실험 기법을 표준화하려는 노력이 국제 연구 협력체인 컨소시엄(consortium) 수준에서 지속되고 있으며, 가장 대표적인 ENCODE에서는 관련 실험 기법은 물론, 해당 실험 결과를 표준화된 형태로 전처리 및 분석해주는 컴퓨터 실험 기법까지도 제공하고 있습니다.
Perspectives on ENCODE
ENCODE 관련 정보가 담긴 홈페이지
이처럼 이미 존잘님들이 온갖 기법을 개발해두고 분석할 수 있는 길을 마련해둔 상태입니다. 그러니 우리 같은 인민들은 괜히 새로 개발하겠다고 까불고 다닐 일이 아니라, 그저 남들이 잘 갖춰둔 형태로 따라 쓰고 베껴 쓰면 될 일입니다.
혹시 지금까지 컴퓨터 배워보겠다고 파이썬이나 R 문법부터 배우려고 시도했던 경험을 갖고 계실지 모르겠습니다. 확실하게 말씀 드리자면, 처음부터 그렇게 가면 거의 90%는 자빠질 겁니다. 코딩하지 마세요. 컴퓨터로 실험하고 분석하는 법을 가장 빠르게 배우는 길은 코딩하겠다고 까불지 않는 것입니다. 존잘님을 따르십시오. 우리 같은 똥멍청이들이 생각할 만한 건 이미 존잘님들이 다 해뒀으니, 그저 어딘가에 계실 존잘님들을 향해 큰절 세 번 하고 복붙하며 베껴 씁시다.
다음 자료에서는 남의 분석법을 베껴쓰는 데 절실하게 필요한 Linux 커맨드라인 사용법에 대해 배워보겠습니다.
이번에는 리눅스에서 작업하는 방법에 대해 배워볼 것입니다. "왜 윈도우즈를 안 쓰고 이상한 걸 배워야 하나요?"라고 질문할 수도 있습니다. 이유는 단순한데, 직접 개발 안 하고 남들이 개발한 거 가져다 쓰려면 어쩔 수가 없습니다. 거의 대부분의 계산 도구들이 리눅스 기반으로 개발되어 있거든요. 우리 목적은 쉽게 쉽게, 직접 개발하기는 커녕 코딩도 안 하고 분석하는 방법을 익히는 것이니, 쉽고 재미있게 리눅스에서 작업하는 법을 배워봅시다. 이번주부터는 본격적인 실습과 숙제가 시작됩니다.
리눅스에서 작업하는 법을 배우려면 가장 먼저 리눅스 기반의 컴퓨터나 서버를 구해야 합니다. 대부분은 없을 테지만 걱정하지 않으셔도 됩니다. 위대하신 빛혜식 선생님(@hyeshik)께서 대안을 찾아주셨거든요. 항상 충성충성입니다.
먼저 다음 링크로 들어가 회원 가입을 하십시다. 국가수리과학연구소에서 제공하고 있는 사용하시기 전에 대전을 향해 세 번 절하고 쓰시면 좋을 것 같습니다. 매주 월요일인가 리셋된다고들 하니 한번에 몰아서 숙제하시면 되겠습니다.
회원가입 후 로그인하면 왼쪽 위에 Jupyter라고 적힌 창이 뜰 텐데, 오른쪽 위, 로그아웃 아래 쪽에 있는 New 클릭하시면 됩니다. 그 뒤 가장 아래에 있는 Terminal을 클릭하시면 리눅스 연습할 수 있는 화면이 뜹니다. 처음 연결하시면 까만 바탕에
어쩌구@저쩌구:~/work$라고 뜨고, 그 옆에 하얀색으로 커서가 뜰 것입니다. 접속 끝입니다.
이제 리눅스에서 작업할 준비가 끝났습니다. 겁먹을 필요는 없습니다. 제가 지금까지 석사 과정, 박사 과정, 교수급 등 다양한 분들에게 가르쳐봤는데, 이 단계는 모두 쉽게 통과했습니다. 완전 집중하고 실습한다고 가정했을 때, 보통 약 2-4시간 정도 걸리는 것 같습니다.
먼저 가장 기본이 되는 몇 가지 명령어에 대해 배워봅시다. 평소에 윈도우즈 쓰는 것과 엄청나게 다르진 않은데, 마우스 대신 커맨드 라인이라고 부르는 키보드 기반 환경에서 작업한다는 게 가장 큰 차이점입니다. 참고로 직접 서버를 세팅하시거나 맥 등을 이용하신다면 마우스로도 똑같이 작업할 수 있으니 나중에 한번 해보세요. 일단은 커맨드 라인에 익숙해져 봅시다.
커맨드 라인이라는 이름 그대로, 명령어를 한 줄 한 줄 입력해서 컴퓨터에게 일을 시킬 수 있습니다. 가장 간단한 명령어를 몇 가지 직접 쳐봅시다. 먼저 ls를 치고 엔터키를 누르시면 됩니다.
lsls는 list의 약자로, 현재 폴더에 어떤 것들이 들어있는지 보여주라는 명령어입니다.
nano testnano는 메모장을 여는 명령어입니다. 한 칸을 띄운 뒤 특정 단어를 적으면 해당 이름을 지닌 파일을 열어줍니다. 만약 그 이름을 가진 파일이 없다면 새로운 파일을 생성하게 됩니다. ctrl+x를 누르시면 빠져나올 수 있습니다. 윈도우즈와 다른 점이 하나 있는데, 확장자를 굳이 쓰지 않아도 된다는 것입니다. 윈도우즈처럼 .txt 등으로 끝나지 않아도 텍스트 파일입니다.
mkdir test
cd testmkdir은 make a new directory, 즉 새로운 폴더를 만들어주는 명령어입니다. 마찬가지로 한 칸을 띄운 뒤 특정 단어를 적어주면 그 이름을 지닌 폴더를 생성해줍니다. cd는 change your working directory, 즉 특정 폴더로 옮겨주는 기능을 지니고 있습니다. 마찬가지로 한 칸을 띄운 뒤 특정 단어를 입력하면 그 이름을 지닌 폴더로 이동해줍니다. 아무것도 입력하지 않으면 홈(home) 폴더로 이동시킵니다. ~ 이 물결 표시가 여러분의 홈 폴더(/home/id)의 약자입니다.
여기서 중요한 개념을 하나 배워야 하는데, 상대경로(relative path)와 절대경로(absolute path)입니다. 상대경로는 A가 B라는 집을 소개할 때 "B는 우리집에서 동쪽 방향의 옆집에 살지."라고 표현하는 것과 비슷합니다. A의 집을 안다면 손쉽게 찾아갈 수 있다는 장점을 지니고 있습니다. 다만 "B는 동쪽으로 두 건물, 북쪽으로 두 건물, 그 건물에서 3개 층 올라가면 나오는 집에 살고 있다" 등으로 너무 복잡하다면 "서울시 관악구 관악로 1 서울대학교"라고 표현하는 것이 훨씬 더 편할 수도 있습니다. 이렇게 접근하는 방식을 절대경로라고 부릅니다.
먼저 커맨드 라인에 cd를 입력해보십시오. $ 표시 다음에 뜬 명령어만 치시면 됩니다.
foo@bar:~/work/test$ cd # 현재 폴더인 ~/work/test 폴더에서
foo@bar:~$ # 홈 폴더로 이동함
foo@bar:~$ ls # 현재 홈 폴더에는 더 이상 test 폴더가 존재하지 않고, work 폴더만 존재함만약 이 홈 폴더에서 cd test를 친다고 하면 다음과 같은 오류 메시지가 나타날 것입니다.
bash: cd: test: No such file or directory왜냐면 현재 폴더 안에는 test라는 폴더가 존재하지 않기 때문입니다. 이처럼 현재 폴더를 기준으로 폴더 위치를 찾아나가는 방법을 상대경로라고 부릅니다. 현재 폴더보다 아래에 있는 폴더는 그냥 폴더 이름만 치면 된다는 장점이 있습니다. 다음 명령어들을 따라서 쳐봅시다.
foo@bar:~$ # 홈 폴더
foo@bar:~$ cd work # 현재 폴더인 ~ 폴더에서 ~/work 폴더로 이동
foo@bar:~/work$ cd # 다시 홈 폴더로 이동
foo@bar:~$ cd work/test # 현재 폴더에서 두 단계 아래에 있는 폴더로 이동
foo@bar:~/work/test$ cd .. # 현재 폴더보다 한 단계 위에 있는 폴더로 이동
foo@bar:~/work$ cd ../.. # 현재 폴더보다 두 단계 위에 있는 폴더로 이동
foo@bar:/home$그런데 몇몇 경우에는 상대경로를 이용하는 게 굉장히 귀찮을 때가 있습니다. 예컨대 현재 폴더와 이동해야 할 폴더 사이의 관계를 모른다거나, 그 거리가 너무 멀어 ..과 폴더 이름을 매우 여러 번 쳐야 할 때가 그런 예입니다. 이럴 때는 절대경로를 이용해 보다 편하게 이동할 수 있습니다.
절대경로는 현재 폴더가 어디이건 관계 없이 항상 일정합니다. 절대경로는 항상 / 에서 시작하는데, 이 시작점을 루트(root)라고 부릅니다. 다음 경로로 한 번 이동해봅시다.
foo@bar:/home$ cd / # 루트로 이동
foo@bar:/$ cd # 홈 폴더로 이동
foo@bar:~$ # 홈 폴더
foo@bar:~$ cd /opt/conda/bin # 절대경로를 통해 특정 폴더로 이동
foo@bar:/opt/conda/bin$ # 홈 폴더 아래에 있는 폴더가 아니라서 ~로 시작하지 않음윈도우즈에서 쓰는 파일 시스템과 비슷하죠? 마우스 클릭 대신 cd로 이동한다는 것이 다를 뿐입니다. 참고로, 키보드에 있는 Tab 키를 이용하면 이동이 더 쉬워집니다.
foo@bar:/opt/conda/bin$ cd # 홈 폴더로 이동
foo@bar:~$ cd w # w만 치고 Tab 키 누르기
foo@bar:~$ cd work/ # 자동완성됨이외에도 굉장히 많은 명령어가 있습니다만, 다 가르쳐드리면 실력이 늘질 않으니 익혀야 하는 리스트를 드리겠습니다. 대부분의 명령어는 ls --help처럼 명령어를 치고 한 칸 띄운 뒤 --help를 치면 간략한 설명서를 띄워줍니다. 여기서도 안 나온다면 (물론 그럴 일은 거의 없습니다) 구글에 열심히 검색하면서 찾아보시면 되겠습니다. 대부분이 똑같은 얘길 하면서 알려줄 것입니다. 바로 RTFM, "Read the f****** manual"의 약자입니다. 항상 매뉴얼 먼저 읽고 따라하는 습관을 생활화합시다. 물론 매뉴얼이 이해가 안 될 가능성도 매우 높습니다만, 그럴 땐 구글에 검색해서 찾아봅시다(리눅스 ls 사용법 또는 linux ls usage 등등). 여러분과 같은 중생들이 도처에 널려 있고, 깨달음을 얻은 부처들이 욕은 하면서도 결국 여러분을 구제해줄 것입니다.
명령어를 그냥 입력하면 지루하게 배우겠죠? 텍스트를 읽으면서 좀 더 재밌게 명령어를 익힐 수 있는 게임이 있으니 그걸 이용해봅시다. 이름은 보물 찾기(Treasure Hunt)로, UCSF에 계시는 Stephan J. Sanders 박사가 개발한 것입니다. 제가 번역 초안을 맡았고 시바님이 감수해주셨습니다. 사용법은 다음과 같습니다. 일단 홈 폴더로 이동한 뒤 다음 명령어를 치고 엔터를 눌러주세요.
mkdir treasure && cd treasure
git clone https://github.com/JunKim-Worms/unix_tutorial.git
perl unix_tutorial/treasureHunt_v2_kor.pl
rm -rf unix_tutorial/
각 줄(line)이 어떤 의미인지는 차츰 알게 될 겁니다. 예컨대 세 번째 줄에 있는 treasureHunt_v2_kor.pl 파일이 보물 찾기 게임을 실행해주는 파일입니다. 이 게임은 Perl이라는 프로그래밍 언어로 짜였습니다. .pl 이라는 확장자로 끝나면 perl file.pl 등으로 실행할 수 있습니다.
파일들이 잘 생성됐다면 이제 nano 또는 cat을 이용해 단서를 열고 읽어봅시다.
cat Clue01_S.txt # or nano Clue01_S.txt원래 보물찾기가 끝나면 선물을 달라고 하시면 되는데, COVID-19이 끝날 즈음 서울대로 찾아오시면 차라도 한 잔 대접하겠습니다.
보물찾기가 끝나면 다음 링크를 클릭해서 들어가서 시키는 대로 하시면 됩니다.
1주차 숙제
숙제 파일은 다음 명령어를 커맨드 라인에 입력하면 다운 받을 수 있습니다.
git clone https://github.com/JunKim-Worms/sequencing_for_genetics.git
cd sequencing_for_genetics/week1/참고로 저는 처음에 리눅스 컴퓨터 세팅하는 것부터 배웠는데, 고통스러웠지만 상당히 재밌는 작업이었습니다. 이번 튜토리얼을 다 끝낸 뒤 본격적으로 작업해보고 싶으신 분들은 어딘가 남는 컴퓨터 하나 있다면 싹 밀고 우분투 18.04 설치 등을 검색해서 설치해보셔도 좋을 것 같습니다. 저는 익숙해지는 데 3-6개월 정도 걸렸습니다만, 좀 더 빠르게 해보실 수도 있을 겁니다. 연구실에 직접 세팅하고 쓰실 분들은 업체에 설치 및 관리 등을 맡기시는 것도 좋습니다. 선충 분석용 컴퓨터는 약 500-700만 원 정도 들 텐데, 그 정도면 제가 냈던 Genome Research 논문에 들어간 분석은 아주 빠르게 할 수 있는 수준입니다.
일단 숙제는 다 하셨겠죠? 안 했다면 다시 위로 올라가 숙제부터 하고 오도록 합시다. 숙제를 하지 않으면 눈곱만큼도 늘지 않습니다. 그러면 전혀 쓸데없고 재미도 없는 글만 읽으며 시간 낭비하는 꼴밖에 되지 않으니 차라리 잠이나 잡시다.
다양한 기본 명령어를 익히는 게 첫 숙제였는데, 다들 한번씩 직접 쳐보셨길 바랍니다. 모두 제가 매일 같이 쓰는 명령어들이고 굉장히 유용하기 때문에 자주 사용하면서 익숙해지도록 노력해봅시다.
몇 가지 헷갈리거나 차이가 뭔지 알기 어려운 명령어들만 간단하게 설명해드리겠습니다. 일단 ctrl + c와 ctrl + z는 정지와 일시정지의 차이입니다. ctrl +z를 눌러서 일시정지 된 명령어는 fg를 눌러 재개할 수 있습니다. 또는 bg를 쳐도 되는데, 이렇게 하면 해당 명령어가 백그라운드(background)에서 실행되기 때문에 계속해서 새로운 명령어를 입력하는 작업을 할 수 있습니다. 뭔 말인지 궁금하다면 직접 한번 해보세요. 참고로 ctrl + r은 최근에 썼던 명령어를 검색할 때 쓸 수 있는데, history를 쳐서 최근 명령어 리스트를 쭉 부르는 것보다 좀 더 간편합니다.
ls *에서 *은 **와일드카드(wildcard)**라고 부르는 굉장히 유용한 문자입니다. 이 *은 어떤 문자로도 대체할 수 있는데, 예컨대 sequencing_for_genetics/week1 폴더로 이동해봅시다. 그 뒤 ls *fa를 치면, 무슨 파일이든 fa로만 끝나면 전부 리스트로 나타난다는 것을 확인할 수 있을 것입니다.
ls
# 2.0.bed cb4856.substr.fa hw1.md lsp41_scaffolds.substr.fa
# 4개 파일 모두 나타남
ls *fa
# cb4856.substr.fa lsp41_scaffolds.substr.fa
# 2개 파일만 나타남이 와일드카드는 다양한 명령어 뒤에 쓰일 수 있는데, 예컨대 fa로 끝나는 파일만 골라서 지우고 싶다면
rm cb4856.substr.fa
rm lsp41_scaffolds.substr.fa이렇게 두 번 쳐도 되지만, 한번에 지울 수도 있습니다.
rm *fa보물찾기 숙제에서도 말씀드렸습니다만, rm은 함부로 쓰지 않도록 주의하도록 합시다. rm * 등 치다가 엔터 잘못 누르기라도 하면 정신도, 졸업도 날아가는 수가 생깁니다.
모든 문자를 대체하는 와일드카드 말고 숫자나 특정 알파벳만 대체하는 와일드카드도 만들 수 있습니다. 예컨대 ls [0-9]*이라고 치면 숫자로 시작하는 파일, 즉 2.0.bed 파일만 뜨게 될 것이고, ls *[2-4]*을 치면 뭘로 시작하고 끝나든 숫자 2부터 4까지 들어가는 파일, 즉 2.0.bed cb4856.substr.fa lsp41_scaffolds.substr.fa 세 개 파일 이름이 출력될 겁니다. ls [a-c]*이라고 치면 알파벳 a, b, c로 시작하는 파일 이름이 출력될 것입니다. 이런 것들에 관심이 간다면 와일드카드 또는 정규표현식(regular expression) 등을 검색해보시면 좋을 것 같습니다.
텍스트 파일을 여는 방법도 여러 가지가 있습니다. 찾아서 편한 걸 쓰시면 되는데, 저는 새로 배우기 너무 귀찮아서 그냥 파일 크기가 작으면 nano로 열고, 파일 크기가 크면 less로 열곤 합니다. 컴공과 친구들이 다른 에디터 쓰라고 몇 번이나 강요하고 조롱하고 비난하였으나 너무나 귀찮아서 그냥 쓰고 있습니다.
파이프(pipe)와 리디렉션(redirection) 또한 매우 유용합니다. 스탠다드인풋을 받아 다음 명령어가 처리하도록 넘기거나 파일로 저장하는 용도라고 생각하면 좀 쉬운데, 예컨대
cat cb4856.substr.fa만 치면 256줄짜리 매우 긴 파일이 화면에 다 뜨게 됩니다. 때로는 이 중 첫 5개 염색체/컨티그/리드의 이름만 보고 싶을 때도 있겠죠? 이럴 때 파이프와 리디렉션을 적절하게 사용해주면 편합니다.
grep ">" cb4856.substr.fa | head -5
# 화면에 다음 다섯 줄이 출력됨.
# >tig00000001|quiver|quiver|pilon|pilon
# >tig00000022|quiver|quiver|pilon|pilon
# >tig00000031|quiver|quiver|pilon|pilon
# >tig00000033|quiver|quiver|pilon|pilon
# >tig00000042|quiver|quiver|pilon|pilon
grep ">" cb4856.substr.fa | head -5 > name_list_5
# 파일로 저장됨grep, awk, sed 등은 정말 정말 많이 쓰는 명령어입니다. 이들은 모두 행렬을 다뤄주는 도구들인데, 이외에도 cut 등을 사용하면 행렬에 들어있는 문자열을 검색하거나 연산하는 것이 아주 손쉽게 가능해집니다. 마찬가지로 필요하다면 검색해서 한번 살펴보도록 합시다.
보물찾기도 잘하셨으리라 믿습니다. 다 하신 분은, 그리고 저랑 차 한 잔 하실 분은 서울대 오실 때 연락주세요. 간단하게라도 대접하겠습니다.
남은 숙제도 살펴보면서 중요한 마음가짐을 새기도록 합시다.
컨티그(contig) 수를 세려면 > 로 시작하는 줄(line) 수를 세면 됩니다. 세는 방법은 정말 많고, 다음 명령어들은 모두 같은 결과를 보여줍니다.
grep -c ">" lsp41_scaffolds.substr.fa
grep ">" lsp41_scaffolds.substr.fa | wc -l
grep ">" lsp41_scaffolds.substr.fa > tmp; wc -l tmp
grep ">" lsp41_scaffolds.substr.fa # 화면에 뜨는 거 워드 파일로 복사해서 줄 수 세기중요한 건 어떤 명령어를 치든 상관이 없다는 사실입니다. 심지어 nano lsp41_scaffolds.substr.fa를 친 다음 직접 세도 아무 문제가 되지 않습니다. 정확하게만 할 수 있다면, 몸은 좀 힘들어도 됩니다. 코딩 잘하게 되면 시작할 생각하지 말고, 그냥 하세요. 그냥 결과를 낼 수 있고, 논문에 쓸 데이터를 만들어낼 수 있는 게 중요합니다. 처음에는 삽질 좀 해도 상관 없습니다. 제가 컴퓨터 앞에 앉아서 명령어 하나 끝나면 다음 명령어 엔터 치는 꼴을 컴공과 친구가 본 적이 있는데, 그때 그 친구가 안쓰럽게 쳐다봤지만 결국 저는 그걸로 논문 냈으니 된 거 아니겠어요? 그리고 그 과정에서 엄청나게 빠르게 늘었습니다. 이제는 훨씬 더 빠르고 정확하게 해낼 수 있지만, 그때는 별 수 없었던 거예요. 100줄이든 1000줄이든 그냥 하나하나 엔터 치면서 정성스레 하셔도 됩니다. 하다 보면 늘 거예요.
컨티그 이름에서 |quiver|quiver|pilon|pilon를 제거하는 숙제는 sed를 이용해 해결할 수 있습니다. 물론 nano로 파일 열고 직접 지우셔도 됩니다만, 좀 더 편하게 가보도록 합시다. 한 줄이면 됩니다. sed 's/|quiver|quiver|pilon|pilon//' 어떤 경우에는 sed 's/지울 단어/바꿀 단어/g' 로 끝에 g가 들어가는 경우도 있는데, 이건 같은 줄에 지울 단어가 연달아 있을 때 싹 다 바꾸라는 뜻입니다. 예컨대 acacacac라는 문자열에 sed 's/ac/gt/'를 적용하면 gtacacac로 바뀌지만, sed 's/ac/gt/g'를 치면 acacacac로 바뀌게 됩니다.
echo acacacac | sed 's/ac/gt/'
# gtacacac
echo acacacac | sed 's/ac/gt/g'
# gtgtgtgt여기서 |quiver|quiver|pilon|pilon라는 단어를 없애거나 다른 걸로 치환하는 이유는 | 이 파이프 기호가 이름에 들어가 있을 때 해당 이름을 처리하지 못하는 프로그램이 꽤 많기 때문입니다.
2번 염색체를 잡아내는 방법도 마찬가지로 여러 가지가 있습니다. grep "II"를 치면 2번 염색체에 대한 정보 뿐만 아니라 3번 염색체(chromosome III)에 대한 정보도 함께 딸려나오기 때문에, 이것만 제거해주면 됩니다.
grep "II" 2.0.bed | grep -v "III" # 2번 들어가는 거 찾은 뒤 3번 제거
grep -w "II" 2.0.bed # 정확하게 II 라는 단어가 들어가는 줄 찾기
grep -P "^II\t" 2.0.bed # II 라는 단어로 시작(^)하는 줄을 찾고, 곧바로 탭(\t)이 들어가는 줄 찾기어떻게든 결과만 내면 되겠죠?
파일을 열어보시면 알겠지만 Insertion 등 SV 유형에 대한 정보는 일곱 번째 컬럼에 들어있습니다. awk '{print $7}' | sort | uniq를 쳐서 어떤 유형이 존재하는지 살펴보도록 합시다. SV 크기에 대한 정보는 다섯 번째 컬럼에 들어있습니다. sort -k5rV | head -5를 치면 가장 큰 SV 5개를 확인할 수 있습니다.
참 쉽죠? 이제 본격적으로 염기서열 정보에 대해 살펴보도록 합시다.
컴퓨터로 일하는 과정 상당 부분은 특정한 자료 포맷(format)을 다루는 일에 가깝습니다. 다루려는 자료를 인풋으로 프로그램에 던져넣으면 그 프로그램이 특정한 포맷을 지닌 아웃풋으로 내놓게 됩니다. 워낙 생물학하는 사람들이 많다 보니 정말 많은 포맷이 있습니다만, 염기서열과 관련해서는 몇 가지 정해진 포맷이 있어 기억해둘 필요가 있습니다. 먼저 시퀀싱 방법에 대해 알아보고 그로부터 만들어지는 염기서열 결과 파일에 대해 살펴보도록 합시다.
세상에는 정말 많은 시퀀싱 기법이 있지만, 두 번째로 많이 쓰는 방법은 아마 일루미나(Illumina) 사에서 제공하는 페어드 엔드 시퀀싱(paired-end sequencing) 기법일 겁니다. 이는 각 DNA 분자의 양끝을 읽어내 동시에 사용할 수 있도록 해줍니다.
생물 안에 들어있는 DNA는 정말 긴 분자입니다. 이 긴 DNA를 한번에 정확하게 읽어내기란 아직까진 불가능하고, 아마 한동안 어려울 것 같습니다. 이 때문에 애초에 싹 읽는 걸 포기하고, 한번에 길게는 150 bp 정도만 읽되 정확도를 높인 기법이 널리 쓰이고 있습니다. 염색체처럼 기다란 DNA를 300-500 bp 정도의 적당한 길이로 잘라내고, 이 DNA 분자의 양끝을 읽어내는 것이죠. 이렇게 양끝을 모두 읽어내는 방식을 페어드 엔드 시퀀싱이라 부릅니다. 한쪽만 읽는 싱글 엔드 시퀀싱도 쓰이긴 하는데, 고생물 다루지 않는 이상 거의 쓸 일 없을 겁니다.
페어드 엔드 시퀀싱 기법을 사용하면 DNA 한 분자당 150 bp씩 2개의 염기서열 정보가 생산되게 됩니다. 이 두 염기서열 정보는 굉장히 가까운 위치에 있었던 만큼, 이 두 개를 한번에 활용할 수 있다면 150 bp짜리가 아니라 300-500 bp짜리 정보를 통째로 사용할 수 있겠죠? 그래서 두 정보를 한번에 사용할 수 있도록, 두 파일을 하나의 데이터로 사용하는 방식을 택하고 있습니다. 실제 시퀀싱 결과를 받아보면, 보통 실험 하나를 맡겨도 2개 파일이 동시에 옵니다. 예컨대 nematode_1.fastq와 nematode_2.fastq 등으로 말이죠. 그리고 이 2개 파일, _1과 _2가 들어가는 두 파일이 한 쌍의 데이터입니다. 때로는 .1과 .2로 짝을 이루기도 하고, f와 r로 짝을 이루기도 합니다. 시퀀싱 과정이 짝지어져있었듯(paired), 시퀀싱 결과 파일도 짝지어져 있는 셈입니다.
각 파일의 첫 줄을 읽어보면 흥미로운 점을 확인할 수 있습니다. 두 파일의 첫 번째 줄이 똑같다는 점입니다.
head -1 nematode_1.fastq
# @SRR1234567.11111
head -1 nematode_2.fastq
# @SRR1234567.11111둘째 줄부터 넷째 줄까지는 또 두 파일이 서로 다른 정보를 담고 있지만, 다섯째 줄도 서로 같습니다. 첫째 줄부터 시작해 4의 배수로 두 파일은 같은 정보를 담고 있는데, 이 정보는 한번에 읽어낸 염기서열인 리드(read)의 이름을 가리킵니다. 한 파일 안에는 같은 이름이 존재하지 않지만, 두 파일에는 같은 이름이 존재하고 이 같은 이름은 두 염기서열 정보가 1개의 DNA 분자 조각에서 나왔다는 것을 뜻합니다.
이제 한 파일 안에 들어있는 정보에 대해 살펴보도록 합시다. 4줄씩 1개 리드의 정보를 담고 있는데, 각 줄은 다음과 같은 정보를 포함하고 있습니다.
@SRR3441150.2241699 2241699 length=90 # 리드의 이름. 항상 @로 시작한다.
CATTATAAGAATGTAGATATGCTTTGGGAATTGTTTCAATAAACCTTAAGCAATCAGTGCGACTGATTGAGAACAAACCCTGTTCGCATT # 리드의 염기서열. DNA라면 대개 ATGC로 이뤄져 있다.
+ # 보통 + 하나만 기록돼 있음. 어떤 경우엔 리드 이름 정보가 똑같이 기록되기도 함.
BCCFDFFFHDFHH4EFHIGGIJIIGIJJIGIIJDFGHGIHGIEGIEHJIJJJJJJJJGHIGJIEGHHEGCHIIHFCHFFFFEDEEDBBDD # 해당 리드가 지닌 각 염기의 품질값(quality score).이 중 리드 염기서열 부분과 품질값도 한 쌍을 이룬다는 것을 기억해두면 좋습니다. 각 품질값 문자열 1개는 같은 순서에 자리 잡은 염기가 얼마나 정확한지 추정할 수 있는 지표가 됩니다. 이 값은 -10 log P_error이라고 생각하시면 되는데, 예컨대 처음 3개 염기서열인 CAT가 얼마나 정확한지를 확인하려면 같은 자리에 있는 BCC를 살펴보면 됩니다. 여기서 B와 C는 문자가 아니라 품질값 33과 34를 가리키는 것이며, 이는 염기 C가 틀릴 가능성은 10^-3.3 이고 그 뒤에 나오는 A와 T는 틀릴 가능성이 10^-3.4 라는 것을 의미합니다. 품질값은 보통 0에서 40까지 나타내는데, 무조건 틀릴 경우나 0.01%로 틀릴 경우 사이에 들어간다는 뜻이 됩니다.
같은 짝을 이루는 두 번째 파일, 즉 _2.fastq 파일도 거의 똑같이 생겼습니다만, DNA의 다른 부분을 읽어낸 것이다 보니 리드 염기서열 정보와 품질값 정보는 값이 다릅니다.
@SRR3441150.2241699 2241699 length=90
CATTATAAGAATGTAGATATGCTTTGGGAATTGTTTCAATAAACCTTAAGCAATCAGTGCGACTGATTGAGAACAAACCCTGTTCGCATT
+
@CCFDFFFHDFHH4EFHIGGIJIIGIJJIGIIJDFGHGIHGIEGIEHJIJJJJJJJJGHIGJIEGHHEGCHIIHFCHFFFFEDEEDBBDD이처럼 4줄로 리드의 이름, 염기서열 정보, 품질값 등을 표기하는 포맷을 FASTQ 포맷이라 부릅니다. 패스트큐 등으로 읽고, 보통 쉽게 알아볼 수 있도록 .fq 또는 .fastq 등을 파일 끝에 붙여둡니다. 만약 특정 염기의 퀄리티가 낮다면 그 부분을 없애야 오히려 더 정확한 정보를 얻을 수 있겠죠? 이처럼 시퀀싱이 개판으로 된 부분을 제거하는 용도로 가장 많이 쓰이는 포맷입니다.
어떤 때는 퀄리티 정보따윈 필요없고 이름과 염기서열 정보만 필요한 경우도 있습니다. 이전에 살펴봤던 cb4856.substr.fa나 표준 유전체 같은 경우가 대표적입니다. 다시 한 번 해당 파일의 앞부분만 head -4로 살펴봅시다. FASTQ 포맷과 달리, @로 시작하지 않는다는 것을 확인할 수 있을 겁니다.
>tig00000001|quiver|quiver|pilon|pilon
cagagtcacttatttggtgaacggcc...
>tig00000022|quiver|quiver|pilon|pilon
aaacctaatgagaccccaaaatactg...이처럼 @ 대신 >를 이름 구분자로 사용하고, 그 외 나머지 줄에는 모조리 염기서열 정보를 때려넣는 포맷을 FASTA 포맷이라 부릅니다. 패스트에이, 파스타 등으로 읽고 .fa 또는 .fasta 등을 파일 끝에 붙여둡니다. 어떤 경우에는 한 줄에 염기서열 60개만 들어가도록 엔터 쳐서 좀 더 보기 편하게 바꿔놓은 경우도 있는데, 이렇게 엔터를 쳐놔도 다음 > 기호가 나오기 전까지는 모두 1개의 연속된 염기서열 정보를 가리킨다는 것을 기억해두시면 좋습니다.
본격적인 분석을 해보기 전에, 자그마한 데이터를 이용해 장난질을 한번 해봅시다. 일단 뭐라도 장난을 쳐보려면 FASTA나 FASTQ를 다룰 줄 알아야 합니다. 둘 다 특정 구분자로 구별되고, 각 줄에 정보를 담고 있으니까 이전 숙제에서 써봤던 awk, sed, grep 등등을 이용하면 행렬로 접근해서 살펴볼 수 있지 않을까요? 물론 가능합니다. 그러나 세상에는 이미 존잘들이 만들어둔 정말 유용한 프로그램이 많습니다. 이쯤에서 RTFM만큼이나 위대하고 중요한 격언을 마음에 새로 새깁시다. Don't reinvent the wheel. 이미 위대하신 분들이 다 짜놨는데, 우리 같은 똥멍청이들이 새로 짠다? 머리 쓰지 말고 산소를 아끼면서 기후위기를 막는 데 동참합시다. 코딩을 하지 말고, 존잘들이 만들어둔 프로그램을 갖다 쓰는 것이야말로 빠르게 성장하는 길입니다.
이번 숙제에서는 FASTA와 FASTQ 파일을 다룰 수 있는 몇 가지 프로그램을 배우게 될 겁니다. samtools, bioawk, seqtk 등 온갖 도구들을 만들어둔 Heng Li 선생님이 계신 미국을 향해 큰 절 세 번 올리고 시작하도록 합시다.
이번 숙제를 하시려면 https://jupyter.nims.re.kr/ 에 접속하신 뒤, 터미널로 다시 들어갑시다. 접속 후에는 다음 명령어를 한 줄씩 입력해주세요. 필요한 프로그램과 다뤄볼 자료가 자동으로 설치될 겁니다.
conda init
source ~/.bashrc
conda create -n seq
conda activate seq
conda install -c bioconda assembly-stats samtools bioawk seqtk hisat2
git clone https://github.com/JunKim-Worms/sequencing_for_genetics.git
cd sequencing_for_genetics/week2/2주차 숙제로 들어가 설명을 읽고 따라해주시면 됩니다.
숙제 다 하셨겠죠? 숙제를 하면 반드시 늘게 되어 있습니다. 숙제를 꼭 해주시기 바랍니다.
2주차 숙제는 FASTA와 FASTQ를 다룰 수 있는 다양한 프로그램을 직접 써보는 것이었습니다.
먼저 리드 수, 평균값, N50, 가장 긴 리드 길이 등의 정보는 assembly-stats를 사용해 쉽게 뽑아낼 수 있습니다. 이 프로그램을 쓰려면 반드시 압축을 해제해야 한다는 걸 기억해두시면 좋습니다.
gzip -dk cb4856.pacbio.subseq.fq.gz # 또는 gunzip -k cb4856.pacbio.subseq.fq.gz
assembly-stats cb4856.pacbio.subseq.fq그러면 다음과 같은 창이 뜰 겁니다.
sum = 859272, n = 100, ave = 8592.72, largest = 22593
N50 = 11696, n = 29
N60 = 10824, n = 36
N70 = 9134, n = 45
N80 = 7058, n = 56
N90 = 4886, n = 70
N100 = 261, n = 100
N_count = 0
Gaps = 0여기서 sum은 총 리드 길이의 합을 가리키며 n은 개수를 나타냅니다. ave는 평균 리드 길이, largest는 가장 긴 리드 길이 등을 나타냅니다.
몇 가지는 bioawk를 이용해서도 해결할 수 있습니다.
bioawk -c fastx 'END{print NR}' cb4856.pacbio.subseq.fq.gz # 총 리드 개수
bioawk -c fastx '{sum+=length($seq); a+=1}END{print sum, a, sum/a}' cb4856.pacbio.subseq.fq.gz
# 순서대로 총합, 총 리드 개수, 평균값
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -1 # 가장 긴 리드 길이N50 같이 귀찮은 metric은 직접 계산하는 것도 가능합니다. 처음 논문 쓸 때만 해도 N50 구하는 거 찾기가 너무 귀찮아서 훨씬 귀찮게 직접 계산한 적도 있습니다. 예컨대 리드 길이를 쭉 뽑아놓고, 위쪽에서부터 몇 번째까지 더했을 때 전체 길이의 절반을 넘는지 확인하는 식으로 하나하나 일일이 쳤습니다.
# 총 길이 859272니까 절반은 429636
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -50 | awk '{sum+=$1}END{print sum}'
# 647651 한참 넘침
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -20 | awk '{sum+=$1}END{print sum}'
# 330685 모자람
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -30 | awk '{sum+=$1}END{print sum}'
# 452128 넘침
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -28 | awk '{sum+=$1}END{print sum}'
# 428998 모자람
bioawk -c fastx '{print length($seq)}' cb4856.pacbio.subseq.fq.gz | sort -rV | head -29 | awk '{sum+=$1}END{print sum}'
# 440694 처음으로 절반 넘었으니 29번째 값이 N50이 됨어떻게든 답만 도출할 수 있으면 됩니다.
특정 길이 이상의 리드에 대해서만 살펴보는 것도 여러 방법을 통해 실험해볼 수 있습니다. 제일 쉬운 건 특정 길이 이상의 리드만 따로 추출해서 살펴보는 거겠죠. awk와 마찬가지로 bioawk에서도 조건을 잡아줄 수 있습니다.
bioawk -c fastx 'length($seq)>9999{print "@"$name; print $seq; print "+"; print $qual}' > 10k.fq # 10 kb 이상의 리드 추출
assembly-stats 10k.fq숏리드 시퀀싱 결과도 마찬가지로 다뤄볼 수 있습니다.
gzip -dk srr3440952.sub.1.fq.gz
assembly-stats srr3440952.sub.1.fq숏리드 시퀀싱은 대개 특정 길이를 넘어가면 정확도가 엄청나게 떨어지기 때문에 리드 길이를 제한합니다. 여기서는 100 bp로 정해놨나 보네요. 거의 모든 리드가 100 bp인 걸 볼 수 있습니다.
시퀀싱 결과는 파일 크기가 워낙 커서 매번 모조리 다루기는 어려운 경우가 있습니다. 그래서 처음에 간단한 실험들을 해보거나 새로운 코드를 짜볼 때는 일부 리드만 뽑아내서 다뤄보기도 합니다. 예컨대 지금 2주차 숙제에서 다루고 있는 롱리드 데이터와 숏리드 데이터 모두 제가 특정 리드만 추출해둔 겁니다. 원본 파일 크기가 3 GB도 넘는 것들이라 어디 올려두기도 힘들거든요. 이런 일도 다양한 방식을 통해 해볼 수 있는데, 예컨대 seqtk에 포함된 sample을 이용해서 추출할 수도 있습니다.
seqtk sample # 요것만 치면 사용법인 usage와 사용 가능한 옵션이 다옵니다.
seqtk sample srr3440952.sub.1.fq.gz 100 > n100.1.fq
seqtk sample srr3440952.sub.2.fq.gz 100 > n100.2.fq참고로 가능한 옵션 중 -s가 있는데요, 설명에는 RNG seed [11]이라고 적혀 있습니다. 여기서 RNG는 random number generation의 약자로, 난수표라고 생각하시면 됩니다. 컴퓨터 안에서 완벽한 랜덤을 구현하기는 매우매우 어려워서 랜덤하게 뽑아둔 숫자가 적힌 난수표를 이용해 대체하는 경우가 많습니다. 이 중 seed는 다양한 난수표 중 어떤 걸 쓸 건지를 정해달라는 뜻이고, [11]은 그 중 11번째 난수표를 기본값으로 사용한다는 뜻이라고 대충 알아두심 됩니다. 참고로 srr3440952.sub.1.fq.gz과 srr3440952.sub.2.fq.gz은 짝지어진(paired) 데이터이므로, 똑같은 RNG seed를 줘야만 추출된 리드들도 짝지어져 나오게 됩니다. 예컨대
seqtk sample -s 100 srr3440952.sub.1.fq.gz 100 > n100.1.fq
seqtk sample -s 101 srr3440952.sub.2.fq.gz 100 > n100.2.fq이렇게 하면 두 파일이 짝지어진 상태로 나오지 않아 문제가 생깁니다. 또 논문 쓸 때 seed 값을 명확하게 표기하는 경우도 종종 있는데, 몇몇 양아치 과학자들이 자기들이 원하는 결과 얻으려고 자기들한테 필요한 일부만 뽑아놓고는 "응 랜덤~"이라고 속이는 경우가 있기 때문입니다. 그래서 seed 값을 표기해서 재현 가능한 결과로 만들어내기도 합니다. "동작그만, 밑장빼기냐? 패 건들지 말어! 손모가지 날아가붕께!" 요새 컴퓨터가 워낙 빨라서 리뷰어 선에서 바로 재현 안 된다는 쌍욕 날아올 수도 있습니다. 게으른 주작은 하지도 맙시다.
물론 훨씬 더 귀찮게 코드를 짜는 법도 있습니다.
paste srr3440952.sub.1.fq srr3440952.sub.2.fq | \ # 두 FASTQ 파일을 컬럼 2개인 1개 파일로 합침
awk -F "\t" '{ printf("%s",$0); n++; \ # 모든 컬럼을 문자열로 받고, n은 1씩 커짐. 인풋 field seperator는 탭
if(n%4==0) {printf("\n")} else {printf("\t")} }' | \ # 나머지가 4이면 엔터, 그 외에는 탭을 쳐서 paired end 리드 정보를 한 줄로 모아줌
shuf -n $NUM | \ # 랜덤하게 섞고 $NUM 숫자만큼의 line을 추출
awk -F "\t" '{ OFS="\n"; \ # 아웃풋 field seperator는 엔터
print $1,$3,$5,$7 >> "sub.1.fastq"; \
print $2,$4,$6,$8 >> "sub.2.fastq" }'이런 걸 이해하느니 seqtk sample 사용법을 잘 숙지하는 게 낫겠죠? 바퀴를 재발명하는 일은 가능하면 하지 맙시다.
마이토콘드리아 지놈도 한번 살펴봅시다. 크기는 assembly-stats 등을 이용하면 쉽게 확인할 수 있습니다. 시퀀스 대소문자를 통째로 바꾸는 건 bioawk에서도 가능합니다.
bioawk -c fastx '{print ">"$name; print tolower($seq)}' wbcel235.mt.fa > wbcel235.mt.lower.fa
bioawk -c fastx '{print ">"$name; print toupper($seq)}' wbcel235.mt.fa > wbcel235.mt.upper.faawk에도 tolower와 toupper가 존재하는데, bioawk와 달리 컨티그 이름까지도 대소문자를 바꿔버리니 주의하시기 바랍니다.
awk '{print tolower($1)}' wbcel235.mt.fa > lower.fa
awk '{print toupper($1)}' wbcel235.mt.fa > upper.fa
# MtDNA라는 이름이 mtdna 또는 MTDNA로 바뀌게 됨특정 구간의 시퀀스만 뽑아낼 때는 samtools가 유용합니다. samtools에 딸린 하위 프로그램 중 faidx를 쓰면 좋은데, FASTA 파일을 다루는 데 쓰입니다.
samtools faidx wbcel235.mt.fa MtDNA:100-200만약 awk로 컨티그 이름까지 바꿔버렸다면 MtDNA 대신 mtdna 또는 MTDNA 등으로 바꿔주셔야 합니다. 대소문자 등을 다 신경 쓰기 때문에 안 돌아갈 거나 다른 시퀀스 정보를 뽑아낼 겁니다. 이름도 중요하니 항상 주의합시다.
hisat2를 이용해 리드 매핑하는 것도 간단하게 해볼 수 있습니다. 마찬가지로
hisat2-build
hisat2이렇게 명령어만 치시면 설명이 쭉 나옵니다. 일단 hisat2-build는 리드 매핑을 훨씬 더 빠르게 하기 위해서, 표준 유전체의 인덱스(index)를 만들어줍니다. 예컨대 사람 유전체 크기는 3 Gb인데, 리드 하나를 집어다 "음 얘는 어디에 붙을까?"라고 하면서 처음부터 끝까지 쭉 맞춰보는 짓을 했다간 몇 년이 지나도 졸업논문을 못 쓸 겁니다. 인덱스를 만들어주면 빠른 검색을 할 수가 있어서 속도를 엄청나게 높여주고 졸업 가능한 수준으로 시간을 단축시켜 줍니다. 궁금하시면 원 논문을 살펴보시면 좋습니다.
여튼 hisat2-build를 실행하고 나면 mt라는 접두어(prefix)가 달린 파일이 쭉 만들어집니다. 그러고 나면 hisat2를 이용해서 실제 리드 매핑을 수행할 수 있습니다. 설명서에 나와있듯 -x 다음에는 앞서 만든 인덱스 파일의 접두어(여기서는 mt)를 입력해주면 되고, -1과 -2 다음에는 리드를 입력해주면 됩니다.
이렇게 매핑 과정을 거친 뒤에는 정렬해줘야 합니다. 이렇게 정렬하는 이유도 마찬가지로 매핑이 끝난 뒤 리드들을 빠르게 검색할 수 있도록 하기 위함이라고 생각하심 됩니다. 정렬은 samtools sort를 이용해 할 수 있습니다. 보통 이렇게 매핑되거나 얼라인된 리드 정보를 저장하기 위해 SAM이라는 포맷을 이용합니다.
매핑된 리드 수나 리드 비율은 samtools flagstat 등으로 확인할 수 있습니다만, hisat2를 돌리면 간단한 통계값을 보여주기 때문에 그걸 그냥 살펴보셔도 좋습니다.