혼동을 줄이기 위하여 염기는 b, 바이트는 B로 표시하겠다. 컴퓨터 세상에서는 1 KB = 1024 B(byte)이지만 나머지 세상에서는 1 K = 1000이다. 좀 더 솔직하게 말하자면 컴퓨터 세상에서도 1 K = 1,000, 1 M = 1,000 K로 표시할 때가 있다. 바로 HDD 용량을 말할 때. 그래야 조금이라도 더 크게 보이기 때문이다.
pair end sequencing 결과 파일이 있을때 압축을 하지 않은 상태로 파일의 크기는 얼마인지, 포함된 read와 염기쌍의 수는 얼마인지, 마지막으로 target genome의 추정 크기에 대한 시퀀싱 배수는 얼마인지에 대한 감을 잡아보도록 하자. 우선 오래전에 생산한 HiSeq 101 x 2 결과물을 하나 꺼내어 보자. fastq 파일의 수치를 뽑는 도구는 여러가지가 있겠지만 khmer 패키지의 readstats.py가 요즘 마음에 든다.
$ du -sh *fastq
5.8G BL21TKR_1.fastq
5.8G BL21TKR_2.fastq
파일 크기가 제법 크다.
$ readstats.py BL21TKR_1.fastq BL21TKR_2.fastq
(화면 표시 과정 중간 생략)
... BL21TKR_2.fastq 23400000
... BL21TKR_2.fastq 23500000
... found 2379737962 bps / 23561762 seqs; 101.0 average length -- BL21TKR_2.fastq
---------------
2379737962 bp / 23561762 seqs; 101.0 average length -- BL21TKR_1.fastq
2379737962 bp / 23561762 seqs; 101.0 average length -- BL21TKR_2.fastq
---------------
4759475924 bp / 47123524 seqs; 101.0 average length -- total
파일 하나에 대하여 2.38 Gb를 수록하고 있다. 파일 쌍 전체로는 4.76 Gb, 2천4백만 read pair이다. 이 시퀀싱 결과물은 대장균 BL21에서 유래한 것이니 sequencing coverage는 1,000x 정도가 되겠다.
따라서 단일 fastq 파일에 대하여 5.8 GB(기가바이트)의 크기라면 2.4 Gb(기가베이스페어)라는 뜻이다. 요즘의 미생물 시퀀싱 서비스에서는 보통 1개 샘플에 대해서 1 Gb 생산, 즉 5 Mb 게놈에 대해 200X를 목표로 하고 있으니 지나치게 많은 분량을 시퀀싱한 셈이 된다.
다음으로는 이 데이터를 가지고 k-mer abundance 분포를 그리는 방법에 대해서 알아보겠다. 게놈 서열 자체가 아니라 충분히(과하게!) 오버샘플링이 된 read임을 염두에 두도록 한다. 메타게놈 데이터를 이용하는 것이 더욱 바람직하겠지만, 갖고 있는 데이터는 너무 분량이 커서 분석에 시간이 많이 걸린다. 사용하는 프로그램 패키지는 요즘 집중적으로 공부하고있는 khmer이다. 우선 위에서 다룬 fastq 파일 쌍을 하나의 interleaved file로 바꾸어보자. 예전에는 velvet distribution에 포함된 shuffleSequences_fastq.pl을 주로 사용했었다.
$ interleave-reads.py -o BL21TKR.fastq.pe BL21TKR_1.fastq BL21TKR_2.fastq
역시 이런 종류의 작업은 시간이 많이 걸린다. 솔직하게 말하자면 khmer의 유용성과 개발 철학에는 동의하지만 소요시간이 너무 길다는 것이 불만이다. 물론 비슷한 부류의 다른 프로그램을 전부 테스트해보고 하는 말은 아니지만... unique k-mer의 수는 322919167였다.
$ load-into-counting.py -k 20 -N 4 -x 4e9 BL21.count.ct BL21TKR.fastq.pe
$ abundance-dist.py BL21.count.ct BL21TKR.fastq.pe BL21.count.hist
만들어진 k-mer abundance histogram을 한번 열어보자.
$ head -5 BL21.count.hist
abundance,count,cumulative,cumulative_fraction
0,0,0,0.0
1,268607656,268607656,0.832
2,33381122,301988778,0.935
3,9908255,311897033,0.966
첫번째 컬럼은 abundance이다. 즉 주어진 read에 대해서 단 1번 출현한 길이 20-mer가 무려 2억6860개나 된다는 뜻이다. 5 Mb 게놈을 20-mer로 나누는(1-bp shift) 경우를 상상해 보자. 모든 위치에서 얻어지는 20-mer가 전부 다 다르다고 가정하고(실제 repeat에서 유래하는 것도 있으니 그럴 수는 없다), complementary sequence까지 감안을 해도 천만개 정도에 불과하다. 그런데 단 1번 출현하는 20-mer가 2억하고도 거의 7천만개? 이는 바로 일루미나 특유의 시퀀싱 에러에서 기인하는 것이다. 즉 실제 genome에는 존재하지 않는 서열인 것이다.
자, 그럼 이걸 가지고 gnuplot에서 그림을 그려보자. 필드 구분자가 공백이 아닌 콤마이니 gnuplot 내에서 set datafile separator ','를 실행해야 한다.
gnuplot> set datafile separator ','
gnuplot> plot 'BL21.count.hist' using 1:2 with lines
플롯 창이 열렸지만 아무것도 안보인다. 왜? 극단적인 값들 때문에 그렇다. 각 스케일을 로그로 변환하여 다시 그리자.
gnuplot> set logscale x
gnuplot> replot
바로 왼쪽의 파랑색 상자로 둘러싼 low abundance 영역 내의 read가 바로 error에서 기인한 것이다. abundance ~750 근처에 나타난 피크가 바로 실제 시퀀싱 커버리지가 된다. 그보다 더 빈도가 높게 나타나는 피크는 repeat에서 유래한 것들이다. 만약 다른 피크가 더 있다면 타 샘플로 오염이 일어났음을 증명하는 셈이 되겠다. 만약 시퀀싱 리드가 아니라 완성본 게놈 서열을 가지고 에러가 전혀 없는 시퀀싱 리드를 시뮬레이션하여 k-mer 분포를 조사하면 저 왼쪽 영역의 수치는 나오지 않을 것이다. 내가 가끔 사용한 NGS read simulator ART는 너무 충실하게 데이터를 생성해주는 바람에 error-free data를 만들지는 못하는 것 같다.
다음으로 생각해 볼 것은 히스토그램의 끝부분에는 어떤 데이터가 있느냐 하는 것이다. 이 히스토그램 파일의 길이는 무려 65537 라인이나 된다. abundance가 수천을 넘는 것은 손가락에 꼽을 수준으로 매우 적게 존재한다. 그런데 히스토그램의 마지막 줄을 보니 65535회 존재하는 20-mer가 53종류나 있다. 이건 뭘까? 이건 틀림없이 어댑터 서열에서 유래한 20-mer일 것이다.
그런데 이상의 분석은 kmerspectrumanalyzer로 전부 가능한 것이 아닐까?
댓글 없음:
댓글 쓰기