[QIIME2 튜토리얼] “Moving Pictures” (1)

2024.04.18.업데이트

 

 

 

마이크로바이옴을 공부하면 아마 가장 먼저 배우게 되는 것이 이 QIIME2의 사용법입니다. Moving pictures tutorial을 참고하여 각 단계별로 세세하게 알아봅시다.🙉  

 

 

QIIME2  pipeline

 

분석 데이터 관찰하기

- QIIME tutorial 홈페이지: https://docs.qiime2.org/2024.2/tutorials/moving-pictures/

- 관련 영상: https://www.youtube.com/watch?v=RcdTZE8VbJg&list=PLOCEVoX6zu2Ii8RD7i9Oi7Pbot_5WF08n

 

QIIME2의 moving-picture tutorial에서 사용된 데이터는 사람의 마이크로바이옴 데이터입니다. 이 데이터는 항생제 사용에 관하여 두 명의 대상자(Subject)에서 몸에 4개의 부위(장, 혀, 양 손바닥)의 마이크로바이옴을 시간에 따라 관찰한 데이터입니다.(논문 : Caporaso et al., 2011). 본 논문에서 일부 데이터만 추출해서 사용되었습니다. 

이 데이터는 Earth Microbiome Project(EMP) 중 일부로, Illumina Hiseq 시퀀서를 이용하여 16S rRNA의 V4(길이는 약 254 bp) 부분이 사용되었습니다.

 

 

 

먼저, 분석이 수행되는 위치를 정해봅시다.  

mkdir qiime2-moving-pictures-tutorial # 새로운 폴더 생성
cd qiime2-moving-pictures-tutorial    # 새로운 폴더로 이동

 

 

 

1. 샘플의 메타데이터

metadata란 샘플의 정보가 담긴 표를 말합니다. QIIME2에 사용되는 metadata 형식은 아래 요소를 포함해야 합니다(실제 분석 시 없는 내용은 NA 처리해도 됨).

 

위 예제데이터의 metadata는 아래와 같습니다. 각 열의 이름과, 열 데이터의 속성값(numeric, categorical) 1, 2행에 표시되어 있습니다. 

metadata를 살펴보면, 이 연구는 항생제 투약(reported-antibiotic-usage) 이후 0일(투약 당일), 84일, 112일, 140일 뒤의 마이크로바이옴을 관찰하였습니다. 또한 2명의 subject에서 샘플을 추출하였습니다. 

 

튜토리얼에 사용된 메타데이터는 아래 링크에서 볼 수 있습니다(다운로드도 가능!)

-링크:  https://docs.google.com/spreadsheets/d/1fAxItWc_JXVuXEhBiBNFMi3sLwjgiC9wUtIFn6q2W9g/edit#gid=0 

 

 

샘플의 메타데이터 다운로드

wget -O "sample-metadata.tsv" "https://data.qiime2.org/2021.2/tutorials/moving-pictures/sample_metadata.tsv"

 

 

2-1. Obtaining data

분석에 사용되는 서열을 다운로드하여 봅시다. 참고로 single end 서열이 사용되었습니다.  

 

- 바코드 서열은 시퀀싱 시 한 번에 수행된 여러 서열들을 각샘플로 다시 분류할 때 사용됩니다. 이 단계를 demultiplexing이라고 합니다. 

- 실제 회사에 수주를 맡기면, demultiplexing이후 샘플을 받을 수 있습니다. 

mkdir emp-single-end-sequences # 시퀀스를 보관할 장소 만들기
#바코드 서열 다운
wget -O "emp-single-end-sequences/barcodes.fastq.gz" "https://data.qiime2.org/2021.2/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz"

#시퀀스 서열 다운
wget -O "emp-single-end-sequences/sequences.fastq.gz" "https://data.qiime2.org/2021.2/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz"

 

 

 

2-2. Importing data

다운로드한 fastq서열을 QIIME2에서 사용되는 데이터 포맷(qiime artifact : qza, qzv)으로 변환해 줍니다. 

이는 qiime tools import 함수를 사용하여 진행됩니다. 

qiime tools import \
  --type EMPSingleEndSequences \
  --input-path emp-single-end-sequences \
  --output-path emp-single-end-sequences.qza

 

 

3. Demultiplexing sequences

- Demultiplexing을 위해서는, 각 샘플이 어떤 바코드 서열을 가지고 있는지 표시된 metadata를 필요로 합니다. 

- 바코드 서열은 metadata 파일에 존재합니다.

qiime demux emp-single \
  --i-seqs emp-single-end-sequences.qza \  # input file
  --m-barcodes-file sample-metadata.tsv \   # metadata file
  --m-barcodes-column barcode-sequence \    # metadata안에 있는 barcode sequence가 있는 column의 열이름을 입력
  --o-per-sample-sequences demux.qza \    # output file
  --o-error-correction-details demux-details.qza

 

- 위 코드의 실행 결과로 생기는 demux.qza는  demultiplex 된 서열들의 정보를 담고 있습니다.

- demux-details.qza은 Golay error correction details을 보여주며 이 튜토리얼에서 자세하게 다루지는 않습니다..

 

 

 

 

리눅스에서 "head demux.qza" 코드를 입력하면 , 위 그림처럼 깨진 파일이 보일 겁니다.

완성된 demux.qza 파일은 리눅스 환경에서 볼 수 없습니다. 이유는 qiime형식으로 된 압축파일이기 때문입니다.

 

 

즉, 압축을 푼 다음에 볼 수 있습니다.  혹은 qzv로 변환한 후 볼 수 있습니다.

먼저 압축을 풀어서 보는 방법을 알아봅시다. 

qiime tools export \
--input-path [file_name].qza \
--output-path [file_name]

 

 

 

하지만 압축을 풀면 fastq 파일뿐임으로, 요약된 정보를 얻기 위해서 demux.qza파일을 demux.qzv파일로 바꾸어 보자.

qiime demux summarize \
  --i-data demux.qza  \
  --o-visualization demux.qzv  # output file (v:visualization)

 

이후 qzv 파일은  https://view.qiime2.org/에서 확인 가능합니다!

 

 

https://view.qiime2.org/

 

위 사이트에서 demux.qzv을 상세히 살펴봅시다. 

 

 

[Overview]

1) Demultiplexed sequence counts summary

전체 샘플의 요약된 정보를 보여줍니다. 각 샘플에서 read개수의 최소, 최대, 평균, 중앙값을 나타내고 있습니다. 

전체적으로 26만 read가 생성 되었으며, 샘플당 평균적으로 7700개의 read개 생성되었다는 것을 볼 수 있다.

 

2) Reads Frequency Histogram

single-end 서열이기 때문에 forward만 나타나고 있습니다. 

이를 보면 약 2,600개의 read를 가진 샘플이 가장 많고, 그다음으로는 약 100,000개의  read를 가진 샘플이 9개 정도로 보입니다. 

 

위의 표를 막대그래프로 옮겨놓은 것, 이 그래프를 보면 sequence가 2500개 정도인 샘플이 12개, 10000개정도인 샘플이 10개인것을 알 수 있음

 

3) Per-sample sequence counts

각 샘플당 실제적인 read의 개수를 보여줍니다. L4S137 샘플이 18,787개의 read를 가진 것으로 보입니다.

 

각 샘플당 몇개의 sequence가 생성 되었는지를 나타낸다

 

[Interactive Quality Plot]

전체 샘플의 퀄리티를 나타냅니다.  110번째 염기 이후에는 QC가 매우 나빠지는 것을 볼 수 있습니다. 

퀄리티를 위해서는 80~90번대에서 잘라내는 게 좋다는 생각이 들지만, 이는 여러 연구와 연구자의 고민으로 결정해야 합니다.

 

특히 paired-end 서열일 때, 각 forward와 reverse서열을 최소 겹치는 basepair가 12자리가 존재해야 합니다. 이를 반드시 유의하여 잘라내야 합니다. 

- 만약 분석 영역이 V3V4라면, 약 forward와 reverse서열을 겹쳤을 때, 약 460~470bp의 길이를 가져야 합니다.

- 이를 위해 일부 샘플을 사용하여 여러 영역에 따른 퀄리티를 계산하는 것도 좋은 방법입니다(참고)

 

 

4. Sequence quality control and feature table construction

QIIME2에서 사용 가능한 quality control methods로는 DADA2, Deblur, and basic quality-score-based filtering가 있습니다. 위 튜토리얼에서는 DADA2와 Deblur만 다루고 있습니다.

 

전통적인 방법인 Operational Taxonomic Units (OTU)는 각 read를 97% 유사도를 기준으로 묶는 방법입니다. 반대로 Amplicon Sequence Variant(ASV) read에서 에러로 보이는 노이즈를 제거해 가면서 유니크한 서열을 골라내는 것입니다. 

 

마이크로바이옴 연구에서 최근 많이 사용되는 것은 ASV방식입니다. 이는 OTU 방식으로 만들어지는 단위가 ASV방식보다 더 많기 때문입니다. 이는 노이즈를 잘 삭제하지 못하고 하나의 OTU로 만드는 경우가 종종 있기 때문입니다. 

 

4-1) DADA2

실제 분석에서는 Trimmomatic, Cutadapt을 사용하여 primer서열을 제거하고, QC를 진행합니다. 

그러나 이 예제 분석 과정에서는 포함되어있지 않습니다. 

 

물론, DADA2의denoise-single은 quality score filtering과정도 포함되어 있습니다. 

 

qiime dada2 denoise-single \# forward, reverse를 모두 고려하는 것이 아닌 forward만 고려
  --i-demultiplexed-seqs demux.qza \
  --p-trim-left 0 \# QC filtering
  --p-trunc-len 120 \
  --o-representative-sequences rep-seqs-dada2.qza \
  --o-table table-dada2.qza \
  --o-denoising-stats stats-dada2.qza

 

위 예제데이터는 각 서열을 120bp에서 잘라내고, QC filtering은 수행하지 않았습니다.

 

 

이후 분석에서 dada2 결과물을 사용하고 싶다면, 5단계로 넘어가기 전에 이름 바꾸어 줍니다.

mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza

 

 

DADA2 분석 결과 결과도 확인해 봅시다.  

qiime tools export --input-path rep-seqs.qza  --output-path rep-seqs # 압축해제
cd rep-seqs
less -S  dna-sequences.fasta # 일부 확인하기

 

ASV가 fasta형식으로 저장되어 있습니다. 

fasta 데이터는 '>' 다음이 각 feature의 이름이며, 이는 qiime2에서 랜덤으로 부여합니다.

 

 

 

Q. 그럼 Dada2 분석 이후 남은 서열의 개수를 보려면 어떻게 해야 하는가?

A. '>'의 수를 계산하면 구할 수 있다

grep '>' -c aligned-dna-sequences.fasta
#770

 

이는 현재 770개의 feature(=ASV)가 존재한다는 것을 의미합니다.

 

 

리눅스에서 FASTQ파일 다루기

위 글을 참고해 주세요. 

 

4-2) Deblur

Deblur은 DADA2와 달리 전에 QC filtering과정 따로 필요하며, 각 샘플마다 read길이를 맞춰야 합니다.

즉, 필터링 과정에서 일부 정보가 손실됩니다.

 

(1) QC 체크

qiime quality-filter q-score \
 --i-demux demux.qza \
 --o-filtered-sequences demux-filtered.qza \
 --o-filter-stats demux-filter-stats.qza

(2)  denoise

qiime deblur denoise-16S \
  --i-demultiplexed-seqs demux-filtered.qza \
  --p-trim-length 120 \
  --o-representative-sequences rep-seqs-deblur.qza \
  --o-table table-deblur.qza \
  --p-sample-stats \
  --o-stats deblur-stats.qza

파일을 보기 위해  이전과 마찬가지로 metadata tabulate를 사용해 줍니다.

 

qiime metadata tabulate \
  --m-input-file demux-filter-stats.qza \
  --o-visualization demux-filter-stats.qzv
qiime deblur visualize-stats \
  --i-deblur-stats deblur-stats.qza \
  --o-visualization deblur-stats.qzv

이 중 deblur-stats.qzv를 보면 아래와 같습니다.

 

 

이후 분석도 Deblur를 사용하고 싶다면, 5단계로 넘어가기 전에 이름 바꾸어 줍니다. 

mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza

 

 

 

5. FeatureTable and FeatureData summaries

이제 각 ASV에 따라 각 샘플의 read가 몇 개가 존재하는지 확인해 봅시다. 이러한 정보가 담긴 표를 기존에는 OTU table, feature table이라고 불렀으며, ASV 방식을 사용하였으면 ASV table이라고 부릅니다(그러나 OTU table이라는 명사가 굳어져서, ASV를 사용해도 그냥 OTU table이라고 많이 부릅니다...).

 

qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file sample-metadata.tsv
  
qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv

 

 

완성된 table.qzv을 qiime view 홈페이지에서 확인해 봅시다. 

 

[Overview]

1) Table summary

샘플의 총 개수와, 전체 ASV개수를 보여주고 있습니다.

[Tools]	[Dada2]	[Deblur]
metric	sample	sample
Number of samples	34	34
Number of feature	770	485
Total frequency	153,807	102,557

 

2) Frequency per sample

각 샘플당 몇 개의 read count가 존재하는지 보여줍니다. 

[Tools]	[Dada2]	[Deblur]
	frequency	frequency
Minimum frequency	897.0	514.0
1st quartile	1,838.25	1,367.5
Median frequency	4,010.5	2,581.0
3rd quartile	7,013.0	4,952.25
Maximum frequency	9,820.0	6,770.0
Mean frequency	4,523.735294117647	3,016.3823529411766

 

3) Frequency per feature

각 샘플당 몇 개의 ASV가 존재하는지 보여줍니다.  평균적으로 한 샘플당 약 200 여개가 존재합니다.

[Tools]	[Dada2]	[Deblur]
	frequency	frequency
Minimum frequency	2.0	10.0
1st quartile	8.0	16.0
Median frequency	23.0	33.0
3rd quartile	81.75	88.0
Maximum frequency	11,373.0	8,224.0
Mean frequency	199.7493...	211.45773...

비교해 본 결과 Deblur가 Dada2보다 조금 더 까다롭게 체크하여 read를 더 많이 버리는 것을 알 수 있다. 

 

 

[Interactive Sample Detail]

두 번째 페이지인 interactive Sample Detail를 보면 각 샘플에 포함된 read의 수를 볼 수 있습니다.

여기서 Sampling Depth란 한 샘플이 포함하고 있는 read 수의 합을 말합니다.

 

Sampling Depth 바를 움직이면, cut-off에 따라 미달되는 샘플을 분홍색으로 표시해 줍니다.

 

각 연구 분야마다 적절한 sampling depth가 요구됩니다. 만약 장내 미생물의 경우 적어도 1만~2만 depth는 나와야 시퀀싱이 충분히 되었다고 할 수 있습니다. 그러나 피부의 경우 상대적으로 바이오매스가 적기 때문에 2,000~5,000 depth를 기준으로 잡는 연구도 많습니다. 

 

 

[Feature Detail]

세 번째 페이지인 Feature Detail를 보면 각 FeatureID의 정보를 볼 수 있다.

• 각 FeatureID에 해당하는 read의 수는 전체 **"Frequency"**만큼 존재
• 해당 FeatureID를 가진 샘플이 전체에서 **"Of Samples Observed In"**개가 존재

 

각 FeatureID에 해당하는 은 rep-seqs.qzv 보면 알 수 있습니다.

 

 

 

 

각 ASV 아이디에 해당하는 실제 서열을 보여줍니다. 우리가 서열의 길이를 120으로 잘랐기 때문에, 모든 ASV의 길이는 120bp입니다. 

 

 

 

 

 

 

 

 

다음 단계는 2편에서 보도록 합시다!

[QIIME2 튜토리얼] “Moving Pictures” (2)

 

 

 

 

Reference

- https://docs.qiime2.org/2024.2/tutorials/moving-pictures/

- http://compbio.ucsd.edu/wp-content/uploads/2018/07/20180621_oslo_university_microbiome_analysis_with_qiime2_tutorial.pdf

- [KGOL] 마이크로바이옴분석(경희대 이재형)