여기에서는 주로 QIIME 1.9.1을 이용한 16S rRNA 서열 분석 방법에 대하여 다루고 있다. 시퀀싱 데이터 타입에 대해서는 QIIME 2에서 더욱 명확하게 정의하고 있어서 첫 섹션에 소개하였다. 인터넷에 돌아다니는 QIIME 1용 실습 자료는 주로 454 혹은 일루미나의 multiplexed single-ended fastq file(barcode reads가 별도 파일로 주어지는)에 적합하게 만들어진 것 같다.

[QIIME 2 docs] Importing data에서 각 포맷에 대한 샘플 데이터 파일 링크를 제공하고 있다. 그러나 16S rRNA sequence analysis 자체를 위한 샘플 데이터로도 적당한지는 알 수 없다. 파일의 형식을 이해하는 정도로만 활용하는 것이 좋을 것이다.

• “EMP protocol” multiplexed single-end fastq
• “EMP protocol” multiplexed paired-end fastq
• Casava 1.8 single-end demultiplexed fastq
• Casava 1.8 paired-end demultiplexed fastq
• Fastq manifest formats: SingleEndFastqManifestPhred33V2, SingleEndFastqManifestPhred64V2, PairedEndFastqManifestPhred33V2, PairedEndFastqManifestPhred64V2
• FASTA (sequences without quality information)

EMP SE에 해당한다. QIIME 1과 2의 실습 사이트에서 같은 데이터가 쓰였다. 2019년 현재의 시퀀싱 서비스 수준을 생각한다면 더 이상 이런 형태의 데이터가 얻어질 것 같지는 않다. Demultiplexed paired-end sequencing 자료가 제공될 것이 거의 100% 확실하다.

• QIIME 1 Illumina Overview Tutorial
• QIIME 2 docs

EDAMAME 2016 워크샵wiki에서 쓰인 자료. Demultiplexed paired-end fastq에 해당한다.

• GitHub 메인 페이지
• 실제 데이터 다운로드와 처리방법은 여기에서 얻을 수 있다.

• Gut microbiome modulates response to anti–PD-1 immunotherapy in melanoma patients. Gopalakrishnan et al., Science 359, (2018)
• European Nucelotide Archive under accession numbers PRJEB22894 (fecal 16S), PRJEB22874 (oral 16S), PRJEB22895 (murine 16S), and PRJEB22893 (fecal WGS). ⇐ 링크한 URL은 실제 데이터가 보이는 “Experiment”에 해당함.

여기서 다운로드한 시퀀싱 파일은 서열 ID 라인에 샘플 정보를 담고 있지만 하나의 파일에 모든 샘플의 것이 수록된 상태이다. 따라서 ID를 참조하여 demultiplexing을 해야 된다. 그러나 QIIME 스크립트로는 전환할 수가 없으므로 AWK를 이용한 fasta file splitter에서 소개한 방법을 사용하라.

• The commensal microbiome is associated with anti-PD-1 efficacy in metastatic melanoma patients. Matson et al., Science 359, 104–108 (2018) 5 January 2018 논문 링크
• https://www.ncbi.nlm.nih.gov/sra/SRP116709 - 172 SRA experiments
• Custom code and additional processed data used in this study are publicly available on GitHub at https://github.com/cribioinfo/sci2017_analysis.

마우스 데이터는 amplicon이 91건, 사람 데이터는 amplicon과 WGS가 각각 42건과 39건이다. 한 건의 experiment는 하나의 샘플에서 유래한 paired read 데이터에 해당한다. Run Selector에서 원하는 데이터의 Accession List를 받은 다음 SRA toolkit으로 한꺼번에 다운로드하는 것을 권장한다.

인터넷에 널린 자료를 바탕으로 검토해 본 결과 이러한 타입의 데이터는 QIIME v1에게는 상당히 성가신 입력물이 아닌가 한다. EDAMAME workshop의 방법을 변형한 내 고유의 방법을 쓰거나 아예 VSEARCH pipeline을 쓰는 것이 가장 바람직한 방법일 것으로 믿는다.

INPUT_DIR 디렉토리에 demultiplexed paired-end fastq 파일이 있다고 가정한다.

QIIME의 pair end 연결용 스크립트는 기본적으로 ea-utils에 포함된 fastq-join을 사용한다. 시퀀싱 업체에서 제공한 fastq 파일을 사용하는 경우에는 --read1_indicator _R1_ --read2_indicator _R2_로 지정한다. 별로 친절한 스크립트가 아니라서 join 성공률이나 진행 과정을 보여주지도 않고, 성공적으로 연결된 서열 파일의 이름을 자동적으로 샘플 이름에 맞게 고쳐주지도 않는다.

$ multiple_join_paired_ends.py -i INPUT_DIR/ --read1_indicator _1.fastq --read2_indicator _2.fastq -o output_jpe
# join에 실패한 read는 삭제함
$ find output_jpe -name fastqjoin.un* -exec rm {} \;
$ cd output_jpe
# 디렉토리 이름을 단순화하기: YUN1242.3_44_L001_R1_001 => YUN1242.3
$ for f in $(find . -type d ! -name .) 
> do
> mv $f $(cut -d'_' -f1 <<< "$f")
> done
$ cd ..

$PWD/output_jpe에 각 샘플에 대한 디렉토리가 생기고 각각의 하위에 fastqjoin.join.fastq 파일이 생긴다.

워크샵에서 소개한 Merged_Reads_Script.sh를 수정한 merge_paired_fastq.sh를 사용한다. 내부적으로는 QIIME 1의 join_paired_ends.py를 사용한다. shell script 내에서 fwd 및 rev sequence file을 정확히 가리키도록 수정해야 한다. 오리지널 Merged_Reads_Script.sh에서는 list.txt 파일이 필요하고, join이 완료된 fasta 파일만을 남기고 fastq 파일은 제거한다.

$ cd INPUT_DIR
$ merge_paired_fastq.sh
$ cd ..

$PWD/INPUT_DIR/merged_dir/에 (sample)_merged.fastq 및 (sample)_merged.fasta 파일이 생긴다.

Torbjørn Rognes의 VSEARCH pipeline을 수정한 vsearch_pipeline_from_PE_fastq.sh를 사용한다. Paired-end read merging, quality filtering, chimera removal 및 OTU clustering이 한꺼번에 수행된다. 실행 속도는 매우 빠르며, 모든 결과물은 INPUT_DIR에 만들어진다. 유익한 정보가 많이 출력되므로 파일로 리다이렉션하여 두기를 권장한다.

$ vsearch_pipeline_from_PE_fastq.sh INPUT_DIR 2>&1 | tee vsearch_log_`date +%Y-%m-%d`.txt

$ multiple_split_libraries_fastq.py -i output_jpe -o output_mslf --include_input_dir_path --remove_filepath_in_name

실행이 완료되면 $PWD/output_mslf/seqs.fna가 만들어진다. 샘플에 따라서는 간혹 “skbio.parse.sequences._exception.FastqParseError: Failed qual conversion for seq id: 2. This may be because you passed an incorrect value for phred_offset.” 에러가 발생한다. split_libraries_fastq.py를 실행할 때에는 --phred_score 33 옵션을 주어 해결했으나 multiple_split_libraries_fastq.py에서는 어떻게 해야 좋을지 모르겠다.

혹은 다음의 트릭을 사용하여 quality filtering이 가능한 split_libraries_fastq.py로 sequence join을 할 수도 있다. -i와 --sample_ids에 제공할 문자열 묶음(콤마로 연결됨)을 만들어내는 것이 관건이다. 맨 마지막의 unset 명령어는 만약 코드를 복사하거나 옮겨적어서 실행하다가 실수가 있어서 다시 실행할 때를 대비하여 변수 sl과 bl을 초기화하는 것이다. 이 명령이 없으면 sl과 bl의 값은 계속 길어진다. split_libraries_fastq.py를 사용할 때 원래 phred encoding을 자동으로 검출하지만 간혹 에러가 발생할 수 있으니 명시적으로 --phred_score 33 옵션을 주는 것이 바람직하다.

$ cd output_jpe
$ for f in $(find . -name *join.fastq)
> do
> f=$(sed -e "s/^\.\///" <<< "$f")
> sl+=",$f"
> b=$(cut -d'/' -f1 <<< "$f")
> bl+=",$b"
> echo "$f ===> Sample ID: $b"
> done
$ sl=$(sed -e "s/^,//" <<< "$sl")
$ bl=$(sed -e "s/^,//" <<< "$bl")
$ split_libraries_fastq.py -i $sl --sample_ids $bl -o ../output_slf_q20 -q 19 --barcode_type 'not-barcoded' --phred_offset 33
$ cd ..
$ unset sl bl
$ count_seq.py -i output_slf_q20/seqs.fna

출력물은 $PWD/output_slf_q20/seqs.fna이다.

Sequence join이 성공적으로 끝났다면 단일 디렉토리인 $PWD/INPUR_DIR/merged_dir/ 아래에 (sample)_merged.fastx 파일이 생겼을 것이다. fasta(.fna) 파일을 사용하려면 metadata mapping file을 적당히 만든 다음 add_qiime_labels.py를 실행한다. merged.fasta 파일의 이름이 mapping file의 InputFastaFileName 컬럼에 있어야 한다. 이 방법에서는 (multiple_)split_libraries_fastq.py을 사용하지 않는다.

$ add_qiime_labels.py -i INPUT_DIR/merged_dir/ -m mapping.txt -c InputFastaFileName -n 1

결과물은 $PWD/INPUT_DIR/mergied_dir/combined_seqs.fna이다.

$PWD/INPUT_DIR/merged_dir/(sample)_merged.fastq를 사용하려면 위에서 소개한 방법을 응용하여 다음과 같이 split_libraries_fastq.py 스크립트를 돌린다.

$ cd INPUT_DIR/merged_dir
$ for f in $(ls *merged.fastq)
> do
> sl+=",$f"
> b=$(cut -d'_' -f1 <<< "$f")
> bl+=",$b"
> echo "$f ===> Sample ID: $b"
> done
$ sl=$(sed -e "s/^,//" <<< "$sl")
$ bl=$(sed -e "s/^,//" <<< "$bl")
$ split_libraries_fastq.py -i $sl --sample_ids $bl -o ../../output_slf_q20 -q 19 --barcode_type 'not-barcoded' --phred_offset 33
$ cd ../..
$ unset sl bl

출력물은 $PWD/output_slf_q20/seqs.fna이다.