약 일주일에 걸쳐서 QIIME 1을 집중적으로 공부하면서 느낀 것은 sequencing data 파일로부터 seqs.fna를 만드는 과정이 전체 수고의 1/3 정도는 되는 것 같다는 점이다. 더군다나 아직 실제 16S rRNA 데이터를 생산한 경험이 없어서 SRA에서 자료를 받아서 연습을 하니 그 과정이 더욱 어렵게 느껴졌다. 왜냐하면 서열 ID 라인의 형식이 실제 MiSeq에서 생산된 직후의 fastq 파일 모습을 그대로 따르고 있지 않으며, 메타데이터 역시 제대로 마련되지 않았기 때문이다. 하지만 그 덕분에 SRA-Toolkit의 prefetch와 fastq-dump의 기능을 제대로 이해하게 되었으며, NCBI 웹사이트에서 run accession의 목록 파일을 손쉽게 입수하는 방법을 알게 되는 성과는 있었다.
QIIME 1에서는 Roche/454 pyrosequencing data 또는 barcode read가 별도로 주어지는 Illumina multiplexed single end data가 가장 다루기 편한 것 같다. 바코드 서열과 샘플의 연결 관계는 metadata mapping file(map.tsv)에서 주어진다. 이는 간단히 mapping file이라고도 부른다. 그러나 2019년 9월 초순 현재 국내 시퀀싱 서비스 업체에서는 16S rRNA 시퀀싱 데이터를 이러한 형태로 고객에게 전해주는 것 같지는 않다. 결과 파일인 seqs.fna에서 각 염기서열의 ID는 (샘플이름)_(일련번호)의 형식을 갖는다.
$ split_libraries_fastq.py -o slout/ -i forward_reads.fastq.gz -b barcodes.fastq.gz -m map.tsv
Paired end read라면
join_paired_ends.py 스크립트로 연결한 뒤 추가적으로 바코드 정보를 업데이트해야 하므로 살짝 작업이 번거로워진다. 바코드 업데이트가 필요한 이유는 모든 read pair가 전부 성공적으로 연결되는 것은 아니기 때문이다. 지난주에 천랩에 고객에게 제공되는 16S rRNA amplicon sequencing data의 타입에 대해 문의하니 MiSeq 한 레인 단위의 복수의 샘플을 받아서 시퀀싱을 한 다음 전부 demultiplexing을 해서 고객에게 보내준다고 하였다. 이러한 데이터는 QIIME 1에서 오히려 다루기가 더 불편하다. QIIME 2의 문법으로 말하자면 이는 Casava 1.8 paired-end demultiplexed fastq 파일에 해당한다.
Demultiplexed paired end 데이터를 이용한 QIIME 분석의 좋은 사례는 2016년도
EDAMAME 워크샵(
위키 및
GitHub)에서 제공한 교육 자료에 있다. 이 자료를 만든 사람들 중에는 한국 사람이 틀림없는 이름(Sang-Hoon Lee)이 있어 반갑게 느껴지기도 했다.
여기에서 샘플 파일을 다운로드하여 압축을 풀면 V4 region을 증폭하여 MiSeq으로 양쪽에서 150 bp를 읽은 파일 54쌍이 나온다(C01D01F_sub.fastq, C01D01R_sub.fastq...). QIIME 1 script인 join_paired_ends.py는 한번에 한 쌍의 파일만 처리하므로, 파일이름 베이스를 수록한 list.txt를 만들어서 각 파일에 대하여 순차적으로 join_paired_ends.py을 실행한 다음 연결이 된 sequence만 가져다가 사용하는 스크립트인
Merged_Reads_Script.sh를 사용하는 방법을 소개하였다. 이 스크립트를 실행하면 몇 단계를 거쳐서 Merged_Reads라는 디렉토리에 52개의 fasta 파일이 만들어진다. 다음 단계에서는
add_qiime_labels.py를 사용하여 combined_seqs.fna라는 최종 결과물을 만든다. add_qiime_labels.py는 mapping file을 필요로 한다. 원래 이 파일에는 BarcodeSequence 컬럼이 있어서 각 read에 부가된 바코드로부터 sequence를 구분하는 중요한 역할을 하지만, 이 경우에는 이미 demultiplexing이 되어 있어서 이는 중요하지 않다. 단 InputFileName 컬럼에 있는 fasta file 이름이 각 sequence 파일을 분류하는 기준이 된다. 이 경우에는 필수 컬럼인 BarcodeSequence 컬럼은 공백으로 비워져 있어도 된다.
나는 Merged_Reads_Script.sh를 응용한
merge_paired_fasta.sh 스크립트를 만들어 보았다. 별도의 list.txt 파일이 없이도 쌍을 이룬 fastq 파일이 있으면 작동하며, 최종 결과물로서 fasta + fastq를 같이 만들어내는 것이다. Merged_Reads_Script.sh에서는 merged fastq를 지우도록 만들어 놓았다.
combined_seqs.fna 또는 seqs.fna 파일 내에서 각 서열의 번호는 새로운 샘플에 대해서 무조건 0또는 1부터 시작한다고 생각했는데 그건 착각이었다. 위에서 설명한 방법으로 만들어진 combined_seqs.fna의 서열 일련번호는 _1에서 시작하여(_0으로 시작하게 만들 수도 있음; 전산인들이 좋아하는 방법) 샘플이 바뀌어도 계속 증가한다. 사실 각 read가 구별만 되면 충분한 것이다. 만약 샘플마다 1번부터 시작하게 만들면 일련번호만으로 sequence를 구별하기 어려워지므로 이는 당연한 것이다.
Demultiplexed paird end data를 다루는 방법은 이상의 것이 가장 권장할 만한 것으로 보인다. 이것 외에도 오늘 오후를 종일 투자하여 QIIME 스크립트 자체에 충실한 방법을 만들어내기는 했는데 하도 신경쓸 것이 많아서 과연 이것을 블로그에 기록하는 것이 현명한 일인지 아닌지 잘 판단이 서지 않는다. 간단히 소개하자면 여기에서는
multiple_join_paired_ends.py와
multiple_split_libraries.py를 사용한다. 하나의 디렉토리 안에 샘플별로 이름이 붙은 read pair가 전부 들어있는 경우(아마 대부분의 경우일 것임)와, 샘플별로 구분된 디렉토리 안에 한 쌍씩의 파일이 들어있고 이것 전체가 하나의 인풋 디렉토리에 들어있는 경우에 따라서 옵션이 약간 달라진다. 후자의 경우라면 fwd 및 rev read fastq file이 같은 이름을 갖고 있을 수도 있다. 샘플별로 이름이 붙은 디렉토리에 들어있으니 당장은 문제가 되지 않는다.
EDAMAME 데이터를 예로 들어서 설명해 본다. 54쌍의 read는 전부 샘플 이름을 기준으로 이름이 지어졌으니 하나의 인풋 디렉토리(예: Fastq)에 넣어서 작업을 해도 문제가 되지 않는다. --include_input_dir_path, --include_input_dir_path --remove_filepath_in_name 옵션을 주고 실행한 다음 결과 디렉토리가 어떻게 다르게 만들어지는지 한번 확인해 보라. --remove_filepath_in_name는 단독으로 쓰이지 못하고 항상 --include_input_dir_path와 같이 쓰여야 한다.
$ multiple_join_paired_ends.py -i Fastq --read1_indicator F_sub --read2_indicator R_sub -o output_folder
$ ls output_folder
C01D01F_sub C01D02F_sub C01D03F_sub C02D01F_sub C02D02F_sub log_20190904164417.txt
$ ls output_folder/C01D01F_sub/
fastqjoin.join.fastq fastqjoin.un1.fastq fastqjoin.un2.fastq
각 샘플에 대하여 결과물에 대한 서브디렉토리가 생겼고 그 하위에 연결(join)된 것과 그렇지 않은 것들에 대한 fastq file이 생겼다. 여기서부터 주의를 해야 한다. 연결되지 않은 파일은 일단 제거를 해야 하고, 결과 디렉토리며에서 _sub는 없애는 것이 여러모로 편하다. 이것은 알아서 해결을 했다고 가정하고, 다음 단계인 multple_split_libraries_fastq.py로 넘어가자.
$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final
$ head -n 1 final/seqs.fna
>fastqjoin.join.fastq_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final2 --include_input_dir_path
$ head -n 1 final2/seqs.fna
>C01D01fastqjoin.join.fastq_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final3 --include_input_dir_path --remove_filepath_in_name
$ head -n 1 final3/seqs.fna
>C01D01_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
--include_input_dir_path --remove_filepath_in_name 옵션을 전부 준 세번째 것이 가장 바람직한 서열 ID를 만들어 내었음을 알 수 있다. 이것으로 pick_open_reference_otus.py를 아무런 에러 없이 실행하였다.
split_libraries_fastq.py를 무식하게 실행하여 소기의 목적을 달성할 수도 있다. 다음은 QIIME의
스크립트 설명에 나오는 가장 마지막 예제이다.
$ split_libraries_fastq.py -i lane1_read1.fastq.gz,lane2_read1.fastq.gz --sample_ids my.sample.1,my.sample.2 -o slout_not_multiplexed_q20/ -q 19 --barcode_type 'not-barcoded'
Fastq 형태로 join을 마친 파일을 한 곳에 가져다 놓고 그 이름을 콤마로 연결하여 -i file1,file2,file3...으로, 그리고 순서를 맞춘 샘플 이름을 --sample_ids sample1,sample2,sample3...으로 제공하면 될 것이다. Shell script를 이용하여 이는 어렵지 않게 구현해 보았다. 단, 바코드는 의미가 없으므로 --barcode_type 'not_barcoded' 옵션을 준다. 상당히 번거로왔지만 테스트는 잘 되었다.
마지막으로 가장 무식한 방법을 생각해 보았다. fastq 파일에서 barcode 판독 정보를 뽑아내는 것이다(
multiple_extract_barcodes.py). 그 다음에 이들을 concatenation하여 도로 multiplexed fastq로 전환한 후 split_libraries_fastq.py를 실행하는 것. 원본 파일에 가까운 것을 재구성할 수는 있지만 너무 미련한 방법이다. 그리고 이미 demultiplexing이 된 파일을 입수해 보면 바코드 정보가 fastq 파일 내에 남아있지 않은 경우가 많다. 이론적으로는 안 될 이유가 없지만 시도하지는 말자.
2019년 9월 12일 업데이트
나의 위키 사이트(genoglobe.kr/kribb)에 일루미나 16S rRNA 시퀀싱 데이터를 처리하는 방법을 체계적으로 설명한 글을 올렸다.
