데이터분석/Bioinformatics

[한주현님강의] DNA 분석 파이프라인(3)

묘걍 2023. 10. 24. 21:26

Genome Analysis Tutorial - 툴 체크와 셋팅

에서 시작!

👉🏻Samtools

나는 첫 번째 영상 보고 따라해서

이 위치에 samtools가 깔려 있다.

그래서 강사님과 달리

이렇게 입력했고

이렇게 결과를 볼 수 있다

이 자체를 실행해보면

다음과 같은 결과가 나온다

 

👉🏻 BWA

BWA 부터는 이 튜토리얼을 따라 설치했기 때문에 강사님과 같은 명령어를 입력하면 됐다

다시 이 자체를 실행해서

이렇게 결과가 나오면 된다

 

👉🏻 GATK4

- jar 파일은 자바로 실행

자바가 잘 깔려 있는지 확인

잘 깔려 있다

- Java버전이 1.8이상이어야한다

자바 버전 체크

 

8보다 높아서 괜찮음

GATK 실행해보기

이렇게 나오면 잘 된 것!

 

👉🏻환경설정

(필수는 아님)

samtools 사용할 때 마다 해당 경로를 써 줄 수는 없으니

- path를 지정하거나

- 환경 설정 변수에 samtools를 집어넣는 방법이 있음

환경 설정 변수에 집어넣기

(영상 설명이 정확하지 않아 내가 이해한대로 한거라서 야매일 수 있음)

를 입력하고 엔터 치면 어떤 내용이 나온다

↓키를 눌러 라인을 이동하며 아래와 같은 내용을 찾는다

강사님과 똑같은 위치에 똑같은 내용을 입력해준다

입력이 안돼서 따로 찾아보니

~/.bashrc 파일은 사용자의 bash 셸 설정이 담긴 파일이다. 여기에 환경 변수를 추가하면 셸을 시작할 때마다 해당 환경 변수가 설정된다.
vi는 텍스트 에디터 중 하나다. 기본적인 키 바인딩으로 텍스트 수정이 가능하다.
vi ~/.bashrc 명령을 사용하면 .bashrc 파일을 편집하기 시작한다. 이제 파일 내용을 수정하기 위한 기본적인 vi 명령들을 알아보자:

입력 모드로 전환하기:
i 키를 누르면 현재 커서 위치에서 텍스트 입력이 가능하다.
a 키를 누르면 현재 커서 위치의 다음에서 텍스트 입력이 가능하다.
o 키를 누르면 현재 커서 아래에 새 줄이 추가되고, 그곳에서 텍스트 입력이 가능하다.

텍스트 저장 및 종료:
입력 모드가 아닐 때, :w 를 입력하면 파일 저장이 된다.
:q 를 입력하면 vi가 종료된다.
수정한 내용을 저장하고 종료하려면 :wq 를 입력한다.
수정 없이 종료하려면 :q! 를 입력한다.

환경 변수 설정을 예로 들면, SAMTOOLS의 경로를 환경 변수로 추가하려면 아래와 같이 입력하면 된다:
export SAMTOOLS=/your/samtools/path​

(여기서 /your/samtools/path 부분은 실제 SAMTOOLS의 경로로 변경하면 된다.)

결론적으로, vi ~/.bashrc를 실행하고 파일의 맨 아래로 이동 (G 키로 맨 아래로 이동) 후 환경 변수를 추가하고, :wq를 입력하여 저장하고 종료하면 된다.

또 한가지!

일단 X해서 나가고 다시 들어가면 될줄 알고 다시 들어왔더니 경고창이 나왔다

그러면서 다음 옵션 중 하나를 골라야 했다

이 알림은 vi나 vim에서 .bashrc 파일을 편집하는 중에 에디터가 제대로 종료되지 않았을 때 생기는 .swp 파일 때문에 나타난다. .swp 파일은 편집 중인 파일의 내용을 임시로 저장하는 용도로 사용되며, 에디터가 비정상적으로 종료될 경우 데이터 손실을 방지하기 위해 존재한다.

이런 경고 메시지가 나타나면, 사용자는 다음과 같은 선택을 할 수 있다:

[O]pen Read-Only: 파일을 읽기 전용 모드로 연다.
(E)dit anyway: .swp 파일을 무시하고 파일 편집을 계속한다.
(R)ecover: .swp 파일에서 변경 내용을 복구한다.
(D)elete it: .swp 파일을 삭제한다.
(Q)uit: 편집을 종료한다.
(A)bort: 작업을 취소한다.

만약 .bashrc 파일 편집을 정상적으로 마치지 않았다면, "(R)ecover" 옵션을 사용해 .swp 파일에서 마지막으로 변경한 내용을 복구할 수 있다. 그 후, :wq 명령으로 파일을 저장하고 종료하면 된다.

하지만 이전의 변경 내용이 중요하지 않다면, "(D)elete it" 옵션을 선택해 .swp 파일을 삭제하고 원래 파일을 다시 편집하면 된다.

가장 안전한 방법은 "(R)ecover" 옵션으로 이전 변경사항을 확인한 후 필요한 부분만 저장하고, .swp 파일은 삭제하는 것이다.

나는 이전에 임시 저장된 파일이 필요 없어서 D를 선택했다

 

다시 해당 위치로 돌아가서

위 설명대로 i를 누르니

INSERT가 나타났다

여기에

여기서 또 한가지!

입력모드에서 :wq 하면 저장및 종료 되지 않고 그냥 ':wq'라는 문자가 입력되는 것이다

저장 및 종료를 원한다면 esc를 눌러 명령 모드로 이동한 다음 :wq를 입력해주어야 저장 및 종료가 된다

 

 

Samtools 인덱스 명령어

이제 인덱스 명령어를 다음과 같이 입력할 수 있다

samtools는 생물정보학 툴로서, 다양한 기능을 제공한다.
index는 그 중 하나의 명령어로, BAM 파일에 대한 인덱스를 생성한다. 이 인덱스는 대량의 시퀀싱 데이터를 효과적으로 검색하고 액세스하기 위해 필요하다.

 

 

Genome Analysis Tutorial - Reference 서열에 Mapping

이미지 출처: 강사님 화면
이미지 출처: 강사님 화면

 

ll을 통해 확인

data 폴더로 들어가기

ll로 확인해보기

그런데 강사님 결과랑 조금 다른 것 같다

나의 결과
강사님 결과

(일단은 진행)

 

gz 파일 들여다보기

gz 파일은 파일 내용을 압축한 형태이다. 이는 파일의 크기를 줄이기 위해 사용되는 기술 중 하나이다.

zless는 Linux와 같은 UNIX 계열 운영체제에서 gz 압축된 텍스트 파일의 내용을 페이지 단위로 볼 수 있게 해주는 명령어이다. 일반적인 less 명령어는 텍스트 파일을 페이지 단위로 읽을 때 사용되는데, zless는 이 기능에 압축된 파일도 읽을 수 있는 기능이 추가된 것이다. 따라서 gz 파일을 압축 해제하지 않고도 그 내용을 확인할 수 있다.

구성 요소

- FASTQ는 네 줄이 하나의 리드

  1. 헤더 (Header):
    • 첫 번째 줄은 헤더다. 보통 '@' 문자로 시작한다.
    • 예시: @SRR1518133.318
    • 이 줄은 각각의 DNA 리드에 대한 고유한 식별자나 설명을 포함하고 있다.
  2. 서열 (Sequence):
    • 두 번째 줄은 DNA 서열을 나타낸다.
    • 예시: CCTAAACTGAGTCCAGCTGGCTAACTCTAAATATATGTGTATCTTTTCAGCATAAAAAAAATAATGTTTTTCATAA
    • 이 줄은 주로 A, T, C, G와 같은 네 가지 문자로 구성된다. 이는 DNA의 네 가지 염기를 나타낸다.
  3. 구분자 (Separator):
    • 세 번째 줄은 '+'로 시작하는 구분자다.
    • 예시: +
    • 이 줄에 추가적인 설명이 포함될 수도 있지만 대부분 의미 없는 구분자로만 사용된다.
  4. 퀄리티 점수 (Quality Score):
    • 네 번째 줄은 각 염기에 대한 퀄리티 점수를 나타낸다.
    • 예시: CCCFFFEDDFDDFHGBGII<EHIJGIIIIJIIEIJ@F?CFFEGIGIJGIEDEGIE@HGIHDCG@CEHEEHFBCFED
    • 이 점수는 ASCII 문자로 표현되며, DNA 시퀀싱 과정에서 해당 염기가 얼마나 확실하게 읽혔는지를 나타낸다. 높은 퀄리티 점수는 해당 위치의 염기가 더 믿을 만하다는 것을 의미한다.

이 네 줄이 반복되면서 FASTQ 파일 내의 모든 리드를 나타낸다.

 

https://github.com/KennethJHan/GenomeAnalysisTutorial#3-reference-%EC%84%9C%EC%97%B4%EC%97%90-mapping

 

GitHub - KennethJHan/GenomeAnalysisTutorial: Genome Analysis Pipeline for Lecture

Genome Analysis Pipeline for Lecture. Contribute to KennethJHan/GenomeAnalysisTutorial development by creating an account on GitHub.

github.com

$BWA2 mem -t 1 -R "@RG\tID:sample\tSM:sample\tPL:platform" ../resource/reference/hg38.chr21.fa sample_1.fastq.gz sample_2.fastq.gz > sample.mapped.sam

붙여 넣기

실패

출처: 강사님 화면

강사님이 보여주신 실패 케이스랑 내 케이스랑 다르다.

BWA 실행파일(?)을 계속 못 찾는 것 같은데.. 아무리 확인해도 제대로 설치되어 있다.

❗환경 변수 PATH 설정을 해줘야 했다!

강사님이 안 알려주신건지, 내가 놓친건지ㅠㅠ

더보기

PATH 설정이 필요한 이유?

  1. 환경 변수:
    • 환경 변수는 운영 체제에서 사용되는 변수다. 이 변수는 특정 값을 저장하고, 필요할 때 그 값을 참조할 수 있게 도와준다. 예로, $BWA2 환경 변수는 bwa-mem2의 위치를 저장하고 있다. 이런 환경 변수를 통해 다른 스크립트나 프로그램에서 $BWA2 값을 사용해 그 경로를 참조할 수 있다.
  2. PATH 환경 변수:
    • PATH는 특별한 환경 변수로, 실행 가능한 파일의 경로를 저장한다. 운영 체제는 사용자가 터미널에서 명령어를 입력하면 PATH에 저장된 디렉토리들을 검색하여 해당 명령어를 찾아 실행한다.
    • 어떤 프로그램을 터미널에서 전체 경로를 지정하지 않고 실행하려면, 그 프로그램이 위치한 디렉토리를 PATH에 추가해야 한다.

환경 변수 설정은 특정 값(예: 경로, 설정 값, 버전)을 저장하고 쉽게 참조하기 위한 것이다. 그러나 PATH 설정은 사용자가 터미널에서 특정 명령어를 입력할 때 그 명령어의 실행 파일을 찾기 위한 경로를 지정하는 것이다. 따라서 $BWA2 환경 변수만 설정하면 bwa-mem2의 위치는 알 수 있지만, 터미널에서 바로 bwa-mem2를 실행하려면 해당 프로그램의 위치를 PATH에 추가해야 한다.

아무튼

vi ~/.bashrc로 들어가서

export 하고 환경변수 설정해준 곳 밑에

export PATH=$PATH:~/tool/bwa-mem2-2.0pre2_x64-linux/

를 입력해주니 강사님 화면과 같은 화면이 나왔다

안되는 이유

파일이 없다

파일을 만들어서 주긴 했는데 사양에 따라 다르게 나올 수 있다 (가상환경을 조그만 컴퓨팅 사양에서 만들어서)

파일 새로 만들어야함

나와서 리소스 디렉토리로 들어간다

ll로 확인해보면

reference가 있다. reference로 들어간다

ll로 확인하면

다음과 같은 명령어를 입력해 index 작업을 해준다.

$BWA2 index hg38.chr21.fa 명령어는 DNA 서열 데이터를 빠르게 검색할 수 있도록 BWA-MEM2라는 도구로 참조 서열 파일인 **hg38.chr21.fa**에 대한 인덱스를 만드는 작업이다.

상세하게 보면:
$BWA2: 환경 변수로서 BWA-MEM2의 위치를 나타낸다. 이 변수는 BWA-MEM2 프로그램을 실행하는 데 사용되는 경로이다.
index: BWA-MEM2의 여러 기능 중 하나로, 참조 서열 데이터 파일을 이용해 빠른 검색을 위한 인덱스를 만드는 기능이다.
hg38.chr21.fa: DNA 참조 서열 파일의 이름이다. 이 파일은 인간의 21번 염색체의 서열 정보를 담고 있다.

이 명령어를 실행하면, BWA-MEM2는 hg38.chr21.fa 파일의 DNA 서열 데이터에 대한 인덱스를 만든다. 이렇게 만들어진 인덱스는 후에 다른 DNA 서열 데이터를 이 참조 서열과 비교할 때 빠른 작업을 위해 사용된다.

간단히 말하자면, 이 명령어는 DNA 참조 서열 파일을 효과적으로 정리하여 후에 DNA 샘플 데이터와의 비교 작업을 빠르게 진행할 수 있게 도와주는 작업이다.

완료되고나서 ll로 확인해보면

다음과 같은 파일이 생성된 것을 볼 수 있다.

cd .. 두 번으로 Genome...까지 나간다

다시 data 디렉토리로 들어가서

$BWA2 mem -t 1 -R "@RG\tID:sample\tSM:sample\tPL:platform" ../resource/reference/hg38.chr21.fa sample_1.fastq.gz sample_2.fastq.gz > sample.mapped.sam

다시 위 명령어를 붙여넣는다

레퍼런스 서열에 리드를 매핑하는 과정이다

→ FASTQ 파일에서 매핑된 BAM파일이 나오게 된다

 

sam과 bam의 차이

- sam: 그냥 텍스트 파일

- bam: 압축하여 binary 형태로 만든 것

- 압축 하는 이유? 지금은 작은 데이터이지만 실제 업무나 연구에서 사용하게 되면 무지무지하게 커져 단순히 sam으로 쓸 수 없고 bam으로 바꾸게 된다.

잘 만들어진 것을 확인할 수 있다.

 

매핑된 파일 확인하기

강사님이 | 뒷 명령어 알려주시기 전에 무작정 명령어를 입력했더니

이런 화면이 우르르르르 내려간다.. 멈추는 방법도 모르겠다

다른분들은 꼭 아래 명령어로 입력하세요

그러면 다음과 같은 결과가 나온다

아까 FASTQ 파일에서 서열 - 퀄리티 있던 내용을 chr 버전에 맞게 매핑이 된 것

 

sam파일을 bam파일로 바꾸기

$SAMTOOLS view -Sb sample.mapped.sam > sample.mapped.bam

ll로 확인해보면

sam파일과 bam파일을 확인할 수 있다.

ll에서 -h 옵션을 주면 파일 사이즈가 나온다

42에서 12로 줄어든 것을 볼 수 있다.

 

 

 

 

 

 

출처: https://youtu.be/CNTr_FNcUcI?si=RnHehlrmU_Ks50qb

https://youtu.be/n7RQS8DB1Qg?si=QSyNQTp3VpXdKTeZ