2014년 11월 3일 월요일

쉬운 cutadapt 사용법

조금 전에 올린 포스트에서 어댑터 서열 제거와 관련한 배경 지식을 살펴보았다. One Tip Per Day라는 블로그를 보면 어댑터/프라이머 서열 제거를 위한 세가지 방법을 소개하고 있다. Paired-end sequencing에서는 어느 한쪽 read에서 어댑터가 발견되는 경우 반대편 read도 같이 제거해야 한다. Short insert read는 paired end read 중 어느것도 쓸모가 없기 때문이다.
  1. FastqMfc
  2. Trimmomatic
  3. HTSeq
오늘의 포스팅에서는 좀 더 간단한 방법인 cutadapt에 대해서 알아보도록 한다. 독일 도르트문트 대학의 Marcel Martin이 발표한 도구로서 논문 사이트는 여기에 있다.

먼저 pip를 설치하자. 설치 방법은 여기를 참조하라.

pip install cutadapt

일루미나 TruSeq 어댑터가 쓰인 경우 다음과 같이 하면 된다.

1. read 1에서 TruSeq indexed adapter 및 이후의 서열을 제거.

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -o trimmed.1.fastq.gz reads.1.fastq.gz

2. read 2에서 TruSeq universal adapter의 reverse complement를 제거.

cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o trimmed.2.fastq.gz reads.2.fastq.gz

-a ADAPTER는 cutadapt의 가장 간단한(그리고 가장 일반적인 조건에서의) 실행 방법이다. 즉 read의 3'에 어댑터 서열이 나오게 되고 어댑터와 그 이후의 서열을 전부 제거할 때의 실행 모드라고 보면 된다.

혹은 좀 더 간단한 방법을 쓸 수 있다. 두 어댑터 서열이 모두 공통적인 접두어(AGATCGGAAGAGC)를 갖고 있음에 착안하는 것이다.

cutadapt -a AGATCGGAAGAGC -o trimmed.1.fastq.gz reads.1.fastq.gz
cutadapt -a AGATCGGAAGAGC -o trimmed.2.fastq.gz reads.2.fastq.gz

입력 파일이나 출력 파일 전부 read의 수는 같다. 즉, read 전체가 다 trim되어 제거가 되어도 서열 헤더와 길이 '0'짜리 read가 출력 파일에 남아 있다는 뜻이 된다.


Paired-end read의 처리

paired-end read를 처리할 때에는 trim 작업 결과 길이가 어느 기준 이하로 줄어들었다면 그 read를 제거함은 물론 반대편 read(mate)도 없애야 한다. cutadapt에서는 이를 위해서 아직까지는 두 번의 작업을 거쳐야 한다. 실제 작업은 다음과 같이 세 단계를 거친다(홈페이지에 있는 예제를 가만히 들여다 보니 생각보다 어렵다. -p 옵션의 의미를 도대체 모르겠다. 한참을 더 공부하여 겨우 확인하였다. 어쩌면 SolexaQA 혹은 SolexaQA++가 더 이해하기 쉬울지도 모른다. 그러나 SolexaQA는 내가 알기로 어댑터 서열 제거는 하지 못한다).

(1) cutadapt -q 10 -a ADAPTER_FWD --minimum-length 20 -o tmp.1.fastq -p tmp.2.fastq reads.1.fastq reads.2.fastq

자, 조금 어렵다. -p FILE은 무슨 의미인가? cutadapt --help를 해 보면 "Write reads from the paired-end input to FILE."라고 되어 있다. 이게 도대체 무슨 말인지를 모르겠다. 1st read file에 trim이 적용되고 나서 -o로 지정된 파일에 씌여질 것은 당연한데, -p로 지정된 파일에는 도대체 뭐가 기록된단 말인가? tmp.1.fastq는 길이가 제각각이지만 tmp.2.fastq는 그렇지 않다. 그러나 두 파일에 기록된 read의 수는 같다. 그렇다면 이렇게 해석을 할 수 있다. -o로 지정된 tmp.1.fastq에는 실제 어댑터 제거와 트리밍까지 이루어져 살아남은 read만이 기록되고, -p로 지정된 tmp.2.fastq에는 이와 짝을 맞춘 second read가 기록된다. 매우 독특한 작동 방식이다.

(2) cutadapt -q 15 -a ADAPTER_REV --minimum-length 20 -o trimmed.2.fastq -p trimmed.1.fastq tmp.2.fastq tmp.1.fastq

그렇다면 두번째 단계에서 할 일은 명백하다. tmp.2.fastq와 tmp.1.fastq가 순서대로 인풋 파일이 된다. 최종적으로 쓰일 파일은 -o trimmed.1.fastq -p trimmed.1.fastq로 지정된 2개이다. 이제서야 어렴풋하게 이해를 하겠다.

(3-임시 파일 제거) rm tmp.1.fastq tmp.2.fastq

아직 이해가 완벽하게 된 것은 아니다. 예제에서는 왜 first와 second read에 대해서 -q 값을 다르게 적용했는가 하는 점이다. second read의 quality가 first read보다는 더 나쁘므로(per base quality plot을 해 보면 쉽게 알 수 있다), 일부러 더 stringent하게 trim을 유도하는 것인지도 모르겠다.

first 및 second read에 똑같은 옵션(-q 20 --minimum-length 50 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC)을 적용해 보았다. 어댑터 서열을 지나치게 길게 설정한 것이 보기는 좋지 않다. 어찌되었든 15,833,979개의 read 중에서 first file에서는 15,720,585개, second file에서는 15,467,766개가 살아남았다. 따라서 같은 조건을 적용하면 당연히 second file에서 더 적은 read가 살아남는다.


만약 일루미나 TruSeq 방법으로 시퀀싱 라이브러리가 만들어졌다면, NGS raw data 등록 시 어댑터 서열을 공급하지 않아도 대부분 같은 원리가 적용됨을 알 수 있지 않을까? NGS에 사용되는 거의 대부분의 어댑터/프라이머 서열이 여기에 수록되어 있으니 참조하도록 하자.

Trimming과 관련한 최근 논문이 있어 소개한다.

An Extensive Evaluation of Read Trimming Effects on Illumina NGS Data Analysis
오랜만에 효모의 RNA-seq 데이터를 받아들었다. 시료는 총 4개. Sample에 대해 기술한 파일을 열어보니 다음과 같은 어댑터 정보가 있다.

Adapter,AGATCGGAAGAGCACACGTCTGAACTCCAGTCA,,,,,,
AdapterRead2,AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT,,,,,,

자, 이게 뭔 말이냐? 일루미나의 MiSeq Sample Sheet Quick Reference Guide 2012년 11월판 10-11쪽을 보자.

Adapter. Specify the 5' portion of the adapter sequence to prevent reporting sequence beyond the sample DNA.
AdapterRead2. Specify the 5' portion of the Read 2 adapter sequence to prevent eporting sequence beyond the sample DNA

즉 insert(sample)의 길이가 sequencing read보다 짧을 경우, 샘플 DNA 서열을 지나서 나타나게 되는 어댑터의 서열을 read 1과 read 2에 대해서 표시한 것이다. 서열을 살펴보니 매우 낯이 익다. 이 포스팅 윗부분의 cutadapt 사용례에서 보인 것 그대로이다. Read 1에서 sample DNA를 넘어서 나올 수 있는 adapter는 TruSeq indexed adapter(P5)이고, read 2에서는 TruSeq universal adapter(P7)의 reverse complement이다. 즉 이번 데이터는 지극히 보편적인 어댑터를 사용한 것이다.

구글링 과정에서 잘 정리된 문서를 또 하나 발견하였다. 


댓글 없음: