먼저 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개이다. 이제서야 어렴풋하게 이해를 하겠다.
그렇다면 두번째 단계에서 할 일은 명백하다. 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에 사용되는 거의 대부분의 어댑터/프라이머 서열이 여기에 수록되어 있으니 참조하도록 하자.
아직 이해가 완벽하게 된 것은 아니다. 예제에서는 왜 first와 second read에 대해서 -q
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이다. 즉 이번 데이터는 지극히 보편적인 어댑터를 사용한 것이다.
구글링 과정에서 잘 정리된 문서를 또 하나 발견하였다.
Peter N. Robinson. Next-generation sequencing: Methods and computational analysis
댓글 2개:
오랜만에 찾아왔어요. 좋은 정보 얻고 갑니다~^^ (김병용 올림)
감사합니다! 많은 도움이 되었어요!
댓글 쓰기