김준네 연구실 튜토리얼

목록

• 튜토리얼을 시작하시는 걸 환영합니다
• 대용량 시퀀싱 데이터에 대한 기본적인 이해
• 까만 바탕에 하얀 글씨-리눅스에서 작업 하기
• 행렬 데이터 다루기
• Conda를 이용한 프로그램 설치
• 기본적인 프로그램 다뤄보기
• 간단한 변이 추출 따라해보기
• 변이가 주는 영향을 이해해보기
• 유전체를 조립하기
• 구조 변이 분석하기

튜토리얼을 시작하시는 걸 환영합니다

• 이 자료는 충남대학교 생명시스템과학대학 생명정보융합학과에서 일하고 있는 김준이 만들었습니다.
• 일단은 제 연구실에 관심 있는 학생들에게 컴퓨터를 이용해 어떻게 생명현상을 이해할 수 있는지 가르쳐주기 위해 만든 자료입니다.
• 주로 롱리드 시퀀싱(long-read sequencing) 데이터과 이를 활용한 다양한 생물학 연구를 다룹니다.
• 교육용 목적으로는 어떤 식으로 이용해도 좋으니 편하게 가져다 쓰시길 바랍니다.
• 궁금한 점은 위쪽에 보이는 Issues에 남기시면 됩니다.
• 각 실습 자료는 읽는 데 1시간 이내, 실습하는 데 3-4시간 정도를 들이면 충분하도록 짰습니다만, 익숙해지는 정도에 따라 시간 차가 있을 순 있습니다.
• 중요한 건 하다 보면 금방 익숙해질 거라는 점입니다. 사람마다 시간 차가 있을지언정, 충분한 시간만 들이면 누구나 익숙해질 수 있습니다.
• 기본적으로는 코딩하는 법을 다루지 않습니다. 좋은 프로그램이 이미 정말 많이 개발돼 있기 때문에 코딩할 필요가 거의 없기 때문입니다.
• 그렇기 때문에 주로 다루는 일은 이런 프로그램을 갖다 쓰는 방법입니다. 코딩보다는 스크립팅에 가깝습니다.
• 현존하는 다양한 프로그램만 잘 활용할 줄 알아도 수많은 생물학 연구를 해낼 수 있습니다.
• 이 튜토리얼을 다 끝낸 뒤에는 본격적인 연구를 시작할 수 있는 수준일 겁니다. 그렇게 되길 기대합니다.

대용량 시퀀싱 데이터에 대한 기본적인 이해

• 본 튜토리얼에서는 대용량 시퀀싱 데이터에 대해 이해하고, 고성능 컴퓨터를 이용해 이를 처리하는 방식에 대해 배웁니다.

시퀀싱이란?

• 그 첫 번째 단계로 시퀀싱 데이터가 무엇인지, 왜 고성능 컴퓨터를 활용해야 하는지 알아보고자 합니다.
• 시퀀싱이란 DNA, RNA, 단백질 등 실처럼 생긴 선형(linear) 분자 안에 담긴 단량체(monomer)의 정보를 읽어내는 과정을 가리킵니다.
• 이들 DNA, RNA, 단백질은 생물학의 중심 원리(central dogma)에 등장하는 가장 핵심 생체 분자입니다. 그러니 이 분자들에 담긴 정보를 읽어내 사람마다 어떤 DNA/RNA/단백질이 다르기에 어떤 사람은 특정한 병에 더 잘 걸리는데 다른 사람은 그렇지 않은지 확인하고 싶은 연구자가 많습니다. 이런 표현형 차이를 만들어내는 유전형 등의 차이를 이해할 수 있다면 진단이나 약물 개발과 치료 등 다양한 시도를 해볼 수 있을 테니까요.
• 우리가 맨눈으로도 이런 선형 분자에 담긴 정보를 읽어낼 수 있다면 좋겠지만, 현실적으로는 불가능합니다. 그러니 우리 눈을 대신할 수 있는 간접적인 방식을 통해 선형 분자 안에 담긴 단량체 정보를 확인하는 일이 필요합니다. 그리고 이렇게 확인한 정보를 우리가 이해할 수 있는 문자 형태의 분자로 바꾸면 좋겠죠. 이 과정이 시퀀싱 과정인 셈입니다. 시퀀싱이란 말 그대로, 기다란 선형 분자 안에 담긴 단량체의 정보를 순서대로 읽어내는 것입니다.

DNA 시퀀싱

• DNA는 이러한 주요 선형 분자 중에서도 가장 효과적으로 시퀀싱 데이터를 대량 생산할 수 있는 주요 생체 분자입니다. 생체의 DNA는 주로 아데닌, 타이민, 구아닌, 사이토신이라는 4개 분자로 이뤄져 있습니다. 그리고 이 4개의 분자가 각각 양손을 내민 채 옆 분자를 붙들고 있는 게 DNA라는 선형 분자의 생김새입니다.
• 예컨대 사람의 DNA는 약 30억 개의 단량체로 이루어져 있는데요, 이를 30억 염기쌍(basepair; bp)이라고 부릅니다. 이들 중 거의 대부분은 23개의 염색체에 담겨있으며, 이러한 염색체는 그 자체로 나름 안정화된 구조를 이룹니다. 참고로 DNA는 2개의 실이 서로 마주 보고 뱀처럼 휘감긴 형태를 띠고 있는데, 이중나선이라는 표현은 여기에서 나온 것입니다.
• DNA를 시퀀싱한다는 것은, 이처럼 사람의 DNA에 있는 단량체 정보를 순서대로 읽어내 우리가 이해할 수 있는 문자로 변환하는 것입니다. 여기서 아데닌은 A, 타이민은 T, 구아닌은 G, 사이토신은 C로 표현합니다. 시퀀싱이 끝난 데이터는 예를 들면 ATGCTAGCTTTAGCCTAGCTACGT 같은 문자열이 됩니다.
• DNA를 시퀀싱할 수 있는 기법은 상당히 많고 지금도 계속해서 개발되고 있습니다만, 본 튜토리얼에서는 이를 크게 둘로 나누고자 합니다. 하나는 숏리드(short-read) 시퀀싱 기법이며, 다른 하나는 롱리드 시퀀싱 기법입니다. 이는 한번에 읽어내는 DNA 분자의 길이에 따라 나누는 구분입니다. 숏리드 시퀀싱이란 대개 300염기쌍 미만의 DNA 단량체 서열을 한번에 읽어내는 기법을 가리키며, 롱리드 시퀀싱이란 보통 대부분의 서열 정보가 10,000염기쌍 이상의 DNA 단량체 서열을 한번에 읽어내는 기법을 가리킵니다. 이렇게 한번에 읽어낸 DNA 서열 정보를 리드라고 부릅니다.
• DNA 시퀀싱 실험 과정을 이해함으로써 리드에 대해 조금 더 정확하게 알아봅시다.
• DNA 시퀀싱 과정은 실험마다 차이가 있긴 합니다만, 보통은 DNA를 자르는 데에서 시작합니다. DNA를 자르는 이유는 단순합니다. 한없이 긴 DNA라면, 그 안에 담긴 단량체 정보를 처음부터 끝까지 한번에 읽어낼 수 있는 기술이 없는 것이나 다름 없기 때문입니다(있긴 합니다). 그렇기 때문에 DNA를 시퀀싱하는 경우에는 대부분 DNA를 잘라서 읽기 편한 길이로 맞춰줍니다. 숏리드 시퀀싱 기법이라면 500염기쌍에서 700염기쌍 정도로, 롱리드 시퀀싱 기법이라면 20,000염기쌍이나 그보다 긴 길이 정도로 맞춰줍니다. 그리고 이 조각난 DNA에 담긴 단량체 정보를 A, T, G, C 등으로 바꿔주는 것입니다. 다시 한 번 강조하자면, 이렇게 한번에 읽은 DNA 조각을 리드라고 부릅니다.
• 리드가 길면 하나의 분자에서 훨씬 더 많은 정보를 얻을 수 있으며, 그로 인해 리드가 짧을 때는 분석할 수 없는 생명현상을 들여다볼 수 있게 됩니다. 이는 뒷부분에서 기술할 예정입니다.
• 보통 이런 DNA 조각, 즉 리드를 한두 개만 읽진 않습니다. 사람의 유전체 정보를 분석한다고 하면 전체를 적어도 한 번은 읽어야 할 테니 30억 염기쌍은 읽어야겠죠. 대충 150염기쌍을 읽어내는 숏리드 시퀀싱 기법으로는 이런 리드를 적어도 2천만 개를 읽어야 합니다. 대충 15,000염기쌍을 읽어내는 롱리드 시퀀싱 기법이라면 20만 개면 거뜬하죠. 물론 실제로는 이러한 숫자의 30배 정도를 생산하긴 합니다.

DNA 시퀀싱 데이터 양식

• DNA에 담긴 정보는 여러 형태로 저장될 수 있습니다만, 가장 표준화된 데이터 양식(format)은 크게 두 가지가 있습니다. 하나는 FASTA, 다른 하나는 FASTQ입니다. 파스타와 파스트큐, 패스트에이와 패스트큐 등등 다양하게 불립니다. 대충 아무렇게나 불러도 다 알아듣는다는 뜻입니다.
• 참고로 DNA에만 쓰이는 양식은 아니며, 선형 분자라면 어떤 것이든 표현할 수 있습니다.
• FASTA 포맷은 DNA 분자의 이름과 단량체 서열 정보를 담고 있습니다. 예를 들면 다음과 같은 형태를 띱니다.

>read_00001
ATGCTACGTACGTACGTAGCTACGATCGATCGATCGTAGTACACAATTTCGCGCCGCTCG
AAATCTCGGGCTCCCATCACGTACGAGCTAGCTAGCTAGCTACGATCGATCGATCGATGT

• FASTA 포맷에서 DNA 분자의 이름은 >로 시작합니다. 여기서는 read_00001이라는 것이 해당 DNA의 이름이 됩니다. 그리고 ATGCT......GATGT에 해당하는 부분이 이 read_00001이라는 DNA 분자에 담겨있던 단량체의 서열 정보가 됩니다.
• FASTA는 한 파일에 여러 개의 DNA 서열 정보가 담기기도 합니다. 다음과 같이 여러 줄로 DNA 정보가 표현되는 것이죠.

>read_00001
ATGCTACGTACGTACGTAGCTACGATCGATCGATCGTAGTACACAATTTCGCGCCGCTCG
AAATCTCGGGCTCCCATCACGTACGAGCTAGCTAGCTAGCTACGATCGATCGATCGATGT
>read_00002
ATGATAAATCGTACGATTTTTTTCCCCCCCGATGCTAGTCGATCGAAAACGTACGTAGCT
TAGCTAGCTTTACACACGGGGAAAGTGCAAAGTCGTACGGGCACATGTCGGACTAGCTAG
GTAGCTAGCTAACAGTCGCTGGCATGCATCATTGCTAGCTAGCCGCTAGCTAGTCGATTT

• 이때 한 DNA의 정보가 시작되고 끝나는 위치는 DNA의 이름을 알려주는 > 문자를 통해 확인할 수 있습니다. 위 예제에서처럼 한 DNA 분자는 120염기쌍 정보를 담고 있고, 다른 DNA 분자는 180염기쌍을 담고 있으면 길이가 달라서 어디서부터 어디까지가 이름인지 알기 어렵겠죠? 그렇기 때문에 길이가 아닌 > 문자가 시작하면 새로운 DNA 분자에 대한 정보가 나온다는 약속을 해둔 것입니다. DNA 길이는 엄청나게 길어질 수도 있습니다. 다음 > 문자가 나오기 전까지는 아무리 길어도 하나의 DNA 정보라고 생각하면 됩니다.
• FASTQ는 FASTA와 정말 비슷하지만, 중요한 차이점을 몇 가지 지니고 있습니다. FASTA는 DNA의 이름과 서열 정보만 담고 있지만, FASTQ는 이름과 서열 외에도 부가적인 옵션 정보와 각 염기의 정확도를 알려주는 퀄리티(quality) 정보가 함께 담겨 있다는 것이 가장 중요한 차이점입니다. 예를 들면 다음과 같은 형태를 띠고 있습니다.

@read_00001
ATGCTACGTACGTACGTAGCTACGATCGATCGATCGTAGTACACAATTTCGCGCCGCTCG
AAATCTCGGGCTCCCATCACGTACGAGCTAGCTAGCTAGCTACGATCGATCGATCGATGT
+
DBD0@GHIIHHHHIGEHFH@G1@GCGHHFHIE1GHIIIHIIHHHHHIIG@GFCF@D@DHC
HHHHIEHI1@@CGHI@GH1G1D0@CC@CE@CCFH@1@CEHHHHFCHHHCHHIHEHHEHEH

• FASTQ에서 DNA의 이름은 @ 문자 다음에 나옵니다. 그 다음 줄부터 서열 정보가 시작되며, + 기호가 나오면 DNA 서열 정보가 끝났다는 뜻이 됩니다. 그리고 마지막 줄에 퀄리티 정보가 담기게 되죠. 마찬가지로 이 품질 정보는 다음 @ 기호가 나오면 끝나게 됩니다.
• 현존하는 시퀀싱 기법은 대부분 오류를 만들어낼 가능성이 있기 때문에 이러한 퀄리티 정보가 포함된 FASTQ 양식이 널리 쓰이고 있습니다.
• 여기서 퀄리티 정보의 길이는 DNA의 서열 정보와 같은데요, 이는 퀄리티 줄의 각 문자 하나가 짝지어진 DNA 서열의 정확도를 가리키기 때문입니다.
• 예를 들면 첫 번째 DNA 서열 정보인 A, 그 다음의 T, 그 다음의 G는 각각 "D", "B", "D"라는 퀄리티를 지니고 있다는 뜻이 됩니다.
• 이 퀄리티는 기본적으로는 -10×log(Error rate)입니다.
• 예를 들면 시퀀싱 결과 첫 번째 단량체가 A라고 판정했는데, 실험적 오류로 인해 10% 정도는 T일 것 같다고 쳐봅시다. 그러면 이때 퀄리티 값은 -10×log(0.1) = 10 이 되겠죠.
• 그런데 이렇게 십진법으로 나타내면 퀄리티 값이 2개의 문자로 표현하게 될 가능성이 매우 높다는 문제가 생깁니다. 서열 정보와 짝 지어질 수 없게 되는 게 문제입니다. 예를 들면 ATGC라는 서열 1021121이라는 퀄리티 정보를 지닌다고 했을 때, 각 단량체 정보의 정확도가 10, 21, 12, 1인지, 아니면 10, 2, 11, 21인지 등등 알기가 어렵겠죠. 서열 문자와 퀄리티 정보가 모두 1개의 문자로 표현되는 게 필요한 겁니다.
• 이런 문제를 해결하기 위해 퀄리티 값을 1개의 문자만으로 표현할 수 있도록 하는 장치가 쓰이고 있습니다. 십진법으로는 두 자리인 숫자도 모두 1개의 문자로 표현할 수 있는 방식이죠.
• 예컨대 read_00001의 첫 3개 정보인 A, T, G는 각각 D, B, D라는 퀄리티 값을 지니고 있는데요, 이는 각각 35, 33, 35를 가리킵니다. 오류일 가능성은 각각 10^-3.5, 10^-3.3, 10^-3.5인 셈이니 상당히 믿을만한 정확도라고 볼 수 있겠죠.

시퀀싱 플랫폼별 차이점

• 현재 가장 널리 쓰이는 숏리드 시퀀싱 기법은 Illumina, BGI 등이 제공하고 있으며, 가장 널리 쓰이는 롱리드 시퀀싱 기법은 PacBio와 ONT 등이 제공하고 있습니다.
• 숏리드 시퀀싱은 리드 길이가 짧은 대신 정확도가 높고 값이 매우 저렴하다는 장점을 지닙니다. 인간 유전체의 30배(30×라고도 부름) 정도를 생산한다고 치면 900억 염기쌍, 다른 말로는 90 Gb (90×10⁹염기쌍)를 얻는 수준인데요, 한국 돈으로 대략 100만 원 정도 소요됩니다(업체 기준).
• 롱리드 시퀀싱은 업체별로 기술의 장점이 다릅니다만, 둘 다 비싸다는 단점은 똑같습니다.
• PacBio만이 제공하고 있는 롱리드 시퀀싱 기법은 HiFi 시퀀싱 기법이라고 불립니다. 이는 다음과 같은 특징을 지닙니다.
• HiFI 시퀀싱 기법은 정확도가 매우 높은 롱리드 시퀀싱 기법입니다. 2023년 현재 기준으로는 숏리드 시퀀싱 기법의 정확도를 뛰어넘어 Q40 (퀄리티 40, 즉 오류 가능성 10^-4 수준)을 넘는 리드를 제공해줍니다. DNA 시퀀싱의 표준이 되고 있습니다.
• HiFi 시퀀싱 기법의 리드 길이는 약 20,000염기쌍, 즉 20 kb 수준입니다. 이는 시퀀싱 화학 반응의 특징 때문인데, HiFi 시퀀싱 기법은 하나의 DNA 이중나선을 원형(circular)으로 만든 뒤 여러 차례 읽어 오류를 줄입니다. 이처럼 원형을 만든 뒤 여러 번 읽어내려면 DNA 길이가 현재로서는 20 kb 수준이어야 가장 적당합니다.
• HiFi 시퀀싱 기법은 DNA 메틸화(methylation) 등을 확인할 수 있게 해줍니다.
• HiFi 시퀀싱 기법을 이용해 30× 인간 유전체 데이터를 생산한다면, 현재 한국 돈으로 대략 500만 원 정도 소요됩니다.
• ONT에서 제공하는 시퀀싱 기법은 나노포어 단백질을 활용합니다. DNA나 RNA 같은 선형 분자가 이 나노포어 단백질을 지날 때 나타나는 전류 차이를 이용해 단량체 서열 정보를 확보할 수 있습니다.
• ONT에서 제공하는 가장 두드러진 시퀀싱 데이터 유형은 울트라 롱리드(ultra-long-read) 시퀀싱 데이터입니다. 이는 상당수의 DNA 정보가 200 kb 수준이고 가장 긴 DNA 분자는 Mb 수준으로 읽힌다는 강점을 지니고 있습니다. 유전체 지도 작성시 보조 데이터로 많이 활용됩니다.
• ONT 울트라 롱리드 시퀀싱 데이터는 생산해주는 업체가 많지 않지만, 30× 인간 유전체 데이터 기준 대략 800만 원 정도 듭니다.
• ONT는 현존하는 시퀀싱 기법 중 유일하게 DNA는 물론 RNA에 있는 변형(modification)까지도 확인할 수 있다는 강점을 지닙니다. 관련해서는 다양한 실험 기법과 분석 프로그램이 계속해서 개발되고 있습니다.
• ONT의 정확도는 그렇게 높지는 않으며, 대략 Q20 내외라고 생각합니다.
• 이외에도 Dovetail이나 Arima Genomics 사의 HiC 시퀀싱 기법(physical mapping/chromatin interaction 등 관찰 가능), Bionano 사의 옵티컬 매핑(optical mapping) 기법 등도 널리 쓰입니다.
• 이러한 시퀀싱 데이터는 대부분 컴퓨터에 저장된 채로 활용됩니다. 거의 대부분은 고성능(high-performance) 컴퓨터를 이용해 시퀀싱 데이터를 처리하게 됩니다.

왜 고성능 컴퓨터를 활용하나?

• 시퀀싱 데이터를 활용하려면 고성능 컴퓨터가 거의 필수적입니다. 데이터 자체가 엄청난 양으로 생산되다 보니, 이러한 대용량 데이터를 다루려면 일반 컴퓨터로는 매우 오랜 시간이 걸리기 때문입니다. 데이터를 빠르게 처리해서 좋은 논문으로 낸 뒤 졸업하려면 고성능 컴퓨터를 다룰 수 있어야만 합니다.
• 고성능 컴퓨터의 운영체제에는 대개 리눅스(Linux) 기반 운영체제가 쓰입니다. 참고로 윈도우즈 같은 게 운영체제입니다.
• 예컨대 제가 처음으로 썼던 서버는 우분투 16.04라는 운영체제가 설치돼 있었습니다. 여기서 "우분투"는 "윈도우즈"와 대응하는 개념이며, "우분투 16.04"는 "윈도우즈 10"과 비슷한 개념입니다. 리눅스 기반 운영체제인 우분투의 16.04 버전을 썼다는 뜻입니다.
• 이러한 리눅스 기반 서버에 익숙해지는 것이 다음 단계의 목표입니다.
• 정확하게는 커맨드라인 인터페이스(command-line interface; CLI라고도 부름)에 익숙해지는 과정이 필요합니다. 이는 많은 사람들이 익숙한 그래픽 사용자 인터페이스(graphical user interface; GUI)와 대응하는 개념입니다.
• GUI는 쉽게 말하면 화면에 보이는 아이콘 등을 마우스로 클릭하는 것을 가리킵니다.
• CLI는 화면에 아이콘 같은 그림은 아예 없고 문자만 나오는 것을 가리키며, 키보드만 가지고 모든 작업을 할 수 있도록 설계돼 있습니다.
• 대부분의 서버를 활용하려면 이러한 커맨드라인에 익숙해지는 과정이 필요합니다. 리눅스 운영체제도 GUI를 활용하는 경우가 매우 많습니다만, 일반적인 서버는 대개 커맨드라인을 기반으로 운용됩니다. 그러니 이 커맨드라인에 익숙해지는 연습을 꾸준히 반복해서 컴퓨터와 친해지도록 합시다.

까만 바탕에 하얀 글씨-리눅스에서 작업 하기

리눅스 터미널 열기

• 리눅스 운영체제에서 작업하는 법을 익히려면, 리눅스 운영체제가 깔린 컴퓨터가 있어야 합니다.
• 제 수업 듣는 학생이나 연구실 사람이라면 제가 운영하고 있는 서버에 접속하시면 됩니다. 그전에 먼저 서버 접속에 필요한 프로그램을 설치해야 합니다. 맥을 쓰신다면 터미널 열고 접속하시면 되겠습니다만, 윈도우즈를 쓰신다면 MobaXterm을 설치하시길 추천 드립니다. Home Edition 클릭한 뒤 Installer edition 또는 Portable edition 클릭해서 다운 받고 실행하시면 되겠습니다. 둘 다 차이는 없으니 아무거나 쓰면 됩니다. 이후 접속에 필요한 정보를 입력하시면 됩니다.
• 참고로 비밀번호 입력 시 키보드를 쳐도 커서가 안 움직이고 입력이 되지 않는 것처럼 나올 수도 있는데요, 화면에 아무 변화가 없다고 해도 입력되고 있는 상태입니다. 비밀번호 끝까지 친 다음 엔터 치시면 됩니다.
• 다른 분들은 본인 운영체제에 맞춰 다음과 같이 진행하시면 됩니다.
• 맥 계열 운영체제를 쓰고 있다면 바로 터미널을 열고 작업을 시작하면 됩니다.
• 윈도우즈 계열 운영체제를 쓰고 있다면 WSL을 이용해 리눅스를 사용할 수 있습니다. 다만 윈도우즈 10 또는 윈도우즈 11이어야 활용 가능합니다. 설치는 간단하게 진행할 수 있는데요, Power Shell을 관리자 권한으로 실행한 뒤, wsl --install을 입력합니다. 그 뒤 Microsoft Store에 들어가서 Ubuntu 다운로드하고 컴퓨터 재시작하면 설치 완료입니다. 그 뒤 윈도우즈키 누르고 "Ubuntu" 검색하시면 됩니다.
• 이제 검은 배경화면에 하얀 글씨가 뜬 터미널이 열렸을 겁니다. 다음과 같은 줄이 보이면 성공입니다.

어쩌구@저쩌구:~$

• 이제 본격적인 학습을 시작해봅시다.

터미널에서 작업하기 위한 기본 명령어

• 시작하기에 앞서 한 가지 유의할 점을 미리 알려드리겠습니다. 여기에 나와있는 명령어들은 마이크로소프트 워드 같은 데 저장하시면 매우 높은 확률로 명령어가 이상해질 가능성이 높습니다. 워드에 내장된 자동 변환 기능이 특정한 문자를 다른 문자로 바꿔버리기 때문에 나타나는 현상인데요, 대표적인 예로 '는 ‘ 및 ’와는 완전히 다른 문자로 인식됩니다. 마찬가지로 " 또한 “ 및 ”와는 완전히 다른 문자입니다. 반드시 주의하시기 바랍니다. 메모장 같이 자동 변환 기능이 없는 텍스트 에디터를 활용하시길 추천 드립니다.
• 터미널에서 작업한다 해도, 컴퓨터를 사용하는 기본적인 방식은 아이콘 클릭하는 것과 다를 게 없습니다. 마우스로 클릭하는 대신 키보드를 활용해 다양한 명령어(command)를 입력하는 게 다를 뿐이죠. 여전히 우리는 특정한 디렉토리(폴더)에서 작업하게 되며, 디렉토리와 디렉토리를 이동하며 다양한 파일들을 읽고 쓰게 실행하게 될 겁니다. 키보드로 입력하는 다양한 명령어를 활용해서 말이죠.
• 가장 기본적인 명령어를 입력해봅시다. 그래픽 기반으로 작동하는 경우에는 디렉토리를 더블클릭해서 들어가면 자동으로 현재 디렉토리 안에 있는 다양한 파일과 하부 디렉토리 등을 보여줄 겁니다. 그렇지만 지금 터미널에는 아무것도 나와있지 않고 본인 아이디와 위치 정도만 나와있을 겁니다.

어쩌구@저쩌구:~$

• 현재 디렉토리에 들어있는 파일이 무엇인지 확인하고 싶다면 ls를 적고 엔터를 쳐봅시다. ls란 list의 약자로, 현재 디렉토리 안에 있는 내용물의 목록을 뽑아주는 명령어입니다.

어쩌구@저쩌구:~$ ls
directory1 directory2             # 처음 접속했을 때는 홈 디렉토리에 아무런 파일이 보이지 않을 겁니다. 아무것도 안 나와도 이상 없으니 진행하시면 됩니다.
어쩌구@저쩌구:~$

• 현재 디렉토리 안에 뭐라도 들어있다면 위에 보이는 것처럼, ls 이후 엔터를 쳤을 때 그 내용물이 뭔지를 보여줄 겁니다. 그래픽 기반처럼 자동으로 내용물을 보여주진 않지만, 아주 간단하게 내용물을 볼 수 있는 셈이죠. 참고로 제 서버에서는 기본적으로 l만 쳐도 ls를 친 것과 거의 동일한 결과를 보여줍니다. (QUIZ) ls는 파일 이름만 보여주지만, ls -l을 치면 디렉토리 안에 있는 내용물들의 다른 정보들도 제공해줍니다. 직접 한번 이 명령어를 쳐보고 각 줄이 어떤 걸 뜻하는지 검색해서 이해해봅시다. 참고로 제 서버에서 ll만 쳐도 ls -l을 친 것과 거의 동일한 결과를 확인할 수 있습니다.

집중해주세요!

• 참고로, 거의 모든 명령어는 명령어 --help 또는 명령어 -h를 치면 설명서가 나오게 설정돼 있습니다. 이해가 잘 안 된다면 검색창에 검색해서 사용법을 익히도록 합시다.

• 다음으로 가장 간단한 텍스트 에디터를 사용해봅시다. 터미널에 nano test라고 쳐봅시다.

• 참고로 nano와 test는 반드시 그 사이에 빈칸이 하나 들어가야 합니다. 빈칸이 들어가야 nano가 명령어로 인식되고, 뒤의 test는 이 nano를 사용해 열고자 하는 파일 이름이라는 것이 인식됩니다. 이후에 쓰는 다양한 명령어들도 모두 마찬가지입니다. 중요하니 꼭 주의해주세요.

어쩌구@저쩌구:~$ nano test

• 그러면 빈 파일이 하나 열리게 될 겁니다. 가운데 위쪽을 보면 "test"라는 글자가 보일 텐데, 이게 파일 이름을 가리킵니다. 이 상태에서 원한다면 아무 문자나 입력하시면 됩니다. 그러면 메모장에 글을 쓰듯 test라는 파일에 내용이 채워지죠. 이 파일에 Hello world라는 문자열을 입력해봅시다.
• 이 상태에서 나가고 싶다면 아래쪽에 쭉 적혀 있는 내용들을 읽으시면 되는데요, 자세히 보시면 ^X Exit라고 적힌 게 보일 겁니다. 여기서 ^ 문자는 키보드의 컨트롤 키를 가리킵니다. 여기 적힌 것처럼 파일을 수정 중인 메모장에서 나가고(exit) 싶다면 ctrl+x를 치시면 됩니다. 키보드 왼쪽 아래 컨트롤 키와 x 키를 누르면 나가지게 될 겁니다.
• 다만 우리처럼 뭐라도 내용을 채웠다면 저장할 건지 등을 물어봅니다. Save modified buffer?라는 말이 보일 텐데요, 이때 왼쪽 아래에 있던 명령어들이 Yes, No, Cancel 등으로 바뀌어있는 걸 보실 수 있을 겁니다. 우리는 저장해야 하니, Yes 앞에 적힌 대로 Y를 누릅시다. 그러고 나면 File Name to write라는 창으로 바뀌게 될 겁니다. 다른 이름으로 저장하고 싶다면 이름을 수정하면 됩니다. 바꾸지 않고 그냥 엔터를 치면 test라는 파일에 해당 내용이 저장될 겁니다. 지금 실습에서는 파일 이름이 test여야만 계속 진행하실 수 있으니, 일단 그대로 저장해봅시다.
• 여기서 하나 중요한 게 있는데요, 리눅스 운영체제에서는 기본적으로 대부분의 파일 확장자가 중요한 역할을 하지 않습니다. test.txt처럼 이 파일이 텍스트 파일이라는 걸 지정해주지 않더라도 상관이 없다는 뜻입니다. test라는 이름만으로도 파일은 텍스트 파일 역할을 하게 됩니다.
• 그리고 생각보다 이런 텍스트 에디터에서 다시 터미널로 나가지 못해 갇힌 사람들이 많습니다. 혹시 나중에 나가는 방법을 까먹었다면 인터넷 검색창에 "how to exit nano" 등을 검색하시면 됩니다. 실제로 다른 텍스트 에디터 중 하나인 Vim도 탈출하기가 쉽지 않다 보니, How do I exit Vim?과 같은 질문은 해당 사이트에서 5천 개 이상의 추천을 받았을 정도입니다.
• 중요한 건 잘 검색하면 뭐든 다 나온다는 겁니다. 여러분은 이 세계에 처음 발을 들인 사람도 아니고, 여러분보다 먼저 당하고 헤맨 사람은 세상에 정말 많습니다.

신경 써주세요

nano를 끄려고 하다가 ctrl + z를 누르는 일은 삼가도록 합시다. ctrl + z는 프로그램을 끄는 게 아니라, 잠시 멈추는 명령어에 가깝습니다. 이 때문에 화면에 보이지는 않지만 프로그램은 계속 대기를 하게 되는데요, 이게 nano와 같은 텍스트 에디터에 적용되면 안 좋은 결과를 일으킬 수 있습니다. 왜냐면 ctrl + z를 누른 다음에 같은 파일을 다시 nano로 여는 일은, 윈도우즈에서 똑같은 파일을 메모장으로 여러 번 여는 일과 똑같기 때문입니다. 같은 파일을 메모장으로 여러 번 열고 각자 수정하면 언제 어떻게 수정됐는지 추적도 안 되고 이상하게 되겠죠? ```nano``도 이런 식으로 열면 에러를 내보냅니다. 주의하도록 합시다.

• nano라는 텍스트 에디터 말고 다른 명령어를 활용해 파일 내용을 확인하는 법도 배워봅시다. 두 번째로 써볼 명령어는 less입니다. 터미널에 다음과 같이 입력해봅시다.

어쩌구@저쩌구:~$ less test

• 이번에는 파일 내용이 화면에 보이게 될 겁니다. 제대로 입력했다면 Hello world라는 말이 적혀있겠죠. 마찬가지로 나가봅시다. less에서 나가려면 q를 치시면 됩니다. 키보드 왼쪽 위에 있는 "Q" 키입니다.
• nano와 less는 용도가 다릅니다. nano는 텍스트 에디터, 즉 파일 내용을 바꾸고 싶을 때 주로 활용할 수 있습니다. 보통 큰 문제 없이 쓸 수 있습니다만, 파일 크기가 크면 바로 랙 걸린다는 문제를 지니고 있습니다. 파일 크기가 1 Mb를 넘어간다면 nano를 쓰는 순간 바로 랙 걸려서 모든 작업이 멈추게 됩니다. 유의하시면 좋습니다.
• 반면 less는 파일 크기가 커도 문제 없이 잘 작동한다는 강력한 장점을 지니고 있습니다. 대신 파일 내용을 수정하기에는 적합하지 않습니다. 읽기 전용이라고 생각하고 쓰셔도 무방합니다.

주목해주세요!

참고로 이전 명령어를 그대로 쓰고 싶다면 방향키를 활용하시면 좋습니다. 방금 less test라고 쳤던 걸 다시 치려면 할 수야 있지만 좀 귀찮잖아요? 방금 자기가 썼던 내용이 길면 길수록 더 그럴 겁니다. 그럴 때는 방향키 중 위쪽 방향키를 눌러봅시다. 그러면 방금 자기가 친 명령어가 고스란히 나올 겁니다. 계속 위로 올라가면 이전에 쳤던 것들이 쭉 나올 거고요, 다시 아래쪽 방향키를 내리면 보다 최근에 친 명령어로 이동할 수도 있습니다.

• 마지막으로 써볼 명령어는 cat입니다. cat은 "concatenate"의 약자로, 파일의 내용을 연달아 보여주는 명령어입니다. 1개의 파일에 대해서도 쓸 수 있으며, 여러 개의 파일에 대해서도 사용할 수 있습니다. 마찬가지로 터미널에 다음 명령어를 입력해봅시다.

어쩌구@저쩌구:~$ cat test

• 그러면 이전에 nano나 less를 썼을 때와는 달리 파일 내용이 화면에 바로 출력되는 걸 확인할 수 있을 겁니다.

어쩌구@저쩌구:~$ cat test
Hello world
어쩌구@저쩌구:~$

• 이처럼 cat은 파일의 내용을 수정하거나 한 페이지씩 차근차근 보여주는 것이 아니라, 화면에 고스란히 출력한다는 특징을 지니고 있습니다.
• 참고로 이렇게 화면에 출력된 결과를 STDOUT, 표준 출력이라고 부릅니다. 이건 꽤 중요한 개념이니 기억해두시면 좋겠습니다.
• 다음으로 파일 내용이 아니라, 우리가 원하는 문자를 곧바로 표준 출력해주는 명령어인 echo를 사용해봅시다. 다음 명령어를 쳐보시기 바랍니다.

어쩌구@저쩌구:~$ echo "Hello world"

• 그러면 이전에 cat test를 입력했을 때와 마찬가지로 화면에 Hello world라는 문자열이 고스란히 표준 출력으로 나오게 됐을 겁니다.
• 이번에는 이렇게 나온 표준 출력을 새로운 파일에 저장해보도록 하겠습니다. 다음과 같이 명령어를 입력해봅시다.

어쩌구@저쩌구:~$ echo "Hello world!" > test2

• 이전과 달리, 이번에는 화면에 Hello world!라는 문자열이 출력되지 않았을 겁니다. 그 이유는 우리가 >라는 연산자를 활용했기 때문입니다. 이 >는 리디렉션(redirection)이라고 부릅니다.
• 리디렉션을 활용하면, 표준 출력되어야 할 Hello world!라는 문자열이 test2라는 새로운 파일에 저장되게 됩니다. 다시 말하자면, 출력되어야 할 값이 표준 입력으로 변환돼 파일에 저장되는 것입니다.
• 리디렉션이 어떻게 작동했는지 확인해보고 싶다면 다음 명령어를 입력해봅시다.

어쩌구@저쩌구:~$ cat test2

• 이렇게 하면 test2라는 파일에 담긴 내용물이 화면에 출력되겠죠? 제대로 작동했다면, 그 내용인 Hello world!가 화면에 나올 겁니다. 리디렉션을 활용해 화면에 출력될 내용을 test2라는 파일에 저장해뒀기 때문에 일어난 일입니다.
• 마지막으로 cat을 활용해 두 파일의 내용을 동시에, 연달아 확인해봅시다. 다음 명령어를 차례대로 쳐보시기 바랍니다.

어쩌구@저쩌구:~$ cat test test2
어쩌구@저쩌구:~$ cat test2 test

• 그러면 두 파일의 내용이 한번에 화면에 출력될 겁니다. (QUIZ) 위 명령어와 아래 명령어는 다른 결과를 보일 겁니다. 이유는 한번 생각해보시기 바랍니다.
• 파일을 만드는 법을 배워봤으니, 이번에는 디렉토리를 만드는 법을 배워봅시다. 윈도우즈에서 "새로 만들기 - 폴더"를 실행하는 것과 동일합니다. 여기서는 마우스 대신, 키보드를 활용한다는 것만 다르죠. mkdir을 사용해봅시다. mkdir은 "make directory"라는 뜻입니다. 다음 명령어를 입력해주세요.

어쩌구@저쩌구:~$ mkdir sampleDirectory

• 엔터 치면 아무 변화도 일어나지 않는 것처럼 보일 겁니다. 하지만 실제로는 무언가 변화가 일어났습니다. 그 차이를 확인하고 싶다면 현재 디렉토리에 들어있는 내용물을 확인해보면 되겠죠? ls를 쳐봅시다.

어쩌구@저쩌구:~$ ls
directory1 directory2 test test2 sampleDirectory

• 그러면 이전과 달리, 우리가 새로 만든 파일인 test, test2와 함께 sampleDirectory라는 디렉토리가 새롭게 형성돼있는 것을 확인할 수 있을 겁니다.
• 이제 새로 만든 이 디렉토리로 이동해봅시다. 이동에 쓰이는 명령어는 cd로, "change directory"라는 뜻입니다. 다음 명령어를 입력해주세요.

어쩌구@저쩌구:~$ cd sampleDirectory

(QUIZ) 참고로 cd s까지만 치고 키보드에 있는 Tab 키를 눌러보시면 좋습니다. 이름을 다 적지 않은 상태에서 Tab을 치면 어떤 일이 벌어지는지는 직접 확인해보세요. 또 ls t까지만 친 다음에도 Tab을 쳐보시길 바랍니다. 이때는 Tab을 한번만 치는 게 아니라 두세 번 정도 쳐보시면 좋습니다. 어떤 차이가 있는지 한번 생각해보시고, 인터넷에 검색해서 확인해봅시다.

• 위에서처럼 이동하고 나면, 터미널이 뭔가 바뀌어있는 걸 확인할 수 있을 겁니다.

어쩌구@저쩌구:~$ cd sampleDirectory
어쩌구@저쩌구:~/sampleDirectory$

• 명령어를 입력하던 부분이 ~$에서 ~/sampleDirectory$로 바뀌어있을 겁니다. 디렉토리를 이동했으니, 지금 어디에 있는지 그 위치를 표시해주는 것입니다. (주의 사항) 만약에 cd sampleDirectory를 아무리 쳐도 계속 그런 파일이 없다고 나온다면 다음과 같은 내용을 확인해주세요. (1) mkdir sampleDirectory를 입력할 때 오타를 낸 건 아닌지 확인해주세요. 오타가 있다면 오타를 수정해서 다시 디렉토리를 만드시면 됩니다. (2) 현재 위치가 어쩌구@저쩌구:~$이 맞는지 확인해주세요. 만약 아니라면 터미널에 cd만 입력하시면 됩니다. 그러면 어쩌구@저쩌구:~$로 돌아와있을 겁니다. 다시 ls를 쳐서 sampleDirectory가 있는지 확인해주시면 됩니다.

• 다시 원래 있던 디렉토리로 이동해봅시다. 다음 명령어를 입력하시면 됩니다.

어쩌구@저쩌구:~/sampleDirectory$ cd ..

• 마찬가지로 cd는 디렉토리를 이동하라는 명령어이며, ..은 상위 디렉토리라는 뜻입니다. 다시 말해 cd ..이란, sampleDirectory라는 하위 디렉토리에 들어오기 전에 작업하고 있던 곳으로 돌아가라는 뜻이 되죠.

어쩌구@저쩌구:~/sampleDirectory$ cd ..
어쩌구@저쩌구:~$

• 그러면 위에 보이듯 이전 디렉토리인 ~로 돌아오게 될 겁니다.
• 그러면 sampleDirectory에 들어오기 전까지 작업하고 있던 곳, ~, 이 물결 표시는 대체 어느 위치였던 걸까요? 이런 걸 확인하고 싶다면, pwd 명령어를 활용하시면 됩니다. 이는 "print working directory"를 뜻하는 명령어로, 현재 디렉토리의 위치 정보를 보여줍니다.

어쩌구@저쩌구:~$ pwd
/home/어쩌구

• 이렇게 pwd를 치면, 화면에 현재 위치 정보가 출력될 겁니다. 보통은 /home/본인ID라고 나올 겁니다. ~란 /home/본인ID/를 줄여둔 표현이었던 것입니다.
• 이제 이러한 위치 개념을 조금 더 깊게 배워보겠습니다. 바로 상대 경로(relative path)와 절대 경로(absolute path)라는 개념입니다.
• 상대 경로는 내가 현재 작업하고 있는 위치가 기준이 됩니다. 현 위치로부터 하위 디렉토리를 검색하거나, 상위 디렉토리로 이동하는 식으로 작동하죠. 예를 들면 이전까지 했던 cd sampleDirectory와 cd .. 같은 명령어는 모두 상대 경로를 기반으로 작동하는 것입니다.
• 상대 경로가 작동하는 방식을 알아보고 싶다면 이게 제대로 작동하지 않는 경우를 보면 좋겠죠. 홈 디렉토리에서 다시 한 번 상위 디렉토리로 이동해봅시다.

어쩌구@저쩌구:~$ cd ..
어쩌구@저쩌구:/home$

• 그 뒤 우리가 만들었던 sampleDirectory로 이동해봅시다.

어쩌구@저쩌구:/home$ cd sampleDirectory
-bash: cd: sampleDirectory: No such file or directory

• 그러면 위에 보이듯, 그런 파일이나 디렉토리가 없다는 문구가 출력될 겁니다. 왜냐면 이 sampleDirectory는 현재 디렉토리에서는 검색할 수 없는 곳에 있기 때문입니다. 현재 디렉토리에 있는 내용물을 다시 한 번 확인해봅시다.

어쩌구@저쩌구:/home$ ls
어쩌구

• 그러면 위와 같이 본인 아이디가 나올 겁니다. 현재 디렉토리에서는 sampleDirectory가 검색되지 않습니다. sampleDirectory는 /home/어쩌구 디렉토리의 하위 디렉토리이지, /home 디렉토리에서 바로 검색 가능한 하위 디렉토리는 아니기 때문입니다.
• 이처럼 상대 경로를 이용하려면, 현재 위치에서 내가 이동하려는 위치 또는 사용하려는 파일이 있는 위치를 정확하게 지시해줘야 합니다. 현재 상황에서 sampleDirectory로 곧장 이동하고 싶다면, 내 홈 위치를 거쳐서 가야만 하는 겁니다. 예를 들면 1번, 2번, 3번 방이 순서대로 있다고 쳤을 때, 1번 방에서 3번 방으로 이동하려면 2번 방을 반드시 거쳐야 하는 것처럼 말이죠.
• 구체적으로는 다음과 같이 이동하면 됩니다.

어쩌구@저쩌구:/home$ cd 어쩌구/sampleDirectory
어쩌구@저쩌구:~/sampleDirectory$

• 여기서 어쩌구 대신에 ls를 쳤을 때 나온 본인 디렉토리 이름을 넣어주시면, 위에 보이듯 한번에 이동하게 될 겁니다. 이번에는 내가 이동하려던 sampleDirectory가 현재 디렉토리에서 바로 접근하는 것이 아닌, 어쩌구라는 디렉토리를 거친 뒤에 갈 수 있는 곳이라고 컴퓨터에 알려줬기 때문입니다. 이때 디렉토리 이름과 디렉토리 이름 사이에는 반드시 /가 들어가야 합니다. 그래야 컴퓨터가 어쩌구sampleDirectory라는 이름을 가진 바로 하위에 있는 디렉토리를 검색하지 않고, 어쩌구라는 하위 디렉토리 안에 있는 sampleDirectory를 찾아서 들어가줍니다.
• 마찬가지로, 상위 폴더의 상위 폴더로 이동하고 싶을 때도 여러 방식을 쓸 수 있습니다.

어쩌구@저쩌구:~/sampleDirectory$ cd ..
어쩌구@저쩌구:~$ cd ..
어쩌구@저쩌구:/home$

• 또는, 한번에 이동할 수도 있죠.

어쩌구@저쩌구:/home$ cd 어쩌구/sampleDirectory   # 안쪽 디렉토리로 이동한 뒤,
어쩌구@저쩌구:~/sampleDirectory$ cd ../..       # 상위 디렉토리의 상위 디렉토리로 한번에 이동
어쩌구@저쩌구:/home$

• 이런 상대경로 방식을 활용해서, 내가 작업하고 있는 디렉토리에서 자유롭게 이동할 수 있습니다.
• 물론 상대경로 방식이 항상 유용한 것은 아닙니다. 예를 들면 내가 디렉토리의 너무 깊숙한 곳까지 이동했다면 어떻게 해야 할까요? 다음과 같이 명령어를 입력해봅시다.

어쩌구@저쩌구:/home$ cd 어쩌구/sampleDirectory                                 # 안쪽 디렉토리로 이동한 뒤,
어쩌구@저쩌구:~/sampleDirectory$ mkdir -p test/test/test/test/test/test/test/  # 하위 폴더 아래에 계속해서 하위 폴더 만들기
어쩌구@저쩌구:~/sampleDirectory$ cd test/test/test/test/test/test/test/
어쩌구@저쩌구:~/sampleDirectory/test/test/test/test/test/test/test$

• 여기서 다시 홈 디렉토리(~)까지 이동하려면 ../를 한참이나 입력해야겠죠. 아주 귀찮은 일입니다. 이렇게 한참이나 멀리 떨어진 경로를 이동해야 할 때, 절대 경로가 아주 유용하게 쓰일 수 있습니다.
• 컴퓨터에 있는 모든 디렉토리와 파일은 저마다 정해진 위치 정보가 있습니다. 예를 들면 우리가 위에서 pwd를 쳤을 때 확인했던 /home/어쩌구 같은 게 그 예죠. 위에서 /는 디렉토리와 디렉토리를 나눠주는 구분 기호라고 말씀 드렸는데요, 여기서 보시면 절대 경로는 항상 가정 첫 /에서부터 시작한다는 걸 알 수 있습니다. 바꿔 말하면 모든 파일과 디렉토리는 가장 첫 번째 /의 하위에 위치하는 것입니다. 이렇게 맨 처음 등장하는 /를 루트(root)라고 부릅니다. 마치 뿌리에서 뻗어나오는 가지처럼, 컴퓨터 안의 모든 파일과 디렉토리는 이 루트에서부터 뻗어나간 뒤 그로부터 제 위치를 배정 받는 것입니다.
• 이번에는 이 절대 경로를 이용해서 홈 디렉토리로 돌아가는 방법을 익혀봅시다. 홈 디렉토리로 한번에 이동하고 싶다면 다음과 같이 입력하면 됩니다.

어쩌구@저쩌구:~/sampleDirectory/test/test/test/test/test/test/test$ cd /home/어쩌구
어쩌구@저쩌구:~$

• 여기서 지정해준 /home/어쩌구(또는 본인 ID)가 바로 내 홈 디렉토리의 절대 경로가 됩니다. 그러니 그 깊숙한 곳에서 뿌리에 더 가까운 홈까지 한번에 이동할 수 있었던 것이죠.
• 마찬가지로 홈 디렉토리에서 깊숙하게 있는 안쪽 디렉토리로도 한번에 이동할 수 있습니다. cd를 입력한 뒤 이동하고자 하는 위치의 절대 경로만 입력하면 한번에 이동하게 됩니다.

어쩌구@저쩌구:~$ cd /home/어쩌구/sampleDirectory/test/test/test/test/test/test/test
어쩌구@저쩌구:~/sampleDirectory/test/test/test/test/test/test/test$ 

• 그러면 위 예시처럼 위치가 한번에 바뀌게 되죠. cd를 친 다음 단계별로 하나씩 칠 필요가 없는 것입니다.
• 참고로, 홈 디렉토리는 별 다른 위치 정보를 주지 않고도 한번에 이동할 수 있습니다. 아무런 위치 정보를 주지 않고 cd를 치기만 하면 됩니다.

어쩌구@저쩌구:~/sampleDirectory/test/test/test/test/test/test/test$ cd
어쩌구@저쩌구:~$

커맨드라인 작업에 익숙해지기 위한 보물 찾기

• 이제 가장 기본적인 위치 이동과 파일 내용 확인 관련 명령어에 대해서는 다 알려드렸습니다.
• 하지만 어떤 것이든, 아는 것보다 중요한 건 익숙하게, 능숙하게 사용할 수 있느냐 하는 점이겠죠. 이런 커맨드라인에 익숙해질 수 있도록 보물 찾기 게임을 하나 준비했습니다.
• 이런 게임을 준비한 이유는 그냥 하면 지겨워질 수 있기 때문입니다. 그러니 조금이나마 즐기면서 익숙해질 수 있도록 게임을 준비한 거죠. 마치 예전에 마우스가 보편적이지 않던 시절 지뢰찾기 게임을 통해 마우스 사용법에 익숙해질 수 있었던 것처럼, 커맨드라인 연습용 보물 찾기 게임을 진행하면서 커맨드라인과 친해지시길 바라겠습니다.
• 참고로 이 게임은 제가 만든 건 아니고요, 미국 해양연구소(Marine Biology Labs, 메사추세츠 우즈홀)에서 진행된 실습 과정에서 작성된 게임입니다. 정확히는 캘리포니아 대학 샌프란시스코에 계신 Stephan J Sanders 교수님께서 유닉스 강의를 하기 위해 고안하셨다고 하네요. 당시에 조교로 일하다가 지금은 고려대학교에서 일하고 계신 안준용 교수님께서 한국으로 들여와주셨습니다. 다만 영어가 너무 복잡하고 어렵다 보니 한국어로 제가 번역했고, 조금씩 수정을 가했습니다.
• 게임을 설치하고 싶다면 다음 명령어를 복사하고 터미널에 붙여넣기 하신 뒤 엔터 치시면 됩니다.

mkdir 01_treasureHunt && cd 01_treasureHunt
wget https://raw.githubusercontent.com/JunKimCNU/JunKimLabTutorial/main/task01_linux_tutorial/treasureHunt_v2_kor.pl
perl treasureHunt_v2_kor.pl

• 그 뒤 ls를 쳐보면 비어있던 디렉토리에 새로운 파일과 새로운 디렉토리가 잔뜩 만들어진 걸 확인할 수 있을 겁니다.
• 참고로 파일 끝에 붙은 .pl은 "Perl"이라는 프로그래밍 언어로 짜인 프로그램이라는 걸 가리키는 것입니다. 이런 파일은 위에서 perl treasureHunt_v2_kor.pl이라고 한 것처럼, perl 파일.pl 이런 식으로 실행할 수 있습니다. 보물 찾기 중간에 이런 문제가 하나 있으니 기억해두세요.
• 이제 보물 찾기를 시작해봅시다. 보물을 찾아내려면 보물 지도나 단서가 있어야겠죠? 이 게임에서는 "Clue"로 시작하는 파일에 그 단서가 저장돼 있습니다. 첫 번째 단서를 열어봅시다.

cat Clue01_S.txt

• 앞쪽에 헛소리(...)가 잔뜩 나오고 뒤에 실제로 힌트가 적혀 있습니다. 첫 번째 힌트를 예로 들자면 물의 탑(water tower)으로 향하라고 적혀 있는데요, 이 말은 디렉토리 중 waterTower 디렉토리로 이동하라는 뜻입니다. 힌트에 적힌 "서쪽으로 세 걸음, 남쪽으로 다섯 걸음"이란 말은, waterTower 디렉토리 하위에 있는 디렉토리 중 westThreeSteps로 이동한 뒤, 다시 southFiveSteps로 이동하라는 뜻입니다. 힌트가 대부분 이런 식이나 잘 따라해보시길 바랍니다.
• 중간중간 다양한 명령어를 활용할 수 있도록 예제가 포함돼 있습니다. 반드시 다양한 명령어를 검색하시면서 어떻게 쓰는 것인지 살펴보시길 추천 드립니다. 알려드리지 않은 명령어도 많으니 직접 찾아보셔야 할 겁니다.
• "보물 찾기 종료. 수고하셨습니다!"라는 말이 나와야 끝난 겁니다.
• 보물은 실재합니다. 충남대학교 오시면 제가 커피 한 잔 대접하겠습니다.
• 보물 찾기가 끝나고 나면 첫 번째 과제 링크로 들어가셔서 거기 적혀있는 다양한 명령어도 찾아보고 써먹어봅시다.
• 보물 찾기 즐겁게 하시길 바라겠습니다!

행렬 데이터 다루기

• 보물 찾기 반복 연습 많이 하셨나요? 두세 번쯤 반복하셨다면 충분히 익숙해졌을 겁니다. 이제 커맨드라인에서 쓸 수 있는 다양한 명령어를 좀 더 본격적으로 활용해봅시다.
• 이전 단계인 보물 찾기의 목표는 cd를 비롯한 디렉토리 간 이동에 익숙해지는 것이었습니다.
• 이번 단계의 목표는 awk, grep, sed, sort 등 데이터를 다루는 데 필요한 다양한 명령어에 좀 더 익숙해지는 것입니다. 실제 유전체 및 생물학 데이터를 다뤄보면서 조금 더 익숙해져보도록 합시다.
• 데이터 분석을 시작하려면 예제 파일을 다운 받아야 합니다. 다음 명령어를 실행해주시기 바랍니다.

mkdir 02_parsing_matrix
cd 02_parsing_matrix
wget https://raw.githubusercontent.com/JunKimCNU/JunKimLabTutorial/main/task02_parsing/Assemblytics_structural_variants.bed

• 다운로드가 잘 됐다면 화면에 ‘Assemblytics_structural_variants.bed’ saved [363093/363093] 같은 문구가 뜰 겁니다.
• 이제 이 파일을 한번 들여다봅시다. 저는 보통 파일을 확인할 때는 ls -l을 이용해 파일 크기를 체크하고, 너무 크면 아예 nano는 쓸 생각을 하지 않습니다. 그 다음에는 head로 먼저 파일의 첫 10개 줄을 화면에 출력해 살펴봅니다. 더 많은 데이터를 보고 싶다면 less를 이용해 파일을 열어보고요. 순서대로 한번 쭉 따라하면서, Assemblytics_structural_variants.bed 파일을 들여다봅시다. cat, less, head 등으로 파일 내용을 간략하게 살펴보고 다음 단계 진행하시면 됩니다.
• 살펴보셨나요? 들여다보면 알겠지만, 이 파일은 행렬(matrix) 형태를 이루고 있습니다. 각 열(column)과 열 사이를 나누는 기호는 탭(tab)입니다. 이처럼 탭으로 나뉜 파일 형식을 TSV (tab-separated values)라고 부릅니다.
• 이 파일은 BED라는 생물학에서 굉장히 널리 쓰이는 포맷을 따르고 있기도 합니다. 이 BED 포맷은 다양한 정보를 담을 수 있는데요, 변이(variant)에 대한 정보를 담아낼 때 아주 유용하게 쓰입니다. 기본적으로는 첫 번째 열에는 염색체 이름이, 두 번째 열에는 그 염색체에서 시작되는 위치가, 세 번째 열에는 그 염색체에서 끝나는 위치가, 그리고 나머지 열에는 다양한 정보가 기록되게 되어있고요, 이를 이용해 다양한 변이에 대한 정보를 효율적으로 담아낼 수 있습니다.
• DNA가 잘려나가거나 끼어들어간 결손과 삽입 돌연변이 같은 경우를 예로 들어봅시다. 염색체 1번의 300 bp 지점부터 350 bp 지점까지 DNA가 잘렸고 염색체 2번은 500 bp 지점에 100 bp의 DNA가 추가로 끼어들어갔다고 쳐봅시다. 이는 예를 들면 다음과 같은 BED 포맷으로 저장될 수 있습니다. Assemblytics_structural_variants.bed 파일에 담긴 내용보다는 훨씬 단순하지만, 이것도 엄연한 BED 포맷이라고 볼 수 있습니다.

#chromosome	start	end	type	size
chr1	300	350	deletion	50
chr2	500	500	insertion	100

• 이 행렬은 첫 번째 열에 염색체 정보, 두 번째 열에 변이가 시작하는 위치, 세 번째 열에 변이가 끝나는 위치, 네 번째 열에 변이의 유형, 다섯 번째 열에 그 변이의 크기가 담겨있는 데이터가 됩니다. 그러면 변이에 대한 다양한 정보를 담은 자료가 만들어지는 셈이죠. 이처럼 BED 포맷은 다양한 생물학 분석에 쓰이기 때문에, 간단하게나마 알아두는 것이 좋습니다.
• 다시 Assemblytics_structural_variants.bed 파일로 돌아가봅시다. 이 파일은 구조 변이(structural variant)라고 부르는 특정한 형태의 변이 정보를 담고 있는 파일입니다. 앞에서 언급한 삽입, 결손과 같은 돌연변이를 포함해 DNA가 크게 바뀐 변이를 가리킵니다. 현재는 50 bp 이상 바뀐 것을 구조 변이라고 부르고 있습니다.
• 또 구조 변이란 2개 이상의 유전체를 비교해서 얻어내는 정보입니다. 보통은 참조 유전체(reference genome)에 비교 대상이 되는 유전체(query genome)를 비교함으로써, 참조 유전체 기준으로 염색체의 어떤 지역에 삽입/결손과 같은 돌연변이가 존재하는지를 확인할 수 있게 해줍니다.
• 이제부터는 이 구조 변이 정보를 담은 BED 파일을 자세하게 들여다보는 연습을 해보겠습니다.
• 먼저 파일 형식에 대해 이해해봐야겠죠? 이 파일의 맨 첫 번째 줄에는 위 예제와 마찬가지로 각 열의 정보가 담겨있습니다. 이처럼 값이 아닌 정보를 담고 있는 행(row)을 헤더(header)라고 부릅니다. 헤더를 뽑는 명령어는 다음과 같이 쓸 수 있습니다.

어쩌구@저쩌구:~/02_parsing_matrix$ head -1 Assemblytics_structural_variants.bed
reference	ref_start	ref_stop	ID	size	strand	type	ref_gap_size	query_gap_size	query_coordinates	method

• 각 구조 변이에 대한 정보는 각 행에 담기게 되는데, 한 행의 각 열에는 다음과 같은 정보가 담기게 됩니다. 하나의 구조 변이에 대해, 첫 번째 열에는 참조 유전체 기준 염색체 정보, 두 번째 열에는 그 염색체에서 구조 변이가 시작되는 위치, 세 번째 열에는 그 구조 변이가 끝나는 위치 등이 기본적으로 담기죠. 그리고 다섯 번째 열에는 그 구조 변이에 부여한 이름(ID)이, 여섯 번째 열에는 그 구조 변이의 크기(size)가 담기며, 방향(strand)과 유형(type) 같은 정보도 담기게 됩니다.
• 이런 정보를 이용해 수많은 구조 변이 중 우리가 원하는 것만 쏙쏙 뽑아 들여다볼 수 있습니다. (QUIZ) 먼저 이 파일에 몇 개의 구조 변이 정보가 담겨있는지를 살펴봅시다. 이전 과제를 잘 해결했다면, 제가 넣어둔 명령어 목록 중에서 wc를 잘 활용하면 이 정보를 알아낼 수 있다는 걸 이해하실 수 있을 겁니다.
• 정답은 3,398개입니다. 헤더는 구조 변이 정보를 담고 있지 않으니 1개 빼야 한다는 걸 신경 써주세요.
• 이 파일에 담겨 있는 구조 변이의 유형(type)은 뭐가 있는지도 알아보면 좋겠죠? 일곱 번째 열에 그 정보가 담겨 있을 겁니다. 그러면 일곱 번째 열을 뽑는 명령어를 활용해봅시다. 다양한 명령어가 있습니다만, 저는 주로 awk를 활용하고 있습니다만, cut을 이용해도 좋습니다. 3천 줄이나 화면에 출력되면 귀찮을 테니, head를 이용해 앞에서 10개만 뽑아봅시다.

어쩌구@저쩌구:~/02_parsing_matrix$ awk '{print $7}' Assemblytics_structural_variants.bed | head
어쩌구@저쩌구:~/02_parsing_matrix$ cut -f 7 Assemblytics_structural_variants.bed | head

• 위 두 명령어는 동일한 결과를 줄 겁니다. "type", "Insertion", "Deletion"이 마구 뒤섞인 결과를 말이죠. 이 열에 뒤섞여 있는 정보를 정렬하고, 중복되는 정보를 제거하면 전체 구조 변이 유형이 몇 종류나 있는지 알기 좋을 겁니다. 이럴 때는 sort와 uniq를 이용하면 좋습니다. 이때는 head는 빼시면 되겠습니다.

어쩌구@저쩌구:~/02_parsing_matrix$ awk '{print $7}' Assemblytics_structural_variants.bed | sort | uniq
어쩌구@저쩌구:~/02_parsing_matrix$ cut -f 7 Assemblytics_structural_variants.bed | sort | uniq

• 그러면 이 BED 파일에는 총 6가지 유형의 구조 변이가 있다는 걸 알 수 있을 겁니다. Deletion, Insertion, Repeat_contraction, Repeat_expansion, Tandem_contraction, Tandem_expansion이 그것이죠. 이런 식으로 데이터를 확인해볼 수 있는 겁니다.
• 이전 단계의 과제를 열심히 하셨다면, 제가 여기서 열심히 쓴 | 요 기호의 역할에 대해 제대로 이해하고 계실 겁니다. 이 기호는 파이프(pipe)라고 부르는데요, 파이프는 리디렉션과 비슷하게 작동합니다. 다만 차이가 있긴 한데요, 리디렉션을 활용하면 이전 명령어의 표준 출력을 새로운 파일로 저장할 수 있었죠? 파이프를 이용하면 이전 명령어의 표준 출력을 이 다음 명령어의 표준 입력으로 바꿔주게 됩니다.
• 예를 들어 봅시다. cut -f 7 Assemblytics_structural_variants.bed | sort | uniq이라는 명령어는 파이프로 두 번 연결돼 있습니다. 이 명령어는 우리가 7번째 열(column)을 추출한 뒤(cut -f 7 부분), 그 열에 있는 데이터를 정렬하고(sort 부분), 이중 중복을 제거하는 것(uniq)을 목표로 할 때 쓸 수 있습니다. 물론 이를 다음과 같이 순서대로 진행할 수도 있습니다.

어쩌구@저쩌구:~/02_parsing_matrix$ cut -f 7 Assemblytics_structural_variants.bed > tmp1 # 7번 열에 담긴 정보를 tmp1 파일에 저장함
어쩌구@저쩌구:~/02_parsing_matrix$ sort tmp1 > tmp2                                     # 7번 열에 담긴 정보를 정렬해서 tmp2 파일에 저장함
어쩌구@저쩌구:~/02_parsing_matrix$ uniq tmp2                                            # 7번 열에 담긴 정보를 정렬한 뒤 중복 결과를 제거함

• 그렇지만 이렇게 파일에 저장하고 다시 읽고, 새로운 명령어로 처리한 뒤 다시 저장하고 다시 또 읽는 과정은 생각보다 귀찮고 시간도 많이 걸립니다. 이럴 때 파이프를 이용해 3개의 명령어를 하나로 묶어버릴 수 있는 거죠. 이건 꽤 중요한 개념이니 꼭 기억해두시길 바랍니다.

집중해주세요!

• 뜬금 없지만 여기서 정말 중요한 말씀을 드리고 싶은데요, 세상에는 동일한 결과를 내주는 수많은 명령어 조합, 코드 조합이 존재합니다. 그리고 생물학자로서 데이터를 다룰 때는 거의 모든 경우에, 뭔 삽질을 하든, 얼마나 돌아서 가든, 결과만 정확하다면 아무 문제가 없습니다.

• 이는 생물학에서 다루는 수많은 문제들은 처리하는 속도가 그렇게 중요하지 않기 때문입니다. 유전체 데이터는 당연히 엄청나게 크지만, 앞으로 배울 프로그램들을 잘 활용할 수만 있다면 시간이 오래 걸리는 부분은 모두 아주 쉽게 해결됩니다. 위대하신 분들이 이미 복잡한 프로그램을 모두 개발해서 세상에 공짜로 뿌려두셨기 때문입니다. 우리가 할 일은 그 프로그램에서 나온 결과물을 다루는 일이 됩니다. 그리고 이렇게 전처리가 끝난 데이터는 보통 수천 줄에서 수십만 줄 수준밖에 되지 않습니다. 물론 많은 양이고, 직접 눈으로 들여다본다면 몇 시간, 아니 몇 십 시간이 들겠죠. 그렇지만 컴퓨터를 이용하면 아무리 느려도 5초면 끝날 겁니다. 삽질을 수백 번해서 잘못 처리한다고 해도 수천 초, 1-2시간밖에 걸리지 않을 테죠. 게다가 여러분의 실력은 점점 늘 겁니다. 그러면 그렇게 쓰이는 시간은 점점 짧아져 어떤 데이터든 원하는 대로 몇 초만에 처리할 수 있게 됩니다.

• 그러니 이상한 명령어 조합 또는 한참이나 돌아가는 명령어 조합을 쓰는 것을 결코 두려워하지 마세요. 뭔 결과를 내든 여러분이 눈으로 들여다보는 것보다는 빠르고, 그렇게 삽질하다 보면 점점 좋은 방식으로 나아갈 수 있을 겁니다.

• 물론 제 연구실 사람이라면 작성한 명령어 조합 등을 정리해서 보여주시길 바라겠습니다. 좀 더 효율적인 명령어 조합 등을 알려드리겠습니다. 또 박사과정 하다 보면 수십억 줄은 되는 데이터를 다루게 될 때도 있을 텐데요, 그런 데이터라면 제가 코드를 짜드리거나 스스로 짤 수 있도록 도와드릴 테니 걱정 안 해도 됩니다. 성장하는 데 집중하시길 바랍니다.

• 이번에는 참조 유전체(reference genome)에 대해서도 간략하게 살펴봅시다. 참조 유전체란 특정한 생물 종의 DNA 정보를 고품질로 확보해둔 것입니다. 가장 초기에는 생어 시퀀싱 기법 등을 이용해서 제작했고, 한때는 숏리드 시퀀싱 기법을 이용해서 많이들 제작했지만, 요새는 거의 모든 사람들이 롱리드 시퀀싱 기법을 이용해 다양한 생물의 유전체 지도를 작성하고, 이를 참조 유전체로 활용하고 있습니다.

• 먼저 예제 데이터를 다운 받아야 합니다. 다음 명령어를 입력해주시기 바랍니다.

wget -O ceN2.fa.gz https://downloads.wormbase.org/releases/WS289/species/c_elegans/PRJNA13758/c_elegans.PRJNA13758.WS289.genomic.fa.gz
gzip -d ceN2.fa.gz

• 이 데이터는 다세포생물 중 가장 처음으로 유전체 지도가 완성된 생물인 예쁜꼬마선충(Caenorhabditis elegans)의 참조 유전체 데이터입니다. wget은 웹에서 데이터를 다운 받을 때 많이 쓰는 명령어이고, gzip은 압축하거나 압축 풀 때 쓰는 명령어입니다. 각각 -O와 -d가 붙어있는데, 이런 옵션들이 무슨 역할을 하는지는 "help" 페이지를 읽어보거나 해당 내용을 검색해서 공부해봅시다.
• 파일을 다운 받았으니, 마찬가지로 ls -l 등을 이용해서 데이터 크기를 살펴봅시다. 매우 크니까 nano 같은 텍스트 에디터를 이용해서 ceN2.fa라는 파일을 여는 것은 권하지 않습니다. 파일이 워낙 크다 보니 그대로 컴퓨터가 한동안 멈춰버릴 겁니다. less 등으로 파일을 열어보시기 바랍니다.
• 그러면 DNA 정보가 FASTA 형식으로 쭉 적혀 있는 걸 볼 수 있을 겁니다. 참고로 파일 이름이 *.fa로 되어있는데요, 이렇게 .fa로 끝나는 파일은 FASTA 포맷이라는 걸 가리킵니다. FASTQ 포맷은 보통 *.fq 등으로 적혀 있습니다.
• 파일 앞쪽을 head를 이용해 열어 봅시다.

어쩌구@저쩌구:~/02_parsing_matrix$ head ceN2.fa
>I
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
GCCTAAGCCTAAAAAATTGAGATAAGAAAACATTTTACTTTTTCAAAATTGTTTTCATGC
TAAATTCAAAACGTTTTTTTTTTAGTGAAGCTTCTAGATATTTGGCGGGTACCTCTAATT

• 그러면 예쁜꼬마선충의 1번 염색체 왼쪽 끝에 담긴 서열 정보가 나타나게 됩니다. 이전에도 말씀 드렸듯 > 표시는 해당 서열의 이름이 그 줄에 나온다는 뜻이 되고, 여기서는 I이 염색체의 이름이 됩니다. 앞쪽에는 GCCTAA가 반복해서 나타나는데요, 이는 예쁜꼬마선충의 텔로미어 반복서열인 TTAGGC가 reverse complement되어있는 형태입니다. 보통 진핵생물의 유전체 지도가 염색체 수준으로 잘 완성돼 있다면, 이처럼 왼쪽 끝과 오른쪽 끝에 모두 텔로미어 반복서열이 cluster를 이루며 등장하게 됩니다.
• 이번에는 이 파일을 이용해 염색체 정보를 뜯어보는 연습을 좀 더 해봅시다. 먼저 이 참조 유전체에 들어있는 염색체 또는 염색체의 일부 조각 개수가 얼마나 되는지 살펴보도록 합시다. 어렵지 않게 살펴볼 수 있는데요, FASTA 포맷에서는 모든 DNA 정보의 이름에 >가 포함되어있기 때문입니다. 그러니 >의 개수만 세면, 염색체나 기타 등등 DNA 조각이 몇 개나 있는지 알 수 있죠. 다음 명령어를 입력해봅시다.

어쩌구@저쩌구:~/02_parsing_matrix$ grep ">" ceN2.fa
>I
>II
>III
>IV
>V
>X
>MtDNA

• 그러면 위와 같이 DNA의 이름 정보가 쭉 출력될 겁니다. 참고로 예쁜꼬마선충은 상염색체를 5개, 성염색체를 1개, 미토콘드리아 염색체를 1개 지니고 있습니다. 각각 I, II, III, IV, V, X, MtDNA라는 이름이 붙어있고요.

집중해주세요!

• 여기서 한 가지 주의할 사항을 말씀 드리려고 합니다. > 문자를 사용할 때는, 이게 리디렉션으로 인식될 수 있다는 점을 반드시 주의하셔야 합니다. 컴퓨터에서 이 기호가 리디렉션으로 잘못 인식되면 작업하려던 파일이 빈 파일이 되는 참사가 일어날 수 있기 때문입니다.
• 예를 들면 grep ">" ceN2.fa라고 쳐야 하는데, 쌍따옴표를 빼먹어서 grep > ceN2.fa라고 쳤다고 해봅시다. 이러면 >가 grep의 검색 대상이 아닌 리디렉션으로 인식되게 됩니다. grep의 결과를 ceN2.fa라는 파일 이름으로 저장하라는 명령어가 되는 셈이죠. 그리고 이렇게 되면 ceN2.fa 파일은 곧장 빈 파일이 되어버리게 됩니다.
• 이유는 다음과 같습니다. grep > ceN2.fa라는 명령어를 보면, 마치 앞에서부터 차례대로 실행될 거라고 생각하기 쉽습니다. 다시 말해 grep으로 뭔가를 하려고 한 다음, 그 결과가 STDOUT으로 나오면 이를 STDIN으로 리디렉션해서 ceN2.fa라는 파일에 저장하라는 식으로 말이죠. 이 경우라면 별 문제가 없을지도 모릅니다. 왜냐면 grep은 각 행을 검색해서 우리가 검색하고 싶은 문자열이 있는지를 확인하는 명령어인데, grep > ceN2.fa라는 명령어에서 >가 리디렉션으로 인식됐다면 grep에는 검색할 문자열도, 검색할 파일도 지정이 되지 않았으니 그냥 헬프 페이지를 띄우겠죠.
• 문제는 이 실행 순서가 정반대로 일어난다는 데 있습니다. 요컨대 grep > ceN2.fa라는 명령어가 실행되면, 일단 다짜고짜 ceN2.fa라는 빈 파일을 새로 생성합니다. 그 다음에 grep을 실행하죠. 그러면 짜잔! 분석하려던 데이터가 싹 날아가고 grep은 헬프 페이지를 띄워줍니다.

어쩌구@저쩌구:~/02_parsing_matrix$ grep > ceN2.fa
Usage: grep [OPTION]... PATTERNS [FILE]...
Try 'grep --help' for more information.
어쩌구@저쩌구:~/02_parsing_matrix$ cat ceN2.fa
어쩌구@저쩌구:~/02_parsing_matrix$                  # 파일이 텅 비어서 아무것도 출력하지 않음.

• 분석하려던 파일이 쉽게 복구할 수 있는 것이라면 별 타격이 없겠지만, 몇날며칠 고생해서 정리해둔 데이터라면 진짜 피눈물 납니다. 조심 하시기 바랍니다.

• 이번 단계에서 배워야 할 건 이 정도면 끝난 것 같습니다. 마찬가지로 두 번째 과제를 한번 따라해보시면서 좀 더 익숙해지시길 바라겠습니다.

Conda를 이용한 프로그램 설치

• 지금까지 가장 기본적인 명령어들을 이용해서 커맨드라인과 간단한 행렬 데이터 처리에 익숙해지는 과정을 거쳤습니다. 이번 단계에서는 좀 더 본격적인 프로그램들을 사용하는 방법들에 대해 익혀봅시다.
• 프로그램을 활용하는 데 있어서 가장 큰 걸림돌이 되는 것 중 하나는 설치가 쉽지 않았다는 점입니다. 하지만 이제는 Conda나 Docker, Singularity 등 다양한 기법들이 개발되면서 프로그램 설치가 아주 쉬워졌습니다. 프로그램 설치까지만 되면 본인 데이터만 잘 활용해서 프로그램을 돌리기만 하면 끝나거든요. 우리가 할 일은 이런 기법들을 활용해 원하는 프로그램을 설치하고 그 프로그램을 잘 쓰는 수준이면 충분합니다.
• 이러한 다양한 기법 중에 가장 보편적으로 쓰이는 프로그램 설치법 중 하나는 Conda입니다. Conda를 쓰면 프로그램마다 각자의 환경(environment)을 부여할 수 있는데, 이렇게 하면 프로그램들끼리 충돌하지 않고 쉽게 설치된다는 강력한 장점이 있습니다. 환경이랑 충돌이라는 게 뭔지는 이제부터 설명 드리겠습니다.
• 이전에 프로그램들을 설치하기 어려웠던 이유는, 프로그램들마다 서로 다른 작동 환경을 필요로 했기 때문입니다. 환경이란 프로그램이 작동하기 위해 필요한 다양한 하위 프로그램, 라이브러리 등을 가리킨다고 보시면 됩니다. 2개 이상의 서로 다른 프로그램이 있어야 데이터 분석을 하는데, 그 두 프로그램들이 작동하는 환경이 서로 배타적이라면 2개 프로그램을 동시에 설치할 수 없었던 겁니다.
• 분자생물학 실험으로 비유하자면 제한효소 2개를 써서 DNA를 잘라야 하는데, 이 2개의 제한효소를 작동시킬 수 있는 버퍼는 세상에 존재하지 않는 상황이랑 비슷합니다. 별 수 있나요? 실험을 나눠서 하든지 해서, 두 제한효소가 서로 다른 화학적 환경에서 작동하도록 하는 수밖에 없죠. 프로그램 설치도 이랬습니다.
• 예를 들면 이런 겁니다. 프로그램 1은 라이브러리 A의 버전 1.9가 있어야 하는데, 프로그램 2는 라이브러리 A의 버전이 3.0은 되어야 한다고 생각해봅시다. 그러면 내 컴퓨터에 프로그램 1과 프로그램 2를 동시에 설치하려면, 라이브러리 A는 대체 어떤 버전을 써야 하는 걸까요? 라이브러리 A의 버전을 1.9로 맞추자니 프로그램 2를 쓸 수가 없고, 반대로 3.0으로 맞추면 프로그램 1을 쓸 수가 없는 문제가 생깁니다. 이런 문제 때문에 컴퓨터 안에 엄청나게 다양한 라이브러리를 버전마다 다 따로 설치하고, 프로그램마다 따로따로 라이브러리를 짝지어줘야 했는데, 이 작업이 보통 귀찮은 게 아니었습니다. 그러니 설치가 쉽지 않아 프로그램을 돌려보지도 못하고 포기하는 일이 부지기수였죠.
• Conda는 이런 문제를 손쉽게 해결해줍니다. Conda를 활용하면 프로그램마다 서로 다른 환경을 부여해 각자에게 필요한 라이브러리 등을 짝지어주는 게 아주 쉬워집니다. 예를 들면 Conda를 활용하면, 환경 1을 만들어 라이브러리 A 버전 1.9와 프로그램 1을 설치하고, 동시에 같은 컴퓨터에 환경 2를 만들어 라이브러리 버전 3.0과 프로그램 2를 설치하는 게 명령어 두 줄이면 가능합니다. 엄청나게 편해진 거죠.
• 그럼 지금부터 이 Conda를 서버에 설치하고, 이를 활용해 다양한 생물학 관련 프로그램들을 설치하는 과정을 한번 따라가봅시다. Conda 홈페이지에서 본인 운영체제에 맞는 걸 골라 다운로드 하고 설치해봅시다.
• 저처럼 Ubuntu 운영체제를 사용한다면 다음 단계를 따르시면 됩니다. 다른 운영체제는 제가 안 해봐서 모르겠어요. 찾아보셔야 할 것 같습니다.

cd
mkdir 03_conda && cd 03_conda
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh   # 인스톨러 다운로드
bash Miniconda3-latest-Linux-x86_64.sh                                       # 인스톨러 실행. *.pl을 perl로 실행했던 것처럼, *.sh 파일은 bash로 실행하면 됩니다.

• 이렇게 하면 엔터 치고 약관 읽으라는 말이 나올 겁니다. 엔터를 치면서 내리시면 되고, 동의하시면 yes 치시면 됩니다.
• 그 뒤 어디에 설치할지 설치 경로를 물어보는데, 원하는 곳 지정하시면 되겠습니다.
• 마지막으로 Conda 실행에 필요한 초기 설정을 하겠냐고 물어보는데, yes 고르시면 됩니다.
• 그러면 화면에 뭐가 엄청나게 뜰 겁니다. 터미널로 돌아오고 나면, 다음과 같이 입력하고 난 뒤 터미널이 변화하는지 확인하시면 됩니다.

어쩌구@저쩌구:~/03_conda$ source ~/.bashrc           # Conda 관련된 기본 환경을 작동시켜줌.
(base) 어쩌구@저쩌구:~/03_conda$                   # 왼쪽에 (base)가 뜨면 성공

• 이제 Conda 설치가 끝났습니다. 쉽죠? 이제는 conda라는 명령어를 활용할 수 있게 됐습니다. 마찬가지로 터미널에 다음과 같이 쳐서 헬프 페이지가 잘 나오는지 확인해봅시다.

(base) 어쩌구@저쩌구:~/03_conda$ conda
usage: conda [-h] [-V] command ...
[후략]

• (base)라고 적힌 부분이 현재 작업 중인 Conda 환경의 이름을 보여주는 부분입니다. 현재 환경은 (base)라는 말 그대로 가장 기반이 되는 기본 환경이라는 뜻이 되겠습니다. 이후에 새로운 환경을 만들고 그 환경을 활성화시키면, 예를 들면 (assembly), (rnaseq) 등등으로 (base) 부분이 바뀌게 될 겁니다. 한 가지 작업만 더 해놓고 바로 새로운 환경을 만들어봅시다.
• Conda가 잘 설치됐다면, 이번에는 Mamba를 설치할 차례입니다. Conda는 참 좋은 도구이지만 몇 가지 문제가 있는데요, 그 중 하나가 프로그램을 여러 개 동시에 설치하려고 하면 시간이 엄청나게 오래 걸리고 잘 안 된다는 점입니다. 이런 문제를 해결해주는 게 Mamba입니다. Mamba를 활용하면 다양한 프로그램이 아주 빠르게 설치됩니다.
• 앞으로 더 자세하게 말씀 드리겠습니다만, 앞으로는 conda 명령어 대신 mamba 명령어를 쓰시면 됩니다.
• Mamba를 설치하고 싶으시면 마찬가지로 공식 홈페이지에서 알려주는 설치 방법을 따르시면 됩니다. 2023년 7월 6일 현재에는 다음 명령어를 입력해서 Mamba를 설치할 수 있습니다. 중간에 y 한 번 눌러서 설치 진행해주세요.

(base) 어쩌구@저쩌구:~/03_conda$ conda install -n base --override-channels -c conda-forge mamba 'python_abi=*=*cp*'
[중략]
(base) 어쩌구@저쩌구:~/03_conda$ mamba
usage: mamba [-h] [-V] command ...
[후략]

• 위에 보이는 것처럼 설치가 끝나면 이제 mamba를 치기만 하면 헬프 페이지가 나오게 될 겁니다. 다음 명령어도 쳐줍시다.

(base) 어쩌구@저쩌구:~/03_conda$ mamba init
(base) 어쩌구@저쩌구:~/03_conda$ source ~/.bashrc

• 이제 모든 준비는 끝났습니다. 이제 새로운 환경을 만들어봅시다. 다음과 같이 명령어를 쳐주세요.

(base) 어쩌구@저쩌구:~/03_conda$ mamba create -n basicGenomics

• 마찬가지로 설치하겠냐는 말이 나올 텐데 y 누르고 엔터 치면 됩니다. 환경이 새롭게 만들어졌습니다. 이제 이 환경을 활성화시켜봅시다.

(base) 어쩌구@저쩌구:~/03_conda$ mamba activate basicGenomics
(basicGenomics) 어쩌구@저쩌구:~/03_conda$ 

• 그러면 위에 보이는 것처럼, 왼쪽에 환경 이름을 보여주는 곳이 (basicGenomics)로 바뀌게 될 겁니다.
• 이제 이 환경에서 새로운 프로그램들을 설치할 수 있고요, 설치가 끝나면 자유롭게 쓸 수 있습니다. 다음과 같이 명령어를 입력해줍시다.

(basicGenomics) 어쩌구@저쩌구:~/03_conda$ mamba install -c conda-forge -c bioconda assembly-stats samtools bioawk seqtk hisat2 bedtools bcftools bwa fastqc minimap2 hifiasm svim svim-asm busco
# 참고: conda-forge 또는 bioconda는 다양한 프로그램이 공개돼있는 채널입니다. Mamba는 이러한 채널을 뒤져서 가장 적합한 프로그램들을 죄다 긁어모다 자동으로 내가 제공한 프로그램 리스트를 설치해줍니다.
# 여기서 설치시킨 프로그램은 총 14개입니다(assembly-stats, samtools, bioawk, seqtk, hisat2, bedtools, bcftools, bwa, fastqc, minimap2, hifiasm, svim, svim-asm, busco).

• 이렇게 입력하고 나서 y를 치면 우리가 앞으로 분석에 쓸 프로그램이 한번에 그리고 자동으로 설치됩니다. 참 쉽죠? 제 서버에서는 5.5초만에 설치가 끝났네요.
• 혹시 내가 쓰고 싶은 프로그램이 더 있다면 직접 검색해보시면 됩니다. 검색창에 "(프로그램이름) bioconda" 등으로 검색하면 "anaconda" 또는 "bioconda" 웹페이지가 나올 겁니다. 예를 들면 "bwa bioconda"라고 검색하시면 이런 페이지가 검색될 텐데요, 이 페이지에 적힌 설치 방법을 따라하시면 됩니다. 정확히는 conda install -c bioconda bwa라고 적혀 있는데, 여기서 conda만 mamba로 바꿔주시면 됩니다. mamba install -c bioconda bwa 이런 식으로 말이죠.

집중해주세요!

• Conda의 또 다른 문제점은 base 환경에 프로그램을 설치하다 보면 기존에 깔린 모든 conda 환경이 망가질 가능성이 존재한다는 것입니다. 이렇게 되면 프로그램 싹 다 다시 설치해야 돼서 진짜 열받겠죠? 게다가 논문 쓰던 도중이라면 프로그램 버전이 바뀌지 않도록 신경 써야 하는데, 프로그램을 재설치하다 보면 버전까지 바뀔 수 있다 보니 정말정말 고통스러워질 수 있습니다. 논문 쓸 때는 프로그램 버전도 다 기록을 해줘야 하거든요.

• 그러니 새로운 프로그램을 설치하고 싶다면 반드시 새로운 환경을 먼저 생성합시다. 그 다음에 그 환경을 활성화시킨 뒤 새로운 프로그램을 설치하도록 합시다. 그래야 인생이 피곤해지는 일을 덜 마주칠 수 있습니다.

• 작업이 다 끝났다면 원래 환경으로 돌아가볼까요? 다음과 같이 진행해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/03_conda$ mamba deactivate
(base) 어쩌구@저쩌구:~/03_conda$

• 그러면 다시 원래 환경으로 돌아오게 될 겁니다.
• 이번에는 새로운 환경에 다른 프로그램도 설치해봅시다.

(base) 어쩌구@저쩌구:~/03_conda$ mamba create -n repeat
(base) 어쩌구@저쩌구:~/03_conda$ mamba activate repeat
(repeat) 어쩌구@저쩌구:~/03_conda$ mamba install -c conda-forge -c bioconda repeatmodeler
# 혹시 설치 안 되면,
(repeat) 어쩌구@저쩌구:~/03_conda$ conda update -n base conda
# conda를 업데이트한 뒤 다시 다음 명령어를 실행해주세요.
(repeat) 어쩌구@저쩌구:~/03_conda$ mamba install -c conda-forge -c bioconda repeatmodeler

• 이제 또 다른 프로그램이 설치됐습니다. 이번에만 그런 건지는 모르겠지만, 저는 이 repeatmodeler가 다른 프로그램과는 동시에 설치가 안 되더라고요. 이런 식으로 동시에 설치되지 않는 충돌이 벌어지더라도, 새로운 환경에 설치하면 아무 문제 없습니다. 아주 편한 세상입니다.
• 이번 단계에서는 이 정도 하면 된 것 같습니다. 다음 단계부터는 이번에 설치한 다양한 프로그램들을 하나하나씩 써먹어봅시다.
• 이번 단계 과제는 별 건 없습니다. 다음 질문을 읽고 답을 검색해서 알아둡시다.

기본적인 프로그램 다뤄보기

• 이번 단계에서는 이전에 Conda를 이용해 설치한 다양한 프로그램들을 하나씩 다뤄봅시다.
• 가장 먼저 써볼 프로그램은 bioawk입니다. 시퀀싱과 관련된 다양한 데이터 포맷을 다룰 때 쓸 수 있는 프로그램이고, 유용한 기능이 많습니다. 먼저 다음과 같이 쳐봅시다.

(base) 어쩌구@저쩌구:~$ mkdir 04_basic_programs && cd 04_basic_programs
(base) 어쩌구@저쩌구:~/04_basic_programs$ bioawk
bioawk: command not found

• 위에 보이는 것처럼, (base) 상태에서 bioawk를 검색하면 프로그램이 작동하지 않습니다. 이전에 말씀 드렸듯, bioawk를 설치한 건 (basicGenomics)라는 특정한 환경이기 때문입니다. 그러니 프로그램을 쓰기 전에 그 환경을 먼저 활성화시켜줘야겠죠. 다음과 같이 진행해봅시다.

(base) 어쩌구@저쩌구:~/04_basic_programs$ mamba activate basicGenomics
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk
usage: bioawk [-F fs] [-v var=value] [-c fmt] [-tH] [-f progfile | 'prog'] [file ...]

• 위에 보이는 것처럼 bioawk의 사용법 정보가 나오면 성공입니다. Conda를 통해 제대로 설치가 됐다는 뜻이죠.
• 이제 bioawk를 이용하는 방법을 몇 가지 배워봅시다.

Bioawk 활용하기

• bioawk는 기본적으로는 awk와 거의 똑같이 돌아간다고 보시면 됩니다. 중요한 차이는 awk가 행렬로 이뤄진 파일의 내용을 직접적으로 다룬다면, bioawk는 시퀀싱 관련 데이터를 행렬처럼 바꿔서 다룬다는 것입니다. 조금 더 자세하게 설명 드리도록 하겠습니다.
• 이전에 써본 것처럼, awk는 행렬을 다룹니다. 각 행(row)에 있는 여러 열(column)에 대해 다양한 연산을 처리해서 그 결과를 보여주죠. 예를 들면 awk '{print $1}' file.txt라는 명령어를 입력하면 첫 번째 열의 정보를 화면에 프린트해주고, awk '{print $23}' file.txt라는 명령어를 입력하면 스물세 번째 열의 정보를 출력해주는 식입니다.
• awk는 어떤 파일을 집어넣든, 그 파일의 내용을 행렬로 인식해서 처리해주는데요, 보통 각 열을 구분하는 방식은 필드 구분 문자(field separator)를 중심으로 이뤄집니다. 기본적으로는 빈칸을 구분 문자로 인식하기 때문에 사람이 인식하는 것과 비슷하게 빈칸이 있으면 새로운 열이라고 인식하는 식이죠. 구분 문자를 직접 지정해줄 수도 있습니다. 예를 들면 awk -F "TT" '{print $1,$2,$3}'라는 명령어를 입력하면, 스페이스건 탭이건 빈칸은 전부 무시하고 TT가 나올 때만 열로 구분해준 뒤, 1번 2번 3번 열의 정보를 출력해줍니다. 실제로 해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ echo "ATTGCCTAATTCG" > awk.test1.txt
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ awk -F "TT" '{print $1,$2,$3}' awk.test1.txt
A GCCTAA CG

• awk는 이런 식으로, 우리가 입력해준 구분 문자를 활용해 데이터를 행렬로 변환하고 처리하는 것이죠.
• bioawk도 매우 비슷하게 작동합니다. awk와 다른 점은, bioawk는 지정해준 포맷에 맞춰 데이터를 행렬로 변환하고 처리해준다는 것입니다.
• 예제의 결과를 들여다보면서 bioawk가 작동하는 방식을 이해해봅시다. 다음과 같이 명령어를 입력해주세요.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ echo -e ">testSeqName1\nATTGCCTAATTCG\n>testSeqName2\nAAGTCGATCGATCG" > bioawk.test1.txt  # 간단한 FASTA 파일 생성
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ cat bioawk.test1.txt                                      # 파일 내용 확인
>testSeqName1
ATTGCCTAATTCG
>testSeqName2
AAGTCGATCGATCG
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{print $name}' bioawk.test1.txt         # FASTA 파일을 행렬로 변환해 "이름"에 대한 정보 추출
testSeqName1
testSeqName2
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{print $seq}' bioawk.test1.txt          # FASTA 파일을 행렬로 변환해 "서열"에 대한 정보 추출
ATTGCCTAATTCG
AAGTCGATCGATCG
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{print $qual}' bioawk.test1.txt         # FASTA 파일을 행렬로 변환해 "퀄리티"에 대한 정보 추출


(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$

• 그러면 위와 같은 결과가 나올 겁니다. 기본적으로 bioawk를 활용하려면, 내가 입력하는 파일을 어떤 포맷으로 인식해야 하는지를 bioawk에게 알려줘야 합니다. 그게 -c fastx라고 적힌 부분입니다. 이러면 내가 입력하는 파일 포맷이 FASTA 또는 FASTQ이니, 그에 맞게 처리해달라고 bioawk에게 명령하는 셈이 됩니다. 이제 bioawk는 입력 파일을 FASTA/Q 포맷에 맞춰서 적절하게 행렬로 변환해주죠. 그리고 각 열을 가리키는 것은 name, seq, qual 등이 되게 됩니다. (물론 현재 입력 파일은 FASTA 파일이라 퀄리티 정보는 없으니 해당 값은 처리되지 않습니다)

name seq qual
testSeqName1 ATTGCCTAATTCG 
testSeqName2 AAGTCGATCGATCG 

• 이렇게 행렬이 만들어진 뒤에는 awk와 똑같이 작동하게 되는 겁니다. 몇 번째 컬럼이냐를 다루는 대신, 각 포맷에 해당하는 정보에 맞춰서 컬럼이 지정된다는 게 다를 뿐이죠. 더 많은 정보는 공식 홈페이지나 다른 분들이 만들어둔 매뉴얼을 읽어보시면 확인 가능합니다.
• 그러면 bioawk를 이용해 시퀀싱 데이터를 살짝 다뤄봅시다. 먼저 우리가 확보한 DNA 시퀀스의 길이를 측정해볼까요? 다음 명령어를 입력해주시면 됩니다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{print length($seq)}' bioawk.test1.txt
13
14

• 그러면 위와 같이 시퀀스의 길이를 출력해줄 겁니다. 말씀 드렸듯 -c fastx는 우리가 분석하려고 넣은 인풋 파일이 FASTA/Q 포맷이라는 걸 가리키는 거고요, '{print length($seq)}' 부분이 시퀀스의 길이를 출력하라는 부분이 됩니다. print는 말 그대로 출력하라는 것이며, 구체적으로는 length($seq), 즉 시퀀스($seq)의 길이(length함수)를 재서 그 값을 print하라는 말이라고 보시면 되겠습니다. 참고로 length()는 awk에서도 특정 필드의 문자열 길이를 잴 때 똑같이 쓸 수 있습니다.
• 다음으로는 DNA를 reverse complement해봅시다. 마찬가지로 함수를 이용하면 쉽게 작동하는데요, 다음과 같이 revcomp() 함수를 활용하면 쉽게 결과를 얻을 수 있습니다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{print revcomp($seq)}' bioawk.test1.txt
CGAATTAGGCAAT
CGATCGATCGACTT

• 이번에는 조금 더 어려운 걸 해봅시다. 간단한 조건문을 다뤄볼게요. 조건문이란 우리가 원하는 조건에 해당하는 데이터만 얻어낼 수 있도록 하는 기법이라고 보시면 됩니다. 예를 들면 길이가 시퀀스 길이가 20 bp 이상인 것만 출력하게 한다든지, 5 bp 이상 30 bp 미만만 출력하게 한다든지 할 수 있습니다. 다음 명령어를 입력해서 예제 데이터를 만들어 봅시다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ echo -e ">testSeqName3\nAGTACGTCAGTCAGCTAGCTAGCTAGCATCGATCGATCGATCGATCGATCGAT\n>testSeqName4\nCAGTCGATCGATCAGTCAGTCAGCTAGT\n>testSeqName5\nAGTT" >> bioawk.test1.txt

• 조건문은 다음과 같이 활용할 수 있습니다. 20 bp 이상을 먼저 출력해볼까요.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{ if(length($seq)>=20){print ">"$name; print $seq} }' bioawk.test1.txt
>testSeqName3
AGTACGTCAGTCAGCTAGCTAGCTAGCATCGATCGATCGATCGATCGATCGAT
>testSeqName4
CAGTCGATCGATCAGTCAGTCAGCTAGT
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx 'length($seq)>=20{print ">"$name; print $seq}' bioawk.test1.txt
>testSeqName3
AGTACGTCAGTCAGCTAGCTAGCTAGCATCGATCGATCGATCGATCGATCGAT
>testSeqName4
CAGTCGATCGATCAGTCAGTCAGCTAGT

• 위와 아래에 넣은 조건문은 형태만 다를뿐 둘 다 같은 결과를 출력해줍니다. 여기서는 if()라는 조건문을 활용했는데요, 이는 괄호 안에 들어있는 조건에 부합하는 경우(TRUE인 경우)에만 뒤에 나오는 중괄호 안의 명령을 실행합니다. 여기서는 length($seq)>=20이라는 조건을 줬으니, 시퀀스의 길이가 20 이상인 경우에만 중괄호 안의 명령이 실행되겠죠. 중괄호 안에 들어있는 명령어는 시퀀스 데이터를 FASTA 포맷에 맞춰 출력하라는 명령어입니다. print ">"$name 부분은 먼저 > 문자 다음에 해당 시퀀스의 이름($name)을 출력하라는 말이 되죠. FASTA의 이름 양식을 맞춘 것입니다. 그 뒤에 나오는 ; 표시는 엔터 치라는 뜻이 되고요, print $seq은 시퀀스를 출력하라는 명령어가 됩니다. 즉 시퀀스 길이가 20 bp 이상인 경우에만 >시퀀스 이름을 출력한 뒤 엔터 치고 나서 시퀀스를 출력하게 되는 것입니다. 길이가 짧으면 아무런 명령어도 실행하지 않으니 출력되는 게 없습니다.
• 이번에는 시퀀스 개수와 이 파일의 전체 시퀀스 길이를 확인하는 명령어를 살펴봅시다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ bioawk -c fastx '{number+=1; sum+=length($seq)}END{print number,sum}' bioawk.test1.txt
5	112

• 여기서 number+=1 부분은, 각 행(row)마다 number라는 변수에 숫자 1을 계속 더하라는 뜻입니다. 계속 진행해서 파일 전체를 처리하고 나면 number에는 전체 행의 숫자가 저장되겠죠. sum+=length($seq)도 비슷합니다. 각 행마다 sum이라는 변수에 각 시퀀스의 길이를 더 하라는 뜻이 됩니다. 파일 전체를 처리하고 나면 sum에 전체 시퀀스 길이의 합이 저장될 겁니다. END는 그 앞의 명령어를 먼저 파일 끝까지 처리하라는 의미이며, 다 끝나고 나면 END 뒤에 있는 명령어가 실행되게 됩니다. 뒤에 있는 명령어인 print number,sum은, 이전까지 저장된 number와 sum이라는 변수에 저장된 값을 출력하라는 뜻입니다. 그러니 화면에 결과적으로 시퀀스 개수인 5와 전체 길이의 합인 112 (bp)가 출력되는 것이죠.
• 이외에도 bioawk를 이용해 다양한 데이터 처리를 손쉽게 진행할 수 있습니다. 위에서 언급한 예제 페이지에 들어가서 한번 따라해보셔도 좋습니다. 공식 홈페이지에는 별 내용이 없으니 참고하세요.

Seqtk 활용하기

• 이번에는 다른 유용한 프로그램인 seqtk도 한번 사용해봅시다. 참고로 bioawk와 seqtk는 모두 위대하신 Heng Li 선생님이 개발하셨습니다. 브로드 인스티튜트(Broad Institute) 방향으로 절하고 쓰시면 됩니다.
• 마찬가지로 다음과 같이 명령어를 입력해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ seqtk

Usage:   seqtk <command> <arguments>
Version: 1.4-r122

Command: seq       common transformation of FASTA/Q
         size      report the number sequences and bases
         comp      get the nucleotide composition of FASTA/Q
         sample    subsample sequences
         subseq    extract subsequences from FASTA/Q
         fqchk     fastq QC (base/quality summary)
         mergepe   interleave two PE FASTA/Q files
         split     split one file into multiple smaller files
         trimfq    trim FASTQ using the Phred algorithm
[후략]

• 그러면 위와 같은 설명이 뜰 겁니다. 가장 중요한 건 "Usage" 부분인데요, 여기에 적힌 것처럼 seqtk <command> <arguments>를 입력하면 이 프로그램을 활용할 수 있습니다. 이중 ""는 화면에 출력된 것처럼 seq, size, comp 등을 가리키는데요, 이러한 커맨드는 seqtk에 내장된 하위 분석 모듈이라고 보시면 됩니다. 이중 size를 활용해보면서 연습해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ seqtk size
Usage: seqtk size <in.fq>   # 인풋 파일로 FASTQ 파일을 넣으라는 뜻
(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ seqtk size bioawk.test1.txt
5	112

• 어디서 많이 본 숫자죠? 이전에 bioawk를 이용한 결과와 똑같은 결과가 나온 셈입니다. 지난번에도 말씀 드렸듯, 같은 결과에 도달하기 위한 방법은 무수히 많습니다. 어떻게든 정확한 답에만 도달하시면 대부분 문제 없습니다.

assembly-stats 활용하기

• 마지막으로 assembly-stats를 사용해봅시다. 이 프로그램은 FASTA/Q 파일에 대해 좀 더 자세한 통계수치를 뽑아줍니다.

(basicGenomics) 어쩌구@저쩌구:~/04_basic_programs$ assembly-stats bioawk.test1.txt
stats for bioawk.test1.txt
sum = 112, n = 5, ave = 22.40, largest = 53
N50 = 28, n = 2
N60 = 28, n = 2
N70 = 28, n = 2
N80 = 14, n = 3
N90 = 13, n = 4
N100 = 4, n = 5
N_count = 0
Gaps = 0

• 그러면 위와 같이 정보를 줄 겁니다. 첫 번째 줄은 파일 이름을 알려주고 있네요. 두 번째 줄에 적힌 sum = 112는 이 파일에 총 112 bp의 데이터가 들어있다는 것을 가리키며, n = 5는 전체 시퀀스 개수가 5개임으로, ave = 22.40은 시퀀스의 평균 길이가 22.4 bp임을 보여줍니다. largest = 53은 제일 긴 시퀀스 길이가 53 bp임을 보여주며, 다음 줄에 나오는 N50 = 28, n = 2는 N50이라는 특정한 수치에 해당하는 시퀀스의 길이가 28 bp이며, 전체 중 두 번째로 긴 시퀀스임을 보여주고 있습니다.
• N50은 나름 중요한 수치이니 좀 더 알아봅시다. 이는 유전체 조립(genome assembly) 과정에서 많이 쓰이는 수치인데요, 보통은 N50 길이가 길면 길수록 좋습니다. N50 값을 계산하는 방식은 다음과 같습니다. 먼저 DNA 리드의 겹치는 부분을 최대한 이어붙였을 때 모든 시퀀스를 가장 긴 것부터 나열한 뒤, 차례차례 그 길이를 더합니다. 그리고 그 길이의 합이 전체 길이의 합의 50%를 넘었을 때, 마지막으로 더한 시퀀스의 길이가 N50이 됩니다. 그리고 그 시퀀스가 몇 번째로 긴 시퀀스였는지가 N50 = 28, n = 2에서 n = 부분에 적히게 되는 것입니다. 좀 더 자세한 내용은 다음 두 링크를 살펴보시면 되겠습니다. 위키피디아 및및 그림이 포함된 자료
• 한 가지 주의해야 할 점은, assembly-stats는 압축이 풀린 FASTA/Q 파일에 대해서만 작동한다는 것입니다. 압축이 안 풀려 있으면 작동하지 않습니다. 또한 표준 출력된 결과를 인풋으로 받아 처리하지도 못해서 파이프 마지막에 쓸 수도 없습니다. 뭔가 처리한 FASTA/Q의 통계치를 확인하고 싶다면, 그 결과를 일단 파일로 저장한 뒤 assembly-stats를 사용하시면 되겠습니다.
• 이제 이번 단계에서 배울 건 다 배운 것 같습니다. 마찬가지로 이번 단계 과제를 진행하면서 이러한 프로그램들에 익숙해져 봅시다.

간단한 변이 추출 따라해보기

• 이번 단계부터는 시퀀싱 데이터를 활용해 생물학적 의미를 찾는 방법을 본격적으로 배워보겠습니다. 가장 먼저 해볼 일은 변이 추출(variant calling)입니다.
• 거의 모든 생물은 저마다 서로 다른 DNA를 지니고 있습니다. 사람만 해도 저마다 서로 DNA가 다르고, 이러한 DNA의 차이 중 일부는 심각한 유전병과 같은 결과로 이어지기도 합니다. 이처럼 개체마다 서로 다른 DNA를 유전 변이(genetic variant)라고 부르며, 그 결과 나타나는 생물의 차이를 표현형 변이(phenotypic variation)라고 부릅니다.
• 사람을 대상으로 유전학 연구를 하는 수많은 연구자들은 다양한 표현형 변이의 원인이 되는 유전 변이를 찾아내는 것을 목표로 합니다. 유전자를 특정하는 것뿐만 아니라, 그 유전자에 생길 수 있는 다양한 변이 중 대체 무엇이 인구 집단 내에서 표현형 변이를 일으키는 것인지, 변이마다 얼마나 크게 영향을 주는지 등을 연구하기도 합니다. 가장 대표적인 예는 사람의 키 차이를 연구한 사례죠. 2022년에 발표된 연구에서는 500만 명을 조사해 키에 영향을 주는 1만 2천 여 개의 단일염기다형성(single-nucleotide polymorphisms; SNPs)을 밝혀내기도 했습니다.
• 이런 연구를 해내려면 일단 표현형이 다른 집단을 모으고, 이러한 사람들이 지닌 유전 변이를 찾아내는 일을 먼저 진행해야 합니다. 요즘은 보통 사람의 혈액 등을 수집하고, DNA 시퀀싱을 진행해 막대한 양의 데이터를 생산합니다. 그러곤 그 DNA를 참조 유전체와 비교 분석해 개개인이 지닌 유전 변이를 먼저 찾아냅니다. 마지막으로 이 유전 변이와 표현형 변이의 연관성을 조사함으로써, 표현형이 다른 원인 유전 변이를 확인하는 일을 진행하죠.
• 이번 단계에서는 이 유전 변이를 찾아내는 프로그램을 직접 다뤄보고자 합니다. DNA 분석의 가장 기초적인 단계일뿐만 아니라 거의 모든 DNA 데이터 분석에 쓰이는 단계이니만큼 잘 익혀두시길 기대합니다.
• 변이 추출은 (다양한 필터링 단계가 추가되긴 하지만) 크게 인덱싱(indexing), 매핑(mapping), 콜링(calling)이라는 3단계를 거칩니다. 이에 대한 개념을 간략하게 이해해보고, 실제 프로그램을 돌려봅시다. 숏리드 시퀀싱 데이터 처리를 먼저 익혀보겠습니다.
• 변이 추출은 기본적으로는 참조 유전체와 개개인의 DNA가 다른 부분을 찾아내는 과정입니다. 예를 들면 참조 유전체의 염색체 1번 48만 번째 위치가 A였는데, 어떤 사람은 그 자리가 엄마 염색체도 A, 아빠 염색체는 T라면, 이는 A/T 헤테로(heterozygous)라고 부릅니다. 다른 사람은 A/A 호모(homozygous)일 수도 있으며, 또 다른 사람은 T/T 호모일 수도 있죠. 이러한 정보를 얻어내는 과정이 위 3단계인 셈입니다.
• 이러한 변이 추출을 하려면 가장 먼저 시퀀싱된 리드(read)를 처리해줘야 합니다. 이전에 말씀 드렸듯 리드란 한번에 읽어낸 DNA 정보가 담겨있는 가장 작은 조각을 가리킵니다. 숏리드 시퀀싱이라면 보통 100 bp에서 200 bp 정도 길이의 DNA 서열 정보가 하나의 리드에 담기게 됩니다.
• 보통은 변이 추출의 정확성을 높이기 위해 한 자리가 30번 정도 독립적으로 읽을 수 있도록 시퀀싱을 진행합니다. 사람은 이배체(diploid)다 보니 한 자리가 적어도 2회는 읽혀야 엄마한테 받은 염색체와 아빠한테 받은 염색체 정보를 정확하게 확인해줄 수 있겠죠. 또 유전체가 30억 염기쌍인 걸 감안하면, 우연히 특정한 위치는 다른 지역보다 2-3배 더 많이 읽히고, 또 다른 위치는 다른 지역보다 2-3배 더 적게 읽힐 수도 있습니다. 게다가 시퀀싱 자체에 오류가 생겨 A인 자리가 T, G, C 등으로 읽힐 수도 있을 테니, 오류가 생길 걸 감안해 3배 정도는 읽어줘야 그 오류를 보정할 수 있을 겁니다. 이러한 모든 경우를 감안해서 보통은 평균적으로 각 위치가 30회 정도 읽힐 수 있도록 시퀀싱을 진행합니다. 사람이라면 30억 염기쌍의 30배인 900억 염기쌍을 생산하는 겁니다. 이 정도면 변이 추출이 상당히 정확하게 진행될 수 있다는 것이 알려져 있기도 합니다.
• 그리고 이렇게 많은 데이터가 있다 보니, DNA 길이가 짧다는 게 굉장히 큰 문제가 됩니다. 각각의 100개짜리 DNA 조각이 대체 염색체 몇 번의 어느 위치에서 나온 DNA인지 알 수가 없다는 문제가 생기기 때문입니다.
• 그렇기 때문에 숏리드 시퀀싱 데이터를 분석하는 과정에는 바로 이런 DNA 조각의 제 위치를 찾아주는 일이 굉장히 중요한 단계가 됩니다. 900억 염기쌍에 달하는 수 억 개의 DNA 조각의 위치를 찾아줘야 하는 일이니 엄청난 연산이 들어가야 하기도 합니다. 한 조각 한 조각 비슷한 위치를 검색해주는 일을 사람이 한다면 도저히 할 수 없는 일이겠지만, 컴퓨터를 이용하면 상당히 빠른 속도로 진행할 수 있습니다.
• 게다가 이미 이런 일을 다루기 위한 다양한 프로그램들이 아주 오래 전부터 개발돼서 널리 쓰이고 있습니다. 가장 대표적인 프로그램인 BWA 논문은 인용수가 4만 회에 육박합니다. 많은 사람들이 쓸 정도로 쓰기도 쉽습니다. 그래서 코딩을 할 필요도 없습니다. 우리가 할 일은 프로그램을 설치하고, 인풋 파일과 아웃풋 파일을 지정해 엔터만 쳐주면 됩니다. 설치도 이미 해놨으니 직접 돌려보기만 하면 되겠네요.
• 이제부터 bwa를 직접 사용해봅시다. 다음과 같이 명령어를 입력해주세요.

(basicGenomics) 어쩌구@저쩌구:~$ mkdir 05_variant_calling && cd 05_variant_calling              # 새 폴더 형성
(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ ln -s ../04_basic_programs/wbcel235.mt.fa   # C. elegans 미토콘드리아 DNA에 대한 바로 가기 형성
(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ ln -s ../04_basic_programs/srr3440952.sub.1.fq.gz                # C. elegans 숏리드 시퀀싱 데이터에 대한 바로 가기 형성
(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ ln -s ../04_basic_programs/srr3440952.sub.2.fq.gz                # C. elegans 숏리드 시퀀싱 데이터에 대한 바로 가기 형성
# 참고로 여기서 바로가기가 잘 생성돼 있다면 파일 이름이 파란색 등으로 표시됩니다. ls 등을 이용해 확인해주세요.
# 파일 이름이 빨간색이라면 바로가기가 형성되지 않은 겁니다.
# 새로 바로가기를 형성하려면 이 빨간색 파일들을 rm으로 삭제한 뒤 진행해주세요.
# 혹시 이전 단계에서 srr3440952.sub.1.fq.gz 및 srr3440952.sub.2.fq.gz 파일 압축 풀었다면,
# 이번 기회에 gzip 등을 이용해서 압축하는 걸 진행해보시길 바랍니다. 아니면 bwa의 인풋 파일 이름을 바꿔주셔도 됩니다.
(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ bwa index wbcel235.mt.fa                    # C. elegans 미토콘드리아 DNA를 인덱싱
[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 0.00 seconds elapse.
[bwa_index] Update BWT... 0.00 sec
[bwa_index] Pack forward-only FASTA... 0.00 sec
[bwa_index] Construct SA from BWT and Occ... 0.00 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index wbcel235.mt.fa
[main] Real time: 0.006 sec; CPU: 0.005 sec                                                    # 미토콘드리아 정도는 0.01초 안에 끝남

• 이렇게 입력하면 샘플 데이터를 현재 디렉토리에서 활용할 수 있게 되고, 참조 유전체(여기서는 미토콘드리아 DNA)에 대한 인덱스 정보가 생성됩니다. 생물학하는 사람들에게 인덱스가 뭔지 각 잡고 설명하려면 작은 책 한 권 필요한 수준이라 여기서 자세한 내용을 다루긴 어렵고요, 검색을 빠르게 하기 위한 준비 단계라는 정도만 알고 계시면 됩니다.
• 예를 들면 DNA 100개짜리 리드를 검색한다고 쳤을 때, 30억 염기쌍이나 되는 유전체의 처음부터 일일이 검색하면 속도가 매우 느리겠죠? 심지어 수 억 개의 리드를 매번 그렇게 검색한다고 치면 한도 끝도 없이 오랜 시간이 소요될 겁니다. 그렇기 때문에 다양한 프로그램은 처음부터 하나하나 검색하는 대신, 저마다 특정한 알고리즘을 활용해 이런 검색 과정을 매우 효율적으로 진행합니다. BWA는 버로우즈-휠러 변환이라는 방식을 활용하는데요, 관심 있으시면 읽어보셔도 좋을 것 같습니다.
• 이렇게 인덱싱이 끝나고 나면 본격적인 검색 과정을 진행할 수 있는데요, 이 과정을 매핑(mapping)이라고 부릅니다. 참조 유전체에서 리드가 있어야 할 위치를 검색하고 제 위치를 알려주는 과정입니다. 다음과 같이 입력하면 결과를 얻을 수 있습니다.

(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ bwa mem wbcel235.mt.fa srr3440952.sub.1.fq.gz srr3440952.sub.2.fq.gz > srr3440952.sub_to_wbcel235.mt.sam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 20000 sequences (1956180 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 0, 0, 0)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] skip orientation FR as there are not enough pairs
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 20000 reads in 0.509 CPU sec, 0.509 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem wbcel235.mt.fa srr3440952.sub.1.fq.gz srr3440952.sub.2.fq.gz
[main] Real time: 0.562 sec; CPU: 0.563 sec                                        # 데이터가 작아서 0.5초면 끝남

• 그러면 매핑이 자동으로 끝납니다. 참 쉽죠? 아웃풋 파일은 SAM (Sequence Alignment/Map) 포맷으로 저장돼 있는데요, 이에 대해서 간단하게 알아봅시다.
• 아웃풋 파일인 srr3440952.sub_to_wbcel235.mt.sam을 살펴봅시다. 먼저 wc -l 등으로 행(row) 수를 확인해봅시다. 이전 숙제를 마치셨다면 알겠지만, 이 아웃풋 파일은 인풋 파일들에 들어있는 리드 수보다 2줄이 추가돼 있습니다. @로 시작하는 헤더 부분이 추가돼서 그렇습니다. 이 헤더에는 참조 유전체에 대한 간단한 정보와 실제로 입력한 명령어 등이 담기게 됩니다. 나중에 내가 무슨 버전 썼는지 까먹었더라도 이런 걸 활용하면 관련 정보를 추출할 수 있습니다.
• 각 행에는 각 리드별 매핑 정보가 담겨 있습니다. 간단하게 설명 드리자면 1번 열에는 우리가 인풋으로 넣은 리드 이름이 적혀 있습니다. 보통 참조 유전체에 검색하는 대상을 쿼리(query)라고 불러서, 해당 열을 QNAME이라고 부르기도 합니다. 2번 열에는 FLAG 정보가 담겨 있는데요, 이는 각 리드가 어떻게 매핑되어있는지 알려주는 것입니다. 이는 이진법으로 적혀 있는데요, 예를 들면 77이라는 FLAG는 1+4+8+64로 쪼갤 수 있고, 각각은 다음과 같은 의미를 지닙니다.

1 = 조각 나있다
4 = 이 조각은 매핑이 안 됐다
8 = 다른 조각은 매핑 안 됐다
64 = 이게 첫 번째 조각이다

• 3번 열에는 이 리드가 매핑된 참조 유전체의 염색체(또는 컨티그/스캐폴드)의 이름이 들어가게 되고, 4번 열에는 그 위치가, 5번 열에는 매핑이 얼마나 잘 됐는지 등의 정보가 담기게 됩니다.
• 이외에도 다양한 정보가 담겨 있는데요, 자세한 설명은 공식 홈페이지에 들어가신 뒤 SAM 포맷에 대한 정보를 담고 있는 이 PDF 파일을 열어서 읽어보시면 됩니다. 직접 프로그램을 개발하실 일이 있다면 이런 포맷에 대한 정보는 꼭 숙지해두시는 게 좋습니다.
• 그런데 이 아웃풋 파일을 열어보시면 알겠지만, 매핑 안 된 리드가 맨 위에 나올 정도로 정돈이 안 돼 있다는 게 문제입니다. 이런 파일을 정렬해주는 프로그램도 이미 잘 갖춰져 있습니다. 다음과 같이 입력해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ samtools sort srr3440952.sub_to_wbcel235.mt.sam -o srr3440952.sub_to_wbcel235.mt.sorted.sam
(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$

• samtools는 마찬가지로 4만 회 이상 인용된 정말 범용적인 프로그램입니다. 이전에 사용했던 seqtk처럼 다양한 하위 프로그램 모듈이 내장돼 있으며, samtools sort는 그 중 하나입니다. 이는 인풋 SAM 파일을 정렬해주며, 아웃풋 이름은 -o [아웃풋 파일 이름]을 이용해서 지정할 수 있습니다.
• 이제 정렬이 끝난 srr3440952.sub_to_wbcel235.mt.sorted.sam 파일을 살펴봅시다. 마찬가지로 행 수를 살펴보면 2줄이 추가된 걸 알 수 있는데요, 정렬 과정에서 헤더가 추가됐기 때문입니다. head 등을 이용해서 내용을 살펴보면, 정렬된 덕분에 결과가 다르게 출력되는 걸 확인할 수 있을 겁니다.

(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ head -5 srr3440952.sub_to_wbcel235.mt.sorted.sam
@HD     VN:1.6  SO:coordinate
@SQ     SN:MtDNA        LN:13794
@PG     ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem wbcel235.mt.fa srr3440952.sub.1.fq.gz srr3440952.sub.2.fq.gz
@PG     ID:samtools     PN:samtools     PP:bwa  VN:1.17 CL:samtools sort -o srr3440952.sub_to_wbcel235.mt.sorted.sam srr3440952.sub_to_wbcel235.mt.sam
SRR3440952.4460953      113     MtDNA   639     60      100M    =       639     0       AGATTATTTTTAAAATTTTCTTATGTTTTAGGGGAAATAATGTTTTTTTATTTTATGTGTTTTTCTGTTATTTCAAGAATCCTGGGTATGGTAGTTATAG    B/FFFFFFFFFFB<<BFBFFFFFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBBBB NM:i:0  MD:Z:100        MC:Z:100M       AS:i:100        XS:i:0

(QUIZ) 이번에는 처음 나오는 리드에 더 많은 정보가 기입돼 있는 걸 확인할 수 있을 겁니다. 어떤 의미인지 한번 생각해봅시다.

• 이렇게 하면 매핑이 모두 끝난 겁니다. 실제로는 다양한 필터링을 거쳐야 하긴 하지만 그건 나중에 차근차근 배우시면 됩니다.
• 이쯤 되면 이런 결과의 의미를 파악하는 게 어려워서 그렇지, 프로그램을 돌리는 것 자체는 매우 쉽다는 걸 알 수 있을 겁니다.
• 수많은 시퀀싱 데이터 분석은 이처럼 표준화된 프로그램들을 기반으로 진행됩니다. 다른 분석도 생각 이상으로 무척 쉽습니다. 다양한 분석에 도전해보시길 바랍니다.
• 이제 마지막 단계인 콜링을 진행해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/05_variant_calling$ bcftools mpileup -f wbcel235.mt.fa srr3440952.sub_to_wbcel235.mt.sorted.sam | bcftools call -mv -Ov -o srr3440952.sub.calls.vcf

• 이것도 이렇게 진행하면 끝입니다. 그러면 VCF 포맷으로 파일이 만들어질 겁니다. 마찬가지로 긴 헤더와 실제 변이 데이터가 나와야 하는데요, 예제 데이터를 대충 만들었더니 결과가 안 나오네요. 다른 데이터 넣으면 마찬가지로 잘 돌아갑니다. (QUIZ) BCFtools도 마찬가지로 널리 쓰이는 프로그램 중 하나입니다. mpileup, -f, call, -mv, -Ov, -o 등에 대해서 검색해서 알아봅시다.
• 마지막으로 이번 단계 과제에 답해보면서 변이 추출 과정에 좀 더 익숙해져 봅시다.
• 다음에는 이러한 변이들이 유전자에 어떻게 영향을 주는지 살펴보도록 하겠습니다.

변이가 주는 영향을 이해해보기

• 많은 유전학 연구자들이 변이를 연구하는 이유는 변이가 유전자에 영향을 미칠 수 있다는 사실 때문입니다. 변이는 유전자에서 만들어지는 단백질의 서열을 바꿔 단백질의 기능 자체를 바꾸거나 망가뜨릴 수 있으며, 때로는 유전자의 조절 부위에 영향을 줘 단백질이 만들어지는 양이나 시공간을 뒤바꾸기도 합니다. 물론 단백질이 아닌 RNA 수준에서 기능하는 다양한 유전자에 영향을 주기도 하죠. 그리고 이러한 다양한 영향들은 세포의 작동 방식을 바꾸기도 하며, 이는 조직이나 기관의 기능과 개체의 활동에까지 영향을 줄 수 있습니다.
• 다시 말해 이러한 유전 변이 중 일부는 표현형 변이에 영향을 줄 수 있습니다. 그렇기에 유전 변이 수준에서 표현형 변이를 이해하고 그 원인을 제거하거나 조정해 우리가 원하는 형태로 생물 안의 반응을 조절하는 것이 가능합니다. 예컨대 특정한 유전병의 원인이 되는 유전 변이를 이해함으로써 그것이 단백질에 미치는 영향을 이해하고, 이를 기반으로 새로운 약물 개발에 뛰어드는 것 같은 사례가 대표적이겠죠.
• 이번 단계에서는 이러한 변이가 단백질 등에 미치는 영향을 조사할 수 있는 방법에 대해 알아보려고 합니다. 이것도 프로그램이 이미 많이 만들어져 있는데요, 그 중에서도 가장 널리 쓰이고 아주 편하게 활용할 수 있는 프로그램 중 하나인 SnpEff를 다뤄보고자 합니다.
• 다음과 같이 명령어를 입력해봅시다.

(basicGenomics) 어쩌구@저쩌구:~$ mkdir 06_snpeff && cd 06_snpeff
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ unzip snpEff_latest_core.zip
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ ls
snpEff/	snpEff_latest_core.zip

• 이렇게 snpEff라는 디렉토리가 형성되면 성공입니다. 이제 이 디렉토리 안으로 들어가서 파일들을 하나씩 살펴봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ cd snpEff/
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff$ ls
LICENSE.md	examples/	galaxy/	snpEff.config	SnpSift.jar	exec/	scripts/	snpEff.jar 

• 버전이 크게 바뀌지 않는 이상, 위와 같은 파일들이 만들어져 있을 겁니다. 이중에서 우리가 바로 써볼 건 snpEff.jar 파일과 examples 디렉토리 안에 들어있는 예제 파일들입니다. 먼저 예제 파일들을 들여다봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff$ cd examples/
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff/examples$ ls

그러면 매우 많은 VCF 포맷의 파일들이 나올 겁니다. 전부 예제 파일인데요, 이중 제일 작은 파일 중 하나인 test.chr22.vcf 파일을 한번 열어서 내용을 살펴봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff/examples$ head test.chr22.vcf
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
22	17071756	.	T	C	.	.	.
22	17072035	.	C	T	.	.	.
22	17072258	.	C	A	.	.	.
22	17072674	.	G	A	.	.	.
22	17072747	.	T	C	.	.	.
22	17072781	.	C	T	.	.	.
22	17073043	.	C	T	.	.	.
22	17073066	.	A	G	.	.	.
22	17073119	.	C	T	.	.	.

• 이건 가장 기본적인 VCF 파일입니다. 첫 번째 행(row)에는 헤더 정보가 적혀 있고요, 그 다음 행부터는 1번 열(column)에는 염색체 이름, 2번 열에는 변이가 존재하는 위치, 3번 열에는 변이의 이름(없을 경우엔 .으로 표시)이 표시됩니다. 그리고 4번 열에는 참조 유전체의 시퀀스(REF)가, 5번 열에는 새로 시퀀싱해서 이 참조 유전체와 비교했을 때 차이가 나는 경우 그 서열(ALT)이 적히게 됩니다(같은 경우에도 적힐 수 있긴 함). 예를 들어 두 번째 행에 담긴 정보는, 참조 유전체의 염색체 22번 17,071,756 bp에는 T가 있는데, 새로 시퀀싱한 사람은 그 자리가 C로 바뀌어있다는 뜻이 되겠습니다.
• VCF 포맷을 활용하면 시퀀싱 데이터 결과를 엄청나게 압축할 수 있습니다. DNA 시퀀싱 데이터로 따지면 한 사람당 900억 염기쌍의 데이터가 생산되지만, VCF 포맷을 활용하면 보통은 이중 변이가 있는 위치만 데이터가 담기게 되기 때문입니다. 예를 들면 보통 사람과 사람을 비교하면 수백만 개의 변이가 존재하는데요, 900억 개를 수백만 개로 줄일 수 있다면 1만 분의 1 수준으로 데이터를 줄일 수 있는 셈이죠. 매우 유용하게 활용할 수 있습니다.
• 2023.08.22 기준, database download가 되질 않습니다. snpEff에서 관리하던 링크가 깨졌네요. 다음과 같이 진행하면 됩니다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff/examples$ cd ../
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff$ wget https://snpeff.blob.core.windows.net/databases/v5_0/snpEff_v5_0_GRCh37.75.zip && unzip snpEff_v5_0_GRCh37.75.zip

• 이제 이러한 변이가 유전자에 끼치는 영향을 실제로 확인해보도록 합시다. 먼저 초기 디렉토리로 이동해봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff/snpEff$ cd ../
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$

• 그 뒤 nano를 활용해 00_snpEff.sh라는 새로운 텍스트 파일을 열어준 뒤, 다음 내용을 입력하고 저장해봅시다.

#!/bin/bash
# java = Programming language
# -Xmx8g = Memory size
# snpEff/snpEff.jar = snpEff program
# -v = Verbose mode
# GRCh37.75 = Human genome data for snpEff
# snpEff/examples/test.chr22.vcf = Example file containing variant information of human chromosome 22
java -Xmx8g -jar snpEff/snpEff.jar -v GRCh37.75 snpEff/examples/test.chr22.vcf > test.chr22.ann.vcf

• ls를 이용해서 00_snpEff.sh 파일이 잘 생성됐는지 확인해보시고, cat을 이용해서 파일 내용이 우리가 원하는 대로 잘 작성됐는지도 확인해봅시다.
• 이 00_snpEff.sh 파일은 bash를 이용해 실행시킬 수 있는 파일입니다. 이전에 perl을 이용해서 *.pl 파일을 실행시켰던 것과 비슷합니다.
• 이런 파일을 작성하면 좋은 이유는, 한번 작성해두면 매번 똑같은 명령어를 작동시킬 수 있기 때문입니다. 예를 들어 bash 00_snpEff.sh라는 명령어를 실행시키면, 그 안에 적힌 내용에 따라 매번 똑같은 프로그램이 작동하게 됩니다. 아주 편해지는 셈이죠.
• 이제 이 00_snpEff.sh 파일의 내용을 한번 살펴봅시다. 위에서부터 7번 줄까지는 전부 문장 맨 앞에 #이 붙어있는데요, 이는 해당 행에 적혀 있는 정보는 명령어가 아니라는 표현이 됩니다. 그러니 컴퓨터에서 해당 문장은 실행되지 않습니다. 보통은 프로그램에 대한 추가적인 해설을 적어놓는 용도로 쓰이며, 주석(annotation)이라고 부릅니다.
• 주석의 첫 번째 행에는 #!/bin/bash가 적혀 있는데, 이는 이 파일에 담긴 내용을 bash를 이용해서 실행해야 한다는 뜻이 됩니다.
• 두 번째 행부터는 SnpEff를 실행하는 데 필요한 데이터와 옵션에 대한 정보가 담겨 있습니다. 실제로 실행되는 명령어는 java -Xmx8g -jar snpEff/snpEff.jar -v GRCh37.75 snpEff/examples/test.chr22.vcf > test.chr22.ann.vcf인데요, 각각에 대한 해설이 적혀 있습니다. 쭉 한번 읽어보시면 되겠습니다.
• 이제 우리가 작성한 00_snpEff.sh 파일을 실행시켜 봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ bash 00_snpEff.sh            # 기본 실행
(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ nohup bash 00_snpEff.sh &    # 추천 실행 방식

• 윗줄과 아랫줄은 동일한 결과를 줍니다만, 후자가 더 편하고 좋습니다. 두 번째 명령어인 nohup bash 00_snpEff.sh &을 치고 나면 화면에 뭐라고 뜰 텐데요, 그냥 엔터 치시면 됩니다.
• 이 nohup은 no hang up의 약자인데요, 문자 그대로 컴퓨터 꺼도 프로그램이 계속 돌아가게 해주는 매우 유용한 명령어입니다. 지금이야 돌리는 프로그램들이 금방금방 끝나지만, 나중에 본격적으로 연구하다 보면 사나흘씩 프로그램을 돌려야 하는 경우도 종종 생깁니다. 그런데 그렇게 오래 돌리다 보면 문제가 생길 수 있습니다. 예를 들어 인터넷이 불안정해서 서버랑 접속이 끊긴다든지, 서버는 멀쩡한데 누가 내 컴퓨터 전원을 꺼버린다든지, 아니면 터미널을 닫아버린다든지 하면 며칠씩 돌아가던 프로그램이 그대로 멈춰버립니다. 그러면 처음부터 다시 돌려야 해요. 몇 번 당하면 진짜 성질 나빠집니다.
• nohup을 이용하면 이런 문제가 애초에 일어나지도 않도록 막아줄 수 있습니다. nohup을 이용해서 프로그램을 실행하면 내 컴퓨터가 꺼지든 서버랑 접속이 끊기든, 프로그램은 서버 안에서 계속 돌아가기 때문입니다. 그러니 정전처럼 서버 본체가 꺼지는 일이 발생하지 않는 이상 마음 편히 프로그램을 돌릴 수 있습니다. 아주 유용하죠.
• 이런 nohup을 제대로 활용하려면, 명령어 맨 끝에 &를 붙여주면 가장 좋습니다. 그러면 이제 백그라운드(background)에서 프로그램이 실행되게 됩니다. 프로그램이 뒤에서 계속 실행되고 있다는 뜻인데요, 이렇게 프로그램이 백그라운드에서 실행되도록 하면 프로그램이 계속 돌아가든 말든 다시 터미널로 곧장 돌아오게 됩니다. 그러니 같은 터미널에서 계속 새로운 작업을 진행하는 게 가능하죠. 매우 유용합니다.
• 이렇게 글을 읽는 사이에 아마 00_snpEff.sh가 전부 돌아가서 결과를 내줬을 겁니다. 제대로 돌아갔는지 확인해보려면 엔터 한 번 쳐보세요. 화면에 [1]+ DONE nohup bash 00_snpEff.sh 같은 내용이 뜰 겁니다. 만약 안 뜨면 아직도 돌아가고 있는 건데요, 이것도 확인하고 싶으면 ps를 쳐봅시다. 화면에 현재 돌아가고 있는 프로그램들을 띄워주는 명령어인데, 만약 백그라운드에서 잘 돌아가고 있다면 여기에 뜰 겁니다.
• 잘 돌아갔다면 4개 파일이 새로 만들어졌을 겁니다. ls를 이용해서 확인해보시면 nohup.out, snpEff_summary.html, test.chr22.ann.vcf, snpEff_genes.txt 이 4개 파일이 만들어진 게 보일 겁니다.
• nohup.out은 우리가 nohup으로 작동시킨 bash 00_snpEff.sh 명령어의 로그(log)를 담고 있습니다. less 등으로 열어보면 snpEff.jar가 실행되면서 남긴 다양한 정보가 담겨있다는 걸 알 수 있을 겁니다. 여기에는 어떤 참조 유전체가 쓰였는지, 거기에 유전자는 몇 개나 있는지 등등의 정보가 담겨 있습니다. 한번 어떤 게 있는지만 대충 살펴보시면 됩니다.
• 다음으로 snpEff_summary.html 파일은 그래픽 기반으로 통계치 등을 보여주는 파일입니다. 좀 더 이해하기 쉬울 테니 이것도 한번 살펴봅시다.
• 먼저 snpEff_summary.html 파일을 서버에서 내 컴퓨터로 옮겨봅시다. MobaXterm을 쓰고 있다면, 화면 맨 왼쪽에 있는 노란색 지구(?) 모양을 클릭하시고, 그 아래쪽에 있는 "Follow terminal folder" 버튼을 클릭해서 체크 표시 해주세요. 그 다음에 터미널을 다시 클릭한 뒤 엔터를 치면 현재 디렉토리에 있는 정보를 보여줄 겁니다. 만약 안 된다면 pwd를 쳐서 현재 디렉토리 위치를 복사하고, 왼쪽 위에 디렉토리 위치 적혀 있는 곳에 붙여넣기 하시면 이동할 겁니다.
• 그러면 현재 디렉토리에 있는 파일 중 snpEff_summary.html 파일이 보일 텐데요, 이걸 바탕 화면 등으로 드래그 하시면 파일이 자동으로 복사됩니다. 쉽죠?
• 맥이나 다른 리눅스 기반 컴퓨터를 써서 서버에 접속하고 있다면 터미널에서 빠져나간 뒤 scp -rP 포트번호 본인서버아이디@서버아이피:디렉토리위치/파일이름 .을 치시면 서버에 있는 파일이 현재 디렉토리로 복사됩니다.
• 그러곤 방금 파일 복사한 디렉토리로 이동한 다음 마우스로 더블 클릭해서 snpEff_summary.html 파일을 여시면 됩니다. 그러면 다양한 통계치와 의미 있는 숫자 같은 것들을 보여줄 거예요. 편하게 보시면 됩니다.
• 이제 다시 터미널로 돌아와서 SnpEff가 내놓은 다른 아웃풋 파일들도 내용을 한번 확인해봅시다.
• 먼저 test.chr22.ann.vcf 파일을 열어봅시다. 이 파일은 기본적으론 인풋으로 넣었던 test.chr22.vcf 파일에 담긴 변이 정보를 거의 동일하게 담고 있는데요, 마지막 열(column)인 INFO에 새로운 정보가 추가돼 있습니다. 새롭게 추가된 이 정보는 각 변이가 다양한 유전자에 끼치는 영향에 대해 정리해놓은 것입니다. 이를 활용해서 다양한 분석을 할 수 있습니다만, 내용이 너무 많아서 쉽게 다루기는 어려울 겁니다. 이 파일에는 어떤 내용이 들어있는지 한번 훑어보는 정도로만 살펴보시면 됩니다.
• 마지막 파일인 snpEff_genes.txt는 꽤 유용하고 직관적입니다. 이 파일에는 개별 유전자가 다양한 변이에 어떻게 영향을 받고 있는지에 대한 정보가 담겨 있습니다. 다음 명령어를 활용해 살펴봅시다.

(basicGenomics) 어쩌구@저쩌구:~/06_snpeff$ head -3 snpEff_genes.txt
# The following table is formatted as tab separated values. 
#GeneName	GeneId	TranscriptId	BioType	variants_impact_HIGH	variants_impact_LOW	variants_impact_MODERATE	variants_impact_MODIFIER	variants_effect_3_prime_UTR_variant	variants_effect_5_prime_UTR_premature_start_codon_gain_variant	variants_effect_5_prime_UTR_variant	variants_effect_conservative_inframe_deletion	variants_effect_conservative_inframe_insertion	variants_effect_disruptive_inframe_deletion	variants_effect_disruptive_inframe_insertion	variants_effect_downstream_gene_variant	variants_effect_frameshift_variant	variants_effect_intron_variant	variants_effect_missense_variantvariants_effect_non_coding_transcript_exon_variant	variants_effect_splice_acceptor_variant	variants_effect_splice_donor_variant	variants_effect_splice_region_variant	variants_effect_start_lost	variants_effect_stop_gained	variants_effect_stop_lost	variants_effect_stop_retained_variant	variants_effect_synonymous_variant	variants_effect_upstream_gene_variant
A4GALT	ENSG00000128274	ENST00000249005	protein_coding	1	10	9	32	1	1	0	0	0	0	0	1	0	90	0	0	0	0	0	0	0	9	0

• 그러면 위와 같이 엄청 긴 정보가 나오는데요, 첫 번째 행(row)에 적혀 있듯 기본적으로 이 파일은 TSV 포맷으로 구성돼 있다는 걸 알 수 있습니다. 두 번째 행에는 각 열(column)이 어떤 영향을 주는 변이에 대한 값인지 그 의미가 적혀 있습니다. 그리고 세 번째 행부터는 실제 유전자에 대해 각 값을 보여주고 있죠.
• 이중 가장 많이 쓰이고 유용한 열은 1번부터 7번 열까지 담겨있는 정보들입니다. 이는 유전자의 이름뿐만 아니라 각 유전자에 영향을 주는 변이가 얼마나 강하게 영향을 주고 있는지, 그에 해당하는 변이의 개수는 몇 개나 되는지에 대한 정보를 담고 있죠. 여기서는 그 영향을 Impact라고 표현하고 있으며, 그 정도에 따라 high, low, moderate, modifier 등으로 분류하고 있습니다. 각각에 해당하는 영향을 좀 더 자세히 살펴보고 싶으시다면 공식 홈페이지를 들여다보시면 좋을 것 같습니다.
• 이처럼 유전자별로 영향 받은 정도를 활용하여, 유전자가 개개인에게서 얼마나 망가져있을지 등을 추정할 수 있습니다.
• 예를 들어 유전병 환자와 비-환자 집단을 수백 명씩 모으고 시퀀싱했다고 생각해봅시다. 이렇게 시퀀싱한 개개인이 지닌 변이는 한 사람만 수백만 개일 테니 이런 과정을 통해서 확보한 유전 변이만 해도 수천만에서 수억 개에 이를 수 있겠죠.
• 이렇게 많은 것 중 유전병이 원인이 되는 변이가 대체 무엇인지 어떻게 알 수 있을까요? SnpEff처럼 단백질이 영향 받은 정도를 추정할 수 있다면, 이러한 변이의 수를 극도로 줄일 수 있습니다. 예를 들면 high impact처럼 단백질을 크게 뒤바꾸는 변이에 집중할 수 있기 때문입니다. 이렇게 필터링을 거치면 분석이 더 단순해질 수 있는 겁니다.
• 이번 단계의 과제는 이러한 SnpEff 결과를 좀 더 다뤄보는 것을 목표로 합니다. 마지막으로 과제를 진행하면서 복습하고, 좀 더 유연하게 활용할 수 있는 방안들을 익혀봅시다.

유전체를 조립하기

• 이번 단계부터는 본격적으로 롱리드 시퀀싱 데이터를 다뤄봅시다.

롱리드 시퀀싱이 지닌 장점

• 숏리드 시퀀싱 데이터 분석과 비교했을 때 롱리드 시퀀싱 데이터를 다루는 일은, 할 수 있는 일이 좀 더 추가되고 프로그램이 조금 다르다는 정도 말곤 큰 차이가 없습니다. 아직까지는 상당히 비싸다 보니 소규모 데이터가 많고, 그래서 이해하기 좀 더 편하다는 장점도 있습니다. 숏리드 시퀀싱 데이터로 1만 명씩 시퀀싱 해서 논문 내는 일은, 돈도 돈이지만 배경 지식을 이해하는 데에만 몇 년은 걸리거든요. 롱리드 시퀀싱 데이터 분석은 두세 명 분석으로도 (아직까지는) 논문이 되기 때문에 고도의 지적 능력이나 통계 기법을 요구하진 않습니다.
• 이번 단계부터 다루는 프로그램들만 잘 다뤄도 작은 논문은 꾸준히 낼 수 있는 정도의 분석 능력을 갖출 수 있습니다. 물론 데이터 생산 비용이 매우 비싸다 보니 실제 논문을 쓰려면 돈을 많이 써야 하는 문제가 있긴 합니다. 그래도 앞으로 점점 저렴해질 테니 잘 배워서 많이 써먹으시길 바랍니다.
• 먼저 유전체 조립(genome assembly)를 왜 하는지, 그리고 유전체 조립 과정이 무엇인지 등등 알아보도록 합시다.
• 숏리드 시퀀싱은 아주 저렴한 비용만으로도 생물의 유전체 정보를 밝혀주는 정말 유용한 기법입니다만, 몇 가지 한계를 지니고 있습니다. 가장 대표적인 한계는 반복 서열(repetitive sequence) 영역의 정보를 제대로 해석하기 어렵다는 것인데요, 이로 인해 상당수 진핵생물에서 문제가 생깁니다. 진핵생물은 전체 유전체 중 반복 서열의 비율이 30% 이상인 경우가 많기 때문입니다.
• 예를 들면 1번 염색체와 2번 염색체에 동일한 150 bp짜리 반복 서열이 각각 10개와 30개 있다고 생각해봅시다. 이런 유전체 정보를 시퀀싱을 통해 150 bp 조각으로 엄청나게 읽고 나면, 이 리드를 참조 유전체에 매핑해서 제 위치를 찾아야겠죠. 그런데 1번 염색체에 있는 10개짜리 클러스터와 2번 염색체의 30개짜리 클러스터 중, 대체 어디에 이 150 bp짜리 리드를 매핑해야 할까요? 전부 똑같다면 이 리드는 40개 어디에 가든 상관이 없을 텐데 말입니다. 제 위치를 찾는 게 불가능한 겁니다.
• 돌연변이가 생기면 제 위치를 찾기가 어렵다는 것이 더욱 문제가 됩니다. 예를 들어 염색체 1번에 있는 반복 서열 클러스터가 10개짜리가 복제돼 염색체 3번으로 끼어들어갔다고 생각해봅시다. 숏리드 시퀀싱 기법을 이용해 이를 분석하면 이 결과가 제대로 분석되지 않습니다. 늘어난 10개가 어디서 왔는지 알 수가 없기 때문입니다. 늘어난 건 아마 알 수도 있겠지만, 염색체 1번에 있던 10개짜리가 20개로 늘어난 건지, 염색체 2번에 있던 30개짜리가 40개로 늘어난 건지, 반반 나눠서 늘어난 건지, 아예 다른 염색체로 일부 또는 전체가 복제된 건지를 알 방도가 없습니다. 150 bp라는 짧은 길이로 살펴보면, 이 반복 서열은 어디서 왔든 똑같은 150개짜리로 보일 테니까요. 만약 이중 3번 염색체로 옮겨갔을 때만 질환이 생겨난다면 숏리드 시퀀싱 기법으로는 분석할 방법이 없는 겁니다.
• 특히 염색체의 안정성에 중요한 텔로미어(telomere) 인근 지역과 센트로미어(centromere) 지역은 반복서열이 엄청나게 많이 분포해 숏리드 시퀀싱 기법으로는 제대로 분석할 방법이 없다시피 했습니다. 150 bp짜리 리드뿐만 아니라 1 kb 정도되는 리드로도 정확한 분석이 불가능할 정도로 반복서열이 길고, 많고, 복잡하기 때문입니다.
• 이로 인해 2000년도 초반에 초기 결과가 공개된 인간 유전체 사업(Human Genome Project) 결과도 완벽하지 않았습니다. 특히 센트로미어와 텔로미어가 매우 가까이 위치한 아크로센트릭 염색체(acrocentric chromosome)들은 해당 지역의 시퀀스 정보가 아예 없다시피 했습니다. 이러한 기술적 한계로 인해 30억 염기쌍의 인간 유전체 지도 중 약 1억 2천만 염기쌍 정도는 모르는 채 남아있었습니다.
• 이런 문제는 롱리드 시퀀싱 기법이 빠르게 발전하면서 모두 해결되고 있습니다. 수십 kb의 염기쌍이 한번에 읽힌 리드 정보가 확보되면서 상당히 긴 반복 서열 클러스터도 한번에 읽힐 수 있게 됐기 때문입니다.
• 또한 반복 서열뿐만 아니라, 구조 변이(structural variant)라고 부르는 커다란 변이에 대해서도 정확하게 확인하는 것이 가능해졌습니다.
• 구조 변이란, 지금 시점에서는 보통 50 bp 이상의 커다란 결손/삽입 등의 변이를 가리킵니다. 이전 단계에서 숏리드 시퀀싱을 이용해 SNP 등 단일염기(1 bp) 수준에서 나타나는 변이를 확인해봤는데요, 실제로 사람과 사람은 그보다 훨씬 큰 변이도 아주 많이 지니고 있습니다. SNP처럼 수백만 개 정도가 있는 건 아니지만, DNA가 수십 개씩 결손되거나 삽입되는 경우는 한 사람당 수만 개 정도 있고 수만 개에서 수십만 개 수준으로 DNA가 바뀌는 변이도 수십 개씩 존재합니다.
• 이런 구조 변이는 숏리드 시퀀싱으로는 분석하는 게 매우 어렵거나 거의 불가능했습니다. DNA를 읽을 수 있는 길이가 대략 100 bp 수준인데, 수십 수백 bp가 결손되거나 삽입되는 걸 분석하긴 쉽지 않았기 때문입니다.
• 롱리드 시퀀싱은 이런 구조 변이를 확인하기에 아주 적합합니다. 애초에 읽어내는 DNA 길이가 수만 bp에 이르니 상당수 구조 변이는 처음부터 끝까지 덮을 수 있기 때문입니다.
• 이에 더해, 유전체 조립(genome assembly) 과정을 거치면 수십만 bp에 이르는 구조 변이 뿐만 아니라, DNA 손상 회복 과정에서 염색체가 아예 뒤섞이는 전좌 등도 엄밀하게 분석하는 것이 가능해졌습니다. 실질적으로는 개체가 지니고 있는 모든 유전 변이를 확인하는 것이 가능해진 것입니다.

유전체 조립(genome assembly)

• 유전체 조립 과정은 기본적으로는 리드들 사이에 존재하는 겹치는 부분을 이어붙여 더 긴 덩어리의 DNA 조각으로 만드는 과정입니다. 이렇게 만든 더 큰 덩어리를 보통은 컨티그(contig)라고 부릅니다.
• 현존하는 롱리드 시퀀싱 데이터를 활용하면 대개 염색체 또는 염색체에 준하는 수준, 몇 Mb 이상되는 거대한 DNA를 재구성할 수 있습니다. 유전체 조립 과정이 잘 진행된다면 때로는 거의 염색체에 준하는 거대한 덩어리로 DNA를 재구성하는 게 가능하고, 다양한 롱리드 시퀀싱 데이터를 조합하면 텔로미어에서 텔로미어까지 이어진 염색체 수준으로 재구성하는 것도 이제는 가능한 수준입니다.
• 실제로 이런 장점을 살려 30억 염기쌍이 온전히 풀린 최초의 인간 유전체 지도가 공개되기도 했습니다. (실제론 못 푼 지역이 극히 일부 남아있긴 함)
• 참고로 겹치는 부분이 존재하는 이유는, 유전체 조립을 할 때에 유전체 크기의 20배 이상되는 많은 양의 데이터를 활용하기 때문입니다. PCR 등 증폭 과정을 거치지 않고, 같은 자리를 20번 이상 읽다 보니 리드와 리드 사이에는 많은 경우 겹치는 시퀀스가 존재하게 됩니다. 그러니 이렇게 겹치는 시퀀스들을 통합해서 더 길게 이어붙일 수 있는 거죠.
• 현재 유전체 조립 분석에서 가장 표준이 되고 있는 데이터는 PacBio의 HiFi 기법입니다. 오류가 거의 없는 20 kb 가량의 리드를 생산해주기 때문에 반복 서열이 매우 길게 나타나는 지역을 제외하고는 거의 모든 지역을 풀어줍니다. 제가 분석해본 경험으로도 유전체 크기가 수백 Mb 정도라면 거의 염색체 수준의 유전체 조립이 가능하며, 사람 유전체도 N50이 적어도 10 Mb 이상 나올 정도로 분석이 잘 됩니다.
• 이렇게 거대한 덩어리로 DNA를 재구성할 수 있다 보니, 결손이건 삽입이건 전좌건 어떤 형태의 변이든 그 크기와 관계 없이 정확한 검출이 가능해진 것입니다.
• 게다가 이런 롱리드 시퀀싱 기법의 가격이 점차 저렴해지면서, 한 사람의 유전체 정보를 정확하게 확보하는 수준이 아닌 수십 수백 명에 이르는 인구 집단의 유전체 정보를 고품질로 얻어내려는 시도가 계속되고 있습니다.
• 이런 시도를 범유전체(pangenome) 사업이라고 부르는데요, 현재 전세계의 다양한 인구 집단을 반영한 47명의 고품질 유전체 지도 및 중국 내 인구 집단을 반영한 36명의 고품질 유전체 지도가 공개돼 있는 상황입니다. 게다가 빠른 시일 내에 이 숫자는 점점 가파르게 늘어날 겁니다. 숏리드 시퀀싱 기법도 10년 15년 전에는 돈 천만 원씩 들 때는 몇 명씩만 분석되다가, 값 싸지니까 분석 대상이 된 사람 숫자가 수백, 수만, 수백만 명 수준으로 늘게 된 거거든요. 롱리드 시퀀싱도 그렇게 될 겁니다.
• 그럼 이제부터 실제 데이터를 활용해서 유전체 조립 과정을 진행해봅시다.

공공 데이터 다운로드

• 먼저 예제 데이터를 다운 받아볼까요? 제가 이번에 쓰려고 하는 데이터는 특정한 돌연변이 예쁜꼬마선충의 HiFi 리드 데이터입니다. 이런 공공 데이터를 활용하려면 해당 데이터가 보관돼있는 곳으로 가서 데이터를 다운 받아야 합니다. 이러한 정보는 보통 논문의 Data availability 섹션에 담겨 있습니다.
• 이 논문에는 NCBI BioProject 또는 KoNA에 데이터가 저장돼 있다고 나오네요. NCBI는 미국에서, KoNA는 한국생명공학연구원 내 KOBIC에서 운영하고 있는 데이터 저장소입니다. 많은 경우 생산한 시퀀싱 데이터를 이용해 논문을 내려면 데이터를 이런 데이터 저장소에 업로드해야만 합니다. 둘 다 알아두고, 데이터 업로드 하는 걸 논문 낼 때쯤 되면 꼭 진행하시길 바랍니다.
• 한국 데이터 저장소인 KoNA도 좋지만, 아무래도 역사가 깊은 NCBI 데이터 저장소에 더 많은 데이터가 저장돼 있긴 합니다. 이 데이터를 다운 받는 방법을 실습해봅시다.
• 먼저 논문에 적힌 대로, NCBI BioProject 홈페이지에 들어간 뒤, "PRJNA896647"를 검색해봅시다. 그러면 다음과 같은 페이지가 나올 텐데요, 이중 "SRA Experiments"라고 적힌 부분의 오른쪽에 적힌 숫자 링크를 눌러주세요. 그러면 이 BioProject 안에는 2개의 데이터가 포함돼있다는 걸 알 수 있을 겁니다. 이런 페이지는 한번씩 쭉 읽어보면서 내용을 살펴보도록 합시다.
• 데이터를 다운 받는 방법은 여러 가지가 있습니다만, SRA Explorer를 이용하면 쉽게 데이터를 다운 받을 수 있습니다. 홈페이지에 들어간 뒤, 검색창에 BioProject 번호인 "PRJNA896647"를 입력해봅시다. 그러면 아까 BioProject 홈페이지에 검색했던 것과 동일한 데이터 2개가 나올 겁니다. 이 2개의 왼쪽에 있는 박스에 체크 표시를 하면 이 데이터를 다운 받을 수 있는 링크를 얻을 수 있습니다. 저희는 이중 "SRR22137523", 즉 ALT1의 데이터만 다운 받도록 하겠습니다. 해당 파일의 왼쪽에 있는 박스를 클릭하고, 검색 결과 창 오른쪽 위에 있는 "Add 1 to collection"을 눌러줍시다. 그러면 화면 맨 오른쪽 위에 있는 장바구니가 하나 채워지고, "1 saved datasets"라는 문자로 바뀔 겁니다. 이제 이 장바구니를 클릭해줍시다.
• 그러면 이제 데이터를 다운 받을 수 있는 다양한 링크가 나올 겁니다. 저는 주로 "Raw FastQ Download URLs"을 클릭한 뒤, 거기 있는 링크를 복사해서 데이터를 다운 받습니다. 다음과 같이 명령어를 입력해봅시다.

(basicGenomics) 어쩌구@저쩌구:~$ mkdir 07_assembly && cd 07_assembly
(basicGenomics) 어쩌구@저쩌구:~/07_assembly$ wget -O ALT1.hifi.fq.gz 복사한링크

• 위와 같이 입력하면, 복사한 링크에 담겨있는 데이터가 ALT1.hifi.fq.gz라는 이름으로 저장될 겁니다. -O ALT1.hifi.fq.gz를 빼면 기본 파일 이름인 SRR22137523_subreads.fastq.gz로 저장될 겁니다. 그러면 이제 ALT1 샘플의 HiFi 데이터가 다운로드될 겁니다.
• 이 HiFi 리드 데이터를 활용해서 실제 유전체 조립을 실행해봅시다. 이번 단계부터는 이미 다양한 프로그램을 잘 활용할 수 있다는 걸 전제로 하고 진행하도록 하겠습니다. 제가 써준 스크립트 그대로 써도 좋긴 하지만, 좀 더 다양한 프로그램을 활용하려면 계속 검색하고 찾아보시면서 쓰는 게 더 좋을 거예요.
• 이제 바로 이번 단계의 과제로 넘어가서 실습을 마저 진행해보도록 합시다.

구조 변이 분석하기

• 이번 단계에서는 구조 변이에 대해 이해해보고 이를 분석하는 프로그램들을 활용해보도록 합시다.
• 구조 변이란 현재로서는 50 bp 이상의 DNA가 바뀐 경우를 가리킵니다. 50 bp로 기준을 정한 이유는 엄청 명확한 이유가 있는 건 아니고요, 숏리드 시퀀싱 기법으로 이 정도 이상은 잘 분석이 안 되기 때문에 50 bp라는 기준을 잡았다고 보시면 됩니다. 그래서 기술 발전하면 구조 변이라는 말의 정의가 바뀔 수도 있어요.
• 가장 대표적인 구조 변이는 결손(deletion), 삽입(insertion), 치환(substitution), 재배열(rearrangement), 전좌(translocation) 등이 있습니다. 이런 구조 변이가 형성될 수 있는 기작은 다양한데요, 여러 리뷰 논문에서 이런 내용을 설명하고 있으니 한번 살펴보셔도 좋겠습니다. 다양한 구조 변이를 포괄적으로 다룬 리뷰 논문도 있고, 상동(homology)이 있는 경우에 작동하는 경우를 다룬 논문도 있습니다.
• 결손은 참조 유전체와 비교했을 때 내가 넣은 쿼리 유전체에 특정한 DNA가 없는 경우를 가리킵니다. 예를 들어 DNA 이중나선이 통째로 잘리는 DNA 손상이 일어났을 때, 노출된 양끝이 바로 붙지 못하면 끝부분이 깎인 다음에 붙을 수가 있습니다. 이러면 원래 중간에 깎여나간 DNA가 사라지니 그만큼 결손 변이가 생기게 됩니다.
• 삽입은 참조 유전체와 비교했을 때 내가 넣은 쿼리 유전체에 특정한 DNA가 없는 경우를 가리킵니다. 예를 들어 DNA 이중나선이 통째로 잘린 뒤에 회복을 할 때 비슷한 상동 서열을 지닌 DNA를 복제할 수도 있는데요, 이렇게 상동 서열을 찾는 과정에서 잘못 미끄러지면 근처에 있는 DNA를 잘못 복제할 수 있습니다. 이러면 원래 DNA보다 더 긴 DNA가 끼어들어가는 삽입 변이가 생기게 됩니다. 물론 반대로 참조 유전체에서 결손이 일어나도 쿼리 유전체에 삽입 변이가 있는 것처럼 잡히기도 합니다.
• 치환은 참조 유전체와 비교했을 때 내가 넣은 쿼리 유전체의 DNA가 다른 지역의 DNA로 바뀐 경우를 가리킵니다. 예를 들어 삽입과 비슷한 과정을 거쳐 잘린 DNA를 회복할 때 끝부분이 잘려나간 다음에 다시 새로운 서열이 삽입되면 치환 변이가 발생할 수 있습니다.
• 재배열은 참조 유전체와 비교했을 때 내가 넣은 쿼리 유전체의 DNA가 마구잡이로 뒤섞인 경우를 가리킵니다. 예를 들어 DNA 손상이 회복되는 과정에서 복제를 진행하는데, 쪼끔 복제하다가 다시 다른 지역의 DNA를 복제하고, 다시 다른 지역 복제하는 식으로 다양한 지역의 DNA가 복제되면 재배열이 일어날 수 있습니다. DNA가 여러 번 쪼개진 뒤에 회복되는 과정에서 조립이 잘못 일어나면 재배열이 일어날 수도 있죠.
• 전좌는 참조 유전체와 비교했을 때 내가 넣은 쿼리 유전체 내 한 염색체의 매우 거대한 조각이 다른 염색체와 이어져 있는 경우라고 볼 수 있습니다. DNA가 손상된 뒤 회복되는 과정에서 교차(crossover)가 일어나면 염색체의 거대한 조각이 서로 뒤바뀔 수 있습니다. 또는 텔로미어 지역처럼 DNA 손상이 심각하게 일어날 수 있는 지역이라면, DNA가 아예 잘려나간 뒤 서로 다른 염색체가 융합한 경우에도 전좌처럼 보일 수 있습니다.
• 이외에도 구조 변이는 다양하게 정의할 수 있고, 그와 관련된 기작은 생각 이상으로 다양합니다. 몇몇 기작들은 저마다 독특한 DNA 회복 시그니처를 남기는데요, 이를 활용해 구조 변이가 어떤 기작을 통해 회복됐는지 추정할 수도 있습니다.
• 구조 변이를 확인하는 방법은 여러 가지로 분류할 수 있습니다만, 저는 크게 리드 기반의 구조 변이 분석과 유전체 지도 기반의 구조 변이 분석으로 나눠서 이야기 드리고자 합니다.
• 리드 기반의 구조 변이 분석은 유전체 조립을 거치지 않고, 리드를 참조 유전체에 매핑한 뒤 곧바로 구조 변이를 검출하는 방식입니다. 장점은 데이터가 적을 때 쓰기 좋다는 것입니다. 유전체 조립을 하려면 유전체 크기의 20배 정도 되는 데이터가 필요한데, 리드 기반으로 구조 변이를 분석할 때는 5배 정도만 있어도 됩니다. 단점은 정확도가 떨어지며, 리드 길이가 수만 개 정도밖에 되지 않다 보니 크기가 매우 큰 구조 변이 분석은 어렵다는 것입니다.
• 여기서 정확도는 두 가지 정도의 수준에서 이야기할 수 있습니다. 하나는 구조 변이의 시작과 끝 위치를 정확히 확인할 수 있는가 하는 것입니다. 예를 들면 염색체 1번의 305에서 503까지가 실제 구조 변이인데, 정확도가 떨어지면 320에서 480까지가 구조 변이라고 검출할 수도 있는 것이죠. 다른 하나는 구조 변이의 진위성과 관련된 것입니다. 예를 들면 염색체 2번의 3-73 사이는 분명 구조 변이가 있는데 이런 구조 변이를 확인하지 못한다든지(false negative), 반대로 염색체 3번의 50-105는 분명 구조 변이가 없는데 이 지역에 구조 변이가 있다고 확인한다든지(false positive)하는 문제가 생길 수 있습니다. 모두 중요한 오류이니 기억해두시면 좋겠습니다.
• 유전체 지도 기반 구조 변이 분석은 리드를 곧바로 쓰는 게 아니라, 유전체 조립을 거친 뒤 형성된 컨티그를 활용해 구조 변이를 검출하는 방식입니다. 장단점은 리드 기반과 반대라고 생각하시면 되는데요, 두 가지 정확도가 모두 높고 조립만 잘 된다면 변이 크기에 관계 없이 모든 변이를 검출할 수 있다는 게 장점입니다. 대신 더 많은 데이터를 생산해야 하다 보니 비용이 더 많이 필요하다는 게 단점이죠. 가격이 점차 싸지다 보면 이런 단점은 금방 무시될 수 있을 것 같습니다.
• 이제 실습을 통해 이런 구조 변이를 직접 검출해보도록 합시다. 이번 과제를 따라하시면 되겠습니다.