RNA 시퀀싱(RNA sequencing)은 차세대 유전체 분석 기술(Next Generation Sequencing, NGS)을 이용해 DNA로부터 전사된 모든 RNA 전사체(transcriptome)의 염기서열을 분석하는 기술입니다.

RNA 시퀀싱의 주요 목적 중 하나는 두 가지 이상의 생물학적 조건 간 차등적으로 발현되는 유전자(DEG, Differentially Expressed Gene) 또는 분자 경로를 식별하여 유전자 발현을 프로파일링하는 것입니다.

DNA가 RNA로 전사되고 RNA가 단백질로 번역된다는 기본 원리에 따라, 특정 유전자의 RNA 양이 많을수록 해당 유전자의 발현 수준이 높다고 판단할 수 있습니다. RNA 염기서열 분석을 통해 특정 샘플의 유전자 발현량과 발현 양상을 파악할 수 있습니다.

RNA 시퀀싱은 염기서열 변이, 대체 RNA 스플라이싱(alternative RNA splicing), 유전자 융합(gene fusion), 단일 염기 다형성(Single Nucleotide Polymorphism, SNP)과 같은 정보를 탐지할 수 있습니다. 이를 통해 RNA 전사체 수준에서 발생하는 변이와 발현 양상을 분석하며, 특히 대체 스플라이싱 이소폼과 유전자 융합 이벤트를 확인하는 데 유용합니다. 또한, miRNA, tRNA와 같은 작은 RNA의 발현과 엑손 및 인트론의 경계를 분석하여 전사체 구조와 스플라이싱 패턴을 재구성할 수 있습니다. 다만, SNP 탐지는 RNA 편집이나 기술적 오류를 주의해야 하며, 작은 RNA 분석은 해당 전용 라이브러리 준비 방법이 필요합니다.

01. QC (Quality Control)

NGS는 기술적 한계와 실험적 원인에 의해 다양한 오류(error)가 발생할 가능성이 있습니다.

원시 데이터(raw data)에서는 염기(base)별로 추정된 오류 확률을 수치로 나타내며, 이를 Phred quality score라고 부릅니다. Phred score는 각 염기의 품질을 나타내는 지표로 활용됩니다.

Phred score 20(Q20)은 해당 염기서열 결과가 오류일 확률이 1%임을 의미하며, Q30은 0.1%의 오류 확률을 의미합니다. 일반적으로 Q30 이상의 Phred score를 가진 염기는 시퀀싱 품질이 우수하다고 판단되어 분석에 활용할 수 있습니다.

01-1. Quality Control

FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)

FastQC 프로그램은 NGS 결과로 생성된 FASTQ 파일을 입력받아 데이터의 시퀀싱 품질을 평가하는 소프트웨어로, 염기별 시퀀싱 품질, GC 콘텐츠, 어댑터 서열 등의 여러 품질 지표를 그래프 형태로 제공하여 원시 데이터의 품질을 직관적으로 확인할 수 있게 합니다.

주요 활용되는 지표로 Per Base Sequence Quality, Per Sequence Quality Scores, Per Base Sequence Content, Per Base GC Content 등이 있습니다.

• Input
  • 각 샘플에 대한 Raw data file (forward, reverse read)
  • SRR26528168.sra_1.fastq, SRR26528168.sra_2.fastq
• Practice ### 01. FASTQC #/path/fastqc -o /path/output/dir [INPUT FILE] /root/bin/fastqc \ -t 8 \ --outdir ${RESULT_DIR} \ SRR26528168.sra_1.fastq /root/bin/fastqc \ -t 8 \ --outdir ${RESULT_DIR} \ SRR26528168.sra_2.fastq
• Option description -o / --outdir [DIR] fastQC 결과 파일을 생성할 디렉터리의 경로 입력 -t / --threads [NUM] 분석 수행을 위하여 사용하는 threads의 수 입력
• Output
  • 각 read 에 대한 fastqc 의 html 및 zip 파일 생성됨
  • SRR26528168.sra_1_fastqc.html
  • SRR26528168.sra_1_fastqc.zip
  • SRR26528168.sra_2_fastqc.html
  • SRR26528168.sra_2_fastqc.zip
• 특정 디렉터리 내 모든 fastq 파일에 대한 fastQC 수행 시 아래의 코드를 참고 for x in $(ls ${working_path}/00.raw_data/*.fastq) do /root/bin/fastqc -t 8 -o ${RESULT_DIR} $x done
• 결과 해석 (아래의 PDF 참고)
  FastQC_Manual.pdf

02. Trimming

Trimming은 시퀀싱 데이터 분석에서 품질을 높이기 위한 데이터 정제 과정입니다. 이 과정에서는 각 read 말단에서 나타나는 정확도가 낮은 부분을 잘라내고, read의 정확도 평균 점수가 너무 낮거나 염기서열이 정확히 판별되지 않아 N으로 기록된 위치가 많은 경우 해당 read 전체를 제거합니다.

또한, 라이브러리 제작 과정에서 내부 DNA 서열 길이가 짧은 경우, 시퀀싱 장비가 내부 서열을 모두 읽은 후 NGS 어댑터 서열을 읽게 됩니다. 이 어댑터 서열은 원래 존재하지 않는 서열이므로 분석에 오류를 발생시킬 수 있습니다. 따라서, read 말단에서 어댑터 서열이 인식되면 해당 서열을 제거하는 과정도 필요합니다.

출처 : https://www.celemics.com/ko/resources/blogs/bioinformatics-standard-ngs-data-analysis-pipeline-2/

02-1. Trimming

Trimmomatic (https://github.com/usadellab/Trimmomatic)

어댑터 서열 제거 및 품질 낮은 서열을 정제하는 과정은 원시 데이터를 개선하기 위해 수행됩니다. 이 과정에서는 최소 품질 임계값(threshold)을 설정하고, 해당 임계값 이하의 염기를 제거하여 고품질 데이터를 확보합니다.

Cutadapt는 Python 기반 도구로 어댑터 서열을 효과적으로 제거하는 데 사용됩니다. 반면, Trimmomatic은 Java 기반으로 개발된 도구이며, Sickle과 마찬가지로 슬라이딩 윈도우(sliding window) 알고리즘을 사용하여 고품질 서열만을 남기는 방식을 채택하고 있습니다. 이를 통해 분석의 신뢰도를 높이고, 데이터 품질을 보장합니다.

• Input
  • 각 샘플에 대한 Raw data file (forward, reverse read)
  • SRR26528168.sra_1.fastq, SRR26528168.sra_2.fastq
• Practice ### 02.Trimmomatic #java -jar /path/trimmomatic-[VERSION].jar [SE | PE] [-phred33 | -phred64] [INPUT FILE R1] [INPUT FILE R2] [OUTPUT FILE 1P] [OUTPUT FILE 1U] [OUTPUT FILE 2P] [OUTPUT FILE 2U] [trimmer 1] java -jar ${Tool_path}/trimmomatic-0.39.jar \ PE \ -threads 16 \ -phred33 \ ${RAW_DIR}/SRR26528168.sra_1.fastq \ ${RAW_DIR}/SRR26528168.sra_2.fastq \ ${RESULT_DIR}/SRR26528168.R1.Trim.Paired.fastq.gz \ ${RESULT_DIR}/SRR26528168.R1.Trim.UnPaired.fastq.gz \ ${RESULT_DIR}/SRR26528168.R2.Trim.Paired.fastq.gz \ ${RESULT_DIR}/SRR26528168.R2.Trim.UnPaired.fastq.gz \ ILLUMINACLIP: ${ADAPTER_DIR}/TruSeq3-PE.fa:2:30:10 \ LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
• Option description [SE | PE] 데이터가 single end면 SE, paired end면 PE 선택 -phred33 Phred score 33 ILLUMINACLIP:[FILE]:[NUM]:[NUM]:[NUM] {Adapter 정보가 있는 파일}: {maximum mismatch counts}: {paired end에서 adapter reads간의 일치도}: {single end에서 adapter reads 간의 일치도} LEADING:[NUM] Read의 시작 부분에서 base를 제거할 최소 품질 기준 TRAILING:[NUM] Read의 끝 부분에서 base를 제거할 최소 품질 기준 SLIDING WINDOW:[NUM] {window size}:{평균 quality}, window를 이동시켜 평균 qulity 이하로 내려가는 read 제거 MINLEN:[NUM] 유지되어야 할 최소 reads 길이
• Output
  • 각 샘플별 trimming 이후 결과 파일 (R1, R2에 대한 Paired, UnPaired read)
  • SRR26528168.R1.Trim.Paired.fastq.gz
  • SRR26528168.R1.Trim.UnPaired.fastq.gz
  • SRR26528168.R2.Trim.Paired.fastq.gz
  • SRR26528168.R2.Trim.UnPaired.fastq.gz

03. Alignment

전처리가 완료된 Sequence Read 데이터는 일반적으로 참조 유전체 또는 참조 전사체에 정렬(alignment) 및 매핑(mapping)하는 과정을 거칩니다.

이 과정에서 read를 매핑할 적절한 참조 서열을 선택하고 STAR, HISAT2, Bowtie2, Kallisto, Salmon와 같은 도구를 사용하여 정렬을 수행합니다. 도구 선택은 연구자의 선호도 및 사용 가능한 계산 자원에 따라 달라질 수 있습니다.

이 튜토리얼에서는 STAR를 사용하여 FASTQ 파일을 참조 게놈에 매핑하고, 결과로 정렬된 read 정보를 담은 BAM 형태의 파일 생성합니다.

STAR 매뉴얼: https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf

※ 참조 유전체(Reference genome)는 특정 생물 종의 대표적인 유전체 서열로 해당 생물의 유전자와 유전체 구조를 이해하는 데 기본이 되며, 다양한 유전체 연구 및 비교 연구에서 기준점으로 사용됩니다.

※ 매핑 전, 참조 서열(.fa)과 주석(.gff, .gtf) 파일을 기반으로 인덱싱하는 과정이 필요합니다.

03-1. STAR indexing

• Input
  • 샘플에 해당되는 Reference genome sequence(.fa) and annotation (.gff/.gtf), genome.fa, genes.gtf (00-2. 참조 서열 다운로드 참고)
• Practice ### 03-1. STAR indexing #/path/STAR STAR --runMode genomeGenerate --runThreadN 16 \ --genomeDir ${STAR_INDEX_DIR} \ --genomeFastaFiles genome.fa \ --sjdbGTFfile genes.gtf \ --sjdbOverhang 75
• Option description --runMode genomeGenerate 게놈 인덱싱 생성모드 --runThreadN [NUM] 분석 시 사용할 thread의 개수 --genomeDir [DIR] 인덱싱 결과 폴더 --genomeFastaFiles [FILE] FASTA 형식의 인덱싱할 참조서열 파일 --sjdbGTFfile [FILE] GTF 형식의 참조서열 주석 파일 --sjdbOverhang [NUM] 주석이 달린 접합부 주변의 게놈 서열 길이로 [ReadLength - 1]가 이상적이며, defualt는 100임
  --sjdbOverhang의 숫자는 FastQC 결과 리포트의 Basic Statistics에서 [Sequence length -1] (76 -1 = 75) 값을 입력하면 됨
• Output
  • /path/STAR_INDEX_DIR 내에 인덱싱 결과 파일이 생성됨

03-2. STAR mapping

• Input
  • 각 샘플에 대한 Trimmed raw data file (forward, reverse read) : SRR26528168.R1.Trim.Paired.fastq.gz, SRR26528168.R2.Trim.Paired.fastq.gz
  • 참조 게놈 인덱싱 결과 폴더 : /path/STAR_INDEX_DIR
• Practice ### 03-2.STAR mapping #/path/STAR STAR --runThreadN 16 \ --genomeDir ${STAR_INDEX_DIR} \ --readFilesIn SRR26528168.R1.Trim.Paired.fastq.gz SRR26528168.R2.Trim.Paired.fastq.gz \ --readFilesCommand zcat \ --quantMode TranscriptomeSAM \ --outSAMtype BAM SortedByCoordinate \ --outFileNamePrefix ${RESULT_DIR}/SRR26528168. \ --twopassMode Basic
• Option description --runThreadN [NUM] [NUM]개의 thread를 사용 --genomeDir [DIR] 인덱싱 결과 폴더 --readFilesIn [FILE] FASTQ 형식의 입력 파일 (여러 개 가능) --readFilesCommand zcat 입력 파일이 압축되어있을 경우 사용 --quantMode TranscriptomeSAM 정량된 SAM/BAM 파일을 별도의 파일로 출력 --outSAMtype BAM SortedByCoordinate 좌표기준으로 정렬된 BAM 파일을 출력 --outFilesNamePrefix [PRIXED FILEANME] 매핑 결과 폴더 및 pre-fixed 파일명
• Output
  • SRR26528168.Aligned.sortedByCoord.out.bam
  • SRR26528168.Aligned.toTranscriptome.out.bam
  • SRR26528168.Log.final.out
  • SRR26528168.Log.out
  • SRR26528168.Log.progress.out
• 결과 해석

  기본적으로 Log 파일은 분석 시간/속도 및 input 정보와 같이 분석 매핑 결과를 제공하며, 최종 매핑 결과에 대한 통계 수치는 Log.final.out 파일에서 확인할 수 있습니다.

  Uniquely mapped reads % 는 매핑의 품질을 판단하는 중요한 지표로, 값이 80% 이상이면 매핑 결과가 좋다고 판단하며, 값이 50% 미만이면 시퀀싱 또는 참조 서열에 문제가 있다고 의심할 수 있습니다.

  출처 : https://pmc.ncbi.nlm.nih.gov/articles/PMC4631051/

04. Quantification

Quantification은 매핑된 RNA 시퀀싱 데이터를 기반으로 유전자의 발현 수준을 정량화하는 중요한 과정입니다. 이 과정은 유전자 발현량을 평가하여, 서로 다른 조건에서의 발현 차이를 비교 분석하고, 다양한 생물학적 질문에 답을 찾는 데 기여합니다.

정량화를 수행하기 위해 일반적으로 RSEM 또는 HTSeq와 같은 도구가 사용되며, 본 튜토리얼에서는 RSEM을 활용한 정량화 과정을 다룹니다.

RSEM은 Expectation-Maximization (EM) 알고리즘을 사용하여 각 read가 유전자 또는 전사체에 매핑될 확률을 기반으로 발현량을 추정합니다. 결과는 TPM (Transcripts Per Million) 또는 FPKM (Fragments Per Kilobase Million) 등의 형식으로 제공되어, 조건 간 발현 비교가 가능하도록 정규화된 값을 제공합니다.

※ 정량 전, 참조 서열(.fa)과 주석(.gff, .gtf) 파일을 기반으로 인덱싱하는 과정이 필요합니다.

04-1. RSEM indexing

• Input
  • STAR 분석 후 생성된 Aligned BAM file : SRR26528168.Aligned.toTranscriptome.out.bam
  • 참조 서열 인덱싱 결과 폴더 : /path/RSEM_INDEX_DIR/{ref_prefix}
  • 샘플에 해당되는 Reference genome sequence(.fa) and annotation (.gff/.gtf), genome.fa, genes.gtf (00-2. 참조 서열 다운로드 참고)
• Practice ### 04-1. RSEM indexing #/path/RSEM/rsem-prepare-reference [OPTION] [REFERENCE FASTA FILE] [OUTPUT PREFIX NAME] ${Tool_path}/rsem-prepare-reference \ --num-threads 16 \ --bowtie2 \ --bowtie2-path ${bowtie2_path} \ --gtf reference.gtf \ reference.fa \ ${RSEM_INDEX_DIR}/${ref_prefixed name}
• Option description --num-threads [NUM] 분석 시 사용할 thread의 개수 --gtf [FILE] GTF 형식의 참조서열 주석 파일 --bowtie2 Aligner 도구 --bowtie2-path Aligner가 설치된 위치 ※ Aligner로 bowtie, bowtie2, star, hisat2-hca가 존재하며, 분석 방법에 따라 변경 가능
• Output
  • /path/RSEM_INDEX_DIR 내에 인덱싱 결과 파일이 생성됨

04-2. RSEM quantification

• Input
  • STAR 분석 후 생성된 Aligned BAM file : SRR26528168.Aligned.toTranscriptome.out.bam
  • 참조 게놈 인덱싱 결과 폴더 : /path/RSEM_INDEX_DIR/{ref_prefix}
• Practice ### 04-2.RSEM quantification #/path/RSEM/rsem-calculate-expression [OPTION] [INDEX PREFIX NAME] [INPUT FILE] [OUTPUT PREFIX NAME] ${Tool_path}/rsem-calculate-expression \ --bam --no-bam-output -p 16 --paired-end \ ${RSEM_INDEX_DIR}/${ref_prefixed name} \ SRR26528168.Aligned.toTranscriptome.out.bam\ ${RESULT_DIR}/SRR26528168
• Option description --bam Inputs 파일이 BAM format일 때 사용 --no-bam-output 사용 시 여러 BAM 파일을 output 하지 않음 -p / --num-threads [NUM] 분석 시 사용할 threads 개수 --paired-end 데이터가 paired-end reads 기반일 때 사용
• Output
  • SRR26528168.genes.results
  • SRR26528168.isoforms.results
  Output 예시

  

  • gene_id : 유전자 고유 식별자로 Ensembl ID 또는 RefSeq ID가 사용됩니다.
  • transcript_id(s) : 해당 유전자로부터 유래한 전사체 ID 목록으로 하나의 유전자에 여러 전사체를 가질 수 있기 때문에 전사체 ID가 쉼표로 구분되어 나열됩니다.
  • length : 유전자의 총 길이로 주로 엑손(exon) 부분의 길이로 계산됩니다.
  • effective_length : 유전자 발현량을 추정할 때 사용되는 유전자의 실질적인 길이입니다. 이 값은 유전자 길이에서 정렬 가능한 리드의 길이 및 시퀀싱 방법 등을 고려하여 조정된 값입니다.
  • expected_count : 해당 유전자에 정렬된 리드의 예상 개수입니다. 단순한 리드 수가 아니라, 모델을 통해 추정된 값입니다.
  • TPM (Transcripts Per Million) : 백만 개의 전사체 당 특정 유전자의 전사체 수를 나타냅니다. TPM은 샘플 간 발현량을 비교하기 위한 표준화된 척도입니다. TPM은 유전자 길이와 샘플의 총 리드 수를 고려하여 계산되므로 다른 샘플 간의 발현량 비교가 가능합니다.
  • FPKM (Fragments Per Kilobase of transcript per Million mapped reads) : 유전자 발현량을 추정할 때 사용되는 유전자의 실질적인 길이입니다. 이 값은 유전자 길이에서 정렬 가능한 리드의 길이 및 시퀀싱 방법 등을 고려하여 조정된 값입니다.

04-3. Make matrix file (read_count, TPM)

샘플별 분석된 결과를 하나의 파일로 합치는 작업으로 이번 과정부터는 R 프로그램을 기반으로 진행됩니다.

• Input
  • RSEM 분석 후 생성된 genes.results file
  • SRR26528163.genes.results
  • SRR26528164.genes.results
  • SRR26528165.genes.results
  • SRR26528166.genes.results
  • SRR26528167.genes.results
  • SRR26528168.genes.results
• Practice #R library(dplyr) ### 폴더 경로 설정 및 파일 리스트 dir_path <- '/path/to/your/directory' file_list <- list.files(path = dir_path, pattern = '*.genes.results', full.names = TRUE) ### Load expected_count file (1, 5 column) extracted_data <- lapply(file_list, function(file) { data <- read.table(file) # 파일 불러오기 if(ncol(data) >= 5) { # 데이터에 최소 5개의 열이 있는지 확인 selected_data <- data %>% select(1, 5) # 1열과 5열 선택 colnames(selected_data) <- c("key", sub("\\..*$", "", basename(file))) # 열 이름 변경 return(selected_data) } else { warning(paste("File", file, "does not have at least 5 columns. Skipping...")) return(NULL) } }) # NULL 값을 제거 extracted_data <- Filter(Negate(is.null), extracted_data) # 초기 데이터 프레임 설정 combined_data <- extracted_data[[1]] # NULL 값을 제거 for (i in 2:length(extracted_data)) { combined_data <- full_join(combined_data, extracted_data[[i]], by = "key") } # 첫 번째 행 제거 및 쳇 번째 열을 열 이름으로 설정 readcount_data <- combined_data %>% slice((-1)) rownames(readcount_data) <- readcount_data$key readcount_data <- readcount_data[,-1] # 파일명을 샘플명으로 변경 new_colnames <- c("PCL_siSlu7_3", "PCL_siSlu7_2", "PCL_siSlu7_1", "PCL_siGL_3", "PCL_siGL_2", "PCL_siGL_1") colnames(readcount_data) <- new_colnames readcount_data <- readcount_data[,rev(colnames(readcount_data))] # 결과 확인 print(readcount_data) write.table(readcount_data, file="00_readcount.txt", sep="\t", quote=FALSE) rm(extracted_data) rm(combined_data) rm(new_colnames) ### Load TPM file (1, 6 column) extracted_data <- lapply(file_list, function(file) { data <- read.table(file) # 파일 불러오기 if(ncol(data) >= 6) { # 데이터에 최소 6개의 열이 있는지 확인 selected_data <- data %>% select(1, 6) # 1열과 6열 선택 colnames(selected_data) <- c("key", sub("\\..*$", "", basename(file))) # 열 이름 변경 return(selected_data) } else { warning(paste("File", file, "does not have at least 6 columns. Skipping...")) return(NULL) } }) # NULL 값을 제거 extracted_data <- Filter(Negate(is.null), extracted_data) # 초기 데이터 프레임 설정 combined_data <- extracted_data[[1]] # NULL 값을 제거 for (i in 2:length(extracted_data)) { combined_data <- full_join(combined_data, extracted_data[[i]], by = "key") } # 첫 번째 행 제거 및 쳇 번째 열을 열 이름으로 설정 tpm_data <- combined_data %>% slice((-1)) rownames(tpm_data) <- tpm_data$key tpm_data <- tpm_data[,-1] # 파일명을 샘플명으로 변경 new_colnames <- c("PCL_siSlu7_3", "PCL_siSlu7_2", "PCL_siSlu7_1", "PCL_siGL_3", "PCL_siGL_2", "PCL_siGL_1") colnames(tpm_data) <- new_colnames tpm_data <- tpm_data[,rev(colnames(tpm_data))] # 결과 확인 print(tpm_data) write.table(tpm_data, file="00_tpm.txt", sep="\t", quote=FALSE) rm(extracted_data) rm(combined_data) rm(new_colnames) # 본 분석에서는 TPM과 expected_count 만 사용하므로, FPKM을 추출하지 않음 # 필요 시 위의 코드에서 5 또는 6 column을 7로 변경하여 코드를 돌리면 해당 파일을 얻을 수 있음
• Output
  • 00_readcount.txt
  • 00_tpm.txt

05. Expression Profiling

데이터의 정확성과 비교 가능성을 위하여 데이터를 normalization하고 PCA을 적용합니다. Normalization은 각 샘플 간의 리드 수의 차이를 보정하여 유전자의 발현량을 정확하게 비교할 수 있게 합니다. 이는 실험 조건이나 시료 처리로 인해 발생할 수 있는 기술적 변동을 보정함으로써 데이터의 일관성을 유지하고 통계적 검증 및 해석을 정확하게 하기 위한 필수적인 과정입니다. PCA은 고차원의 데이터를 낮은 차원의 데이터로 변환하여 샘플 간의 유사성을 확인하는 기법입니다.

05-1. Normalization & PCA

• Input
  • 00_tpm.txt
• Practice #R library(ggplot2) library(preprocessCore) ### Load a file tpm_data <- read.table('00_tpm.txt', row.names=1, header = T, sep="\t")) ### Normalization data.nor <- normalize.quantiles(as.matrix(log2(tpm_data + 1)), copy = TRUE) dimnames(data.nor) = dimnames(tpm_data) # Figure of normalization colors = c(rep('purple',3), rep('orange',3)) # 설명) 그룹이 구분되도록 color 지정 (그룹의 색(gray), 그룹의 샘플 수(4)) # 예시) colors = c(rep('gray', 4), rep('red', 4), rep('orange',4)) 그룹 3개, 그룹 당 샘플 수가 4개인 경우 pdf("01_boxplot_TPM_Nor_before.pdf") boxplot(log2(tpm_data + 1), boxfill=colors, main="Before Normalization", outline=F, las=2, ylab="log2(TPM + 1)", cex.axis = 0.7, cex.ayis = 0.7) dev.off() pdf("02_boxplot_TPM_Nor_after.pdf") boxplot(data.nor, boxfill=colors, main="After Normalization", outline=F, las=2, ylab="log2(TPM + 1)", cex.axis = 0.7, cex.ayis = 0.7) dev.off() ### PCA library(ggfortify) PCAt = prcomp(t(data.nor)) # Figure of PCA condition = c(rep('PCL_siGL',3), rep('PCL_siSlu7',3)) # 설명) 그룹이 구분되도록 샘플(그룹명) 지정 (그룹명(NC), 그룹의 샘플 수(4)) # 예시) condition = c(rep('NC', 4), rep('T1', 4), rep('T3', 4)) 그룹 3개, 그룹 당 샘플 수가 4개인 경우 PCAb <- cbind(t(data.nor), condition) pdf('03_PCA_TPM.pdf') autoplot(PCAt, data = PCAb, colour = 'condition', size = 4) + ggtitle("PCA") + theme_classic() + theme(plot.title = element_text(hjust = 0.5, face="bold", size=15)) dev.off()
• R function description
  • pcromp : 주성분 분석을 수행
  • autoplot : ggplot2 object 중 하나로 PCA 결과를 시각화 함
    • ggtitle : 그래프의 제목을 설정
    • theme : 그래프의 제목, 범례, 축의 위치 또는 글자 모양 등을 설정
    • theme_classic() : 배경은 흰색으로 설정, 축의 눈금과 레이블은 기본 스타일로 표시하며 제목, 범례, 주석 등의 텍스트는 기본 폰트와 크기로 표시
• Output
  • 01_boxplot_TPM_Nor_before.pdf
  • 02_boxplot_TPM_Nor_after.pdf
  • 03_PCA_TPM.pdf
  Output 예시

  • 01_boxplot_TPM_Nor_before.pdf
    
  • 02_boxplot_TPM_Nor_after.pdf
    
  • 03_PCA_TPM.pdf
    

06. DEG analysis

DEGs (Differential expressed genes)는 특정 조건에서 통계적으로 유의한 수준으로 발현이 달라지는 유전자들을 가리키는 용어입니다. 예를 들어, 암세포와 정상 세포 간에 유전자 발현 차이를 비교하여 암세포에서 높은 수준으로 발현되는 유전자들을 찾는 경우가 있을 수 있습니다.

이러한 DEGs 분석은 유전자 발현이 변화하는 생리학적 현상 및 기전 또는 질병의 원인과 발생 메커니즘 등을 이해하는 데 도움을 줄 수 있으며, 새로운 약물을 개발하기 위한 초석이 될 수 있습니다.

발현량 차이를 규명하기 위하여 목적에 따라 CPM, TPM, FPKM/RPKM, DESeq2, edgeR, limma 등이 사용되며, 이러한 도구는 각 유전자에 대한 일반화 선형 모델을 사용하여 차이를 식별하거나 통계적으로 발현에 대한 차이를 규정하고 DEG를 정의합니다.

06-1. DEG extract

• Input
  • 00_readcount.txt
• Practice #R library("preprocessCore") library("edgeR") # Load a file readcount_data <- read.table("00_readcount.txt", row.names=1, header=T, sep="\t") ### make a group condition = c(rep('PCL_siGL',3), rep('PCL_siSlu7',3)) # 설명) 그룹이 구분되도록 샘플(그룹명) 지정 (그룹명(NC), 그룹의 샘플 수(4)) # 예시) condition = c(rep('NC', 4), rep('T1', 4), rep('T3', 4)) 그룹 3개, 그룹 당 샘플 수가 4개인 경우 ### make a DEGlist readcount_data <- data.frame(sapply(readcount_data, as.numeric)) y <- DGEList(counts = readcount_data, group = condition) ### Filtering low counts keep <- filterByExpr(y) y <- y[keep, , keep.lib.size=FALSE] ### TMM normalize y <- calcNormFactors(y) ### estimate dispersion y <- estimateCommonDisp(y) y <- estimateTagwiseDisp(y) ### making model.matrix design <- model.matrix(~0+group, data=y$samples) ### DEGs - quasi-likelihood F-tests fit <- glmQLFit(y, design) qlf <- glmQLFTest(fit, contrast = c(-1,1)) DEGs <- topTags(qlf, n=Inf,sort.by="logFC") write.table(DEGs, file="04_DEG_list.txt", sep="\t", quote=FALSE)
• Output
  • 04_DEG_list.txt
  Output 예시

  

  • (Gene ID) : 유전자의 고유 식별자로 Ensembl ID, RefSeq ID, 유전자 이름이 될 수 있습니다.
  • logFC (log Fold Change) : 실험 조건 간의 유전자 발현 변화 비율을 로그로 변환한 값입니다. 일반적으로 로그2를 기준으로 변환되며, 값이 1이면 2배 발현 증가, -1이라면 2배 발현 감소를 의미합니다.
  • logCPM (log Counts Per Million) : 유전자의 평균 발현 수준을 로그 변환하여 나타낸 값입니다. CPM은 각 유전자의 read 수를 백만 단위로 표준화 한 값으로 시퀀싱 depth를 고려한 샘플 간의 비교가 가능하도록 합니다.
  • F (F-statistic) : 차등 발현을 검정하기 위한 F-통계량 입니다. F-통계량은 실험 조건 간의 그룹 내 변동과 그룹 간 변동을 비교하여 차이가 유의한지를 평가하는 데 사용됩니다. 값이 클수록 그룹 간의 발현 값의 차이가 크다는 것을 의미합니다.
  • PValue : 유전자가 차등 발현되는지에 대한 통계적 유의성을 나타내는 p-값입니다. p-값이 낮을수록 유전자 발현 차이가 통계적으로 유의미할 가능성이 높습니다. 다만, p-값만으로 유의성을 평가하기에 한계가 있어, FDR 값을 고려합니다.
  • FDR (False Discovery Rate) : Benjamini-Hochberg 절차를 사용하여 다중 검정에 대한 보정을 적용한 p-값으로 발견된 차등 발현 유전자 중 거짓 양성의 비율을 추정합니다. FDR 값이 낮을수록 유전자가 진정으로 차등 발현되었을 가능성이 높으며, 일반적으로 FDR 값이 0.05 이하이면 차등 발현이 유의미한 것으로 간주됩니다.

06-2. Setting cut-off of DEGs

전체 DEGs 리스트 중 | log2FC | 값이 1 이상이며, p-value 나 FDR 이 0.05 이하를 만족하는 DEG를 유의하다고 가정하여 DEGs를 필터링합니다.

하지만, 실험 및 결과에 따라 해당 기준 값은 변경할 수 있으며, 아래와 같이 코드로 필터링하거나 엑셀을 사용하여 필터링할 수 있습니다.

• Input
  • 위의 06-1.의 코드가 이어짐
• Practice #R p_value=0.05 logFC=1 # Filtering DEGs_set <- DEGs$table$FDR <= p_value & abs(DEGs$table$logFC) >= logFC DEGs_set <- DEGs[DEGs_set,] write.table(DEGs_set, file="05_DEG_list_fil.txt", sep="\t", quote=FALSE)
• Output
  • 05_DEG_list_fil.txt

07. Visualization

07-1. Volcano plot

Volcano plot은 mRNA 발현 데이터 분석에서 데이터 시각화를 위해 사용되며, x-축에는 발현의 차이를 나타내는 fold change를, y-축에는 통계적 유의성을 나타내는 -log10(p-value) 값을 표시합니다. Volcano plot에서 발현의 차이가 크고 통계적으로 유의한(p-value가 작은) mRNA는 plot의 양쪽 끝 상단에 위치하며, 이러한 mRNA는 다양한 생물학적 역할을 하는 중요한 유전자들일 가능성이 높다고 판단할 수 있습니다.

• Input
  • 04_DEG_list.txt (위의 06-1.의 코드에서 이어질 경우, 변수 DEGs$table로 입력)
• Practice #R library(EnhancedVolcano) p_value=0.05 logFC=1 DEGs <- read.table("04_DEG_list.txt", row.names=1, header=T, sep="\t") volc <- EnhancedVolcano(DEGs, #DEGs$table lab = rownames(DEGs), x = 'logFC', y = 'PValue', pCutoff = p_value, FCcutoff = logFC, # 그래프를 보고 x축 y축 범위 조정 # x축은 DEGs의 logFC 범위 # y축은 DEGs의 PValue의 지수값으로 계산 xlim = c(-7.5, 5), ylim = c(0, 6), legendPosition = 'bottom', legendLabSize = 10) pdf('06_Volcano_plot.pdf') volc dev.off()
• Output
  • 06_Volcano_plot.pdf
  Output 예시

  

07-2. Heatmap plot

Heatmap은 mRNA 발현 데이터 분석에서 데이터 시각화를 위해 사용되며, 색상으로 유전자 발현 수준의 변화를 나타냅니다. Heatmap에서 유사한 발현 패턴을 가진 유전자들은 클러스터링을 통해 그룹화되며, 이는 샘플 간 또는 유전자 간의 관계를 이해하는 데 유용합니다. 발현 패턴이 유사한 유전자들은 종종 동일한 생물학적 경로에 관여할 가능성이 높으며, 다양한 생물학적 역할을 하는 중요한 유전자들을 발견하는 데 도움을 줍니다.

• Input
  • 00_tpm.txt
  • 05_DEG_list_fil.txt (위의 06-2.의 코드에서 이어질 경우, 변수 DEGs_set 사용))
• Practice #R library(preprocessCore) library(pheatmap) tpm_data = read.table('00_tpm.txt', row.names=1, header = T, sep="\t")) data.nor <- normalize.quantiles(as.matrix(log2(tpm_data + 1)), copy = TRUE) dimnames(data.nor) = dimnames(tpm_data) DEGs_set <- read.table("05_DEG_list_fil.txt", row.names=1, header=T, sep="\t") ## pdf('07_Heatmap.pdf') pheatmap(data.nor[rownames(DEGs_set), ], show_rownames = F, scale = "row") dev.off()
• Output
  • 07_Heatmap.pdf
  Output 예시

  

08. Functional Profiling

08-1. Gene ontology (BP, MF, CC)

Gene Ontology (GO)는 유전자의 기능과 생물학적 과정을 이해하기 위한 분석입니다. 이 분석은 대규모 유전자 데이터 세트에서 유전자들의 기능과 특성을 파악하거나, 유전자와 단백질의 기능과 상호 작용을 이해하는 데 도움을 줍니다.

• Input
  • 05_DEG_list_fil.txt (위의 06-2.의 코드에서 이어질 경우, 변수 DEGs_set 사용))
• Practice #R library(org.Hs.eg.db) library(AnnotationDbi) library(clusterProfiler) DEGs_set <- read.table("05_DEG_list_fil.txt", row.names=1, header=T, sep="\t") DEG_list <- rownames(DEGs_set) go_result <- enrichGO(gene = DEG_list, OrgDb = "org.Hs.eg.db", # 유전자 DB (기본: human) keyType="SYMBOL", # 유전자 ID 형식 ont="all") # 분석할 go의 유형 ("BP", "MF", "CC") pdf('08_GO_dotplot.pdf') dotplot(go_result, split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale="free") dev.off() write.table(go_result, file="09_GO_result.txt", sep="\t", quote=FALSE)
• Output
  • 08_GO_dotplot.pdf
  • 09_GO_result.txt
    Output 예시

    • 08_GO_dotplot.pdf
      
    • 09_GO_result.txt
      
      • Ontology : Ontology 종류
      • ID : GO 항목의 ID (예: "GO:0008150")
      • Description : 해당 GO 항목의 설명 (예: "biological_process")
      • GeneRatio : 해당 GO 항목에 속한 유전자 수와 입력 유전자 리스트의 비율
      • BgRatio : 해당 GO 항목에 속한 전체 유전자 수와 데이터베이스 내 전체 유전자 수의 비율
      • pvalue : 해당 GO 항목의 p-value
      • p.adjust : 다중 검정을 수행한 후의 조정된 p-value
      • qvalue : FDR (False Discovery Rate) 기준의 q-value
      • geneID : 해당 GO 항목에 속하는 유전자 ID
      • Count : 해당 GO 항목에 속하는 유전자 수

08-2. Pathway (KEGG)

KEGG 경로 데이터베이스는 분자 상호 작용, 생화학적 반응, 그리고 생물학적 관계를 기반으로 한 네트워크를 경로 지도 형태로 제공합니다. 이러한 지도는 유전자, 단백질, RNA, 화학 화합물, 글리칸 및 화학 반응뿐만 아니라 질병 유전자와 약물 표적까지 통합할 수 있습니다.

• Input
  • 05_DEG_list_fil.txt (위의 06-2.의 코드에서 이어질 경우, 변수 DEGs_set 사용))
• Practice #R library(org.Hs.eg.db) library(AnnotationDbi) library(clusterProfiler) DEGs_set <- read.table("05_DEG_list_fil.txt", row.names=1, header=T, sep="\t") DEG_list <- rownames(DEGs_set) # entrez id로 바꾸기 DEGs_entrez <- mapIds(org.Hs.eg.db, DEG_list, 'ENTREZID', 'SYMBOL') length(DEGs_entrez) length(na.omit(DEGs_entrez)) ### KEGG kegg_entrez <- enrichKEGG(DEGs_entrez, organism = 'hsa') kegg_entrez pdf('10_KEGG_dotplot.pdf') dotplot(kegg_entrez) dev.off() kegg_result <- setReadable(kegg_entrez, 'org.Hs.eg.db', 'ENTREZID') write.table(kegg_result, file="11_KEGG_result.txt", sep="\t", quote=FALSE)
• Output
  • 10_KEGG_dotplot.pdf
  • 11_KEGG_result.txt
    Output 예시

    • 10_KEGG_dotplot.pdf
      
    • 11_KEGG_result.txt
      
      • Category : KEGG pathway가 속한 상위 카테고리
      • Subcategory : KEGG pathway가 속한 하위 카테고리
      • ID : KEGG 항목의 ID (예: "hsa04668")
      • Description : 해당 KEGG 항목의 설명 (예: "TNF signaling pathway")
      • GeneRatio : 해당 KEGG 항목에 속한 유전자 수와 입력된 유전자 리스트의 비율
      • BgRatio : 해당 KEGG 항목에 속한 유전자 수와 DB 내 전체 유전자 수의 비율
      • pvalue : 해당 GO 항목의 p-value
      • p.adjust : 다중 검정을 수행한 후의 조정된 p-value
      • geneID : 해당 KEGG pathway에 속한 유전자의 리스트.
      • Count : 해당 KEGG pathway와 관련된 유전자 수.

강물이나 바닷물, 흙 같은 자연 환경부터 생물의 장내미생물, 분변에 이르기까지, 세균, 곰팡이, 원생생물, 바이러스 등 다양한 미생물이 존재합니다. 특정 환경이나 숙주 내에 서식하는 모든 미생물의 유전체를 연구하여 어떤 미생물 종이 존재하는지 분석하는 것을 Metagenomics라고 합니다.

Metagenomics는 연구 방법에 따라 targeted metagenomics와 shotgun metagenomics로 구분됩니다.

Targeted metagenomics는 미생물 유전체의 특정 영역을 증폭하여 분석하는 방법입니다. 이 과정에서는 모든 미생물에서 보존되어 있으면서도 종 간의 서열 차이를 식별할 수 있는 유전자 부위를 타겟으로 합니다. 대표적으로 원핵생물에서는 16S 리보솜 DNA(rDNA), 진핵생물에서는 18S rDNA와 internal transcribed spacers(ITS)가 주로 사용됩니다. 이 방법은 적은 시퀀싱 데이터로도 충분히 연구를 수행할 수 있고, 분석이 비교적 간단하다는 장점이 있습니다. 그러나 PCR 과정에서 특정 종의 서열이 과대 또는 과소 증폭되는 편향(bias)이 발생할 수 있으며, 보통 속(Genus) 수준까지의 동정만 가능하다는 한계가 있습니다.

Shotgun metagenomics는 샘플 내 모든 미생물의 전체 유전체 또는 유전체의 일부를 무작위로 시퀀싱하는 방법입니다. 이 방식은 전체 유전체를 분석하기 때문에 종(Species) 수준까지의 정확한 동정이 가능하며, 유전자 구성, 대사 경로 등 군집의 기능적 특성까지 분석할 수 있다는 장점이 있습니다. 또한, PCR 과정을 생략하므로 targeted metagenomics에서 발생할 수 있는 편향이 없습니다. 그러나 상대적으로 비용이 많이 들고, 데이터 분석에 많은 자원이 필요합니다. 특히, 숙주와 관련된 샘플에서는 host DNA가 과도하게 포함되어 분석이 어려운 경우가 있을 수 있습니다.

Targeted metagenomics는 비용 효율성과 분석의 간단함으로 인해 여전히 널리 사용되고 있으며, 미생물 군집 구조를 이해하는 데 유용한 도구로 자리 잡고 있습니다. 이 튜토리얼에서는 targeted metagenomics를 통해 생성된 데이터를 분석하는 방법을 중점적으로 다룹니다. 분석에 널리 사용되는 대표적인 도구로는 Qiime2, Mothur, DADA2, USEARCH/VSEARCH 등이 있으며, 본 튜토리얼에서는 Qiime2를 활용한 분석 과정을 단계별로 자세히 설명합니다.

Qiime2는 python conda 환경에서 실행이 가능합니다. 다음 코드로 Qiime2 분석을 위한 가상 환경을 구축할 수 있습니다. 가상 환경 이름은 qiime2-amplicon-2024.5로 지정했습니다.

# Linux - Install QIIME 2 within a conda environment # we’ll name the environments qiime2-<distro>-2024.5 to indicate what QIIME 2 release is installed conda env create -n qiime2-amplicon-2024.5 --file https://data.qiime2.org/distro/amplicon/qiime2-amplicon-2024.5-py39-linux-conda.yml # Activate the conda environment conda activate qiime2-amplicon-2024.5 # Test your installation qiime --help

Raw data를 해당 프로그램에서 인식할 수 있는 artifact 형태로 변환하는 과정이 필요합니다. 디렉터리를 입력합니다.

Qiime2는 다양한 import 옵션이 존재합니다. 다음 Qiime2 tutorial 페이지를 참조하시길 바랍니다. (https://docs.qiime2.org/2024.10/tutorials/importing/)

• Input
  • Raw data path (emp-single-end-sequences)
• Practice qiime tools import \ --type EMPSingleEndSequences \ --input-path emp-single-end-sequences \ --output-path emp-single-end-sequences.qza
• Option description --type : import 하는 데이터의 유형, 예를 들어 이번에 사용한 예제 데이터는 EMP single-end 데이터임으로 EMPSingleEndSequences를 입력 --input-path : import 하는 파일이나 디렉터리의 위치 --output-path : qiime2 artifact 형태로 변환된 결과 파일 이름
• Output
  • emp-single-end-sequences.qza

예제 데이터인 EMP 데이터와 같이 모든 샘플을 합쳐 시퀀싱하는 방식으로 생산된 데이터는 바코드에 따라 샘플 별로 데이터를 나누는 demultiplexing 과정이 필요합니다.

• Input
  • emp-single-end-sequences.qza
• Practice qiime demux emp-single \ --i-seqs emp-single-end-sequences.qza \ --m-barcodes-file sample_metadata.tsv \ --m-barcodes-column barcode-sequence \ --o-per-sample-sequences demux.qza \ --o-error-correction-details demux-details.qza
• Option description --i-seqs: 입력 시퀀스 파일 (.qza). --m-barcodes-file: 샘플 메타데이터 파일 (TSV 형식). --m-barcodes-column: 메타데이터에서 사용할 바코드 열 이름. --o-per-sample-sequences: 각 샘플의 시퀀스를 포함하는 출력 파일. --o-error-correction-details: 오류 수정 정보 기록 파일
• Output
  • demux.qza
  • demux-details.qza

qza 파일은 Qiime2 tool 내에서만 사용되는 전용 artifact 파일입니다. 때문에 일반적인 방법으로 내용을 확인할 수 없어 qzv 파일 형태로 변환한 뒤 Qiime2 Viewer (https://view.qiime2.org/)에서 열람해야 합니다.

demux.qza를 demux.qzv로 변환하면 시퀀스 데이터에 대해 전체적인 통계 수치를 얻을 수 있습니다.

• Input
  • demux.qza
• Practice qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv
• Option description --i-data : input 파일 이름 --o-visualization : output 파일 이름
• Output
  • demux.qzv

Qiime2는 데이터 품질 관리를 위해 실험 과정에서 발생하는 노이즈 데이터를 제거(denoising)하는 다양한 도구를 제공합니다. 대표적으로 DADA2, Deblur 등이 있습니다.

이후 추출된 서열은 각각 고유한 서열로 간주되어 임의의 feature identifier가 부여됩니다. 이 정보는 rep-seq.qza에 저장됩니다. 그리고 각 고유 서열의 수를 counting한 데이터는 table 형태로 table.qza에 저장됩니다. 이 두 데이터는 이후 Metagenome 분석에 사용되는 기본 데이터인 feature table을 이룹니다.

이 튜토리얼에서는 여러 denoising 방법 중 DADA2 파이프라인을 사용하였습니다.

DADA2는 Illumina 시퀀스 데이터의 품질 관리를 위한 파이프라인으로, raw data에서 phiX 서열 및 키메라 서열과 같은 노이즈 데이터를 제거하는 과정을 포함합니다. 이를 통해 최종적으로 고품질의 시퀀스 데이터만을 확보하여 이후의 분석에 활용할 수 있습니다.

• Input
  • demux.qza
• Practice # DADA2 denoising qiime dada2 denoise-single \ --i-demultiplexed-seqs demux.qza \ --p-trim-left 0 \ --p-trunc-len 120 \ --o-representative-sequences rep-seqs.qza \ --o-table table.qza \ --o-denoising-stats stats.qza
• Option description --i-demultiplexed-seqs: 입력으로 사용할 demultiplex 된 시퀀스 파일 (.qza). --p-trim-left: 서열의 시작에서 자를 염기 수를 지정. 위의 경우 0이므로 자르지 않음. --p-trunc-len: 각 서열을 잘라낼 위치를 지정. 위의 경우 120번째 위치에서 서열을 자름. --o-representative-sequences: denoising 후 생성된 representative sequence를 포함하는 출력 파일 이름 (.qza) --o-table: denoising 결과로 생성된 feature 테이블을 포함하는 출력 파일 이름 (.qza) --o-denoising-stats: denoising 과정에서 생성된 통계 정보를 포함하는 출력 파일 이름 (.qza)
• Output
  • rep-seqs.qza
  • table.qza
  • stats.qza

결과 파일(table.qza, rep-seqs.gza)을 qzv 형식으로 변환하면 Qiime2 view에서 내용을 확인할 수 있습니다.

또 stats.qza는 qzv로 변환하면 전체적인 denoising 결과에 대한 통계 수치를 얻을 수 있습니다.

qzv 파일은 Qiime2 viewer(https://view.qiime2.org/)에서 visualization이 가능합니다.

• Input
  • table.qza
  • rep-seqs.qza
  • stats.qza
  • sample_metadata.tsv
• Practice # table.qza visulization qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file sample_metadata.tsv # rep-seqs.qza visulization qiime feature-table tabulate-seqs \ --i-data rep-seqs.qza \ --o-visualization rep-seqs.qzv # stats.qza visulization qiime metadata tabulate \ --m-input-file stats.qza \ --o-visualization stats.qzv
• Option description --i-table : 입력할 feature table 파일 (.qza) --i-data : 입력할 representative sequence 파일 (.qza) --m-input-file : 입력할 statistic 파일 (.qza) --o-visualization : 시각화된 결과 파일 (.qzv)
• Output
  • table.qzv
  • rep-seqs.qzv
  • stats.qzv

계통발생학적 다양성 분석에는 어떤 공통 조상으로부터 분기한 관계를 보여주는 계통수(rooted tree) 정보가 필수입니다. 이 분석을 위하여 align-to-tree-mafft-fasttree을 사용하여 representative sequence를 처리합니다.

먼저, MAFFT(Multiple Alignment using Fast Fourier Transform)은 Fast Fourier 알고리즘을 기반으로 빠르고 효율적으로 다중 정렬을 수행하는 도구입니다. 이 과정에서 계통학적 정보가 부족하거나 애매하게 정렬된 열은 제거됩니다. 정렬된 데이터는 masked 상태로 표현됩니다. 그 다음, 정렬된 서열을 바탕으로 FastTree를 사용하여 계통수(phylogenetic tree)를 추론합니다. FastTree는 최대 우도(maximum likelihood) 방법을 이용하여 대규모 데이터 세트를 신속하고 효율적으로 처리합니다. 최종적으로 이 과정에서 계통수의 중간 지점에서 공통 조상, 즉 분기점에 대한 데이터도 산출됩니다.

• Input
  • rep-seqs.qza
• Practice qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza
• Option description --i-sequences : 입력할 representative sequence 파일 (.qza) -o-alignment: 정렬된 서열 파일 (.qza) -o-masked-alignment: 정보가 적거나 불분명한 열이 제거된 정렬 서열 파일 (.qza) -o-tree: 생성된 분기 정보가 없는(unrooted) 계통수를 출력 (.qza) -o-rooted-tree: 분기 정보가 있는 (rooted) 계통수를 출력 (.qza)
• Output
  • aligned-rep-seqs.qza
  • masked-aligned-rep-seqs.qza
  • unrooted-tree.qza
  • rooted-tree.qza

샘플 내 혹은 샘플 간에 미생물의 다양성을 수치로 나타낼 수 있으며, 이를 alpha diversity와 beta diversity로 구분합니다.

Alpha diversity에서 단일 샘플 내의 다양성을 수치로 나타낼 수 있는데, 대표적으로 shannon 지수, observed_features 등이 있습니다. 반면 beta diversity는 서로 다른 샘플 간의 다양성을 나타냅니다. 그래서 샘플 조합 별로 수치가 다르게 나옵니다. 대표적인 지표로 bray-curtis, jaccard, unweighted unifrac, weighted unifrac 등이 있습니다.

아래 코드를 통해 이러한 결과를 한 번에 얻을 수 있습니다. 이후 Beta diversity 결과는 emperor이라는 tool에서 주성분 분석(PCoA, Principle Coordinates Analysis)을 통해 시각화할 수 있습니다.

• Input
  • rooted-tree.qza
  • table.qza
• Practice qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 1103 \ --m-metadata-file sample_metadata.tsv \ --output-dir core-metrics-results
• Option description --i-phylogeny : rooted 계통수 파일 (.qza) --i-table : 샘플들에 대한 feature table 파일 (.qza) --p-sampling-depth : 각 샘플에서 무작위로 선택할 read의 수(샘플링의 깊이 설정) --m-metadata-file : Emperor 플롯에 사용할 metadata 파일 (.tsv) --output-dir : 실행 결과를 저장할 디렉터리
• Output
  • alpha diversity result
    • core-metrics-results/shannon_vector.qza
    • core-metrics-results/observed_features_vector.qza
    • core-metrics-results/evenness_vector.qza
    • core-metrics-results/faith_pd_vector.qza
    • core-metrics-results/rarefied_table.qza
  • beta diversity result
    • core-metrics-results/bray_curtis_distance_matrix.qza
    • core-metrics-results/bray_curtis_pcoa_results.qzv
    • core-metrics-results/jaccard_distance_matrix.qza
    • core-metrics-results/jaccard_pcoa_results.qzv
    • core-metrics-results/unweighted_unifrac_distance_matrix.qza
    • core-metrics-results/unweighted_unifrac_pcoa_results.qzv
    • core-metrics-results/weighted_unifrac_distance_matrix.qza
    • core-metrics-results/weighted_unifrac_pcoa_results.qzv
    • core-metrics-results/unweighted_unifrac_emperor.qzv
    • core-metrics-results/jaccard_emperor.qzv
    • core-metrics-results/bray_curtis_emperor.qzv
    • core-metrics-results/weighted_unifrac_emperor.qzv

Alpha diversity는 샘플 내에서 종 풍부도와 다양성을 측정하는 지표로, 샘플 그룹 간의 유의미한 차이를 비교하기 위해서는 샘플 그룹 간의 통계적 검정이 필요합니다. 주로 Kruskal-Wallis test와 같은 비모수 검정을 사용하며, 이들 간의 차이를 평가합니다. 이러한 분석 결과는 boxplot로 시각화되며 p-value 같은 통계적 수치도 함께 표시됩니다. 이는 서로 다른 조건에 따른 생물 다양성의 차이를 평가할 때 유용한 방법입니다.

• Input
  • faith_pd_vector.qza
  • evenness_vector.qza
• Practice qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \ --m-metadata-file sample_metadata.tsv \ --o-visualization core-metrics-results/faith-pd-group-significance.qzv qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/evenness_vector.qza \ --m-metadata-file sample_metadata.tsv \ --o-visualization core-metrics-results/evenness-group-significance.qzv
• Option description --i-alpha-diversity : 입력할 alpha diversity 파일 (.gza) --m-metadata-file : metadata 파일 (.tsv) --o-visualization : 결과 시각화 파일(.qzv)
• Output
  • core-metrics-results/faith-pd-group-significance.qzv
  • core-metrics-results/evenness-group-significance.qzv

샘플마다 어떠한 미생물 분류군이 존재하는지, 그 분류학적 분포를 파악하는 단계입니다. 각 feature의 representative sequence와 참조 데이터베이스의 서열을 비교하여 해당 feature가 어떠한 분류학적 그룹(taxonomy)에 속하는지를 식별합니다. 이번 분석에서는 Qiime2에서 제공한 Greengenes라는 DB 기반의 gg-13-8-99-515-806-nb-classifier.qza 파일을 활용했습니다.

• Input
  • rep-seqs.qza
  • gg-13-8-99-515-806-nb-classifier.qza
• Practice qiime feature-classifier classify-sklearn \ --i-classifier gg-13-8-99-515-806-nb-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza
• Option description --i-classifier : 훈련된 classifier 파일 (.qza) --i-reads : 입력할 representative sequence 파일 (.qza) --o-classification : taxonomy 매칭 결과 파일(.qza)
• Output
  • taxonomy.qza

Qiime2 viewer에서 taxonomy.qza를 볼 수 있도록 qzv 파일로 변환합니다.

• Input
  • taxonomy.qza
• Practice # Making qzv file qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv
• Option description --m-input-file : 입력할 taxonomy 파일 (.qza) --o-visualization : 결과 시각화 파일(.qzv)
• Output
  • taxonomy.qzv

Taxonomy 결과와 기존 feature table 결과를 통합하여 샘플별로 존재하는 taxonomy를 barplot 형태로 볼 수 있도록 시각화합니다. metadata 정보를 입력해 x축, y축을 그룹, 임상 정보 별로 정렬할 수 있습니다.

• Input
  • table.qza
  • taxonomy.qza
  • sample_metadata.tsv
• Practice qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample_metadata.tsv \ --o-visualization taxa-bar-plots.qzv
• Option description --i-table : 입력할 feature table 파일 (.qza) --i-taxonomy : 입력할 taxonomy 파일 (.qza) --m-metadata-file : metadata 파일 (.tsv) --o-visualization : 결과 시각화 파일(.qzv)
• Output
  • taxa-bar-plots.qzv

서로 다른 실험 조건이나 집단 간에 특정 미생물 종이 유의미하게 차이가 나는지를 확인하기 위해 Differential Abundance Testing을 수행합니다. 이는 RNA-seq의 DEG(Differentially Expressed Genes) 분석과 유사한 목적을 가지고 있으며, metagenome 분석에서도 두 그룹 간의 미생물 분류군에서 유의미한 차이를 식별할 수 있도록 도와줍니다. 예를 들어, 건강한 사람과 질병을 앓고 있는 사람의 microbiome 차이를 분석하거나, 환경 변화에 따른 미생물 군집의 변화를 평가할 수 있습니다.

본 튜토리얼에서는 피실험자(subject)의 장내 미생물(gut microbiome) 군집 데이터를 비교 분석하여, 두 피실험자 간에 어떤 미생물이 차별적으로 풍부한지 확인하고자 합니다. 따라서 기존 분석된 데이터에서 gut의 데이터 만을 따로 추출하였습니다.

• Input
  • table.qza
• Practice qiime feature-table filter-samples \ --i-table table.qza \ --m-metadata-file sample_metadata.tsv \ --p-where "[body-site]='gut'" \ --o-filtered-table gut-table.qza
• Option description --i-table : 입력할 feature table 파일 (.qza) --m-metadata-file : 메타데이터 (.tsv) --p-where : table에서 추출한 값의 조건 --o-filtered-table : 필터링된 결과 table 파일(.qzv)
• Output
  • gut-table.qza

ANCOM-BC는 군집 내 미생물의 차등 풍부도 차이를 평가하는 통계적 방법으로, 미생물 군집 데이터의 편향을 교정하여 정확한 결과를 제공합니다. 특히 회귀 모델을 사용하여 각 미생물의 절대 풍부도를 추정하고, 이를 바탕으로 두 집단 간 차이를 평가합니다.

ANCOM-BC를 이용해 두 피실험자 간의 장내 미생물 군집에서 통계적으로 유의미한 차이를 확인하고, 결과를 barplot으로 시각화하였습니다.

• Input
  • gut-table.qza
  • sample_metadata.tsv
• Practice qiime composition ancombc \ --i-table gut-table.qza \ --m-metadata-file sample_metadata.tsv \ --p-formula 'subject' \ --o-differentials ancombc-subject.qza qiime composition da-barplot \ --i-data ancombc-subject.qza \ --p-significance-threshold 0.001 \ --o-visualization da-barplot-subject.qzv
• Option description [ancombc] --i-table : 입력할 feature table 파일 (.qza) --m-metadata-file : metadata 파일 (.tsv) --p-formula : 분석에 사용할 임상정보 열의 이름 --o-differentials : 차등 풍부도 결과 파일 (.qza) [da-barplot] --i-data ancombc-subject.qza : 차등 풍부도(ANCOM-BC) 결과 파일 (.qza) --p-significance-threshold : 유의성 임계값. 위의 경우 p-value 0.001 이하일 때 유의미한 차이로 간주함. --o-visualization : 결과 시각화 파일(.qzv)
• Output
  • ancombc-subject.qza
  • da-barplot-subject.qzv