리눅스 환경에서 NGS 데이터를 이용한 엽록체 조립 및 시각화 가이드 (NOVOPlasty, R)

보통 Illumina로 NGS분석을 진행하면 fastq.gz로 된 파일을 2개 받을 수 있다. 해당 파일을 각각 1.fastq.gz, 2.fastq.gz라고 가정하고 분석을 진행할 수 있다. 모든 것은 리눅스 환경에서 진행한다.
먼저 NOVOPlasty에서 config.txt에서 다음과 같은 값들을 셋팅한다. (괄호 안은 지울 것)

Project:
-----------------------
Project name          = Plant_A (프로젝트 이름)
Type                  = chloro
Genome Range          = 150000-160000 (예상 크기)
K-mer                 = 27 (예시. 조정 필요)
Max memory            = 16
Extended log          = 0
Save assembled reads  = no
Seed Input            = /path/Vn_rbcL.fasta (matK, rbcL 등 seed 서열)
Extend seed directly  = no
Reference sequence    = 
Variance detection    = 
Chloroplast sequence  = /path/chloroplast.fasta (근연종_엽록체_전체_염기서열)

Dataset 1:
-----------------------
Read Length           = NGS Read Length (회사/조건마다 다름)
Insert size           = NGS Insert Size (회사/조건마다 다름)
Platform              = illumina
Single/Paired         = PE
Combined reads        = 
Forward reads         = /path/1.fastq.gz (1.fastq.gz 위치)
Reverse reads         = /path/2.fastq.gz (2.fastq.gz 위치)
Store Hash            =

Heteroplasmy:
-----------------------
MAF                   = 
HP exclude list       = 
PCR-free              = 

Optional:
-----------------------
Insert size auto      = yes
Use Quality Scores    = no
Reduce ambigious N's  = 
Output path           = /path/output/ (output이 생성될 폴더 위치) 
(나머지는 그냥 그대로 두어도 됨)

 

이후 저장하고 리눅스에서 NOVOPlasty로 엽록체 조립을 실행한다.
명령어: perl NOVOPlasty4.3.5.pl -c config.txt

엽록체 조립이 끝나면 Output path 경로에 fasta 파일들이 생긴다. 여러 염기서열이 생기는 경우, k-mer값 조정이나 다른 근연종과의 비교를 통해 어떤 염기서열이 정확한 염기서열이 맞는 건지 확인이 필요하다. 이후 확인이 끝났으면 GeSeq에서 엽록체 지도를 그린다. https://chlorobox.mpimp-golm.mpg.de/geseq.html

조립된 엽록체 지도 예시

이후 Sequencing depth를 확인하는 작업이 필요하다.

Trimmomatic을 통해서 fastq.gz파일에 있는 Illumina 시퀀서들의 adapter들을 제거하고 trimming 과정을 진행한다. Java가 설치되어 있어야 한다.
명령어: java -jar trimmomatic-0.39.jar PE -phred33 /path/1.fastq.gz /path/2.fastq.gz 1_paired.fastq.gz 1_unpaired.fastq.gz 2_paired.fastq.gz 2_unpaired.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

분석을 진행하면 4개의 파일이 생성되는데, 각각 1에서 paired된 것, 1에서 unpaired된 것, 2에서 paired된 것, 2에서 unpaired된 것 등이다. PE 분석에서 쌍이 없는 unpaired된 것은 필요가 없어서 더 이상 분석에 사용되지 않는다. 그다음에 bwa와 samtools를 이용하여 sam파일과 bam파일을 생성한다. 
1단계: bwa로 앞서 조립한 엽록체 염기서열에 인덱스한다.
명령어: bwa index assembled_chloroplast.fasta

2단계: FASTQ 파일을 SAM 파일로 정렬한다.
명령어: bwa mem assembled_chloroplast.fasta 1_paired.fastq 2_paired.fastq > output.sam

3단계: SAM 파일을 BAM 파일로 변환 및 정렬
 SAM 파일 변환: samtools view -S -v output.sam > output.bam
 BAM 파일 정렬: samtools sort output.bam -o output.sorted.bam
 색인 생성: samtools index output.sorted.bam

최종적으로 생성된 output.sorted.bam 파일을 R을 통해서 Sequencing depth 그래프를 그릴 수 있다. 다음은 R 코드이다.

library(Rsamtools)
library(GenomicAlignments)
library(ggplot2)
library(dplyr)
library(scales)

# 1. 데이터 로드 및 전처리
# 파일 경로를 입력하세요
bam_path <- "/path/output.sorted.bam"
cov_data <- coverage(bam_path)

# Coverage 데이터를 데이터프레임으로 변환 (첫 번째 컨티그 기준)
df_cov <- data.frame(
  pos   = seq_along(cov_data[[1]]), 
  count = as.vector(cov_data[[1]]),
  species = "Plant_Species_Name"
)

# 2. 시각화 (Genome Coverage Map)
p <- ggplot(df_cov, aes(x = pos, y = count, color = species)) +
  geom_line(alpha = 0.8) +
  labs(
    title = "Sequencing Depth and Coverage Map",
    x = "Genome Position (bp)",
    y = "Coverage Depth (x)",
    color = "Species"
  ) +
  # X축 설정: 25,000bp 간격
  scale_x_continuous(labels = comma, breaks = seq(0, 150000, by = 25000)) +
  # Y축 설정: 천 단위 구분 쉼표 추가
  scale_y_continuous(labels = comma) + 
  scale_color_manual(values = c("Plant_Species_Name" = "blue")) +
  theme_classic() +
  theme(legend.position = "top")

print(p)

sequencing_depth.R

0.00MB

위와 같은 R 코드를 R 스튜디오에서 실행하면 아래 예시와 같은 그림을 얻을 수 있을 것이다.

엽록체 분석 Sequencing Depth and Coverage Map 예시

엽록체 분석 Sequencing Depth 그래프는 논문에서 자주 사용되지 않지만, Supplementary data로 요구되는 경우가 자주 있어서 분석을 해야 하는 것 같다.