Tutorial – Microbiome Analysis (amplicon sequencing)

Leave a Comment / Microbiome (BDSL) / By

Introduction

Microbimoe은 환경 안에 존재하는 미생물 전체를 의미한다. 미생물 간 상호작용 및 미생물과 환경의 상호작용 연구를 통해, microbiome의 기능을 연구할 수 있다. 특히 사람의 장 혹은 피부 등에 서식하는 미생물이 인체의 건강과 질병에 미치는 영향이 활발하게 연구되고 있다. 이러한 연구를 위해서는 미생물에 대한 데이터를 획득해야 하는데, 두 가지 NGS 실험 방법이 널리 이용된다. Amplicon sequencing은 특정 서열을 가진 DNA를 PCR을 이용해 증폭한 뒤 서열을 결정하는 실험 방법으로, 세균을 대상으로 할 때는 세균의 ribosome을 코딩하는 유전자인 16S rRNA 서열을 증폭하여 분석한다. Whole metagenome sequencing은 샘플(혹은 환경)안에 있는 DNA를 조각낸 후 그 서열을 결정하는 방법으로 특정한 환경에 있는 세균이 가진 DNA 서열을 통해 유전자의 종류와 양을 분석하여, 특정한 분자적 기능의 활성화/저하를 분석하는 방법이다.

Prerequisite

Tutorial 진행에 앞서 아래의 내용을 숙지한다.

Summary

Tutorial 전체 과정 간단 요약

  • Input: amplicon sequencing data (FASTQ format file, paired end sequencing)
  • Intermediate data
    • Representative sequences (sequence of unique amplicons)
    • Feature table (microbial profile of individuals)
  • Output
    • Taxonomic profile
    • Functional profile
    • Host factor association

Tutorial

Summary

(목적) Amplicon sequencing 실험 결과 얻어진 FASTQ format 데이터를 이용하여 미생물의 profile을 구하고, 기초적인 분석 방법의 종류를 이해한다.

(방법) QIIME2 suite를 이용하여 샘플 별 FASTQ format 형태로 정리된 amplicon sequencing 데이터를 샘플에 존재하는 (1) DNA 서열과 (2) 샘플 별 해당 DNA 서열의 개수로 구성된 표로 변환하고, 이렇게 변환된 정보를 이용해 다양성 분석을 수행한다. 여기서 DNA 서열은 representative sequence , 샘플 별 DNA 서열 개수표는 feature table이라 한다.

(추가) 이 tutorial은 amplicon sequencing 실험 데이터 분석의 전산적 과정을 표현한 것으로, 각각의 단계에 대한 자세한 내용을 학습 및 이해가 필요하다. 또한 이 tutorial은 기본적인 분석 과정 만을 포함하고 있으므로, 다양한 응용 연구를 위해서는 통계적인 방법론 및 기계학습 방법론을 학습하여 적용할 필요가 있다.

  1. Tutorial 각 과정의 이론적 배경 학습
  2. Feature table- 건강 지표 상호 관련성 분석
  3. Feature table 이용 건강 지표 예측
  4. Representative sequence 서열 기반 feature 변환 (functionality analysis, like picrust2)

(사전질문) Tutorial 진행에 앞서 아래의 질문에 답해보고, 아래 질문에 답하기 어려우면 prerequisite의 주제들을 다시 한번 살펴본다.

  1. FASTQ format 파일은 무엇인가?
  2. 16S rRNA sequencing을 통해 얻은 FASTQ 파일은 어떤 정보를 가지고 있는가?
  3. Forward read와 reverse read는 무엇인가?
  4. V3-V4 region이란 무엇인가?

1. Tutorial 파일 받기

Tutorial 실행을 위해 Visual Studio Code (VSC)를 이용해 서버에 접속하고, 작업 디렉토리를 생성하고, 필요한 파일을 다운 받는다.

1.1. 작업 directory 만들기

작업 directory를 만들고 Visual Studio Code에서 해당 directory를 연다. 

cd ~
mkdir tutorial
cd tutorial
mkdir qiime2_amplicon

1.2. Tutorial 파일 받기

필요한 파일을 다운 받는다.

cd ~/tutorial/qiime2_amplicon
wget https://bdsl.jbnu.ac.kr/blog/wp-content/uploads/2023/05/qiime2_tutorial.zip --no-check-certificate
qiime2_tutorialDownload

2. Install QIIME2

QIIME2 프로그램을 설치한다. 설치는 conda를 이용해서 진행한다. 아래 명령은 2022.2 버전의 QIIME2 suite 설치를 위한 예제이다.

QIIME2는 microbiome 분석에서 자주 사용하는 프로그램을 모아 놓은 program suite이다. 추후 Web page에서 QIIME2에 대한 정보를 추가로 수집하도록 한다 (https://qiime2.org/).

특정 environment가 활성화되어 있다면, 아래 명령어를 이용해 먼저 불활설화 시키고 설치를 진행한다.

conda deactivate
conda clean --all -y

QIIME2 Environment 설치 

conda env create -n qiime2-amplicon-2024.10 --file https://data.qiime2.org/distro/amplicon/qiime2-amplicon-2024.10-py310-linux-conda.yml

3. Preprocess

Microbiome 분석은 여러 샘플을 비교 분석하는 것이 일반적이다. QIIME2는 흩어져 있는 여러 파일을 새로운 파일 하나로 모은다. 이 과정은 qiime tools import 방법을 이용해 진행한다.

(FASTQ 파일 이름) 분석에 활용되는 FASTQ 파일은 NGS 기기를 통해 생성된 원 데이터를 가공한 결과물이다. Bar code 서열을 이용하여 샘플 별로 별도의 FASTQ 파일을 생성하는 것이 일반적이다. 이 과정을 demultiplexing 이라 한다. 또한, paired sequencing이 일반적인데, 이 때 생성되는 forward read에는 _1, reverse read에는 _2와 같은 tag를 붙여 표기하는 것이 일반적이다.

(Primer 서열) Primer 서열은 PCR을 위해 도입한 서열이고, redundancy를 가진 서열을 이용하는 경우가 많이 있기 때문에 이를 제거하는 것이 바람직할 것이다. 따라서 추가 분석에 앞서 primer 서열 제거가 필요하다.

3.1. QIIME2 환경 활성화

QIIME2 실행을 위해서는 설치한 conda 환경을 활성화 시켜야 한다. 

conda activate qiime2-amplicon-2024.10

3.2. MANIFEST 파일 만들기

Manifest 파일은 선언 파일이라고도 하며, 파일에 대한 설명을 포함하는 파일이다. 분석에 필요한 FASQT 파일의 위치와 각 파일이 대응하는 샘플에 대한 정보를 제공해야 한다. 이를 위해 MANIFEST 파일을 작성해야 한다. 이 파일은 샘플 identifier (sample-id), forward read 파일 위치 (forward-absolute-filepath), reverse read 파일 위치 (forward-absolute-filepath) 등의 정보를 포함하며, 아래와 같은 파일을 생성할 수 있다. 

sample-id	forward-absolute-filepath	reverse-absolute-filepath
S001	$PWD/fastq/ERR3472004_1.fastq	$PWD/fastq/ERR3472004_2.fastq
S002	$PWD/fastq/ERR3472005_1.fastq	$PWD/fastq/ERR3472005_2.fastq
S003	$PWD/fastq/ERR3472006_1.fastq	$PWD/fastq/ERR3472006_2.fastq
S004	$PWD/fastq/ERR3472007_1.fastq	$PWD/fastq/ERR3472007_2.fastq
S005	$PWD/fastq/ERR3472008_1.fastq	$PWD/fastq/ERR3472008_2.fastq
S006	$PWD/fastq/ERR3472009_1.fastq	$PWD/fastq/ERR3472009_2.fastq
S007	$PWD/fastq/ERR3472010_1.fastq	$PWD/fastq/ERR3472010_2.fastq
S008	$PWD/fastq/ERR3472011_1.fastq	$PWD/fastq/ERR3472011_2.fastq
S009	$PWD/fastq/ERR3472012_1.fastq	$PWD/fastq/ERR3472012_2.fastq
S010	$PWD/fastq/ERR3472013_1.fastq	$PWD/fastq/ERR3472013_2.fastq

(*) 이 예제에서 $PWD는 현재 디렉토리를 의미하며, 다른 위치에 파일이 존재하는 경우 해당 위치에 대한 절대 경로를 지정해야 한다.

3.3. FASTQ 파일 수집

MANIFEST 파일에 기록되어 있는 FASTQ 파일을 수집하여 하나의 파일(paired-end-demux.qza) 파일로 만든다.

# Import FASTQ format files
mkdir intermediate
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path manifest.tsv \
  --output-path intermediate/paired-end-demux.qza \
  --input-format PairedEndFastqManifestPhred33V2

3.4. PRIMER 서열 제거

Primer 서열을 제거한다. (paired-end-demux-trimmed.qza)

# Trim primer sequences
qiime cutadapt trim-paired \
  --i-demultiplexed-sequences intermediate/paired-end-demux.qza \
  --p-cores 20 \
  --p-front-f CCTACGGGNGGCWGCAG \
  --p-front-r GACTACHVGGGTATCTAATCC \
  --p-minimum-length 100 \
  --p-discard-untrimmed \
  --o-trimmed-sequences intermediate/paired-end-demux-trimmed.qza

단계 4. Denoising (DADA2)

4. Denoising and Amplicon Sequence Variant (ASV) detection

4.1. Denoising (DADA2) and feature-table 만들기

Amplicon sequencing을 통해 읽은 sequence read는 PCR 과정 및 sequencing 과정의 error를 포함한다. 따라서, 최종적으로 읽은 서열에 오류가 존재한다. 이 오류를 위한 방법으로, 서열 유사도에 기반하여 유사 서열을 모아 operational taxonomic unit (OTU)를 구성하는 방법이 널리 이용되었다. 관련 내용을 수집해서 공부하도록 하자. 다른 방법으로 통계적 분석을 통해 오류로 인해 발생하는 서열을 찾아서 제거할 수 있다. 다양한 방법이 있으나 그 중 하나인 DADA2가 널리 활용되고 있다. 여기서는 이 방법을 이용해 서열의 오류를 줄이는 방법을 이용한다.

QIIME2 DADA2 package는 추가로 일련의 데이터 처리 과정을 거치며, 이를 통해 고유한 서열인 representative sequence(representative_sequences.qza) 와 샘플 별 read 개수를 나타내는 표인 feature-table (table.qza)을 생성한다.

# Generate feature table
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs intermediate/paired-end-demux-trimmed.qza \
  --p-trunc-len-f 270 \
  --p-trunc-len-r 220 \
  --o-representative-sequences intermediate/representative_sequences.qza \
  --p-n-threads 20 \
  --o-table intermediate/table.qza \
  --output-dir ./dada2_output

DADA2 논문

Callahan, B. J., et al. (2016). “DADA2: High-resolution sample inference from Illumina amplicon data.” Nat Methods 13(7): 581-583.

Callahan-2016-Dada-high-resolution-sample-inferenDownload

4.2. Rarefy

각 샘플 별로 read 수의 차이가 발생한다. 이러한 차이는 추후 분석 과정에서 문제를 일으킬 수 있다. 따라서 샘플 별 read 숫자를 같게 만들어 주는 과정이 필요하다. 즉, 샘플의 read 개수를 10,000개로 통일 시킬 수 있을 것이다. 이 경우 10,000개 이하의 read를 가진 샘플은 삭제된다. 샘플 당 read 개수가 많을 수록 정보 양은 증가하지만, 분석에 사용되는 샘플 숫자가 줄어들 수도 있다. 


# Rarefy
qiime feature-table rarefy --i-table intermediate/table.qza \
  --p-sampling-depth 10000 \
  --o-rarefied-table intermediate/feature-table_rarefied.qza \
  --output-dir output_rarefaction

5. Microbiome analysis

5.1. Build phylogenetic trees

서열 유사도에 기반하여 ASV 간 유사도 측정을 위한 계통수를 생성한다. 이 계통수는 추후 ASV 유사도 비교가 필요한 분석 과정에 활용한다. 

qiime alignment mafft   --i-sequences intermediate/representative_sequences.qza --o-alignment intermediate/aligned.qza --p-n-threads 20
qiime alignment mask --i-alignment intermediate/aligned.qza --o-masked-alignment intermediate/masked.qza
qiime phylogeny fasttree --i-alignment intermediate/masked.qza --o-tree intermediate/unrooted-tree.qza --p-n-threads 20
qiime phylogeny midpoint-root --i-tree intermediate/unrooted-tree.qza --o-rooted-tree intermediate/rooted-tree.qza

6. Diversity analysis

Microbiome은 미생물의 생태계를 연구하는 과정으로 생각할 수도 있다. 미생물의 다양성이 높을 수록 변화하는 환경에 적응하기 쉬우며, 다양한 영양소를 대사 하는 등의 기능을 가질 수 있다고 여겨진다. 따라서 미생물 다양성 분석이 가장 첫 단계이다. 다양성도 여러 가지 측면에서 고려할 수 있는데, 일반적으로는 하나의 환경에서의 다양성과, 두 환경의 차이에서 오는 다양성 두 가지를 많이 고려한다. 전자를 alpha diversity 후자를 beta diversity라 한다.

6.1. Alpha diversity

Alpha diversity는 특정 환경 안의 미생물 다양성을 나타내는 지표이다. 장내미생물로 생각하면, 한 사람의 장에 얼마나 다양한 미생물이 존재하고 있는지에 대한 측정치라고 할 수 있다. 이는 diversity alpha 모듈을 이용해 계산할 수 있다.

다양성을 측정하는 방법은 여러 가지가 존재한다. 여기서는 Chao1, observed features, Shannon entropy를 계산하였다.

# Create directories for analysis
mkdir diversity
mkdir results 

# alpha diversity
qiime diversity alpha --i-table intermediate/feature-table_rarefied.qza --p-metric chao1 --o-alpha-diversity diversity/chao1_vector.qza
qiime diversity alpha --i-table intermediate/feature-table_rarefied.qza --p-metric observed_features --o-alpha-diversity diversity/observed_features_vector.qza
qiime diversity alpha --i-table intermediate/feature-table_rarefied.qza --p-metric shannon --o-alpha-diversity divers

6.2. Beta diversity

Beta diversity는 두 환경 간 미생물의 차이를 의미한다. 두 환경을 비교한 결과이기 때문에 2차원 array 혹은 table 형태로 데이터를 표현할 수 있다.

두 개체의 차이를 측정하는 방법도 다양할 수 있다. 이 때 서열의 유사성까지 고려해서 차이점을 계산할 수도 있으며, 이 경우 phylogenetic tree를 이용할 수도 있다. 여기서는 Bray-Curtis, UniFrac 값을 계산한다.


# beta diversity
qiime diversity beta --i-table intermediate/feature-table_rarefied.qza --p-metric braycurtis --o-distance-matrix diversity/braycurtis_dmat.qza
qiime diversity beta-phylogenetic --i-table intermediate/feature-table_rarefied.qza --i-phylogeny intermediate/rooted-tree.qza --p-metric weighted_unifrac --o-distance-matrix diversity/weighted_unifrac_dmat.qza

6.3. Diversity와 연관성 분석

Diversity와 인체의 상태 간의 관련성을 분석하기 위해서는 diversity와 여러 지표사이의 관련성을 분석할 수 있다. 이는 diversity alpha-group-significance 방법이나, beta-group-significance방법을 이용해 분석 할 수 있다.

# Significance group and Visualization 
qiime diversity alpha-group-significance \
  --i-alpha-diversity diversity/shannon_vector.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization results/shannon-group-significance.qzv

# Significance group and Visualization 
qiime diversity beta-group-significance \
  --i-distance-matrix diversity/weighted_unifrac_dmat.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column test-grp \
  --o-visualization results/weighted-unifrac-body-site-significance.qzv \
  --p-pairwise

6.4. Principle Coordinate Analysis

Beta diversity 는 개인 간 미생물 유사성/차이성을 나타내는 지표이다. 다만, 여러 샘플간의 차이를 table로만 표현하면 구분이 어렵기 때문에 2차원으로 시각화해서 한 눈에 알아보기 쉽게 표현하는 것이 바람직할 것이다. 이를 위해 Principle Coordinate Analysis (PCoA) 방법을 이용할 수 있다. 

qiime diversity pcoa --i-distance-matrix diversity/weighted_unifrac_dmat.qza --o-pcoa diversity/weighted_unifrac_pcoa_result.qza

qiime emperor plot \
  --i-pcoa diversity/weighted_unifrac_pcoa_result.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization results/weighted-unifrac-pcoa.qzv

(*) emperor plot은 PCoA 결과를 QIIME2 VIEW를 통해 시각화 할 수 있도록 만들어 준다. weighted-unifrac-pcoa.qzv 파일을 QIIME2 VIEW 페이지에 업로드 하여 PCoA 결과를 시각화 할 수 있다. QIIME2 VIEW 페이지의 주소는 https://view.qiime2.org/이다.

7. Taxonomy annotation

Amplicon sequencing에서 미생물은 16S rRNA의 부분 서열로 표현된다. 그리고 이 서열로 부터 미생물의 종류도 유추할 수 있다.

다만, 서열로 부터 유추할 수 있는 미생물의 수준은 차이가 있다. 특정 서열이 한 종의 미생물에서만 관찰된다면, 샘플에서 해당 서열이 관찰되는 경우 샘플에 해당 미생물이 존재한다고 생각할 수 있다. 다만, Illumina sequencing 방법의 한계로 인해 우리는 16S rRNA의 부분 서열 정보만 가지고 있다. 이 경우 여러 종의 미생물의 동일한 서열을 지니는 것도 가능하다. 다만, 종보다 큰 수준인 속(genus)에서 보면 하나의 속 만이 해당 서열을 가지는 경우도 있을 수 있다. 이 경우, DNA 서열에서 속 수준의 정보를 유추할 수 있을 것이다.

DNA 서열로 부터 미생물의 종류를 예측하기 위해서는 미생물의 계통수를 예측하는 모델이 필요하다. 어떻게 예측하는 모델을 만들 수 있을지 생각해 보자.

7.1. taxonmy annotation

서열로 부터 미생물의 계통수를 예측한다.

mkdir taxonomy
qiime feature-classifier classify-sklearn \
  --i-classifier silva-138.1-ssu-nr99-v3v4-classifier.qza \
  --i-reads intermediate/representative_sequences.qza \
  --o-classification taxonomy/taxonomy.qza \
  --p-n-jobs 20

7.2. Taxonomy visualization

QIIME2를 이용해 taxonomy 조성을 시각화 할 수 있다.

Further work

이 tutorial에서는 amplicon sequencing 분석 과정에서 공통적으로 진행하는 내용을 설명하였다. 각 과정의 의미에 대해서 이해하고, 각 과정에서 필요한 전산적 방법이 무엇인 생각해 보고 공부하여 전 과정을 정확히 이해할 필요가 있다.

Views: 393

Leave a Comment

Your email address will not be published. Required fields are marked *

Scroll to Top