[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 메인 페이지
• 실제 데이터 다운로드와 처리방법은 여기에서 얻을 수 있다.

• 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건이다. Run Selector에서 원하는 데이터의 Accession List를 받은 다음 SRA toolkit으로 한꺼번에 다운로드하는 것을 권장한다.

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

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

시퀀싱 업체에서 제공한 fastq 파일을 사용하는 경우에는 –read1_indicator _R1_ –read2_indicator _R2_로 지정한다.

$ 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/INPUR_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에서는 어떻게 해야 좋을지 모르겠다.

다음의 방법은 quality filtering이 가능한 split_libraries_fastq.py를 사용하기 위한 트릭을 소개하고 있다. -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 컬럼에 있어야 한다.

$ 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이다.