Illumina read에서 일부분만을 취하기

일루미나 시퀀싱 결과물을 이용하여 de novo genome assembly를 할 때에는 다다익선(多多益善)이 결코 아니다. 과유불급(過猶不及)이라고나 할까? 대략 100x 이상의 depth(or coverage)는 결코 도움이 되지 않고 해가 됨을 너무나 많은 프로젝트를 통해서 직접 경험하였다. Depth가 너무 높으면 우선 컴퓨터가 힘겨워한다. 많은 메모리를 소요하기 때문이다. 그렇다면 Depth가 증가함에 따라 assembly 성적이 나빠지는 이유는 뭘까? 가장 흔한 원인은 - 좀 더 정확히 말하자면 내가 만져본 데이터에 한해서는 - 오염이었다. 라이브러리 제작 이전 단계에 섞인 non-target genome이 문제를 일으키는 것이다. 오염 수준의 genome은 매우 분량이 적으므로 당연히 자잘한 contig를 양산하게 된다.

2013년도 PLoS One에 small genome assembly를 위한 최적의 sequencing depth를 알아본 논문이 출판되었다. 이런 부류의 논문이 과거에 없었던 것은 아니지만 전부 simulated read에 기반한 반면, 이 논문은 실제 시퀀싱 데이터를 사용했다는 것을 장점으로 내세우고 있다. 이논문의 두번째 페이지 왼쪽 단에는 다음과 같은 단락이 있다.

... read depth of 50X is enough to get a ‘‘good’’ genome assembly and sequencing at a depth greater than 100X does not provide any additional benefits. Additionally, we observed that computational resources required for assembling read data of 20 X–100 X depth is between 6–40 GB, for small genomes of E.coli and S. kudriavzevii.

즉 sequencing depth가 어느 수준 이상을 넘어가면 더 많은 데이터를 투입해도 개선되지 않는다고 하였다. 나는 실제로 결과가 더 나빠지는 사례를 무척 많이 경험했다. 단순한 오염 이외에 어떤 systematic한 error가 문제를 더 악화시킨다는 멋진 결론을 한번 내 보려고 무던히 애를 썼지만 그런 일은 아직까지 벌어지지 않았다.

Fastq read에서 subsample을 취하려면 어떻게 하면 가장 좋을까? 원하는 depth에 해당하는 read(line 수 나누기 4)를 계산하여 head 명령어로 앞부분 일부를 취하면 어떤가? 가장 단순한 방법이지만 손으로 계산을 해야 되고 무엇보다도 read file 전체에서 무작위적으로 read가 골라지지 않을 것임이 당연하다. 오늘의 포스팅에서는 이 문제를 다루어 보자.

먼저 fastq/fasta 파일에서 read의 분량을 세는 가장 편리한 방법부터 소개한다. 나는 khmer 패키지의 readstats.py를 가장 좋아한다. Quality trimming을 거쳐서 read의 길이가 일정하지 않아도 상관이 없다.

CLC Genomic Workbench

NGS Core Tools 메뉴에 "Sample Reads"라는 도구가 있다. 클릭만 해 보면 무엇을 하는 기능인지 알 수 있다. read 수 또는 퍼센트 단위로 잘라내는 것이 가능하다.

MyPro의 scripts 중에서...

Preprocess.py를 쓰면 된다. 사용법은 다음과 같다. Trimming을 겸하는 것이 특징이다.

Preprocess.py -read1 [/path/read1.fastq] -read2 [/path/read2.fastq] -g [genomesize] -d [depth] -l [limit] -r [reduce] -min [min length]
-g genomesize. [default:9000000]
-d Int. depth. [default:100]
-l Float. cutoff value for quality trimming. [default:0.01]
-r Int. step size for reducing target length. [default: 10]
-min Int.  Discard reads below length.(defult 50)
-h Help.

Seqtk

sample 명령을 사용한다. 두 개의 파일에 대해서 짝을 맞춘 read를 뽑아내기 위하여 random seed(-s INT)를 일시적으로 저장한다. 다음은 1만개의 read pair를 추출하는 사례이다.

seqtk sample -s100 read1.fq 10000 > sub1.fq
seqtk sample -s100 read2.fq 10000 > sub2.fq

Oneliner

awk, paste 등의 일반적인 유틸리티를 이용하는 재미난 방법으로 Pierre Lindenbaum이 BioStar 커뮤니티에 제안하였다. Lindenbaum은 BioStar를 설립한 사람이기도 하다(BioStar: An Online Question & Answer Resource for the Bioinformatics Community. PLoS One 2011)

그 얼개를 들여다 보자.

1. 한 쌍의 read 파일을 paste로 붙인다. 이렇게 만들어진 파일의 첫번째 컬럼은 read 1, 두번째 컬럼은 read 2에 해당한다. Field separator로 어떤 문자를 사용했는지 유의하라.
2. shuf 명령으로 무작위로 라인을 추출한다.
3. head 명령으로 원하는 만큼의 라인을 취한다.
4. sed와 awk로 read를 다시 한 쌍의 파일로 나누어 쓴다.

Pipe를 써서 이를 한 줄로 실행하면 다음과 같다. 코멘트가 달려있어서 읽기는 좀 나쁘다. 정말 재미있는 코드 아닌가?

paste f1.fastq f2.fastq |\ #merge the two fastqs
awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t\t");} }' |\ #merge by group of 4 lines
shuf  |\ #shuffle
head |\ #only 10 records
sed 's/\t\t/\n/g' |\ #restore the delimiters
awk '{print $1 > "file1.fastq"; print $2 > "file2.fatsq"}' #split in two files.