본 문서에서 소개한 사례에서는 gzip 압축이 된 Illumina paired fastq file(file_1.fastq.gz & file_2.fastq.gz)을 시작점으로 한다. 여기에서 다루는 방법은 어디까지나 나의 경험에 의한 제안일 뿐, 반드시 이를 지켜야 하는 것은 아니다. 생소한 이름의 shell script는 내가 직접 작성하거나, 인터넷에 공유된 것을 수정한 것이다. 이 페이지는 아직 완성되지 않은 상태이다.

맞다! PacBio RSII이 de novo sequencing의 새로운 강자로 부상하고 있고, Oxford Nanopore Technologies의 nanopore sequencing도 점점 많은 사람들의 관심을 끌고 있다. 그럼에도 불구하고 일루미나 시퀀싱 데이터를 이용한 분석 작업을 비용 대비 throughput이 가장 높은 방법으로서 여전히 중요하다고 생각한다. 본인 역시 3GS를 활발하게 사용하고 있으며, 이를 활용하기 위한 소프트웨어를 익힌는 데에도 열심이다.

최대 2 × 150 bp의 데이터를 생산하는 일루미나 HiSeq 계열(2016년 12월 6일 현재)과는 달리 MiSeq 유래의 “long read”, 즉 250 bp 혹은 그 이상의 길이를 갖는 read를 처리하지 못하는 프로그램들이 있다. 어떤 프로그램은 컴파일을 할 때 설정을 바꾸어서 해결 가능하다. 이에 대한 해법을 앞으로 계속 업데이트할 예정이다.

Read length가 길어진다는 것은 k-mer size도 이에 맞게 커져야 함을 의미한다. 이에 따라서 시스템의 메모리도 많이 요구하게 되는 어려운 상황에 직면하게 됨을 알아두자.

Sequencing read file의 갖가지 조작을 하기 위한 방법은 여러 단계에서 구현 가능하다. awk나 sed와 같은 유틸리티를 코맨드라인 한 줄에서 절묘하게 사용하거나(one-liner), 이들을 shell script 수준에서 불러다 사용하는 방법, perl/python 등으로 스크립트를 직접 짜서 쓰는 방법, 파일 조작용 전용 유틸리티(예: FASTX-toolkit, seqtk)를 설치하여 까는 방법 등이 가능하다. 사용자마다 저마다 스타일이 다르므로 어느 한 방법을 사용하라고 권장하기는 어렵다. 이에 대해서는 아래의 별도 페이지에서 좀 더 체계적으로 다루도록 하겠다.

Fastq 파일의 조작법

간혹 하나의 시료에서 다음과 같이 여러개의 파일을 받는 경우가 있을 것이다. 레인이나 인덱스와 관련한 정보는 따로 기록을 해 두면 되고, 파일을 열어보면 서열 헤더 라인에 바코드 시퀀스는 들어있다. 다음의 여섯개 파일을 _R#_를 기준으로 하여 sample_1.fastq와 sample_2.fastq를 만들면 된다.

lane6_index12_sample_CTTGTA_L006_R1_001.fastq
lane6_index12_sample_CTTGTA_L006_R1_002.fastq
lane6_index12_sample_CTTGTA_L006_R1_003.fastq
lane6_index12_sample_CTTGTA_L006_R2_001.fastq
lane6_index12_sample_CTTGTA_L006_R2_002.fastq
lane6_index12_sample_CTTGTA_L006_R2_003.fastq
(변경 방법) $ cat *_R1_* > sample_1.fastq
(변경 방법 계속) $ cat *_R2_* > sample_2.fastq

물론 이는 대단히 단순하게 표현한 것이다. 001, 002…의 개별적인 파일의 순서는 R1 그룹과 R2 그룹 내에서 정확히 맞아야 한다. 그리고 바꾸어야 할 샘플이 여러개라면 손으로 일일이 명령어를 타이핑할 수도 없는 노릇이다. 일반적인 파일 병합 및 이름 변경 일괄 작업을 다음과 같이 할 것을 제안한다. 파일 이름을 실수로 바꾸면 돌이킬 수 없으므로 임시 디렉토리(temp)에서 작업을 한다. awk가 어떻게 사용되고 있는지 눈여겨보라. 여기에서 유의할 것은 awk의 맨 마지막에 나타나는 $3을 정확히 지정하는 것이다. 기계 세팅에 따라서는 1-3_TATAAT_L001_R1_001.fastq과 같은 파일이 생성되는데, 여기에서는 '_'로 구분되는 필드 중 첫번째 것($1)이 샘플의 이름이 된다. 따라서 상황에 맞게 유연하게 아래의 명령어를 적용해야 한다.

$ ls lane*fastq > fastq_file_names
$ mkdir temp
$ cd temp
$ grep _R1_ ../fastq_file_names | awk -F_ '{printf "cat ../%s >> %s_1.fastq\n", $0, $3}' > file_1.sh
$ grep _R2_ ../fastq_file_names | awk -F_ '{printf "cat ../%s >> %s_2.fastq\n", $0, $3}' > file_2.sh
$ sh file_1.sh
$ sh file_2.sh

2020년 5월 현재에는 이렇게 무식한(?) 방법으로 fastq 파일 이름을 정리하지 않는다. 원본의 파일 이름이 밑줄('_')로 구분된 컬럼을 갖고 있고, 그 중에서 만약 첫번째와 네번째 것을 남기고 싶다면 다음과 같이 하면 된다. 네번째 필드는 보통 판독 방향을 나타내므로(R1 or R2), 이어서 _1과 _2로 바꾸거나 혹은 그냥 놔두어도 보통의 프로그램은 다 알아먹는다.

$ ls *fastq | while read f
> do
> cut -d_ -f1,3 $(<<<$f)
> done

여기서 QC라 함은 원본 시퀀싱 데이터 파일(fastq)에 변화를 가하지 않고 단지 quality에 대한 평가를 하는 일을 뜻힌다. 따라서 수치와 도표로 표현되는 보고서가 주된 결과물이다. 매우 널리 쓰이는 QC 도구인 FastQC를 다음과 같이 사용하면 하나의 입력 파일 file_1.fastq(.gz)에 대해서 file_1_fastq.html 및 file_1_fastq.zip이 생성된다.

$ /usr/local/apps/FastQC/fastqc *fastq.gz

Illumina sequencing read의 주요 용도라면 de novo assembly나 reference mapping이 가장 대표적일 것이다. Metagenomics 혹은 epigenomics를 주업으로 삼는 연구자라면 다른 의견을 가질 수도 있겠으나, read 자체를 입력으로 삼는 매우 구체적인 '소프트웨어'를 기준으로 삼는다면 이렇게 말해도 그다지 틀린 것은 아닐 것이다. Raw sequencing read의 전처리(preprocess)는 매우 광범위한 작업을 포괄하고 있다. 예를 들자면 다음과 같다.

• Adapter sequence 제거
• Low quality 영역 제거
• low abundant k-mer를 갖는 read의 제거
• Error 교정(바로 위의 것과 겹치기도 한다)
• 포맷의 전환

후속 프로그램에 따라서는 적극적인 전처리(low quality region의 trimming)이 결과에 좋은 영향을 미치는 것도 있고 그렇지 않은 것도 있다. 따라서 위에서 나열한 작업을 다 해야 한다거나, 혹은 어느 것들을 골라서 해야 한다고 원칙적으로 만들기는 대단히 어렵다. 본 사례에서는 sequencing read에 오염이 존재할 경우(매우 흔하다) 이를 k-mer abundance 측면에서 필터링하여 좋은 de novo assembly를 하는 방법에 초점을 맞추었다.

Sequencing reads(fastq file)의 전처리기는 위에서 나열한 작업 종류만큼 다양한 것들이 존재한다. 어떤 것은 QC plot의 생성을 겸하는 것도 있고(Prinseq), 일반적인 서열 조작 기능에 충실한 것도 있다(FASTX-toolkit). Quality-based trimming이 주요 기능인 sickle, 어댑터 제거에도 중점을 둔 trimmomatic, 복합적인 기능의 HTQC 등등.

이 페이지에서 중점적으로 다루는 전처리 기법에서는 trimmomatic과 khmer 패키지(website 논문 링크)의 interleave-reads.py 스크립트가 필요하다. 보통 read가 부족한 경우는 없으므로 orphan file은 제거하도록 하였다. 결과물은 file.pe.fq이다. 총 bp 수, read 수 및 평균 read 길이와 같은 기본 정보가 필요하면 khmer 패키지에 포함된 readstats.py 스크립트가 매우 유용하다. myIllu_01_trimPE.sh 쉘 스크립트가 결코 완벽한 것은 아니다. 생각 같아서는 가장 처음 작업으로서 read의 맨 앞부분을 5-6 bp 정도 제거하고, 처리 후에는 필터나 교정 기능을 적용하는 것이 더욱 완벽한 방법일지도 모른다.

$ myIllu_01_trimPE.sh file_1.fastq.gz file_2.fastq.gz # file-trim.pe.fq 생성

원래 k-mer 분석에 관심을 갖기 된 것은 khmer를 접하면서부터였다. k-mer abundance에 따라서 read file을 필터링하지 않고 단지 분석만 할 것이라면 jellyfish를 쓰는 것이 더욱 빠르다. k-mer 기본 크기는 21이다. 출력 파일은 file.kmer21.txt와 file.jf.hist이다. myIllu_02_jellyfish.sh

$ myIllu_02_jellyfish.sh file.pe.fq

jellyfish는 압축된 read 파일을 직접 다루지 못하므로 파이프나 리다이렉션을 경유해야 한다. jellyfish 사용자 지침서 2장의 FAQ를 참조하라. 다음으로는 분석 결과 파일을 이용하여 그림을 그려보자. 결과물로는 gnuplot script(.gp)와 PNG 파일이 생성된다. myIllu_04_plotKmer.sh

$ myIllu_04_plotKmer.sh file.jf.hist

SGA에서는 본 프로그램인 assembly뿐만 아니라 correct와 filter 등 매우 다양한 기능을 제공한다. myIllu_06_sga.sh는 실제로는 file.pe.fq를 이용하여 SGA의 전 과정(scaffolding 포함)을 실행하도록 만든 스크립트이다. 전체적인 처리 과정이 매우 복잡하니 스크립트를 꼼꼼히 검토해 보는 것이 바람직하다. 상당한 수준의 전처리를 거치게 되므로, raw fastq file을 그대로 쓰는 것을 권장한다. 여기에 소개한 스크립트에서는 interleaved file을 제공하는 것을 가정하고 있는데, 실행 과정에서 다시 두 개의 paired file(temp_1.fq, temp_2.fq)로 나누게 된다. 인수로 주어지는 파일이 interleaved인지 아닌지를 판별해서 파일 분리 작업을 할지의 여부를 결정하게 만들면 좋겠으나. 차마 거기까지 생각하여 세심하게 스트립트를 수정하지는 못하였다. 결과 파일 중에서 file.ec.fastq + file.ec.filter.fq(실제로는 fastq임)가 SPAdes나 A5-miseq과 같은 다른 프로그램에 쓰일 수 있다.

$ myIllu_06_sga.sh file.pe.fq

두번째 인수는 k-mer이다. 지정하지 않으면 기본 수치인 20을 사용한다. 출력 파일은 file.kmer##.ct와 file.kh.hist이다. ##는 실제 분석에 쓰인 k-mer의 길이이다. jellyfish에서와 마찬가지로 myIllu_04_plotKmer.sh file.kh.hist라 실행하면 k-mer abundance plot을 그릴 수 있다.

$ source /usr/local/apps/khmer/khmerEnv/bin/activate
$ myIllu_03_khmer.sh file.pe.fq 21

k-mer abundance profile을 검토하여 적절한 컷오프 coverage가 결정되면, 다음 스크립트를 실행하여 필터를 적용한다. 컷오프는 10으로 한다. 출력 파일은 file.pe.fq.abundantfilt.gz + file.filt.pe.fq + file.filt.se.fq의 세가지이다. 필터 완료된 파일(file.pe.fq.abundantfilt.gz)에서 paired read와 orphan을 추출하여 나머지 두 개의 파일을 작성하는 것이다.

$ myIllu_05_khmerFilter.sh file.kmer21.ct file.pe.fq 10

전처리를 거쳐서 interleaved 형태로 전환된 read 파일(file.pe.fq)을 대상으로 하는 다양한 assembler의 실행 명령어를 정리하였다. 이 페이지에서는 quality trimming 및 adapter removal을 거친 interleaved fastq file을 사용하는 것을 기준으로 작성하였다. 어떤 프로그램은 interleaved file을 다시 separate file로 나누기도 한다. GAGE recipe에서는 대표적인 일루미나 de novo assembler를 사용하는 command를 소개하고 있으니 참조하기 바란다. 아래에서 소개한 myIllu_06_sga.sh 역시 이것을 거의 그대로 활용한 것이다.

이것은 GUI program이니 file.pe.fq를 임포트하여 paired file로 전환한 뒤 알아서 조립하라.

$ /usr/local/apps/VelvetOptimiser-2.2.5/VelvetOptimiser.pl -s ${KSTART} -e ${KEND} -t ${THREADS} -f '-shortPaired -fastq file.pe.fq'

• [웹사이트] http://www.bcgsc.ca/platform/bioinfo/software/abyss

• [웹사이트] http://cab.spbu.ru/software/spades/ latest version: 3.9.1, v3.9.1 매뉴얼

현재 가장 활발하게 업데이트가 이루어지는 de Bruijn graph assembler라고 여겨진다. PacBio read를 포함한 hybrid assembly, transcript assembly, single cell 유래 데이터의 assembly, plasmid sequence assembly 등 다양한 기능이 있다.

실행 첫 단계에 interleaved 파일을 분리한다. gzip으로 압축한 파일도 처리 가능하다. -t(–threads) 옵션의 기본 수치는 16이다.

$ /usr/local/apps/SPAdes-3.5.0-Linux/bin/spades.py -k 21,33,55 -t 24 --12 file.pe.fq -o spades_out

• output_dir/에 input_1.00.0_0.cor.fastq.gz, input_2.00.0_0.cor.fastq.gz 및 input__unpaired.00.0_0.cor.fastq.gz 파일이 생성된다.
• 오류 교정을 이미 한 read를 사용하려면 –only-assembler 옵션을 추가한다.
• plasmid 서열만을 조립하려면 –plasmid 옵션을 추가한다(plasmidSPAdes: assembling plasmids from whole genome sequencing data Nuc. Acid. Res 2016). plasmidSPAdes는 원래 별도의 프로그램으로 개발된 것이지만 version 3.8부터는 SPAdes에 포함되었다.
• (구식 포맷) –12 interleavedFile or –1 forwardReadFile –2 reverseReadFile and/or –s unpairedReadFile
• (신식 포맷) –s<#> unpairedReadFile –pe<#>-12 interleavedFile –pe<#>-1 forwardReadFile –pe<#>–2 reverseReadFile … <#> = 1,2,…,9; 상세한 사용법은 매뉴얼을 참조할 것.

dipSPAdes는 highly polymorphic diploid genome을 위한 assembler로서 SPAdes에 기반을 하고 있다(J. Comput. Biol. 2015)

• [웹사이트] https://sourceforge.net/projects/ngopt/

A5-miseq은 전처리 과정을 포함하는 파이프라인이다. 따라서 원본 파일 쌍을 이용하는 것을 direct mode, 앞에서 설명한 방식의 전처리를 거친 file.pe.fq를 투입하는 것을 mixed mode라고 명명하였다. 내가 만든 전처리 과정은 실제로 A5-miseq의 것(trimmomatic + SGA)을 거의 그대로 참조한 것이다. 다만 trimmomatic의 파라미터를 조금 수정하고, error correction용의 SGA를 최신 버전으로 바꾼 것에 해당한다. gzip으로 압축한 read 파일도 입력물로 제공 가능하다.

Read length의 상한선은 500nt: A5-miseq(ngopt) 버전 20140604의 a5_pipeline.pl을 열어보면 input read의 최대 길이가 120 nt 이하이면 idba_ud를, 이를 초과하면 idba_ud500 바이너리를 쓰게 되어있다. read length($maxrdlen 변수)가 500 nt를 넘어가면 사용하지 못한다. idba_ud500에 기본 설정된 k-mer 범위는 20..100이다(–mink 20 –maxk 100, ⇐504). 이 버전의 A5-miseq에 포함된 idba는 1.1.1d이다.

$ /usr/local/apps/ngopt/bin/a5_pipeline.pl --threads=24 file_1.fastq.gz file_2.fastq.gz a5-out

교정을 마친 read file을 a5-out.ec.fastq로 명명해 둔다. 원본 파일(file.pe.fq)과 이름을 맞추지 않아도 된다. mixed mode에서 인수로 주어지는 file.pe.fq(염기 교정 전)은 실제로는 scaffolding에만 쓰인다. 따라서 a5-out.ec.fastq를 인수로 주어도 상관이 없을 것이다.

$ /usr/local/apps/ngopt/bin/a5_pipeline.pl --threads=24 --begin=2 --end=5 file.pe.fq a5-out

위에서 소개한 “SGA를 이용한 서열 전처리” 항목을 참조한다.

• [웹사이트] http://i.cs.hku.hk/~alse/hkubrg/projects/idba_ud/
• [GitHub] https://github.com/loneknightpy/idba - IDBA, IDBA-UD, IDBA-Hybrid and IDBA-Tran은 2012년 10월 18일 이후로 전부 한 패키지에 들어있다.

IDBA는 100 bp 근방의 짧은 read를 쓰도록 설게되어 있다. Illumina의 long read(PacBio read를 말하는 것이 아님!)를 사용하려면 src/sequence/short_sequence.h 파일의 kMaxShortSequence(default = 128) 변수를 변경한 뒤 빌드하여라.

• Hybrid assembly
• 시퀀싱 라이브러리 크기의 신속한 체크