- 수정 2023.06.05 시각화 방법, 시각화 예시자료 추가
- 수정 2024.01.16 NSTI 부가 설명
- 수정 2024.01.19 contribution 추가
🟦 PICRUSt2란?
단백체 분석이나 RNA seq, Microbiome에서 분석하는 기본 데이터는 table이다. 각 샘플이 행으로 배치되어 있고, 각 열의 이름은 단백질 이름, RNA 이름, ASV 이름, taxa로 바꾸어 갈 뿐이지 형식은 거의 동일하다. 이때 우리는 행 값을 feature라고 부르겠다.
이를 통해서 우리는 각 샘플에 어떤 특정한 열에 해당하는 물질 또는 생물이 풍부하게 존재함을 알 수 있다.
microbiome(Amplicon) 데이터에서 feature는 ASV 혹은 OTU이다. 이 feature을 사용해서 계통정보를 확인해 , 객 taxa의 풍부도를 확인할 수 있다. Metagenome(Shotgun) 데이터는 각 생물의 풍부도뿐만 아니라 모든 서열을 읽어냄으로 어떤 기능을 가진 유전자의 풍부도 또한 확인할 수 있다.
그렇다면 16S rRNA를 이용한 Amplicon 데이터에서 기능 유전체 분석이 불가능할까? 이를 가능하게 해주는 것이 바로 PICRUSt2이다. PICRUSt2는 하버드의 The Huttenhower Lab에서 개발한 도구로, Shotgun metagenome 분석처럼 16S rRNA를 이용해 미생물의 whole genome데이터의 기능과 그 기능을 가진 유전체의 풍부도를 "예측"한다. 이때 각 박테리아의 유전체의 copy 수를 고려한다.
그러나 Shotgun 분석 방식은 각 유전자의 변이 정도(SNV)를 볼 수 있지만, PICRUSt2는 서열에 따른 taxa를 유추하여, 그 기능을 예측하는 프로그램일 뿐 변이정도를 보지는 못한다.
- 관련 강의 정리 : 블로그 글 참고
- 어떻게 구동하는지 간략한 설명본 : https://picrust.github.io/picrust/tutorials/algorithm_description.html
PICRUSt2의 작동방식을 간단하게 말하자면
1. picrust2 내의 reference 파일에 우리가 분석하고자 하는 서열을 정렬(align)시킨다. (fas 파일이 사용)
2. 정렬된 데이터를 reference tree에 위치시킨다.
3. ASV의 copy number를 예측한다.
- reference에 존재하는 taxa의 copy 수 정보를 가져온다
4. 각 샘플마다 gene family의 풍부도를 결정한다. (biom 형식의 otu table 파일이 사용)
5. pathway의 풍부도를 추론한다.
단계별 실제 계산 방식을 아래 간략한 표로 나타내었다.
🟦 설치
https://github.com/picrust/picrust2/wiki/Installation#pre-requisites
Installation
Code, unit tests, and tutorials for running PICRUSt2 - picrust/picrust2
github.com
0) mamba 설치하기
1) Conda환경에서 mamba로 설치하기(개발자 추천👍👍)
mamba create -n picrust2 -c bioconda -c conda-forge picrust2=2.5.1
conda activate picrust2
🚨 오류 발생 시 (보통 mamba가 아니라 conda사용 시) 🚨
오류 1 : HTTP 에러
Collecting package metadata (repodata.json): failed
CondaHTTPError: HTTP 000 CONNECTION FAILED for url <https://conda.anaconda.org/r/linux-64/repodata.json>
Elapsed: -
An HTTP error occurred when trying to retrieve this URL.
HTTP errors are often intermittent, and a simple retry will get you on your way.
'https//conda.anaconda.org/r/linux-64'해결방법
# https://next-brain.tistory.com/40
conda config --set ssl_verify no # ssl 안 씀
conda install openssl
오류 2 : Solving environment에서 계속 막혀있음
Retrieving notices: ...working... done
Collecting package metadata (current_repodata.json): done
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json):
Solving environment: - # 계속 반복해결방법
# https://stackoverflow.com/questions/63734508/stuck-at-solving-environment-on-anaconda
conda config --remove channels conda-forge
conda config --add channels conda-forge
conda update conda혹은 conda 재설치
conda update conda -c conda-canary
conda config --set channel_priority false
conda create --name your_env_name
conda activate your_env_name
conda create -n picrust2 -c bioconda -c conda-forge picrust2=2.5.1pytest에러 메시지 없이 잘 실행되는 것을 확인했다.
2) 직접 실행파일 다운로드하여 설치하기
wget https://github.com/picrust/picrust2/archive/v2.5.1.tar.gz
tar xvzf v2.5.1.tar.gz
cd picrust2-2.5.1/condaenv create -f picrust2-env.yaml
conda activate picrust2 # 실행 완료
conda deactivate picrust2 # 환경 종료
pip install --editable .
# -e, --editable은 편집 가능한 모드로 프로젝트를 설치하는 것이다.
# 가볍게 알아보니 본인의 환경에 맞게 좀더 기능이 추가되는것 같다
위와 같이 Successfuuly installed PICRUST2 가 나오면 다운에 성공한 것이다.
그러나 함수가 잘 돌아가는 시 시험 삼아 테스트해 보면 오류가 발생한 것을 볼 수 있다.
pytest
중간에 오류는 /home/ksy/picrust2-2.5.1/picrust2/pathway_pipeline.py 에서 오류가 났으며, pandas 패키지와 관련이 있는 것 같다. 경고를 피하려면 아래와 같은 내용을 drop 하라고 말한다.
이런 자잘한 오류를 피하기 위해 conda환경에서 설치하는 것을 추천한다.
3) 버전 확인하기
print_picrust2_config.py
🟦 튜토리얼
1) 데이터 다운로드하기
분석 환경 만들기
mkdir picrust_tutorial
cd picrust_tutorial
예제 데이터 다운로드
wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust2_tutorial_files/chemerin_16S.zip
unzip chemerin_16S.zip
cd chemerin_16S
다운로드한 데이터 보기
less seqs.fna
biom head -i table.biom
biom summarize-table -i table.biom
2) 전처리
- https://github.com/picrust/picrust2/wiki/PICRUSt2-Tutorial-(v2.3.0-beta)#data-pre-processing
- 전처리 시 유의할 점
- 오직 한 샘플에만 있는 ASV는 noise로 생각되니 삭제할 것, input ASV를 줄여야 계산 속도가 줄어든다
- 위와 동일하게 낮은 read 수를 가진(low depth) 샘플을 확인하고 제거해야 한다.
3) PICRUSt2 pipeline
- 이제 프로그램을 돌려보자!
- 가장 간단하게 할 수 있는 방법은 picrust2_pipeline.py 을 통해서 한 번에 돌리는 것이다. 그러나 이 소프트웨어를 이해하기 위해 각 단계별로 천천히 진행해 보자.
- 차후 단계에서 -p 1은 사용되는 core 수를 나타내는데, 가용한 만큼 설정하는 것이 좋다.
0. 분석 위치 만들기
mkdir picrust2_out_pipeline
cd picrust2_out_pipeline
1. 분석하고자 하는 서열을 reference tree에 맞추어 정렬 📌
- PICRUSt2는 HMMER 프로그램을 이용해 각 ASV의 서열을 referense 16S 서열에 매칭한다(Align).
- 일단 align 되면, 16S 서열을 기준으로 만든 계통수(newick format의 계통수 사용)에 우리가 분석하고자 하는 query서열들이 위치하게 된다. (이때는 EPA-NG 사용)
- Reference 계통수는 IMG database에서 가져온 20,000여 개의 16S 서열로 만들어졌다.
place_seqs.py -s ../seqs.fna -o out.tre -p 1 \
--intermediate intermediate/place_seqs- s ../seqs.fna : ASV의 서열 정보를 담고 있는 FASTA형식의 input 파일
- o out.tre : newick 형식의 계통수 output 파일
- p 1 : 사용할 프로세서 수
--print_cmds : 프로그램이 돌아가는 중간에 실행되는 command를 보여준다
--intermediate intermediate/place_seqs : 중간 파일이 저장되어야 하는 위치. output 파일은 가장 가능성이 높은 배치만 각 ASV에 대해 출력되지만 중간 파일은 트리에서 ASV당 모든 잠재적 배치의 JPLACE 파일
(intermediate/place_seqs/epa_out/epa_result.jplace)을 모두 찾아볼 수 있다.
- 주의할 점은, 이 command를 할 때마다 같은 tree가 output으로 나온다. 그러나 많은 데이터 양에서 일부 변동이 있을 순 있으며, 이는 분석마다 정확한 재현이 어려울 수 있다.
2. Hidden State Prediction (HSP) of gene families 유전자 familes의 숨겨진 상태 예측 📌
- 위 단계에서는 각 ASV에서 genome을 예측한다.
- hsp.py는 maximum parsimony 알고리즘을 이용한다. 이 스크립트는 각 ASV에 대해 누락된 genome 즉 각 ASV의 copy number(한 gene이 몇 개의 사본 유전자를 가지고 있는지) 수을 예측한다. 이는 다음 단계가 이유가 나온다.
- 여러 family level에서 예측이 가능하지만, 컴퓨팅 시간 감소를 위해 Metacyc(https://metacyc.org/)의 Enzyme Classification(EC) 아이디를 가진 데이터 베이스를 사용한다. 물론 KEGG도 사용 가능하다.
hsp.py -i 16S -t out.tre -o marker_predicted_and_nsti.tsv.gz -p 1 -n
# EC
hsp.py -i EC -t out.tre -o EC_predicted.tsv.gz -p 1
# KO
hsp.py -i KO -t out.tre -o KO_predicted.tsv.gz -p 1hsp.py [-h] -t PATH -o PATH [-i {16S,COG,EC,KO,PFAM,TIGRFAM,PHENO}]
[--observed_trait_table PATH] [-e EDGE_EXPONENT] [--chunk_size CHUNK_SIZE]
[-m {mp,emp_prob,pic,scp,subtree_average}] [-n] [--check] [-p PROCESSES]
[--seed SEED] [--verbose] [-v]-i 16S : 검색할 친구
-t out.tre : input 파일
-o marker_predicted_and_nsti.tsv.gz : output 파일
- n : nearest-sequenced taxon index (NSTI) values for each ASV를 추가한다____nearest-sequenced taxon index (NSTI) 값은 각 tree에 위치시킨 ASV에서 가장 가까운 reference 서열까지의 분기 길이를 말한다.
결과 파일을 관찰해 보자
zless -S marker_predicted_and_nsti.tsv.gzzless : zip 파일을 잘 볼 수 있게 해 줌
-S : tsv파일을 각 tab으로 구분해 주어 사람이 보기 편하게 출력
- Sequence : 첫 번째 열, ASV의 이름이다.
- 16S_rRNA_Count : 예측된 ASV의 copy 수를 말한다.
- metadata_NSTI : 각 ASV마다 NSTI값을 담고 있다. 실제 데이터에서 2 점이상일 경우 대부분 noise(오류)로 밝혀졌다. 이 NSTI값을 통해 전체 ASV의 특이값을 확인할 수 있다.
하지만 논문에서도 cutoff 2 NSTI는 juck sequence를 버리는 정도의 역할이지, 퀄리티가 좋다고 볼 수 없다. (참고)
OTU 서열 기준 PICRUSt에서는 대략 0.06 이하(매우 높은 예측 퀄리티)를 권장하고 있다.
실제 picrust2를 사용한 논문의 Reviewer에서도 이 부분을 관련해서 질문을 한다 (참고). 논문 저자는 0.15 기준으로 했을 때와 2를 기준으로 한 결과물의 상관계수는 0.91이지만, 0.15를 기준으로 필터링한다면, 전체 ASV 중 매우 일부만 남게 되어 이는 전체 taxa를 대표한다고 볼 수 없다고 답변했다.
관련 포럼글을 확인해 보았다(https://forum.qiime2.org/t/q2-picrust2-recommended-max-nsti-value/8990). 위 글에서는 NTSI값의 cut-off 가 전체적인 기능 예측에 미치는 영향은 매우 적다고 한다.
다른 파일도 살펴보자.
zless -S marker_predicted_and_nsti.tsv.gz
이 파일은 각 ASV마다 Enzyme Classification(EC)에서 예측된 Copy 수를 보여준다.
sequence EC:1.1.1.1 EC:1.1.1.10 EC:1.1.1.100 ...
20e568023c10eaac834f1c110aacea18 2 0 3 ...
23fe12a325dfefcdb23447f43b6b896e 0 0 1 ...
288c8176059111c4c7fdfb0cd5afce64 1 0 1 ...
...
- 이 파일에서 NSTI값은 이 테이블에 존재하지 않는 이유는 -n 옵션을 쓰지 않았기 때문임.
- EC numbers는 화학적 반응을 기반으로 하는 gene family의 하나이다. 예를 들어 EC 1.1.1.1은 alcohol dehydrogenase과 일치한다.
- 이 튜토리얼에서는 EC의 수는 MetaCyc pathway 단계를 추론하는 데에 있어 사용되기 때문에 여기에 집중해야 한다.
3. Metagenome 예측을 생성 📌
- 이 전 단계에서는 이제 각 ASV에서 gene families에 따른 copy 수를 예측했다.
- 이 단계에서는 단순히 copy 수의 정보에만 의존하는 것이 아니라, 상대적으로 풍부한 ASV에 의해 가중치를 좀 더 부과한 결과를 얻어보자.
# EC
metagenome_pipeline.py -i ../table.biom -m marker_predicted_and_nsti.tsv.gz -f EC_predicted.tsv.gz \
-o EC_metagenome_out --strat_out
# KO
metagenome_pipeline.py -i ../table.biom -m marker_predicted_and_nsti.tsv.gz -f KO_predicted.tsv.gz \
-o KO_metagenome_out --strat_out-i ../table.biom : input 파일, 보통 OTU table, ASV table이라고 불리는 풍부도 정보를 가진 table을 말한다.
-m marker_predicted_and_nsti.tsv.gz : 이전 단계의 output을 마커로 사용한다다
-f EC_predicted.tsv.gz :
-o EC_metagenome_out : output 파일이 들어있는 directory
--strat_out : output 형식을 계층화된 파일로 출력한다.
--max_nsti : 위 값을 통해서 cut off 설정 가능
EC_metagenome_out directory안에 있는 파일을 보자
- pred_metagenome_unstrat.tsv.gz : 샘플당 전반적인 EC 의 풍부도
- pred_metagenome_contrib.tsv.gz : 각 샘플에서 어떻게 ASV 가 gene family의 풍부도에 기여했는지 보여주는 계층화된 표
- seqtab_norm.tsv.gz : 16S copy수로 표준화 된 ASV 풍부도
- weighted_nsti.tsv.gz : 샘플 당 NSTI값의 평균을 말한다. 이 파일은 각 분석 데이터셋의 이상값을 확인하기에 유용하다. PICRUSt1에서는 weighted NSTI 값이 0.06보다 적거나, 0.15보다 크면 너무 값이 낮거나 높다고 할 수 있다. 물론 이 cut off가 모든 데이터 셋에서 적용되는 것은 아니지만 대략적으로 다른 데이터셋과 비교하기 위함으로 사용하기엔 유용하다.
분석된 파일을 살펴보자
zless -S EC_metagenome_out/pred_metagenome_unstrat.tsv.gz
샘플에 따른 기능 유전체의 풍부도를 의미한다.
- column이 각 샘플을 의미하고, 각 row가 대사의 종류이다.
- 이 표은 그냥 보면 ASV의 풍부도 표와 비슷해 보인다. 그러나 이 표는 상대적 풍부도가 아니다.
- 각 데이터를 계층적으로 정리한 output을 보자
zless -S EC_metagenome_out/pred_metagenome_contrib.tsv.gz
각 기능 유전체에 어떤 ASV가 얼마나 기여했는지 볼 수 있다. 만약 질병 A와 정상대조군의 picrust2 결과에서, 질병 A군에서 B 합성 기능이 차등적으로 많이 분포되었다면, 위 contrin.tsv를 통해 B합성 기능에 영향을 주는 ASV를 역으로 추적할 수 있다.
각 column데이터를 살펴보자
- sample : sample의 이름
- function : EC의 ID
- taxon : ASV의 ID 혹은 taxon의 ID
- taxon_abun : ASV의 read count
- taxon_rel_abun : ASV의 relative abundance (각 sample의 총 read수로 각 ASV의 read count를 나눈 것)
- genome_function_count : 예측된 16S의 copy 수
- taxon_function_abun : taxon_abun * genome_function_count
- taxon_rel_function_abun : taxon_rel_abun * genome_function_count
- norm_taxon_function_contrib : taxon_rel_function_abun의 relative abundance
하나 주의할 것은, 기존에 다른 분석 툴은 rare 한 ASV를 제거하기 위해 최소 몇 read이상, 몇 개 샘플 이상에서만 존재하는 ASV를 걸러내는 필터링 단계가 존재한다. 그러나 위의 pred_metagenome_contrib.tsv.gz 은 이러한 방법을 사용하면 안 된다.
4. Pathway 단계 추론
- default값으로는 EC number abundances를 이용한 MetaCyc pathway abundances를 추론한다.
- 이 스크립트는 아래 단계를 포함하고 있고 있다
1) EC 를 MetaCyc 에 맞게 재 그룹화
2) MetaCyc pathways 를 MinPath기반으로 추론
3) abundance of pathways 계산
# EC to METACyc
# abundance file
pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
-o pathways_out -p 1
# contrib file
pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_contrib.tsv.gz \
-o pathways_out -p 1
pathway_output 경로에 있는 파일은 아래와 같다
- path_abun_unstrat.tsv.gz: MetaCyc의 기능 유전체 풍부도
- path_abun_contrib.tsv.gz: MetaCyc에서 taxa가 각 function에 얼마나 기여하는지
zless -S path_abun_contrib.tsv.gz
- 주의할 점은 위 표는 얼마나 많은 각 ASV가 전체 샘플 전반의 pathway의 풍부도에 기여했는지를 나타내는 것이지, 각 ASV에서 예측된 pathway의 풍부도가 아님!
- 만약 각 ASV에서의 기여도를 보고 싶으면 "--per_sequence_contrib"옵션을 추가하면 된다. 그러나 많은 계산 시간이 소요될 것이다.
zless -S path_abun_unstrat.tsv.gz
5. 기능적인 설명 추가(annotation)
- 이 단계는 추가적인 단계로 각 기능에 대한 설명을 추가해 준다
# EC
add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC \
-o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
# MetaCyc
add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC \
-o pathways_out/path_abun_unstrat_descrip.tsv.gz
# KO
add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO \
-o KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz
6. PICRUSt output을 분석하는 tool
- 시각화 : STAMP (R/Python package) https://beikolab.cs.dal.ca/software/STAMP
> 시각화 방법 : https://github.com/picrust/picrust2/wiki/STAMP-example
- 통계분석 : ALDEx2 https://bioconductor.org/packages/release/bioc/html/ALDEx2.html
- Downstream analysis
> http://picrust.github.io/picrust/tutorials/downstream_analysis.html#downstream-analysis-guide
> https://rstudio-pubs-static.s3.amazonaws.com/566863_687400bd7e8742568e73bf167fc42d3d.html
> https://ycl6.github.io/16S-Demo/4_picrust2_tutorial.html
> https://borenstein-lab.github.io/burrito/tutorial_BURRITO.html
> ggpicrust2 packge
주의해야 할 점은 Differential abundance test의 경우 shotgun과 amplicon 데이터 간에는 상당한 차이가 있다. 📌
🟦 시각화 예시
1) heatmap
아마도 lefse분석 이후 LDA score가 2.5 이상, p-value값이 0.05 이하인 kegg pathway를 각 그룹에 따라 시각화함.
위 히트맵은 상관계수에 따라 색을 달리 하였음
2) Relative abundance bar
각 그룹마다 상대적인 분포를 나타냄. 보통 각 그룹사이에 통계적 계산을 통해 유의성을 표시하기도 함(*)
3) LDA score
보통 lefse분석에서 사용된 LDA 점수를 통한 시각화. 그룹에 따라 축의 방향을 달리 하였다.
🟦 Reference
- 튜토리얼 강의 영상(영어) : https://www.youtube.com/watch?v=xZ7yc-GKcSk
- 논문 : PICRUSt2 for prediction of metagenome functions, 2020, nature biotechnology
