이 작업의 이해를 돕기 위한 상위 페이지: Korea BioData Station (K-BDS)

당시 시퀀싱되었던 균주는 KCTC 번호가 부여된 것(대부분 표준균주)이 가장 많으며, 여기에는 DMSZ나 ATCC에서 들여와서 KCTC의 정식 컬렉션이 된 것도 있다. 샘플 ID를 '숫자' 또는 '숫자T'로 표기한 것은 KCTC 자원으로서 숫자는 KCTC 번호에 해당한다. 시퀀싱 대상에는 당시 KCTC 소속 연구자가 연구 과정을 통해 개별적으로 분리·동정한 것 또는 외부에서 입수한 것 소수를 포함한다. 이런 부류의 것은 대부분 공식 KCTC 컬렉션이 아니다.

균주가 재분류되어 공식 명칭이 바뀐 것은 별 문제가 되지 않는다. 그러나 현재 기준으로 분양 불가라서 KCTC 웹사이트에 정보가 공개되지 않은 것을 어떻게 해야 하는지 아직 판단이 잘 서지 않는다. 예를 들어 Bizionia paragorgiae KCTC 12304^T를 KCTC에서 검색하면 다음과 같이 분양불가자원으로 나온다. 그러나 이 일을 진행하면서 KCTC와 연락을 주고받는 도중에 분양이 가능한 것으로 바뀌었다.

그러나 DSMZ(DSM 23842T)에서는 분양에 제한이 없는 것으로 나타난다. 흥미로운 것은 이 균주가 러시아에서 KCTC를 거쳐 DSMZ에 기탁된 것으로 나온다. 이러한 균주가 DSMZ에서는 분양이 가능하지만, 정작 이를 제공한 KCTC에서는 분양 불가라니? 2005년에 발표된 신종 보고 논문에서도 당시 KCTC에 재직하던 두 사람이 저자로 참여하고 있다.

균주를 보관하고는 있지만 다시 자라나지 않을 경우 분양 불가인 것도 있을 것이다. 자원센터에서는 이러한 사연을 전부 공개하는 것을 좋아하지는 않을 것으로 본다. 따라서 분양 불가 사유를 명확히 구별하는 것이 중요하다. 미생물자원은행의 영문 약자에 관한 설명은 RIKEN BRC의 Acronyms of Other Culture Collection을 참고하자.

KCTC에서 현재 분양 불가인 균주라 하더라도 이미 시퀀싱된 자료를 공개하는 것에는 문제가 없을 것이라고 생각한다.

2013년 시퀀싱할 균주를 선정하던 당시에 붙여진 이름이 지금은 더 이상 유효하지 않은 것들이 꽤 있다. 예를 들어 KCTC 22662T는 당시에는 Formosa spongicola였지만 지금은 Xanthomarina spongicola가 되었다. KCTC 웹사이트에는 이전 이름이 synonym임을 보여주고 있고, LPSN에서도 마찬가지이다.

• (LPSN) https://lpsn.dsmz.de/species/xanthomarina-spongicola Taxonomic status: correct name
• (LPSN) https://lpsn.dsmz.de/species/formosa-spongicola Taxonomic status: synonym ← correct name에 대한 링크 포함

재분류의 근거가 된 논문까지 찾아 놓을 필요는 없다. 어차피 LPSN에서 찾을 수 있기 때문이다. Correct name은 우선 KCTC 번호로 조회하고, LPSN에서 재확인한다. KCTC에서 검색이 불가능한 것은 다른 균주 컬렉션 번호를 이용한다. 그 옆에는 phyloFlash의 결과를 놓는다. phyloFlash는 16S rRNA에 기반한 것이므로, SILVA DB의 hit가 표준 균주의 것이 아니라 해도 별 상관은 없다.

• ryzen server(22.04.5 LTS (Jammy Jellyfish))
• BASE_DIR=/data/project/52_KCTC_72_microbial_genomes_2014_Apr/

아무런 전처리를 하지 않은 read에 대하여 사전 점검을 실시한다.

interleaved fastq file로 전환한 뒤 jellyfish(k-mer = 21)에서 분석을 해 둔 히스토그램이 있어서 이를 plot하였다. KAT는 GC content를 구별할 수 있는 plot을 생성해 주는 좋은 점이 있지만, jellyfish의 카운트 결과를 gnuplot으로 그리는 단순한 방법이 나는 더 좋다. 72개의 plot은 다음 명령어를 이용하여 하나로 합쳤다.

$ montage -mode concatenate -tile 7x11 `sort -V <(ls *png)` AllPlots.png

이것만 봐도 도저히 사용하기가 어려운 샘플이 몇 개 보인다. 오염이 명백한 것도 있고(두번째 줄 두번째 샘플), 도대체 뭘 시퀀싱했는지 알 수 없는 것(맨 아래줄 첫번째 샘플)도 있다. 이런 비정상 샘플의 phyloFlash 결과가 기대된다.

Whole-genome shotgun read를 전부 이용하여 de novo assembly를 하지 않고서도 reference DB를 이용, 16S rRNA sequence를 동정해내는 유틸리티인 phyloFlash를 나는 아주 좋아한다. conda pf environment에 설치가 되어 있으니 ${BASE_DIR}/00_DB에 디렉토리에는 pre-formatted database(SILVA 138.1)를 설치하도록 하자. 각 샘플에 대해 다음과 같은 커맨드 라인을 포함하는 스크립트(RUN.sh)를 만들어서 nohup으로 실행하였다. — Haeyoung Jeong 2025/01/07 12:56 오후 3시 32분에 종료되었다.

(작업장) ${BASE_DIR}/04_phyloFlash # 모든 fastq file의 심볼릭 링크를 여기에 생성함 
phyloFlash.pl -dbhome ../00_DB/138.1 -lib 12221T -CPUs 16 -read1 12221T_1.fastq -read2 12221T_2.fastq -almosteverything
phyloFlash.pl -dbhome ../00_DB/138.1 -lib 12304T -CPUs 16 -read1 12304T_1.fastq -read2 12304T_2.fastq -almosteverything
...

결과 파일은 전부 열어서 full length 16S rRNA sequence > 1인 샘플을 전부 찾아서 아래의 ZGA assembly → CheckM 결과와 같이 종합하였다.

Tosten Seemann의 Shovill과 ZGA 중에서 저울질을 하다가 후자를 택하였다. 이는 이미 conda zga environment로 설치가 되어 있으며, 오늘(2025-01-07) 업데이트하였다. Assembler로는 default(unicycler v0.5.0)를 사용하게 하였다. ZGA는 CheckM에 의한 genome completness/contamination/heterogeneity 점검을 해 주므로 매우 유용하다. k-mer analysis로는 오염 정도를 정성적으로 파악할 수 있지만, CheckM은 이를 수치로 나타내어 준다. Update on the proposed minimal standards for the use of genome data for the taxonomy of prokaryotes(IJSEM 2024)에 의하면, completeness > 90%, contamination< 5%를 high quality genome으로 간주한다. zga 옵션 중에서 --minimum-contig-length ###은 genome assembly 결과물이 아니라 annotation으로 진행되는 서열에 대한 한계치이다. 따라서 *.fasta 결과물에는 짧은 contig가 그대로 있으니 주의가 필요하다. 길이로 정렬되어 있는 것은 매우 다행이다.

(작업장) ${BASE_DIR}/05_zga_assembly
awk '{print $2}' ../sampleNameTable.txt | while read f; 
do 
zga -1 ${f}_1.fastq -2 ${f}_2.fastq --calculate-genome-size --threads 16 --minimum-contig-length 200 -o zga_${f}
done

각 샘플에 대한 zga.log 파일을 확인한 결과 다음의 6개는 low quality로 판명되었다. 괄호 안의 수치는 (% completeness, % contamination, % heterogeneity)이고, 뒤쪽의 species name은 phyloFlash에서 확인된 full-length assembled 16S rRNA sequence이다. 굵은 글씨는 샘플에 원래 붙어있던 이름이다.

1. WARNING: 25222T (100.0, 100.69, 0.0) Clostridium subterminale, Succinivibrio dextrinosolvens
2. WARNING: 3520T (100.0, 13.79, 0.0) Brochothrix thermosphacta, Moraxella osloensis
3. WARNING: DSM1535T [Methanobacterium formicicum] (100.0, 96.55, 2.94) Clostridium botulinum, Paraclostridium benzoelyticum
4. WARNING: JCM15547T [(84.58, 20.25, 80.95) Sunxiuqinia faeciviva
5. WARNING: KIM3 (26.06, 2.04, 0.0) Methanobrevibacter smithii
6. WARNING: Strain15 [Geodermatophilus obscurus] (20.97, 0.0, 0.0)
7. NO WARNING: 15166 [Bacillus sp.] (98.28, 8.62, 20.0) Bacillus velezensis, Bacillus gottheilii
8. NO WARNING: K16 [Catabacter hongkongensis] (96.55, 0.0, 0.0) uncultured Anaerotruncus sp., Flavonifractor plautii, Paraclostridium benzoelyticum

Klenkia brasiliensis: Genome-Scale Data Call for a Taxonomic Rearrangement of Geodermatophilaceae 논문(2017)에서 Geodermatophilus brasiliensis를 Klenkia brasiliensis로 재분류해야 한다고 제안했다. 그러나 이 두 가지 species name은 LPSN에서 전부 유효한 이름이라고 나오고, 전자의 경우 여전히 correct name이라고 되어 있다. 결론적으로 말해서 1, 2는 오염, 3은 라벨 실수+오염, 4는 복잡하고, 5와 6은 시퀀싱 분량이 부족하다는 뜻이 된다. 7, 8은 ZGA log에서는 WARNING이 없었지만 phyloFlash에서는 2개 이상의 16S rRNA가 확인된 것이다. k-mer analysis에서 peak의 분포가 이상한 것은 2, 3, 4, 6, 7, 8이다.

Read에 대한 k-mer analysis로는 이렇게 상세한 결과를 알기 어렵다. High coverage peak가 매우 낮은 것은 genome completenss가 낮다. 다시 말해서 sequencing coverage가 부족함을 시사한다. 그러나 k-mer analysis plot으로는 대략적인 느낌을 받을 수 있을 뿐이다. 어쨌거나 read에 대한 점검 없이 그냥 ZGA를 시작했어도 된다는 뜻이 되겠다. 만약 shovill로 조립을 했다면 이만큼의 QC 결과를 얻지 못했을 것이다.

명백히 문제가 있는 샘플의 sequencing raw data를 등록하는 것이 옳은가? 문제점을 포함해서 등록한다면 반면교사의 사례로 쓰일 수 있을 것이다.

각 샘플마다 별도의 디렉토리에 보관된 zga.log 파일을 꺼내서 샘플명을 붙여서 파일이름을 일괄 변경하였다(zga_12221T/zga.log → zga_12111T.log). 나중에 분석을 하기 위함이다. 이를 위해서 다음을 실행하였다.

find . -name zga.log | while read f
do
cp $f $(echo $f | cut -d/ -f 2).log
done

이렇게 만든 72개 파일에 대해 온갖 꼼수를 동원하여 genome completeness/contamination/heterogeneity 수치를 뽑아낸 뒤 하나의 파일(CheckM_summary.txt)에 담았다. 처음 몇 줄을 확인해 보자.

3539T   100.0   2.87    0.0
3572T   98.28   0.0     0.0
3579T   98.28   0.0     0.0
...

Use the NCBI Toolkit to improve data quality: Foreign Contamination Screening (FCS) tool and Assembly QC service.

mamba를 이용하여 gtdbtk-2.1.1 환경을 구축하였다. 설치 요령은 Bioconda 문서를 참고하여 진행하였다. download-db.sh 스크립트로 가져온 reference database는 설치하고자 하였으나 miniconda 디렉토리를 target으로 하고 있어서 충분한 공간이 없었다. 따라서 스크립트를 수정하여 /data/DB/db로 설정한 뒤 GTDBTK_DATA_PATH 변수가 여기를 가리키도록 하였다. db 디렉토리는 이미 존재해서는 안 된다. Reference database의 총 용량은 66G이다. Third party software를 설치하는 것을 잊지 말라.

ZGA가 조립한 FASTA file은 전부 같은 파일명이므로 이를 수정해야 한다. ${BASE_DIR}/06_GTDB-Tk에서 다음을 실행한다.

find ../05_zga_assembly/ -name assembly.fasta | while read f
do 
cp $f $(echo $f | cut -d/ -f 3 | sed s/zga_//).fasta
done

테스트 실행이 잘 되지 않아서 관련 글을 찾아보니 conda 환경을 새로 설치하는 것이 좋다고 하였다. 아예 gtdb라는 계정을 새로 만들어서 conda를 새로 설치하고 gtdb-tk02.1.1 환경을 새로 만들었다. GTDB-Tk 실행에 필요한 3rd party software는 skani(x86-64 Linux statically compiled executable을 구해서 /usr/bin에 직접 복사)를 제외하면 전부 apt를 이용하여 설치하였다.

Fresh conda 환경 설치는 Miniforge에 나온 방법대로 다음과 같이 하였다.

curl -L -O "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh"
bash Miniforge3-$(uname)-$(uname -m).sh # interactive installation

설치 위치는 일반적인 $HOME 하위가 아니라 /data/gtdb/miniforge3로 하였다. 홈 파티션에는 공간이 별로 남지 않았기 때문이다.

(base) gtdb@ryzen-5950x:/data/gtdb/miniforge3$ conda info --envs

# conda environments:
#
base                 * /data/gtdb/miniforge3
gtdbtk-2.1.1           /data/gtdb/miniforge3/envs/gtdbtk-2.1.1

(base) gtdb@ryzen-5950x:/data/gtdb/miniforge3$ mamba info --envs

          mamba version : 1.5.12
# conda environments:  
#
base                 * /data/gtdb/miniforge3
gtdbtk-2.1.1           /data/gtdb/miniforge3/envs/gtdbtk-2.1.1

GTDB-Tk reference data는 conda package에 포함된 download-db.sh를 이용하여 다운로드하고 있는데, r207_v2(파일명: gtdbtk_r207_v2_data.tar.gz, 63G)로서 2025년 1월 13일 현재의 최종 버전(R220)은 아니다. 다운로드 스크립트 실행 뒤 다음과 같은 메시지가 나온다.

[INFO] - The GTDB-Tk database has been successfully downloaded and extracted.
To make your changes take effect please reactivate your environment
[INFO] - Added GTDBTK_DATA_PATH (/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/share/gtdbtk-2.1.1/db) to the GTDB-Tk conda environment.

어느 설정 파일을 건드렸는지 모르겠지만 GTDBTK_DATA_PATH 환경변수가 잘 설정되었다고 한다. 설치가 끝났으므로 Example에 따라서 테스트를 해 보았으나 numpy 관련 에러가 발생한다. numpy 버전 문제인가? STD ERR 메시지를 보자.

(gtdbtk-2.1.1) gtdb@ryzen-5950x:/tmp/gtdbtk/genomes$ gtdbtk identify --genome_dir /tmp/gtdbtk/genomes --out_dir /tmp/gtdbtk/identify --extension gz --cpus 6
[2025-01-13 15:40:27] INFO: GTDB-Tk v2.1.1
[2025-01-13 15:40:27] INFO: gtdbtk identify --genome_dir /tmp/gtdbtk/genomes --out_dir /tmp/gtdbtk/identify --extension gz --cpus 6
[2025-01-13 15:40:27] INFO: Using GTDB-Tk reference data version r207: /data/gtdb/miniforge3/envs/gtdbtk-2.1.1/share/gtdbtk-2.1.1/db
[2025-01-13 15:40:27] INFO: Identifying markers in 2 genomes with 6 threads.
[2025-01-13 15:40:27] TASK: Running Prodigal V2.6.3 to identify genes.
==> Processed 0/2 genomes (0%) |               | [?genome/s, ETA ?]/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/gtdbtk/biolib_lite/prodigal_biolib.py:378: FutureWarning: In the future `np.bool` will be defined as the corresponding NumPy scalar.
coding_base_mask = np.zeros(self.last_coding_base[seq_id], dtype=np.bool)
Process Process-2:1:1:
Traceback (most recent call last):
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/gtdbtk/biolib_lite/parallel.py", line 107, in __producer
    rtn = producer_callback(dataItem)
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/gtdbtk/biolib_lite/prodigal_biolib.py", line 147, in _producer
    prodigalParser = ProdigalGeneFeatureParser(gff_file_tmp)
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/gtdbtk/biolib_lite/prodigal_biolib.py", line 322, in __init__
    self.coding_base_masks[seq_id] = self.__build_coding_base_mask(
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/gtdbtk/biolib_lite/prodigal_biolib.py", line 378, in __build_coding_base_mask
    coding_base_mask = np.zeros(self.last_coding_base[seq_id], dtype=np.bool)
  File "/data/gtdb/miniforge3/envs/gtdbtk-2.1.1/lib/python3.8/site-packages/numpy/__init__.py", line 305, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'bool'.
`np.bool` was a deprecated alias for the builtin `bool`. To avoid this error in existing code, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
..

흠, numpy 버전을 1.20보다 낮은 것으로 바꾸어야 하는 모양이다.

$ pip show numpy
Name: numpy
Version: 1.24.4

GTDB-Tk 설치 관련 문서에 나오듯이 numpy를 1.19.5(numpy >= 1.9.0라고 하지만…)로 내려 보았다. numpy를 지우고 다시 설치하는 데에는 다음과 같이 mamba를 사용하였다.

(gtdbtk-2.1.1) gtdb@ryzen-5950x:/tmp/gtdbtk/genomes$ mamba install -c conda-forge numpy=1.19.5

다시 한 번 Example에 있는 그대로 샘플 파일 2종을 다운로드한 뒤 테스트 실행을 하였다. 성공이다!

(gtdbtk-2.1.1) gtdb@ryzen-5950x:/tmp/gtdbtk/genomes$ gtdbtk identify --genome_dir /tmp/gtdbtk/genomes --out_dir /tmp/gtdbtk/identify --extension gz --cpus 6
[2025-01-13 15:44:30] INFO: GTDB-Tk v2.1.1
[2025-01-13 15:44:30] INFO: gtdbtk identify --genome_dir /tmp/gtdbtk/genomes --out_dir /tmp/gtdbtk/identify --extension gz --cpus 6
[2025-01-13 15:44:30] INFO: Using GTDB-Tk reference data version r207: /data/gtdb/miniforge3/envs/gtdbtk-2.1.1/share/gtdbtk-2.1.1/db
[2025-01-13 15:44:30] INFO: Identifying markers in 2 genomes with 6 threads.
[2025-01-13 15:44:30] TASK: Running Prodigal V2.6.3 to identify genes.
[2025-01-13 15:44:34] INFO: Completed 2 genomes in 4.57 seconds (2.29 seconds/genome).
[2025-01-13 15:44:34] TASK: Identifying TIGRFAM protein families.            
[2025-01-13 15:44:36] INFO: Completed 2 genomes in 1.65 seconds (1.21 genomes/second).
[2025-01-13 15:44:36] TASK: Identifying Pfam protein families.               
[2025-01-13 15:44:36] INFO: Completed 2 genomes in 0.15 seconds (13.63 genomes/second).
[2025-01-13 15:44:36] INFO: Annotations done using HMMER 3.1b2 (February 2015).
[2025-01-13 15:44:36] TASK: Summarising identified marker genes.
[2025-01-13 15:44:36] INFO: Completed 2 genomes in 0.02 seconds (91.82 genomes/second).
[2025-01-13 15:44:36] INFO: Done.

gtdbtk identiy에 이은 나머지 과정(align, classify)도 전부 무사히 종료하였다. 결과물의 디렉토리 체계를 다음에서 확인해 보자. /tmp/gtdbk/genomes는 분석할 FASTA file을 모아 두었던 곳이다.

$ ls /tmp/gtdbk
align  classify  genomes  identify

Docker를 쓰는 경우 reference data는 최신 버전을 쓰면 되는가? 알아보도록 하자. Bioconda를 사용하여 GTDB-Tk를 설치하면 최신 DB와 맞지는 않는다. 이것은 나중에 다시 테스트해 보자.

/data/gtdb/01_20250113/genomes에 모든 fasta file을 복사하여 gzip으로 압축한 뒤 다음의 명령을 차례로 수행하였다.

gtdbtk identify --genome_dir /data/gtdb/01_20250113/genomes --out_dir /data/gtdb/01_20250113/identify --extension gz --cpus 16
gtdbtk align --identify_dir /data/gtdb/01_20250113/identify --out_dir /data/gtdb/01_20250113/align --cpus 16
gtdbtk classify --genome_dir /data/gtdb/01_20250113/genomes --align_dir /data/gtdb/01_20250113/align --out_dir /data/gtdb/01_20250113/classify -x gz --cpus 16

결과 파일은 /data/gtdb/01_20250113/classify 디렉토리의 두 요약 파일(gtdbtk.ar53.summary.tsv, and gtdbtk.bac120.summary.tsv)이다. 요약 파일의 설명은 여기를 참고하라. 특별히 유의하여 볼 컬럼은 2. classification, 14. classification_method, 15. note, 20. warnings이다.

3520T, 15166, 25222T, DSM1535T, JCM15547T, K16, KIM3 및 Strain15의 8개 샘플은 ZGA(CheckM) 점검 결과 completeness나 contamination이 부적합하여 도저히 사용하기 어려운 것으로 판명하였다. JCM15547T를 제외한 contaminated sample은 phyloFlash 분석에서 복수의 complete 16S rRNA sequence가 발견되었다. GTDB-Tk가 'Genome not assigned to closest species as it falls outside its pre-defined ANI radius'라고 경고를 한 것은 새로운 종이라고 볼 수도 있다. 16S rRNA 염기서열 말고는 아무런 정보가 없던 Bacterium NLAE-zl-H470라는 샘플이 Blautia producta에 속한다는 것을 알게 된 것도 성과이다.

그러나 아직 심각한 문제가 남았다. 샘플에 붙어 있었던 이름과 GTDB-Tk가 classify한 것이 판이하게 다른 것이 몇 개 보인다. GTDB reference DB가 correct name을 최신 상태로 유지하고 있는 것은 아니니 조금만 더 신경을 써서 결과를 들여다 보면 일치도를 더욱 높일 수 있다. 그런데 일부 표준균주 샘플에서 불일치가 있다는 것은 문제가 심각하다.

샘플 라벨과 관련한 오류는 두 가지로 나누어 볼 수 있다. 먼저 correctable labeling error이다. 손으로 8을 썼지만 이를 눈으로 확인하고 옮겨 적는 과정(타이핑도 마찬가지)에서 6이나 3이라고 적는 것을 말한다. fixable labeling mistake라고도 할 수 있겠다. 또 하나의 부류는 label-sample mismatch이다. 이것은 복원 불가능하다. 생물자원은행의 신뢰도를 현격하게 떨어뜨릴 수 있는 일이 되므로 주의가 필요하다.

KCTC 3076T(윽, 8076T로 잘못 이름이 붙었다!)는 Gordonia bronchialis라는 correct name이 붙어 있으며, TYGS에는 표준 균주에 대한 유전체 정보가 존재한다. 그런데 phyloFlash 결과로는 Thermoactinomyces vulgaris가 나왔고, GTDB-Tk 결과는 다음과 같다. 두 결과 모두 시퀀싱한 결과는 Thermoactinomyces vulgaris임을 나타내고 있으며, 원래의 라벨인 Gordonia bronchialis와는 아무런 관계가 없다.

• (GTDB-Tk) d__Bacteria; p__Firmicutes; c__Bacilli; o__Thermoactinomycetales; f__Thermoactinomycetaceae; g__Thermoactinomyces; s__Thermoactinomyces vulgaris
• (NCBI taxonomy: TaxID=2054) cellular organisms; Bacteria; Bacillati; Actinomycetota; Actinomycetes; Mycobacteriales; Gordoniaceae; Gordonia

말이 되지 않는다. 가장 쉬운 설명은 엉뚱한 균주를 시퀀싱했다는 것이 된다. KCTC가 10년 전에 보유한 stock이 잘못되었는지, 분양받은 이후에 다른 샘플과 바뀌었는지는 알기 어렵다.

에혀… TYGS로 분석하니 Thermoactinomyces vulgaris DSM 54016T와 ANI 기준으로 동일하며, 이는 KCTC 9076이다. 9076이라 손으로 적은 것을 3076으로, 이를 다시 8076으로 잘못 썼구나! 3843T도 2843(Marinococcus halophilus)T를 잘못 적은 것으로 보인다. 16S와 GTDB-Tk 분석 결과는 시퀀싱 샘플이 Marinococcus halophilus임을 보여준다.

ATCC 25759T라는 표지가 붙은 시퀀싱 결과는 이미 공개된 표준균주의 것과 판이하게 다른 분석 결과(16S, GTDB-Tk)가 나왔다. kevinchjp@gmail.com의 ANI.pl 스크립트를 사용하여 GCF_900103615.1(Romboutsia lituseburensis ATCC 25759T)와 비교해 보니 66.66%에 불과하다. 중간단계 어디에선가 샘플이 뒤바뀐 것인데, 보다 심각하게는 원본 자체를 의심할 수도 있다.

시퀀싱을 해 두고 오랜 기간 동안 결과를 정리하여 발표하지 않아서 외부에서 동일한 표준균주의 유전체를 해독하여 등록한 사례가 많아졌다. 그런데 지금 성적표를 받아들고 생각을 해 보니 만약 시퀀싱 직후에 유전체 조립 결과를 그대로 발표하였다면 잘못 붙은 라벨을 그대로 내보냈을지도 모르는 일이다.

• 72개 샘플로 시퀀싱을 시작하였다. 이 중에서 10개는 KCTC의 공식적인 번호를 받지 못하였다. 공식 KCTC 균주라고 해서 2025년 1월 현재 전부 분양 가능한 것은 아니다. KCTC 균주가 아닌 다음의 3건, 즉 HR7(GCA_000773685.1), HR18(GCA_000773675.1), 32234(M12-1181)는 KRIBB 논문으로 발표되었다. 이 과정에서는 확실히 본 프로젝트의 raw data가 쓰였다. 세번째 논문에서는 표준균주인 M12-1144T의 유전체를 등록한 것으로 되어 있으나(GCA_001895205.1, submitted by Chang, Y.H.), 72 prokaryotic genome 프로젝트와는 무관하다.
• 논문으로 발표된 것 3개 외에도 3810T(GCF_001650325.1, submitted by Lim S.)와 33142T(GCA_001742425.1, submitted by Lim,S. and Kim,B.-C), 43059T(GCA_004916975.1)는 KRIBB에서 유전체 염기서열을 등록한 것으로 되어 있다. 이상의 조립물이 본 프로젝트의 일루미나 시퀀싱 raw data를 사용한 것이었는지는 아직 확인하지 못하였다. 최소한 43059T의 (KRIBB 주저자 논문)에서 마크로젠을 통해 시퀀싱을 했다고 밝혔으므로 72 prokaryotic genomes project와는 무관하다.
• 8개 샘플은 오염 또는 충분하지 않은 sequencing coverage로 인하여 조립 결과물을 등록하기 곤란하다. Sequencing raw data는 잘못된 사례로서 등록할 수도 있겠으나, 이 경우 KRA data에 꼭 수반되어야 하는 BioSample 정보를 입력하기가 곤란하다. 차라리 GeNA로 등록하는 것이 나을 수도 있다. 여기에 해당하는 샘플 ID는 3520T, JCM15547T, 25222T, DSM1535T, 15166, KIM3, L16, Strain15의 8개이다.
• 최소한 세 개의 샘플은 손으로 적은 균주 번호를 잘못 옮겨 적은 것 같다. 3843T(2843T가 맞을 것 같음), 8076T(3076T로 적은 곳도 있으나 9076T가 맞을 것 같음), 8738T(3738T가 맞는 것 같음)가 여기에 해당한다.
• $$$$(ATCC $$$$$$, DSM $$$$$ 병기)는 KCTC에 존재하지 않는 균주임. 16S + GTDB-Tk 분석 결과는 ATCC $$$$$$, DSM $$$$와 일치함
• Irrelevant original labels - 8개 정도의 샘플은 원본 이름과 16S + GTDB-Tk 분석 결과가 너무나 다르다. 이 경우에서도 16S 결과와 GTDB-Tk 결과는 일치한다. 오타 정도가 아니라 아예 샘플 자체가 잘못된 것으로 보인다. 여기에 해당하는 것은 3761T, 3956T, 15258T, DSM16643T, DSM16632T, ATCC25759T, S2-Ab3)이다. 2025년 2월 10일, 여기에 하나가 더 추가될 것으로 보인다. 3449T도 이상하다.
• K10(Bacterium NLAE-zl-H470, 16S rRNA 염기서열과 contig를 비교하면 100% 일치)의 동정이 이번에 제대로 된 것은 그나마 긍정적인 성과이다.
• 신규 해독 유전체로서 의미를 부여할 수 있는 것은 23669, 3829T, 33143, 43239, MAE15, 15385, S2-Ab3, K10, I-14 정도이다. 그러나 일부는 genus 및 species를 특정하기 어려운 문제가 있다.

K-BDS에 'Examples of low quality Illumina sequencing of prokaryotic genomes'라는 바이오프로젝트KAP241424로 등록하였다. NGS sequencing raw data만 등록한다면 KRA 데이터 타입이 맞겠지만, (1) 일부 오염된 샘플 때문에 의미 있는 BioSample을 정의하기 어렵고, (2) fastq file뿐만 아니라 QC 결과와 조립 결과물 등 다양한 파일을 서브디렉토리별로 나누어 담은 구조이기 때문에 기타 데이터 타입으로 등록하는 것이 가장 적당하다고 생각하였다. 등록한 파일에 대한 상세한 설명 및 분석(QC) 과정은 README 파일에 상세하게 설명하였다.

가장 심각한 문제이다. 배양과 시퀀싱을 포함한 실험 과정의 문제가 아니라 균주 stock 자체에 라벨이 잘못 매겨져 있을 가능성이 큼을 시사한다. 해외로부터 잘못된 균주를 기탁받았을 가능성도 배제하기 어렵다. 바로 위에서 소개한 8개 low quality genomes에서도 이런 사례를 포함하고 있다. 공개된 type strain의 유전체와 비교해 보려고 한다. 이를 위해서 다음에 설명한 방법에 따라 2025년 1월 17일에 GenBank로부터 관련 유전체를 전부 다운로드하였다. 표준 균주의 유전체가 여러개인 경우도 전부 받았다.

이미 공개된 유전체 정보를 다운로드하여 72 prokaryotics genomes 프로젝트의 것과 비교해 볼 수 있다. 쉽게 말해서 '어느 것이 더 좋았다'라는 말을 할 수 있으면 좋을 것이다. 그러면 무엇을 비교 지표로 삼을 것인가? 단순한 assembly 완성도 정도? 또는 CheckM 지표? 이에 대한 고민이 필요하다. getSequenceInfo가 제공하는 NucleScore를 써 보는 것은 어떨까? 이 지표는 오로지 genome sequence에서 나오는 metric(GC content, AT/GC ratio, total sequence length)만을 이용하여 만들어진다.

특정 species의 type strain에 대한 유전체 정보를 다운로드하려면 다음과 같이 실행한다. 실제로 파일을 받으려면 --dry-run 옵션을 빼야 한다. -A <something> 옵션의 기능은 매우 독특하다. assembly accession, 콤마로 구분한 assembly accession은 물론이고 한 줄에 하나씩의 assembly accession을 넣은 파일을 주어도 된다. archaea는 별도로 실행해야 한다.

$ ncbi-genome-download [--dry-run] -s genbank --formats fasta --flat-output -o outdir --genera "Winogradskyella thalassocola" -M type bacteria
Considering the following 1 assemblies for download:
GCF_900099995.1	Winogradskyella thalassocola	DSM 15363

본 프로젝트와 관련된 모든 종의 유전체 정보(표준 균주에 대한 중복된 시퀀싱자료 포함)등은 List of public genomes for comparison에 정리하였다.

다음의 일곱 샘플은 표준균주라는 이름을 달고 시퀀싱을 한 것이다. 시간이 많이 흐르면서 다른 연구자들이 표준균주에 대한 유전체 정보를 해독하여 등록하게 되었고, 이것과 비교를 해 보니 애초에 붙어 있던 라벨(이름)을 믿을 수 없는 것들이다. 동일 균주라 해도 Culture collection에서 얼마나 배양을 했는지에 따라 미소한 차이는 있을 수 있지만, 이하의 것은 뭔가 크게 잘못된 것이다. 만약 KCTC에서 균주 번호를 잘못 붙인 상태로 분양한 것인지, 혹은 분양 후에 실험자의 실수가 있었는지는 알기 어렵다. 이 결과는 GenBank에서 original label에 해당하는 종의 표준균주 유전체 염기서열을 다운로드하여 ANI.pl로 단방향 분석한 것이다. 유전체 정보에서 나온 16S rRNA classification(phyloFlash)와 GTDB-Tk의 결과는 당연한 이야기지만 잘 일치한다.

ANI.pl은 다음의 명령어로 계산하였다. Legacy blastall을 사용하는 오래전 펄 스크립트인데 쓰기가 편해서 아직도 애용한다.

ANI.pl --fd formatdb --bl blastall \
       --qr qr/3447T.fasta --sb sb/GCA_000307855.1.fna \
         --od output

사실 ANI를 계산하는 것이 의미가 없다. Genome size, %G+C가 완전히 다르기 때문이다.

• 3449T(Priestia flexa) vs. GCA_001591565.1 (GTDB-Tk: Pseudomonas_E canadensis)
• 3761T(Paenibacillus kobensis) vs. GCA_000525915.1 - ANI: 64.8220233463035 (GTDB-Tk: Aerococcus urinaeequi)
• 3956T(Paenibacillus alkaliterrae) vs. GCA_021532375.1 - ANI: 62.9229519450801 (GTDB-Tk: Niallia sp001076885)
• 15258T(Clostridium kogasensis) - 이 균주에서 등록된 16S rRNA sequence(JQ423944)와도 맞지 않음; identity가 86%에 불과함
• DSM16643T(Methanobrevibacter millerae) vs. GCA_900103415.1 - ANI: 61.1466666666667 (GTDB-Tk: Clostridium_J sp002341865)
• DSM16632T(Methanobrevibacter olleyae) vs. GCA_900114585.1 - ANI: 61.0463157894737 (GTDB-Tk: Clostridium_J sp002341865)
• ATCC25759T(Romboutsia lituseburensis) vs. GCA_900002825.1|GCA_900103615.1 - ANI: 66.7831618759455|66.7831618759455 (GTDB-Tk: Clostridium_F sp001276215)
  • [R] GCA_900002825.1 vs. GCA_900103615.1 - ANI: 99.9876330076003

다음의 샘플도 매우 수상쩍다. 균주 분리 당시 16S rRNA에 의해서 기초적인 동정을 했다 하더라도 Clostridium이 Microbacterium으로 바뀔 수는 없다. 두 종류의 미생물은 Phylum level에서부터 차이가 난다.

• S2-Ab3(Clostridium sp.): 적당한 type strain 없음 (GTDB-Tk: Microbacterium sp001595495)

다음은 시퀀싱 당시에 정확하게 동정되지 않았던 것이 이번 기회에 제 이름을 찾게 된 것으로 보면 된다. 그러나 만일 이것도 중간에 뒤바뀌었거나 라벨을 잘못 한 것이라면? 당시에 붙인 라벨과 유전체 기반 분류에 어느 정도는 일맥상통하는 면이 있어서 온전한 샘플이라고 믿고 싶다.

• YHPB24(Bifidobacterium sp.) vs. GCA_002259605.1 - ANI: 95.5372918918919 (GTDB-Tk: Pseudoscardovia radai)
• K10(Bacterium NLAE-zl-H470) vs. GCA_004340925.1|GCA_034355335.1|GCA_900461125.1 - ANI: 98.4443814844362|98.4360207543393|98.4362849162 (GTDB-Tk: Blautia coccoides)
• I-14(Roseicyclus sp.): Marivivens sp. 적당한 type strain 없음 (GTDB-Tk: genus Marivivens 하의 신종)

이 범주의 것은 염기서열 리포지토리(예: K-BDS)에 등록하기에 적당한 것이다. 괄호 안의 것은 original label이다. 아직 공개되지 아니한 표준균주의 유전체서열은 특히 의미가 깊다. 하지만 고민이 필요하다. 만약 GenBank에 등록한다면 RefSeq등 여러 경로를 통해 계속 전파될 것이다. 그러나 K-BDS에만 등록하면 그렇게 될 가능성이 매우 낮다. 또한 표준균주가 아닌 것은 어디에서 균주를 입수할 수 있는지 명확하게 밝히기가 어렵다.

1. 3829T (Paenibacillus azoreducens) - 아직 공개되지 않은 표준균주의 유전체라서 학술적 가치 있음. 많은 고민 끝에 2025년 2월 10일 NCBI에 제출하였으며(PRJNA1221642 GCA_047779985.1 'Paenibacillus azoreducens KCTC 3829') 2월 15일에 공개되었다. 권고에 따라 SRA에도 raw data file을 등록하였음. 만약 KRA에 등록한 것인 브로커링을 통해 DDBJ로 넘어가면, 결국 데이터 중복이 될 것인가?
2. 33143 (Bacillus solimangrovi GH2-5)
3. 43239 (Paenibacillus dokdonensis JAE2; 표준균주는 JAE5=43059T) - genome 기준으로는 같은 종이 아님
4. K10 (동정되지 않음) - Blautia coccoides
5. 23669 (Formosa sp. KOPRI 21329)- Xanthomarina sp. 극지연구소 KOPRI 21329(극지연구소에서는 더 이상 균주은행 기능을 하지 않음). 균주 입수 가능성이 불투명하므로 46개 균주 목록에서는 일단 제외하였음
6. MAE15 (Acinetobacter sp.) - Acinetobacter halotolerance. KCTC에서 공식적으로 분양하는 균주가 아니므로 46개 균주 목록에서는 일단 제외하였음. 그러나 기탁자의 허가를 통해 K-BDS에 등록하기로 함.
7. 15385 (Clostridium sp. TAN35) - 종을 특정할 수 없음. 공식적으로 분양하는 균주가 아니므로 46개 균주 목록에서는 일단 제외하였음. 그러나 기탁자의 허가를 통해 Clostridium composti KCTC 15385로 공개되었rh, K-BDS에 등록하기로 함.
8. I-14 - (Roseicyclus sp.) - GTDB-Tk에 의해 Marivivens sp.로 동정됨
9. YHPB24 (Bifidobacterium sp.) - GTDB-Tk에 의해 Pseudoscardovia radai(Bifidobacteriaceae)로 동정. KCTC에서 공식적으로 분양하는 균주가 아니므로 46개 균주 목록에서는 일단 제외하였음. 그러나 기탁자의 허락에 의해 유전체 정보를 K-BDS에 등록하기로 함.

KRA(raw sequencing reads, NGS)를 먼저 등록한 뒤 KNA(genome assembly)를 등록한다. 각 데이터의 등록 전에는 이에 해당하는 BioProject와 BioSample이 등록되어 있어야 한다. BioProject는 내용만 간단히 기술하고 과제 정보만 넣으면 되지만, 하나의 BioProject에 부속되는 수십 건의 BioSample을 등록하려면 상당한 노력을 들여 자료를 조사한 뒤 입력에도 요령이 필요하다. 글이 너무 길어질 것을 우려하여 BioSample 메타데이터 파일 작성 및 등록은 다음의 페이지에 별도로 설명하였다.

BioSample submission to K-BDS

이를 통해서 47건의 바이오샘플을 등록하였다(BioSample메타데이터 파일). 다음으로는 우선 KRA에 NGS read data를 등록하였다(KRA 메타데이터 파일). 선택된 47개 샘플은 KCTC에서 2025년 2월 현재 공식으로 분양하는 것 45건과 이미 Genome Announcements 논문으로 발표된 Bordetella trematum HR18, 그리고 16S rRNA sequence를 통해 세상에 알려진 NLAE-zl-H470(입수 경로 모름)이다. 김병찬 박사팀과 같이 작업했던 HR7 및 HR18의 논문 본문에서 시퀀싱을 한 기관명을 Korea Research Institute of Bioscience and Bioengineering이라고 잘못 적었다는 것을 10년이 지난 지금 발견하였다. 아마 KAIST의 학과 이름과 혼동했을 것이 틀림이 없다. 나의 실수다! 허허허…

다음으로는 genome assembly를 KNA 데이터베이스에 등록하는 일이 남았다. 3개 샘플의 fastq file이 내부에서 duplicate된 상태임을 발견하여 2025년 2월 7일 아예 47개 샘플 전체의 ZGA assembly를 다시 진행하였다. 뭔 말인가 하면, cat reads_1.fastq reads_1.fastq > renamed-reads_1.fastq의 상태였다는 뜻이다. K-BDS에 등록한 KRA 데이터를 관리자가 검증하던 중, 앞부분의 read ID가 중간 부분부터 다시 출현하는 것을 발견하여 나에게 통지해 주었다.

BioSample은 KRA와 KNA에 대해서 별도로 등록해야 한다. 따라서 KRA에서 썼던 메타데이터 파일에서 샘플 이름에 일괄적으로 특정 문자를 붙인 파일을 마련하여 등록한다. 혹은 '복제'기능을 사용해 본다.

KNA에 대한 메타데이터는 샘플 수가 많은 관계로 역시 메타데이터 파일로 작성하여 업로드하고 싶지만 BioProject와 BioSample까지 확정이 되어야 비로소 웹 입력 폼 또는 메타데이터파일 다운로드 페이지에 접근 가능하다. KAP240814에 딸린 조립 관련 정보를 참조하도록 하자. 다음 명령어를 사용하여 재조립 후 결과 파일을 이름을 바꾸어 한데 모은다. Conda zga environment에서 수행한다. KNA 메타데이터 파일을 만들기 위한 수치 등 자료는 ZGA log file에서 뽑는다(40 줄짜리 ZGA log 파일의 형태).

cat 47_samples.txt | while read f
do
zga -1 02_renamed_fastq_files/${f}_1.fastq -2 02_renamed_fastq_files/${f}_2.fastq --calculate-genome-size --threads 16 --minimum-contig-length 200 -o 09_47_samples_zga_again/zga_${f}
cp 09_47_samples_zga_again/zga_${f}/assembly/assembly.fasta /work/2025_museum/02_KNA/${f}.fasta
cp 09_47_samples_zga_again/zga_${f}/annotation/genome.gff /work/2025_museum/02_KNA/${f}.gff
cp 09_47_samples_zga_again/zga_${f}/zga.log /work/2025_museum/02_KNA/${f}.zga.log
done

47개 KRA 데이터의 등록 후 검수를 거치는 동안 혹시 파일이 뒤바뀌거나 하지 않았는지 의심이 들어서 재조립 후 GTDB-Tk 분석을 다시 실시하였다. 그러던 중에 KCTC 3449T 샘플의 동정 결과가 이상함을 알게 되었다. KCTC의 자료에 의하면 이는 Firmicutes의 일종인 Priestia flexa(syn. Bacillus flexus)로 되어 있다. 그러나 이번 프로젝트의 초기 작업부터 이 샘플은 phyloFlash와 GTDB-Tk 공히 Pseudomonas 계열인 것으로 확인이 되었다. 따라서 가장 나쁜 사례인 bad label로 보아야 하며, 이를 등록 대상에서 제외하기로 하였다. 이에 따라서 K-BDS 관리자에게 요청하여 BioSample과 KRA에서 3449T를 제거해줄 것을 요청하였다. 관리자는 데이터 전체를 반려 상태로 되돌려 주었다. BioSample은 웹 입력 양식에서 수정하였고, KRA는 새로 만든 메타데이터 파일을 업로드하여 재제출하였다.

지금까지 살아남은 46개의 sequencing raw data는 2025년 2월 10일 오후 5시를 기하여 ZGA pipeline을 아예 다시 돌렸다. 최소 contig 길이를 200 bp로 설정하는 옵션을 주어도 이는 annotation에만 영향을 주고 assembly에는 손을 대지 않는다는 것을 알게 되었다. 따라서 seqkit utility를 사용하여 제출용 assembly를 필터링하였다.

seqkit seq -m 200 IN.FASTA > FILTERED.FASTA

관계자들에게 요청하여 최대한 많은 균주의 유전체 정보를 공개하는 것으로 가닥을 잡아 나가고 있다. 현재 분양 가능성은 중요하지 않다. 10년 전에 균주를 제공하였던 자원은행이 지금 더 이상 서비스를 하지 않는다고 해서 당시 적법하게 구입한 균주의 유전체 연구 성과물을 K-BDS에 등록하지 못하는 것은 말이 되지 않기 때문이다(바이오샘플 56건의 메타데이터 파일, KRA의 메타데이터 파일)

균주를 처음에 제공했던 사람과 KCTC의 협조를 얻어서 56개 정보를 전부 등록하여 공개하는 것으로 합의하였다. 그 과정에서 KCTC에 공식적으로 기탁된 균주 자원이 비공개 상태로 있다가 공개로 전환된 것이 몇 개 있어서 큰 보람을 느낀다.