본문영역 바로가기 하단 바로가기

국내 유일의 클라우드 기반 통합 데이터 분석 서비스

KOBIC 자체 기술로 구축한 빅데이터 플랫폼을 통한 유전체 빅데이터 고속분석

사용방법

바이오익스프레스 웹 서비스 소개

웹 기반 대용량 데이터 분석 서비스

  • 웹 기반 서비스로써 네트워크 가능한 어디서나 이용 가능
  • 공개된 파이프라인을 활용으로 접근성을 높여 초보 연구자의 분석 경험 증가
  • 사용자의 시스템 사용 현황 실시간 확인 가능
  • 사용자간의 커뮤니티 서비스 제공을 통한 정보 공유 및 교류의 장 확대
  • 공개된 프로그램 및 파이프라인 열람을 통해 분석 기초 정보 확인
[Figure 1.1-1]

[Figure 1.1-1]

[Figure 1.1-2]

[Figure 1.1-2]

바이오익스프레스 워크벤치 다운로드

2.1.1 Windows

“ https://www.kobic.re.kr/bioexpress/” 접속한다.

[Figure 2.1.1-1]

[Figure 2.1.1-1]

[Figure 2.1.1-2]

[Figure 2.1.1-2]

메인화면인 [그림 2.1.1-1]의 우측 퀵메뉴를 통하여 Windows용 워크벤치를 선택하여 다운로드한다.

[Figure 2.1.1-3]

[Figure 2.1.1-3]

다운로드된 .zip 파일을 원하는 디렉토리에 압축 해제 후 Bioexpress 를 실행한다.

2.1.2 Linux

[Figure 2.1.2-1]

[Figure 2.1.2-1]

[그림 2.1.2-1]과 같이 wget 명령어로 원하는 디렉토리에 다운로드 한다.

															
$ wget -O Bio-Express.workbench-linux.gtk-x86_64.tar.gz \

"https://www.kobic.re.kr/bioexpress/common/related_file_down?file_name=Bio-Express.workbench-linux.gtk-x86_64.tar.gz"
                                            
[Figure 2.1.2-2]

[Figure 2.1.2-2]

다운로드된 tar.gz 파일을 tar –zxvf 명령어로 압축 파일 해제 후 Bioexpress 를 실행한다.

$ tar -zxvf Bio-Express.workbench-linux.gtk-x86_64.tar.gz
$ cd Bio-Express.workbench-linux.gtk-x86_64
$ ./Bioexpress
                                         

2.1.3 MacOS

“ https://www.kobic.re.kr/bioexpress/” 접속한다.

[Figure 2.1.3-1]

[Figure 2.1.3-1]

[Figure 2.1.3-2]

[Figure 2.1.3-2]

메인화면의 우측 퀵메뉴를 통하여 MacOS용 워크벤치를 선택하여 다운로드한다.

[Figure 2.1.3-3]

[Figure 2.1.3-3]

[Figure 2.1.3-4]

[Figure 2.1.3-4]

다운로드 된 .zip 파일을 압축 해제 후 응용프로그램으로 옮긴 후 Bioexpress를 실행한다.

바이오익스프레스 워크벤치 실행

2.2.1. 바이오익스프레스 워크벤치 로그인

[Figure 2.2.1-1]

[Figure 2.2.1-1]

KOBIC 통합로그인(SSO)에 가입 되어 있는 아이디로 로그인 한다.

아이디가 없으면 [그림 2.2.1-1]의 로그인 화면 왼쪽 하단 Register 버튼 클릭 후 KOBIC 통합로그인(SSO)으로 이동하여 회원가입 후 로그인한다.

2.2.2. 바이오익스프레스 워크벤치 시작화면

[Figure 2.2.2-1]

[Figure 2.2.2-1]

[그림2.2.2-1]과 같이 가운데 영역에 CLOSHA intro 페이지가 활성되고 각 서비스의 카드영역 부분을 누르면 바이오익스프레스 웹 페이지로 연결되어 자세한 소개와 사용법을 익힐 수 있다.

[Figure 2.2.2-2]

[Figure 2.2.2-2]

[그림2.2.2-2]메인 메뉴 [File], [Window], [Help]는 사용자 편의기능과 워크벤치 업데이트가 가능하다.

  • - [File] > [New]: 새로운 워크스페이스, 카테고리, 파이프라인, 프로그램 생성
  • - [File] > [Register]: 파이프라인 등록, 프로그램 등록
  • - [Window] > [Open Perspective]: 클로사 / 지박스 퍼스펙티브 이동
  • - [Window] > [Reset Perspective]: 기본(Default)으로 설정된 퍼스펙티브로 되돌리기
  • - [Window] > [Preference]: 바이오익스프레스 워크벤치 설정 값 세팅
  • - [Help] > [Help Contents]: 바이오익스프레스 사용자 매뉴얼
  • - [Help] > [Community]: 바이오익스프레스 웹 커뮤니티 페이지 이동
  • - [Help] > [Feedback]: 바이오익스프레스 웹 피드백 페이지 이동
  • - [Help] > [About Bio-Express]: 바이오익스프레스 웹 소개 페이지 이동
  • - [Help] > [Update CLOSHA Workbench]: 바이오익스프레스 워크벤치 업데이트
[Figure 2.2.2-3]

[Figure 2.2.2-3]

[그림2.2.2-3]툴바의 가장 우측 메뉴를 통하여 빅데이터 고속 전송 시스템 지박스과 클라우드 기반 통합 분석 워크벤치 클로사의 퍼스펙티브 이동이 가능하다.
툴바를 통해 워크스페이스, 파이프라인, 프로그램, 카테고리의 컨트롤 메뉴를 이용 할 수 있다.

화면 및 기능 소개

GBox 고속 전송 시스템은 가속기 등의 추가적인 하드웨어 구축이나 ActiveX와 같은 별도의 플러그인 설치 없이 소프트웨어 기술만으로 사용자의 대용량 데이터를 고속으로 전송한다. GBox를 통해 업로드 된 사용자 데이터는 CLOSHA에서 사용 가능하다.

[Figure 3.1-1]

[Figure 3.1-1]

GBox는 쉽고 직관적인 사용자 인터페이스로 설계되어 별도의 교육 없이 사용이 가능하다.

빅데이터 고속 전송

[Figure 3.2-1]

[Figure 3.2-1]

업로드 방법 1: [File Explorer] view에서 전송하고자 하는 폴더 또는 파일 선택 하여 [GBox Explorer] view로 Drag & Drop 한다.

[Figure 3.2-2]

[Figure 3.2-2]

업로드 방법 2: [File Explorer] view에서 전송하고자 하는 폴더 또는 파일 선택 하여 컨텍스트 메뉴의 [Upload] 버튼으로 전송한다.

[Figure 3.2-3]

[Figure 3.2-3]

다운로드 방법 1: [GBox Explorer] 에서 전송하고자 하는 폴더 또는 파일을 선택 하여 [File Explorer] view 로 Drag & Drop 한다.

[Figure 3.2-4]

[Figure 3.2-4]

다운로드 방법 2: [GBox Explorer] 에서 전송하고자 하는 폴더 또는 파일을 선택 하여 컨텍스트 메뉴 [Download] 버튼으로 전송한다.

[Figure 3.2-5]

[Figure 3.2-5]

[그림3.2-5]와 같이 파일 전송이 시작되면 어플리케이션 하단의 [Progress] view로 실시간 전송 상황 모니터링이 가능하다.

빅데이터 파일시스템 고속 동기화

[Figure 3.3-1]

[Figure 3.3-1]

GBox → HDFS: [GBox Explorer] view에서 동기화 할 폴더 또는 파일을 선택하여 컨텍스트 메뉴 [Sync] 버튼으로 동기화 한다.

[Figure 3.3-2]

[Figure 3.3-2]

HDFS → GBox : [HDFS Tree Browser] view에서 동기화 할 폴더 또는 파일을 선택하여 컨텍스트 메뉴 [Sync] 버튼으로 동기화 한다.

[Figure 3.3-3]

[Figure 3.3-3]

양방향 파일 동기화 기록을 [Synchronized History] view로 확인할 수 있다.

데이터 파일 관리 기능

[Figure 3.4-1]

[Figure 3.4-1]

[File Explorer] view와 [GBox Explorer] view의 컨텍스트 메뉴로 데이터를 관리 할 수 있다. ([New], [Rename], [Copy], [Paste], [Delete], [Refresh])

[File Explorer] view와 [GBox Explorer] view의 상단의 툴바 메뉴로 디렉토리 이동 및 바로가기 기능을 사용할 수 있다.

화면 및 기능 소개

바이오익스프레스 워크벤치(Bio-Express Workbench)는 대용량 데이터 분석 서비스를 제공하는 클라우드 기반의 프로그램이다. 사용자 친화적인 인터페이스로 구성되어 복잡한 생명정보 분석과정을 손쉽게 파이프라인으로 디자인 할 수 있다.

[Figure 4.1-1]

[Figure 4.1-1]

사용자 정의 기반 프로그램 등록

4.2.1. 사용자 분석 코드 작성

사용자 정의 기반 프로그램을 등록하기 위해 먼저 사용자 분석 코드를 등록해야한다.

[Figure 4.2.1-1]

[Figure 4.2.1-1]

사용자 분석 코드를 등록할 스크립트 폴더를 생성하기 위하여 [Script] View 컨텍스트 메뉴의 [New Folder] 버튼을 클릭한다.

[Figure 4.2.1-2]

[Figure 4.2.1-2]

생성 폴더의 이름을 작성한 후 [OK] 버튼을 클릭하여 폴더를 생성한다.

[Figure 4.2.1-3]

[Figure 4.2.1-3]

[그림 4.2.1-3]과 같이 생성된 스크립트 폴더에서 컨텍스트 메뉴의 [New Script] 버튼을 클릭하여 사용자 분석 코드를 생성한다. 사용자 분석 코드의 파일 형식으로는 Python, Bash, R을 지원한다.

[Figure 4.2.1-4]

[Figure 4.2.1-4]

생성할 스크립트의 이름을 입력한 후 [OK] 버튼을 클릭하여 빈 스크립트를 생성한다.

[Figure 4.2.1-5]

[Figure 4.2.1-5]

#!/bin/bash

INPUT_DIR=$1
OUTPUT_DIR=$2

mkdir -p $OUTPUT_DIR

in_array() {
	local needle array value
	needle="${1}"; shift; array=("${@}")
	for value in ${array[@]};
		do
			[ "${value}" == "${needle}" ] && echo "true" && return;
		done
	echo "false"
}

array=(`find $INPUT_DIR/ -name "*.fastq.gz"`)
for i in "${array[@]}"
do
	FILE=${i##*/}
	IFS='' read -r -a array <<< "$FILE"
	len=(${#array[@]}-1)
	FASTP_IN_FILE_LIST+=("${array[$len]}")
done


for i in "${FASTP_IN_FILE_LIST[@]}"
do
	basename=$i
	FILE1=${basename%1.fastq.gz}

	IFS='' read -r -a file1_array <<< "$FILE1"
	file1_len=(${#file1_array[@]}-1)

	array_check=`in_array ${file1_array[$file1_len]}2.fastq.gz ${FASTP_IN_FILE_LIST[@]}`

	if [ "${array_check}" == "true" ]; then

		IFS='-' read -r -a SEQ_ID <<< "$FILE1"
		 
		FASTP='/BiO/program/fastp/current/fastp'

		time $FASTP -i ${INPUT_DIR}/${file1_array[$file1_len]}1.fastq.gz -I ${INPUT_DIR}/${file1_array[$file1_len]}2.fastq.gz \
		-o ${OUTPUT_DIR}/${file1_array[$file1_len]}1.fastq.gz -O ${OUTPUT_DIR}/${file1_array[$file1_len]}2.fastq.gz \
		--length_required 20 --average_qual 20 --detect_adapter_for_pe --correction \
		-h ${OUTPUT_DIR}/${SEQ_ID[0]}.html -j ${OUTPUT_DIR}/${SEQ_ID[0]}.json
	fi
done
                                            

생성한 사용자 분석 코드는 [그림4.2.1-5]와 같이 [Script] View에서 더블 클릭하여 에디터를 통해 작성 및 편집이 가능하며, 에디터에서 코드 자동 완성 기능을 제공하여 보다 편리한 코드 작성을 지원한다.

4.2.2. 사용자 프로그램 등록

작성한 사용자 분석 코드를 실행하기 위해서 사용자 프로그램을 등록한다.

  • [Figure 4.2.2-1]

    [Figure 4.2.2-1]

  • [Figure 4.2.2-2]

    [Figure 4.2.2-2]

[Program] View에서 컨텍스트 메뉴의 [New Category] 버튼을 클릭하여 사용자 정의 프로그램을 등록할 루트 카테고리와 서브 카테고리를 생성한다.

워크스페이스, 카테고리, 파이프라인, 프로그램 이름은 알파벳, 숫자, 일부특수문자(-_.)만 가능하며 첫 글자는 알파벳으로 지정되어야 한다. 기타 다른 문자나 공백 입력은 불가하다.

루트카데고리에 서브카데고리가 존재하는 경우와 서브카데고리에 프로그램이나 파이프라인이 존재하는 경우 삭제나 수정이 불가능하다.

[Figure 4.2.2-3]

[Figure 4.2.2-3]

컨텍스트 메뉴로 루트 카테고리 하위에 서브 카테고리를 생성한다.

[Figure 4.2.2-4]

[Figure 4.2.2-4]

생성한 서브 카테고리의 컨텍스트 메뉴 [New Program] 버튼을 클릭하여 사용자 프로그램을 등록한다.

[Figure 4.2.2-5]

[Figure 4.2.2-5]

등록할 사용자 프로그램의 기본정보, 추가정보, 인자값을 설정하여 프로그램을 등록한다.

  • 기본 정보 : 이름, 카테고리, 버전, 키워드, 관련 URL, 설명

    키워드는 콤마(,)로 구별하여 여러 개의 키워드를 넣을 수 있다.

  • 추가 정보 : 스크립트 타입, 싱글/멀티 코어, 관련 프로그램
  • 인자값 : Input, Output, Option
    • - Input : 분석에 필요한 파일 또는 경로를 설정
    • - Output : 분석 결과가 출력될 경로 설정
    • - Option : 프로그램에서 지원하는 Option 설정

으로 구성되어 있으며, 설정한 값들은 Default 값이며, 프로그램 실행단계에서 사용자가 새롭게 설정 가능하다.

[Figure 4.2.2-6]

[Figure 4.2.2-6]

인자값 추가 버튼을 클릭하여 fastp 프로그램의 Input 파라미터를 생성한다. 생성된 파라미터는 수정 및 삭제가 가능하다.

[Figure 4.2.2-7]

[Figure 4.2.2-7]

생성할 파라미터의 유형을 지정한 후 Default 값을 입력한다. 파라미터의 유형이 Folder나 File 일 경우 Browser 버튼을 통해 파일 GBox의 파일 및 폴더를 불러올 수 있다.

  • [Figure 4.2.2-8]

    [Figure 4.2.2-8]

    [그림 3.3.2-8]과 같이 생성된 인자값을 확인할 수 있다.

  • [Figure 4.2.2-9]

    [Figure 4.2.2-9]

    Output 및 Option 인자값도 동일한 방법으로 파라미터를 생성할 수 있다.

분석 파이프라인 개발

4.3.1. 워크스페이스 생성

[그림 4.3.1-1]과 같이 [Workspace] View의 컨텍스트 메뉴 [New Workspace] 버튼을 클릭하여 워크스페이스를 생성한다.

[Figure 4.3.1-1]

[Figure 4.3.1-1]

생성할 워크스페이스의 이름과 설명을 작성하고 연관 키워드를 등록한 후 [OK] 버튼을 클릭한다.

[Figure 4.3.1-2]

[Figure 4.3.1-2]

생성한 워크스페이스의 컨텍스트 메뉴 [New Pipeline] 버튼을 클릭하여 분석 파이프라인을 생성한다.

[Figure 4.3.1-3]

[Figure 4.3.1-3]

파이프라인 유형에 따라 사용자 정의 파이프라인을 생성하거나, 기존에 등록되어 있던 파이프라인을 선택하여 인스턴스 파이프라인을 생성할 수도 있다.

[Figure 4.3.1-4]

[Figure 4.3.1-4]

- [Develop a new pipeline] : 사용자 정의 파이프라인 생성

- [Develop an instance pipeline] : 공개된 파이프라인을 선택

[Figure 4.3.1-5]

[Figure 4.3.1-5]

[그림 4.3.1-5]과 같이 [Workspace] View를 통해 생성된 사용자 정의 파이프라인을 확인할 수 있다.

4.3.2. 분석 파이프라인 디자인

4.3.2-1. 사용자 정의 파이프라인 디자인

[Figure 4.3.2-1]

[Figure 4.3.2-1]

분석 파이프라인 디자인에 필요한 분석 프로그램을 에디터의 ‘Palette’에서 찾는다.

[Figure 4.3.2-2]

[Figure 4.3.2-2]

프로그램을 선택한 후 캔버스 위에 마우스로 드래그 앤 드롭 한다.

[Figure 4.3.2-3-1]

[Figure 4.3.2-3-1]

필요한 분석 프로그램으로 디자인한 파이프라인을 저장하고, 각각의 분석 프로그램을 더블 클릭하여 실행에 필요한 옵션을 설정한다.

[Figure 4.3.2-3-2]

[Figure 4.3.2-3-2]

분석 파이프라인 실행 옵션 설정 후, 캔버스에 위치한 분석 프로그램을 연결하여 하나의 파이프라인을 생성한다.

분석 프로그램간의 연결과 옵션 설정의 순서가 바뀌어도 기능상 문제가 없지만, 분석 프로그램을 연결하면 자동으로 앞 단계 프로그램의 Output을 Input으로 설정하게 되므로 프로그램 설정 후 연결하는 방식을 추천하다.

[Figure 4.3.2-4]

[Figure 4.3.2-4]

분석 프로그램의 Output 파라미터 우측 끝에 위치한 원 모양의 연결 영역을 다음 단계 분석 프로그램의 Input 파라미터 좌측 연결 영역으로 드래그 앤 드롭하여 사용자 정의 파이프라인을 디자인한다.

[Figure 4.3.2-5]

[Figure 4.3.2-5]

연결선의 형태는 [그림 4.3.2-5]과 같이 마우스 드래그를 통하여 사용자가 원하는 대로 자유롭게 수정할 수 있다.

4.3.2-2. 인스턴스 파이프라인 디자인

[그림 4.3.2-6]]과 같이 파이프라인 생성 단계에서 인스턴스 파이프라인 유형을 선택하여 기존에 등록되어 있던 파이프라인을 선택할 수 있다..

[Figure 4.3.2-6]

[Figure 4.3.2-6]

인스턴스 파이프라인 역시 사용자 정의 파이프라인과 같이 각각의 프로그램 별 실행에 필요한 설정을 변경할 수 있다.

4.3.3. Share Analytics Pipeline

[Figure 4.3.3-1]

[Figure 4.3.3-1]

디자인된 분석 파이프라인은 KOBIC SSO에 가입한 회원 사이 공유가 가능하다.

공유하고 싶은 파이프라인의 컨텍스트 메뉴에서 Shared Pipeline을 클릭한다.

[Figure 4.3.3-2]

[Figure 4.3.3-2]

공유하려는 파이프라인의 정보와 공유하려는 멤버를 선택할 수 있는 테이블 확인이 가능하다.

  • [Figure 4.3.3-3]

    [Figure 4.3.3-3]

  • [Figure 4.3.3-4]

    [Figure 4.3.3-4]

  • Share Member 탭으로 이동하여 Add를 누르면 공유 가능한 멤버의 목록이 나온다.
  • Search Member에 검색창에 ID로 검색하면 해당 멤버를 빠르게 찾을 수 있다.
  • 공유된 파이프라인은 바이오익스프레스 웹을 통하여 확인 가능하다.

분석 파이프라인 실행

4.4.1. 분석 파이프라인 실행

에디터 상단의 [Run] 버튼을 클릭하여 디자인한 분석 파이프라인을 실행할 수 있다.

[Figure 4.4.1-1]

[Figure 4.4.1-1]

분석이 실행되면 [Pipeline] View, [Log] View와 에디터를 통해 해당 분석 파이프라인의 진행 상태를 확인할 수 있다. 각 View와 에디터의 [Refresh] 버튼을 통해 분석 파이프라인 실시간 진행 상태를 확인할 수 있다. 분석 파이프라인 실행 시 순차적으로 분석 프로그램이 진행된다.

[Figure 4.4.1-2-1]

[Figure 4.4.1-2-1]

[Figure 4.4.1-2-2]

[Figure 4.4.1-2-2]

실행 중인 프로그램은 run 상태이며 에디터에서 녹색으로 표시되고, 실행을 마친 complete 상태의 프로그램은 파란색으로 표시된다.

4.4.2. 단일 분석 프로그램 실행

[Figure 4.4.2-1]

[Figure 4.4.2-1]

각 분석 프로그램의 단일 실행도 가능하다. 캔버스에서 단일 실행할 분석 프로그램에 마우스 오버 시 나타나는 바로 가기 메뉴에서 [Start Node]버튼을 클릭한다.

4.4.3. 데이터 분석 결과 확인

[Figure 4.4.3-1]

[Figure 4.4.3-1]

[Log] View를 통해 각 분석 프로그램의 로그 정보를 확인한다.

[Figure 4.4.3-2]

[Figure 4.4.3-2]

해당 로그를 더블 클릭하면 [Console] View에서 분석 과정 및 오류 등의 로그를 확인할 수 있다.

[Figure 4.4.3-3]

[Figure 4.4.3-3]

[GBoX Browse] View 를 통해 각 분석 프로그램 별로 지정한 Output 경로로 이동하여 결과를 확인할 수 있다.

[Figure 4.4.3-4]

[Figure 4.4.3-4]

[그림 4.4.3-4]는 Bioinformatics 카테고리의 공개 프로그램 fastp의 분석 결과이다.

공개 서비스 등록

4.5.1. 파이프라인 공개 등록

특정 사용자와의 분석 파이프라인 공유를 위하여 분석 파이프라인을 공개 등록이 가능하다.

[Figure 4.5.1-1]

[Figure 4.5.1-1]

[그림 4.5.1-1]과 같이 [Workspace] View에서 원하는 분석 파이프라인을 선택하고, 컨텍스트 메뉴 [Register Pipeline] 버튼을 클릭한다.

[Figure 4.5.1-2]

[Figure 4.5.1-2]

공개 등록이 요청된 분석 파이프라인은 관리자 승인 절차를 거쳐 공개 분석 파이프라인으로 사용이 가능하다.

4.5.2. 프로그램 공개 등록

공개 등록이 요청된 분석 프로그램은 관리자 승인 절차를 거쳐 공개 분석 프로그램으로 등록이 가능하다.

[Figure 4.5.2-1]

[Figure 4.5.2-1]

[그림 4.5.2-1]과 같이 [Program] View에서 원하는 분석 프로그램을 선택하고, 컨텍스트 메뉴 [Register Program] 버튼을 통해 공개 등록할 수 있다.

공개된 카데고리, 프로그램, 파이프라인은 사용자가 삭제나 수정이 불가능하다.

유전체 분석 파이프라인 소개

5.1.1. WGS

파이프라인

[Figure 5.1.1-1]

[Figure 5.1.1-1]

사용 분야

Whole-Genome Sequencing

기능 요약

정상 조직과 암 조직의 DNA Sequence를 입력 받아 Read Map to Reference, Read Merge, CleanUp등의 과정을 거쳐 bam 파일을 만들어주고, GATK4를 이용하여 보다 빠른 SNP & InDel의 somatic variant를 찾아, variant annotate를 수행해 주는 파이프라인.

상세 설명

Fastp, Cutadapt, BWA(mem), SortSam, MarkDuplicates, CountBase, BaseRecalibrator, ApplyBQSR, HaplotypeCaller, somalier 총 10단계로 구성된다. 각 단계에서 분석되는 진행 과정은 다음과 같다. 첫 번째 단계인 fastp의 입력데이터로 pair-end 된 fastq.gz 파일이 존재하는 폴더 경로를 설정한다. fastp는 fastq에 대한 전처리를 수행한다. 또한, Quality control 및 filtering의 결과를 JSON 및 HTML 보고서로 시각화하여 제공한다. 이후 Cutadapt는 fastp에 의해 전처리된 fastq.gz 파일을 입력 받아 high-throughput sequencing reads로 부터 adapter sequences, primers, poly-A tails 및 기타 유형의 unwanted sequence를 찾아 제거한다. 이와 같이 Adaptive Trimming의 결과인 fastq 파일은 다음 분석단계인 BWA(mem)로 전달된다. BWA는 human genome과 같은 large genome reference에 대하여 low-divergent sequence를 mapping한다. BWA(mem)은 genome reference를 기반으로, 70bp~10mbp query sequences를 최대 완전 일치로 seeding한 다음 정렬하여 SAM 파일을 생성한다. 이후 sam파일은 SortSam을 통하여 coordinate, queryname, sam record의 다른 속성별로 정렬하여 BAM 파일로 반환된다. 다음 MarkDuplicates 는 bam 파일 내 중복의 read에 대하여 태그를 지정하고, 이를 CountBase, BaseRecalibrator 그리고 ApplyBQSR로 전달한다. CountBase는 총 base의 수를 count한 출력파일을 나타내어준다. 그리고 BaseRecalibrator는 BQSR(Base Quality Score Recalibration)의 two-stage 중 첫 번째 단계를 수행하는데, read group, reported quality score, machine cycle, nucleotide context 등 다양한 covariates 를 기반으로 재 교정한 table을 생성한다. ApplyBQSR는 BQSR의 두 번째 단계를 수행하며 BaseRecalibrator 에서 재 교정된 table을 기반으로 input reads의 기본품질을 재 교정하고, 재 교정된 BAM 파일을 생성한다. 이후 HaplotypeCaller은 활성영역에 있는 haplotypes의 SNPs 및 InDels를 call 하여 read를 reassemble 한 결과로 VCF 파일을 생성하며, gvcf.gz 형태로 압축하여 다음단계로 전달한다. 마지막 분석단계인 somalier 는 informative sites를 추출하고, sequencing 에 따른 관련성을 평가한 결과를 SOMALIER 파일 형태로 저장한다.

NGS 데이터 분석 파이프라인 설계 과정

5.2.1. NGS 데이터 전처리

전처리 워크플로우

[Figure 5.2.1-1]

[Figure 5.2.1-1]

FASTQ 포맷

[Figure 5.2.1-2]

[Figure 5.2.1-2]

FASTQ 포맷 cont

Phred* quality score Q with base-calling error probability P

Q = -10 log10P

[Figure 5.2.1-3]

[Figure 5.2.1-3]

[Figure 5.2.1-4]

[Figure 5.2.1-4]

FASTA 포맷

  • One or more sequences per file
  • “ > ” denotes beginning of sequence or contig
  • Subsequent lines up to the next “ > ” define sequence
  • Lowercase base denotes repeat masked base
  • Contig ID may have comments delimited by “|”
[Figure 5.2.1-5]

[Figure 5.2.1-5]

SAM 포맷

[Figure 5.2.1-6]

[Figure 5.2.1-6]

SAM 포맷 cont

[Figure 5.2.1-7]

[Figure 5.2.1-7]

[Figure 5.2.1-8]

[Figure 5.2.1-8]

[Figure 5.2.1-9]

[Figure 5.2.1-9]

전처리 단계

  • Quality Check
    • - FastQC is an excellent program for visualizing the overall quality of all reads in a fastq file
    • - Input : FASTQ
    • - Output : QC Report (html)
    • - Software : FastQC ...
[Figure 5.2.1-10]

[Figure 5.2.1-10]

[Figure 5.2.1-11]

[Figure 5.2.1-11]

[Figure 5.2.1-12]

[Figure 5.2.1-12]

Trimming

  • Adapter trimming
    • - May increase mapping rates
    • - Absolutely essential for small RNA
    • - Probably improves de novo assemblies
  • Quality trimming
    • - May increase mapping rates
    • - May also lead to loss of information
  • Input: FASTQ
  • Output: Trimmed FASTQ
  • Software: Sickle, Cutadapt, Fastx Toolkit ...
[Figure 5.2.1-13]

[Figure 5.2.1-13]

Alignment

  • The process of making a sequence by arranging sequences of DNA and RNA
  • Alignment is separated by mapping and assembly
  • Input: FASTQ, Fasta
  • Output: Sam / Bam
  • Software: bowtie, bwa, star, soap ...
[Figure 5.2.1-14]

[Figure 5.2.1-14]

전처리 단계 cont

Alignment – Aligning to a Reference

[Figure 5.2.1-15]

[Figure 5.2.1-15]

[Figure 5.2.1-16]

[Figure 5.2.1-16]

5.2.2. Whole Genome Sequencing Pipeline

파이프라인 요약

  • 파이프라인 모식도
  • [Figure 5.2.2-1]

    [Figure 5.2.2-1]

  • 설명

    정상 조직과 암 조직의 DNA Sequence를 입력 받아 Read Map to Reference, Read Merge, CleanUp등의 과정을 거쳐 bam 파일을 만들어주고, GATK4를 이용하여 보다 빠른 SNV & InDel의 somatic variant를 찾아, variant annotate를 수행해 주는 파이프라인.

  • 사용 프로그램

    Fastp, Cutadapt, BWA, GATK4, Somalier

  • 입력 파일
    • - Fastq 데이터 파일
    • - Reference FASTA file
    • - Reference & Resource VCF filee
  • 최종 결과 파일

    Variant annotate file of Somatic SNP & InDel

등록 파이프라인

[Figure 5.2.2-2]

[Figure 5.2.2-2]

단계별 설명

1) Quality Check - Quality Control & Adaptive Trimming

Fastp

  • - 버전 및 라이선스: 0. 20.1 / MIT License
  • - 설명: Fastq 파일에 대한 all-in-one preprocessing을 수행하여 HTML 기반 보고서를 만들어주는 프로그램.
  • - 입력 인자
    • [input]: Raw data로 pair-ended 된 fastq.gz 파일이 있는 경로 (path of *.fastq.gz)
    • [output]: preprocessing의 결과 보고서가 생성될 경로 (path of *.html, *.json, *.fastq.gz)
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 fastp 프로그램.
    • 입력 값과 출력 값을 폴더로 받고 폴더 내의 fastq.gz 파일 세트를 선별하여 분석.
  •                                                             #!/bin/bash
    
    INPUT_DIR=$1
    OUTPUT_DIR=$2
    
    mkdir -p $OUTPUT_DIR
    
    in_array() {
        local needle array value
        needle="${1}"; shift; array=("${@}")
        for value in ${array[@]};
            do
                [ "${value}" == "${needle}" ] && echo "true" && return;
            done
        echo "false"
    }
    
    array=(`find $INPUT_DIR/ -name "*.fastq.gz"`)
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        FASTP_IN_FILE_LIST+=("${array[$len]}")
    done
    
    
    for i in "${FASTP_IN_FILE_LIST[@]}"
    do
        basename=$i
        FILE1=${basename%1.fastq.gz}
    
        IFS='' read -r -a file1_array <<< "$FILE1"
        file1_len=(${#file1_array[@]}-1)
    
        array_check=`in_array ${file1_array[$file1_len]}2.fastq.gz ${FASTP_IN_FILE_LIST[@]}`
    
        if [ "${array_check}" == "true" ]; then
    
            IFS='-' read -r -a SEQ_ID <<< "$FILE1"
             
            FASTP='/BiO/program/fastp/current/fastp'
    
            time $FASTP -i ${INPUT_DIR}/${file1_array[$file1_len]}1.fastq.gz -I ${INPUT_DIR}/${file1_array[$file1_len]}2.fastq.gz \
    		-o ${OUTPUT_DIR}/${file1_array[$file1_len]}1.fastq.gz -O ${OUTPUT_DIR}/${file1_array[$file1_len]}2.fastq.gz \
    		--length_required 20 --average_qual 20 --detect_adapter_for_pe --correction \
    		-h ${OUTPUT_DIR}/${SEQ_ID[0]}.html -j ${OUTPUT_DIR}/${SEQ_ID[0]}.json
        fi
    done
                                                            
  • - 파이프라인 연결
  • [Figure 5.2.2-3]

    [Figure 5.2.2-3]

2) Trimming - Quality Control & Adaptive Trimming

Cutadapt

  • - 버전 및 라이선스: 3.4 / MIT License
  • - 설명 A program that generates purified fastq files by finding and removing adapter sequences, primers, poly-A tails, and unwanted sequences from high-throughput sequencing reads
  • - 입력 인자
    • [input]: Fastp에 의해 preprocessing된 fastq.gz 파일이 있는 경로
    • [output]: Adaptive Trimming의 결과인 fastq 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 cutadapt 프로그램.
    • 입력 값과 출력 값을 폴더로 받고 폴더 내의 fastq.gz 파일 세트를 선별하여 분석
  •                                                                 
    #!/bin/bash
    
    INPUT_DIR=$1
    OUTPUT_DIR=$2
    
    mkdir -p $OUTPUT_DIR
    
    in_array() {
        local needle array value
        needle="${1}"; shift; array=("${@}")
        for value in ${array[@]};
            do
                [ "${value}" == "${needle}" ] && echo "true" && return;
            done
        echo "false"
    }
    
    array=(`find $INPUT_DIR/ -name "*.fastq.gz"`) 
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        SORT_IN_FILE_LIST+=("${array[$len]}")
    done
    
    index=1
    
    for i in "${CUTADAPT_IN_FILE_LIST[@]}"
    do
        basename=$i
        FILE1=${basename%1.fastq.gz}
    
        IFS='' read -r -a file1_array <<< "$FILE1"
        file1_len=(${#file1_array[@]}-1)
    
        array_check=`in_array ${file1_array[$file1_len]}2.fastq.gz ${CUTADAPT_IN_FILE_LIST[@]}`
    
        if [ "${array_check}" == "true" ]; then
    
            READ_1=$INPUT_DIR/${file1_array[$file1_len]}1.fastq.gz
            READ_2=$INPUT_DIR/${file1_array[$file1_len]}2.fastq.gz
    
            READ_1_FILE_NAME=$(basename $READ_1) 
            READ_2_FILE_NAME=$(basename $READ_2)
    
            IFS='.' read -r -a split_name_1 <<< "$READ_1_FILE_NAME" 
            IFS='.' read -r -a split_name_2 <<< "$READ_2_FILE_NAME" 
    
            READ_1_NAME=${split_name_1[0]} 
            READ_2_NAME=${split_name_2[0]}
    
           	CUTADAPT='/BiO/program/python/current/bin/cutadapt'
           
            time $CUTADAPT -q 30 -j 2 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT --minimum-length 70 --pair-filter=any -o $OUTPUT_DIR/${READ_1_NAME}.fastq -p $OUTPUT_DIR/${READ_2_NAME}.fastq  $READ_1 $READ_2
    
            index=$(($index+1))
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-4]

    [Figure 5.2.2-4]

3) Alignment - Read Map to Reference

BWA mem

  • - 버전 및 라이선스: 0.7.15 / GNU General Public License version 3.0
  • - 설명: 70bp~10mbp query sequences를 최대 완전 일치로 seeding한 다음 정렬하여 sam 파일로 만들어주는 프로그램
  • - 입력 인자
    • [input]: Cutadapt에 의해 Adaptive Trimming 된 fastq 파일이 있는 경로
    • [reference]: BWA index에 의해 index가 생성된 Reference fasta 파일의 index가 생성된 경로
    • [output]: Read Map to Reference 된 sam 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 BWA mem 프로그램.
    • 레퍼런스 값은 폴더를 받고 폴더 내의 fasta파일만 선별하고, 입력 값과 출력 값을 폴더로 받고 폴더 내의 fastaq 파일을 선별하여 분석
  • #!/bin/bash
    
    INPUT_DIR=$1
    REF_DIR=$2
    OUTPUT_DIR=$3 
    
    mkdir -p $OUTPUT_DIR
    
    REF_FILE=(`find $REF_DIR/ -name "*.fasta"`)
    
    in_array() {
        local needle array value
        needle="${1}"; shift; array=("${@}")
        for value in ${array[@]};
            do
                [ "${value}" == "${needle}" ] && echo "true" && return;
            done
        echo "false"
    }
    
    array=(`find $INPUT_DIR/ -name "*.fastq"`) 
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        SORT_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${BWA_IN_FILE_LIST[@]}"
    do
        basename=$i
        FILE1=${basename%1.fastq}
    
        IFS='' read -r -a file1_array <<< "$FILE1"
        file1_len=(${#file1_array[@]}-1)
    
        array_check=`in_array ${file1_array[$file1_len]}2.fastq ${BWA_IN_FILE_LIST[@]}`
    
        if [ "${array_check}" == "true" ]; then
    
            READ_1=$INPUT_DIR/${file1_array[$file1_len]}1.fastq 
            READ_2=$INPUT_DIR/${file1_array[$file1_len]}2.fastq 
    
            INPUT_FILE_NAME=$(basename  $READ_1 .fastq)
    
            IFS='_' read -r -a NAME <<< "$INPUT_FILE_NAME"
            IFS='-' read -r -a FASTQ_2D_ID <<< "${NAME[0]}"
            
            OUTPUT_SAM=$OUTPUT_DIR/${NAME}.sam
    	
            BWA='/BiO/program/bwa/current/bwa'
            $BWA mem \
            -M -R "@RG\tID:${INPUT_DIR}\tPL:Nova6000S4\tLB:Truseq_PCRFREE\tSM:${FASTQ_2D_ID}" -v 1 \
            -t 1 $REF_FILE $READ_1 $READ_2 > $OUTPUT_SAM
            
        fi
    
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-5]

    [Figure 5.2.2-5]

4) Alignment - Read Merge

GATK SortSam

  • - 버전 및 라이선스: 0.7.15 / GATK 4.2.0.0 / MIT License
  • - 설명: sam 또는 bam 파일을 coordinate, query name 또는 sam 레코드의 다른 속성별로 정렬 하는 프로그램
  • - 입력인자
    • [input]: BWA mem에 의해 생성된 sam 파일이 있는 경로
    • [output]: 속성별로 정렬된 bam 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 SortSam 프로그램.
    • 입력 값과 출력 값을 폴더로 받고 폴더 내의 sam 파일을 선별하여 분석
  • #!/bin/bash
    
    INPUT_DIR=$1
    OUTPUT_DIR=$2 
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.sam"`) 
    
    for i in "${array[@]}"
    do
    	FILE=${i##*/}
    	IFS='' read -r -a array <<< "$FILE"
    	len=(${#array[@]}-1)
    	echo "${array[$len]}"
    	SORT_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${SORT_IN_FILE_LIST[@]}" 
    do
        basename=$(basename  $i)
        INPUT_FILE_NAME=$(basename  $i .sam)
        EXTENSION=${basename#*.}
    
        if [ $EXTENSION == 'sam' ]; then
    
            IFS='-' read -r -a FNAME <<< "$INPUT_FILE_NAME"
            IFS='_' read -r -a SEQ_ID <<< "$FNAME"
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.sam
            
        SORTSAM='/BiO/program/gatk/current/gatk SortSam'
    
        time $SORTSAM -I $INPUT_FILE -O $OUTPUT_DIR/${INPUT_FILE_NAME}.bam -SO coordinate
    
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-6]

    [Figure 5.2.2-6]

5) Alignment - CleanUp

GATK MarkDuplicates

  • - 버전 및 라이선스: GATK 4.2.0.0 / MIT License
  • - 설명: 입력 파일의 다섯 가지 주요 위치에 있는 sequence를 비교하여 중복의 read를 찾아 태그를 지정하는 프로그램
  • - 입력인자
    • [input]: Sort 에 의해 생성된 bam 파일이 있는 경로
    • [output]: 중복의 read를 찾아 태그가 지정된 bam 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 MarkDuplicates 프로그램.
    • 입력 값과 출력 값을 폴더로 받고 폴더 내의 bam 파일을 선별하여 분석
  • #!/bin/bash
    INPUT_DIR=$1 
    OUTPUT_DIR=$2
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.bam"`)
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        MARK_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${MARK_IN_FILE_LIST[@]}"
    do
        basename=$(basename  $i)
        INPUT_FILE_NAME=$(basename  $i .bam)
        EXTENSION=${basename#*.}
    
        if [ $EXTENSION == 'bam' ]; then
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.bam
            
        MARKDUPLICATES='/BiO/program/gatk/current/gatk MarkDuplicates'
    
        time $MARKDUPLICATES \
            -I $INPUT_FILE \
            -O $OUTPUT_DIR/${INPUT_FILE_NAME}.bam \
            -M $OUTPUT_DIR/${INPUT_FILE_NAME}.table \
            --REMOVE_SEQUENCING_DUPLICATES false
    
        echo "$MARKDUPLICATES -I $INPUT_FILE -O $OUTPUT_DIR/${INPUT_FILE_NAME}.bam -M $OUTPUT_DIR/${INPUT_FILE_NAME}.table --REMOVE_SEQUENCING_DUPLICATES true"
        
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-7]

    [Figure 5.2.2-7]

6) Alignment - Base Count

GATK CountBases

  • - 버전 및 라이선스: GATK 4.2.0.0 / MIT License
  • - 설명: SAM/BAM/CRAM 파일에 있는 total number of bases를 표준출력으로 count 하는 프로그램
  • - 입력인자
    • [input]: MarkDuplicates 에 의해 생성된 bam 파일이 있는 경로
    • [intervals]: Genomic Interval 파일이 있는 경로
    • [output]: bam 파일의 전체 bases 수가 출력된 count 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 CountBases 프로그램.
    • 레퍼런스값은 파일로 받고, 입력 값과 출력 값은 폴더로 받고 폴더 내의 bam 파일을 선별하여 분석
  • INPUT_DIR=$1
    REF_FILE=$2
    OUTPUT_DIR=$3
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.bam"`)
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        COUNT_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${COUNT_IN_FILE_LIST[@]}"
    do
        basename=$(basename  $i)
        INPUT_FILE_NAME=$(basename  $i .bam)
        EXTENSION=${basename#*.}
    
        if [ $EXTENSION == 'bam' ]; then
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.bam
    
            /samtools/1.6/samtools index $INPUT_FILE
        
        COUNTBASES='/BiO/program/gatk/current/gatk CountBases'
    
        time $COUNTBASES \
            -I $INPUT_FILE \
            -O $OUTPUT_DIR/${INPUT_FILE_NAME}.count \
            -RF NotSecondaryAlignmentReadFilter \
            -L $REF_FILE
        
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-8]

    [Figure 5.2.2-8]

7) Alignment - Base Quality Score Recalibration

GATK BaseRecalibrator

  • - 버전 및 라이선스: GATK 4.2.0.0 / MIT License
  • - 설명: read group, reported quality score, machine cycle, nucleotide context 등의 다양한 covariates 를 기반으로 base quality score 를 재 교정하여 table 을 생성하는 프로그램
  • - 입력인자
    • [input]: MarkDuplicates에 의해 생성된 bam 파일이 있는 경로
    • [reference]: Reference sequence fasta 파일이 있는 경로
    • [known sites1~3]: known polymorphic sites인 vcf 파일이 있는 경로
    • [output]: Base quality score를 재 교정하여 생성된 table 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 BaseRecalibrator프로그램.
    • 레퍼런스값은 파일로 받고, 입력 값과 출력 값은 폴더로 받고 폴더 내의 bam 파일을 선별하여 분석
  • #!/bin/bash
    
    INPUT_DIR=$1
    REF_FILE1=$2
    REF_FILE2=$3
    REF_FILE3=$4
    REF_FILE4=$5
    OUTPUT_DIR=$6
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.bam"`)
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        BASE_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${BASE_IN_FILE_LIST[@]}"
    do
        basename=$(basename  $i)
        INPUT_FILE_NAME=$(basename  $i .bam)
        EXTENSION=${basename#*.}
    
        if [ $EXTENSION == 'bam' ]; then
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.bam
            
        BASERECALIBRATOR='/BiO/program/gatk/current/gatk BaseRecalibrator'
    
        time $BASERECALIBRATOR \
            -I $INPUT_DIR/${INPUT_FILE_NAME}.bam \
            -R ${REF_FILE1} \
            --known-sites ${REF_FILE2} \
            --known-sites ${REF_FILE3} \
            --known-sites ${REF_FILE4} \
            -O $OUTPUT_DIR/${INPUT_FILE_NAME}.table
    
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-9]

    [Figure 5.2.2-9]

8) Alignment - Apply BQSR

GATK ApplyBQSR

  • - 버전 및 라이선스: GATK 4.2.0.0 / MIT License
  • - 설명: BaseRecalibrator 도구로 생성된 재 교정 테이블을 기반으로 입력 read의 기본 품질을 재 교정하여 bam 또는 cram 파일을 생성하는 프로그램
  • - 입력인자
    • [input]: MarkDuplicates 에 의해 생성된 bam 파일이 있는 경로
    • [reference]: Reference sequence fasta 파일이 있는 경로
    • [bqsr recal file]: BaseRecalibrator 에 의해 table 파일이 생성된 경로
    • [output]: BaseRecalibrator 에 의해 생성된 table 파일 기반으로 재 교정된 bam 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 ApplyBQSR프로그램.
    • 레퍼런스값은 파일, 테이블 값은 table파일이 있는 폴더로 받고, 입력 값과 출력 값은 폴더로 받고 폴더 내의 bam 파일을 선별하여 분석
  • #!/bin/bash
    
    INPUT_DIR=$1 
    REF_FILE=$2
    REF_DIR=$3
    OUTPUT_DIR=$4
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.bam"`) 
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        BQSR_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${BQSR_IN_FILE_LIST[@]}"
    do
        basename=$(basename  $i)
        INPUT_FILE_NAME=$(basename  $i .bam)
        EXTENSION=${basename#*.}
    
        if [ $EXTENSION == 'bam' ]; then
     
            IFS='-' read -r -a FNAME <<< "$INPUT_FILE_NAME"
            IFS='_' read -r -a SEQ_ID <<< "$FNAME"
    
            echo $INPUT_FILE_NAME
            echo $FNAME
            
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.bam
    
        APPLYBQSR='/BiO/program/gatk/current/gatk ApplyBQSR'
    
        time $APPLYBQSR \
            -R $REF_FILE \
            -I $INPUT_FILE \
            -bqsr $REF_DIR/${INPUT_FILE_NAME}.table \
            -O $OUTPUT_DIR/${SEQ_ID[0]}_v1.0.0.bam \
            --static-quantized-quals 10 \
            --static-quantized-quals 20 \
            --static-quantized-quals 30 
      
       fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-10]

    [Figure 5.2.2-10]

9) Alignment - Extract&Filter of SNPs, Indels

GATK HaplotypeCaller

  • - 버전 및 라이선스: GATK 4.2.0.0 / MIT License
  • - 설명: active region 에 있는 haplotypes의 SNP와 indels를 call하여 read를 reassemble 한 결과로 vcf 파일을 생성하고 이를 gvcf.gz 로 압축하는 프로그램
  • - 입력인자
    • [input]: ApplyBQSR 에 의해 생성된 bam 파일이 있는 경로
    • [reference]: Reference sequence fasta 파일이 있는 경로
    • [intervals]: Genomic Intervals list 파일이 있는 경로
    • [output]: HaplotypeCaller 를 통해 gvcf.gz 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 HaplotypeCaller프로그램.
    • 레퍼런스값은 파일로 받고, 입력 값과 출력 값은 폴더로 받고 폴더 내의 bam 파일을 선별하여 분석
  • #!/bin/bash
    
    echo "Total Param= $#,  PROG: $0, param1 = $1, param2 = $2, param3 = $3, param4 = $4";
    
    if [ "$#" -lt 4 ]; then
        echo "$# is Illegal number of parameters."
        echo "Usage: $0 input dir, reference file, output dir"
    	exit 1
    fi
    args=("$@")
    
    for (( i=0; i<$#; i++ ))
    do
        echo "$i th parameter = ${args[$i]}";
    done
    
    INPUT_DIR=$1 # BAM 파일
    REF_FILE1=$2
    REF_FILE2=$3
    OUTPUT_DIR=$4 # VCF 파일
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.bam"`)
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        HAPLO_IN_FILE_LIST+=("${array[$len]}")
    done
    
    
    for i in "${HAPLO_IN_FILE_LIST[@]}"
    do
    
            basename=$(basename  $i)
            INPUT_FILE_NAME=${basename%.*}
            EXTENSION=${basename##*.}
    
        if [ $EXTENSION == 'bam' ]; then
    
            IFS='-' read -r -a FNAME <<< "$INPUT_FILE_NAME"
            IFS='_' read -r -a SEQ_ID <<< "$FNAME"
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.bam
    
        HAPLOTYPECALLER='/BiO/program/gatk/current/gatk HaplotypeCaller'
        
        time $HAPLOTYPECALLER -I $INPUT_FILE -R $REF_FILE1 -O $OUTPUT_DIR/${SEQ_ID}_v1.0.0.vcf -OVI true --emit-ref-confidence GVCF -L $REF_FILE2
        echo "$HAPLOTYPECALLER -I $INPUT_FILE -R $REF_FILE1 -O $OUTPUT_DIR/${SEQ_ID}_v1.0.0.vcf -OVI true --emit-ref-confidence GVCF -L $REF_FILE2"
            
        
        mv $OUTPUT_DIR/${SEQ_ID}_v1.0.0.vcf $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf
        echo "mv $OUTPUT_DIR/${SEQ_ID}_v1.0.0.vcf $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf"
    
        BGZIP='/BiO/program/htslib/current/bgzip'
        $BGZIP -c $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf > $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf.gz
        echo "$BGZIP -c $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf > $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf.gz"
        
        TABIX='/BiO/program/htslib/current/tabix'
    	$TABIX -p vcf $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf.gz
        echo "$TABIX -p vcf $OUTPUT_DIR/${SEQ_ID}_v1.0.0.gvcf.gz"
    
    
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-11]

    [Figure 5.2.2-11]

10) Alignment - Variant Annotation

Somalier

  • - 버전 및 라이선스: 0.2.12/ MIT License
  • - 설명: BAM/CRAM/VCF로부터 informative sites를 추출하고, sequencing에 따른 relatedness을 평가해주는 프로그램
  • - 입력인자
    • [input]: HaplotypeCaller를 통해 생성된 gvcf.gz 파일이 있는 경로
    • [reference]: Reference sequence fasta 파일이 있는 경로
    • [sites] : 추출할 variant의 sites인 vcf 파일이 있는 경로
    • [output]: Somalier를 통해 somalier 파일이 생성될 경로
  • - 스크립트 상세
    • 배시 쉘 스크립트로 작성된 Somalier프로그램.
    • 레퍼런스값은 파일로 받고, 입력 값과 출력 값은 폴더로 받고 폴더 내의 gvcf.gz 파일을 선별하여 분석
  • #!/bin/bash
    
    INPUT_DIR=$1
    REF_FILE1=$2
    REF_FILE2=$3
    OUTPUT_DIR=$4
    
    mkdir -p $OUTPUT_DIR
    
    array=(`find $INPUT_DIR/ -name "*.gvcf.gz"`)
    
    for i in "${array[@]}"
    do
        FILE=${i##*/}
        IFS='' read -r -a array <<< "$FILE"
        len=(${#array[@]}-1)
        echo "${array[$len]}"
        SOMAL_IN_FILE_LIST+=("${array[$len]}")
    done
    
    for i in "${SOMAL_IN_FILE_LIST[@]}"
    do
    
        INPUT_FILE_NAME=`basename "$i" | cut -d'.' -f-3`
        EXTENSION=`echo "$i" | cut -d'.' -f4-`
    
        if [ $EXTENSION == 'gvcf.gz' ]; then
    
            INPUT_FILE=$INPUT_DIR/${INPUT_FILE_NAME}.gvcf.gz
    		
            IFS='-' read -r -a FNAME <<< "$INPUT_FILE_NAME"
            IFS='_' read -r -a SEQ_ID <<< "$FNAME"
    
    	SOMALIER='/somalier/0.2.12/somalier'
        
        time $SOMALIER extract -d $OUTPUT_DIR -s ${REF_FILE2} -f ${REF_FILE1} $INPUT_FILE
    
        fi
    done
                                                                
  • - 파이프라인 연결
  • [Figure 5.2.2-12]

    [Figure 5.2.2-12]

결과 화면

  • Fastp 결과 화면
    [Figure 5.2.2-13]

    [Figure 5.2.2-13]

    [Figure 5.2.2-14]

    [Figure 5.2.2-14]

    [Figure 5.2.2-15]

    [Figure 5.2.2-15]

    [Figure 5.2.2-16]

    [Figure 5.2.2-16]

  • GATK_Cutadapt 결과 화면
    @ERR792543.3 HWI-ST333_0255_FC:4:1101:1235:2392/1 CTGCCAATTCCCAATTATATTCTTGGAGAAGCTCATTAGGATAATGGAAATCAATCACTTTGGTTGATCTATCGAAACTT 
    + 
    CBCFFFFFHGHHGIJJIJIIJJJJJJGIIIJJGIGIJIJJJJGIIEIFEIHGIJJJJJJJFGGHIJIEHJGIJGIGHBEC 
    @ERR792543.7 HWI-ST333_0255_FC:4:1101:1228:2438/1 CTCCGGAGTGTAGGCTCCCCATATCTCAGCAAGTGGCCAAAGGTTCGTAATGCCATCTCTGCACNAATCTCCTCCCCCAT 
    + 
    @@@DFFFFCDHF8EHIIIIIIIIIIIIIIIIIGHGHIIIIIII?DDF;BFG<FHHCHEIFGCHH#-5;?DEEEEEBBDDC 
    @ERR792543.8 HWI-ST333_0255_FC:4:1101:1179:2464/1 GCCAGTCATGCCTGCCTGAGACGCCCCGCGGTTGGTGCCCATCTGTAACCCGATCACGTTCTTGNCCTCTTNCAGCTGGT 
    + 
    @CCFFFFFHHHHHJJJJJIJJJIIJJJIIJIDHIJFHIIJEHFHEHHHFFFFDCBDDBBDDDDC#,58??C#+8<<@CCC 
    @ERR792543.10 HWI-ST333_0255_FC:4:1101:1400:2280/1 CACAAACCTTGTTCTAATTTATGTGCACTTCTGAATAAAAATAACACCCTTTTGCTTTAATGCAATATTGATTTTATCAC 
    + 
    @CCDFFDFFHHHHIIGCHIIIIFIGGGGIIHIGIJIIJJIJIIJIJIGGIJIJJJJIJGGDEHGGIJIIHGIJJIHHIEH 
    @ERR792543.12 HWI-ST333_0255_FC:4:1101:1389:2303/1 ATGCGCTCCGCCGCTGTCCTGGCTCTTCTGCTCTGCGCCGGGCAAGTCACTGCGCTCCCTGTGAACAGCCCTATGAATAA 
    + 
    @@CFDFFFFHFD:FGGFHIIIIIIIIIIIGIIGGIIIIIIIFHFFFEACCEECCCC:@B?:@CEDCCCB?@CCCACDEDC
    @ERR792543.14 HWIST333_0255_FC:4:1101:1495:2309/1
    TCTCTCTCCTAATGAAAAAGTTTGTTTCTCTTTTACATTAGTAACGTATGTAACGGGTGGTATGTCTTTCTTTGCTTTGG
    +
    CCCFFFFFHHHDFIJJJJJJJJIJJIJJJJJJIJIGHIJGGIJJJJJJJIFGGIJFG<AH4=EEGGIJHGHHHGFCDFDD
    @ERR792543.15 HWI-ST333_0255_FC:4:1101:1274:2316/1
    GCCCAACTGTCCTGCTGATCTGAAGTAAGACAACTAGCCTTCTACATCAAATAATTCTGTGAGGCAATTTGAGTCAG
    +
    @@@DD?DDAC:3,<CGBFHHHIBCGHFEDEHF;EFEIFII@BFBGGI<DFIIGF@4BF>84<BDEFEEDGIIFFDF@
    @ERR792543.16 HWI-ST333_0255_FC:4:1101:1467:2316/1
    TGCAGTATAGATGGGATATGATGCATAGGCTTGGAGAACCACAGGCAAGGATGAGAGAAGAGAATATGGAAAGGATTG
                                                            
  • BWA_mem 결과 화면
    @SQ	SN:chr3	LN:198295559
    @SQ	SN:chr4	LN:190214555
    @SQ	SN:chr5	LN:181538259
    @SQ	SN:chr6	LN:170805979
    @SQ	SN:chr7	LN:159345973
    @SQ	SN:chr8	LN:145138636
    @SQ	SN:chr9	LN:138394717
    @SQ	SN:chr10	LN:133797422
    @SQ	SN:chr11	LN:135086622
    @SQ	SN:chr12	LN:133275309
    @SQ	SN:chr13	LN:114364328
    @SQ	SN:chr14	LN:107043718
    @SQ	SN:chr15	LN:101991189
    @SQ	SN:chr16	LN:90338345
    @SQ	SN:chr17	LN:83257441
    @SQ	SN:chr18	LN:80373285
    @SQ	SN:chr19	LN:58617616
    @SQ	SN:chr20	LN:64444167
    @SQ	SN:chr21	LN:46709983
    @SQ	SN:chr22	LN:50818468
    @SQ	SN:chrX	LN:156040895
    @SQ	SN:chrY	LN:57227415
    @SQ	SN:chrM	LN:16569
    @SQ	SN:chr1_KI270706v1_random	LN:175055
    @SQ	SN:chr1_KI270707v1_random	LN:32032
    @SQ	SN:chr1_KI270708v1_random	LN:127682
    @SQ	SN:chr1_KI270709v1_random	LN:66860
    @SQ	SN:chr1_KI270710v1_random	LN:40176
    @SQ	SN:chr1_KI270711v1_random	LN:42210
    @SQ	SN:chr1_KI270712v1_random	LN:1760438
                                                            
  • GATK_BaseRecalibrator 결과 화면
    #:GATKTable:Arguments:Recalibration argument collection values used in this run
    Argument                    Value                                                                   
    binary_tag_name             null                                                                    
    covariate                   ReadGroupCovariate,QualityScoreCovariate,ContextCovariate,CycleCovariate
    default_platform            null                                                                    
    deletions_default_quality   45                                                                      
    force_platform              null                                                                    
    indels_context_size         3                                                                       
    insertions_default_quality  45                                                                      
    low_quality_tail            2                                                                       
    maximum_cycle_value         500                                                                     
    mismatches_context_size     2                                                                       
    mismatches_default_quality  -1                                                                      
    no_standard_covs            false                                                                   
    quantizing_levels           16                                                                      
    recalibration_report        null                                                                    
    run_without_dbsnp           false                                                                   
    solid_nocall_strategy       THROW_EXCEPTION                                                         
    solid_recal_mode            SET_Q_ZERO
                                                                  
    #:GATKTable:3:94:%d:%d:%d:;
    #:GATKTable:Quantized:Quality quantization map
    QualityScore  Count     QuantizedScore
               0         0               6
               1         0               6
               2         0               6
               3         0               6
               4         0               6
               5         0               6
               6      3784               6
                                                            
  • GATK4_HaplotypeCaller 결과 화면
    ##fileformat=VCFv4.2
    ##ALT=<ID=NON_REF,Description="Represents any possible alternative allele not already represented at this location by REF and ALT">
    ##FILTER=<ID=LowQual,Description="Low quality">
    ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
    ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
    ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
    ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
    ##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
    ##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
    ##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
    ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
    ##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
    ##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
    ##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller --emit-ref-confidence GVCF --output /bioex/output/haplotype_caller/ERR792543_v1.0.0.vcf --intervals /bioex/reference/references_hg38_v0_wgs_metrics_intervals.interval_list --input /bioex/output/apply_bqsr/ERR792543_v1.0.0.bam --reference /bioex/reference/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta --create-output-variant-index true --use-posteriors-to-calculate-qual false --dont-use-dragstr-priors false --use-new-qual-calculator true --annotate-with-num-discovered-alleles false --heterozygosity 0.001 --indel-heterozygosity 1.25E-4 --
                                                            

somalier 결과 파일

[그림 5.2.2-17]과 같이 somalier는 결과를 [파일명].somalier 파일로 저장한다.

[Figure 5.2.2-17]

[Figure 5.2.2-17]