Abstract
Early selection of grain quality traits in rice (Oryza sativa L.), is essential to improve the yield and quality of this staple crop. We analyzed four key traits—protein content, grain filling rate, height, and panicle length—in 85 Korean cultivars. Through whole- genome resequencing we identified 12,718,879 raw single nucleotide polymorphisms (SNPs); after PLINK-based quality control (bi-allelic selection, call rate≥0.90, MAF≥0.03), ~2.20 million high-quality SNPs remained for machine-learning (ML) pre-screening. To rank the features (without marker-level inference), we applied a liberal univariate PLINK case-control scan using the top and bottom 30% per trait. We also analyzed associations with a linear mixed model (GCTA v1.93.2, MLMA; fixed covariates: ecotype, PC1, PC2; random effect: GRM) to verify calibration under population structure; with n=85, no genome-wide significant hits were detected, and QQ-plots indicated adequate calibration (per-trait effective tests m≈1.54-1.57 million under stricter filters). The random forest feature importance prioritized 26, 51, 19, and 20 core SNPs for the four traits, respectively. Across the algorithms, the best models achieved mean accuracies of 81.8% (protein content), 81.0% (grain filling rate), 73.1% (height), and 94.0% (panicle length). All selected SNPs met the Fluidigm array design requirements, supporting its deployment as a compact, genotype-based panel for early selection and a practical step toward digital breeding in rice.
-
Keywords: rice; SNP marker; machine learning; GWAS; protein content; grain filling rate; height; panicle length
서언
벼(
Oryza sativa L.)는 세계적으로 가장 중요한 식량작물 가운데 하나이며, 수량과 미곡 품질을 좌우하는 조기 선발 역량은 현대 육종의 핵심 과제다(
Ding et al. 2024,
Quille-Mamani et al. 2025). 특히 단백질 함량, 등숙률, 초장(height), 수장(panicle length)은 수량성과 품질⋅재배 안정성과 직결되는 대표 형질로, 정확하고 재현성 있는 선발 지표를 마련할 필요가 있다(
Ding et al. 2024,
Dong et al. 2024,
Liu et al. 2022). 그러나 이들 형질은 환경 변이의 영향이 크고, 패널의 유전적 배경 차이에 의해 신호가 왜곡될 수 있어 기존 접근만으로는 정확한 선발에 한계가 있다(
Liu et al. 2022). 특히 국내 육성 품종처럼 생태형(ecotype; japonica, tongil-type) 이 혼재된 패널에서는 집단분화 신호가 형질 연관 신호로 과대추정될 수 있으므로, 주성분(PC)과 친연도(Genomic Relationship Matrix; GRM)를 반영한 혼합모형 GWAS로 보정하는 것이 필수적이다(
Kang et al. 2010,
Price et al. 2006,
VanRaden 2008,
Yang et al. 2011,
Zhou & Stephens 2012). 한편 전장서열 변이⋅GWAS 신호를 특징 선택(feature selection)으로 추려 기계학습(ML) 분류⋅예측 모델과 결합하는 디지털 육종 프레임워크가 빠르게 도입되고 있으며, 작물 및 원예 작물 전반에서 실용 사례가 축적되고 있다(
Quille-Mamani et al. 2025,
Yu et al. 2021). 이에 본 연구는 국내 육성 벼 85개 품종의 전장 유전체 데이터를 바탕으로, GCTA v1.93.2의 MLMA를 이용해 ecotype (고정효과) 및 PC1-PC2/GRM (임의효과)을 포함한 보정형 GWAS를 수행하고(
Yang et al. 2011) 이를 통한 머신러닝 기반의 예측 모델의 실용 가능성을 평가하였다. 벼(
Oryza sativa L.)는 세계적으로 가장 중요한 식량작물 가운데 하나이며, 수량과 미곡 품질을 좌우하는 조기 선발 역량은 현대 육종의 핵심 과제다(
Quille-Mamani et al. 2025). 특히 단백질 함량, 등숙률, 초장, 수장은 수량성과 품질⋅재배 안정성과 직결되는 대표 형질로, 정확하고 재현성 있는 선발 지표를 마련할 필요가 있다(
Ding et al. 2024,
Liu et al. 2022). 그러나 이들 형질은 환경 변이의 영향이 크고, 패널의 유전적 배경 차이에 의해 신호가 왜곡될 수 있어 기존 접근만으로는 정확한 선발에 한계가 있다(
Liu et al. 2022). 특히 국내 육성 품종처럼 생태형(ecotype; Japonica, Tongil-type)이 혼재된 패널에서는 집단분화 신호가 형질 연관 신호로 과대추정될 수 있으므로, 주성분(PC)과 친연도(Genomic Relationship Matrix; GRM)를 반영한 혼합모형 GWAS 보정이 필수적이다(
Kang et al. 2010,
Price et al. 2006,
VanRaden 2008,
Yang et al. 2011,
Zhou & Stephens 2012). 한편 전장서열 변이⋅GWAS 신호를 특징 선택(feature selection)으로 추려 기계학습 분류⋅예측 모델과 결합하는 디지털 육종 프레임워크가 빠르게 도입되고 있으며, 작물 및 원예 작물 전반에서 실용 사례가 축적되고 있다(
Kim et al. 2025,
Quille-Mamani et al. 2025,
Yu et al. 2021).
이에 본 연구는 국내 육성 벼 85개 품종의 전장 유전체 데이터를 바탕으로, PLINK 단변량(case⋅control) 스캔(
Purcell et al. 2007)을 기계학습용 사전 스크리닝으로 활용하고, 동시에 GCTA v1.93.2 (MLMA)에 ecotype⋅PC1, PC2 (고정효과)와 GRM (무작위효과)를 포함해 집단구조 보정과 캘리브레이션(QQ-plot)을 확인하였다. 형질별 상⋅하위 30% 이분화를 통해 후보 마커를 조합하고, Random Forest 기반 중요도 등으로 선별한 소수의 핵심 SNP 패널을 다양한 ML 모델과 결합하여 유전형 기반 조기 선발의 실용 가능성을 평가하였다.
재료 및 방법
재료 및 식물 재배 조건
본 연구는 농촌진흥청 국립식량과학원에서 분양받은 국내 육성된 대표적인 85개 벼(
Oryza sativa L.) 품종을 대상으로 수행하였다. 각 품종은 Japonica와 통일벼(Tongil-type)를 포함한다(
Table 1). 연구에 사용된 식물체는 유전형 분석의 일관성을 유지하기 위해 같은 재배 환경에서 균일하게 관리하였다. 실험에 사용된 벼는 온실 조건(주간 28±2℃, 야간 23±2℃, 광주기 14시간)에 파종하여 약 3주간 재배하였으며, 충분히 자란 어린 잎 조직을 채취하여 즉시 액체 질소로 동결하여 DNA 추출 전까지 -80℃에 보관하였다.
Genomic DNA는 CTAB (Cetyltrimethyl ammonium bromide) 방법(
Murray & Thompson 1980)을 이용하여 추출하였다. 동결한 잎 조직 약 100 mg을 액체 질소에서 곱게 분쇄한 뒤, 2% CTAB buffer를 첨가하여 65℃에서 60분간 인큐베이션 하였다. 이후 chloroform:isoamyl alcohol (24:1) 혼합액으로 2회 추출한 뒤, isopropanol을 이용하여 DNA를 침전시키고, 70% 에탄올로 세척하고 TE buffer (pH 8.0)에 용해시켰다.
추출된 DNA의 농도와 순도는 NanoDrop Spectrophotometer (ND-1000, Thermo Scientific, USA)를 이용하여 260 nm 및 280 nm 파장에서 측정하여 OD260/280 비율이 1.8-2.0인 DNA만 사용하였다. 이후 Illumina Novaseq 6000 플랫폼을 이용하여 평균 25× 이상의 coverage로 paired-end (PE) 방식으로 전장유전체분석(whole-genome resequencing)을 수행하였다.
변이분석 및 필터링
Resequencing으로 확보한 raw reads는 Trimmomatic v0.39로 어댑터 및 저품질 염기(Quality Score<20)를 제거하였다(
Bolger et al. 2014). 정제된 reads는 BWA-MEM v0.7.17을 이용해 벼 표준 참조유전체(IRGSP v1.0)에 매핑하였고(
Li 2013), SAMtools v1.10로 정렬한 뒤 Picard v2.22.0로 중복 읽기를 처리하였다. GATK v4.1.8.1 파이프라인에서 Base Quality Score Recalibration 후 HaplotypeCaller로 시료별 gVCF를 생성하고, CombineGVCFs/GenotypeGVCFs로 통합 변이를 산출하였다(
De Summa et al. 2017,
Koboldt 2020,
Van der Auwera et al. 2013). 이 과정을 통해 12,718,879개의 초기 SNP를 획득하였다.
품질 필터링은 bi-allelic SNP만 선별한 뒤 genotyping rate≥0.90을 공통 적용하였다. 벼는 자가수분 작물이며 초기 변이 세트에 야생 자원이 포함되어 Hardy-Weinberg 평형(HWE) 가정이 성립하기 어렵다고 판단되어, 본 분석에서는 HWE 필터를 적용하지 않았다. 이후 분석 목적에 맞춰 두 개의 경로로 분기하였다(
Fig. 1).
1. ML 사전 스크리닝(주 분석, PLINK)
머신러닝 입력 차원 축소를 위해 PLINK v1.9 단변량(case⋅control) 스캔을 사전 스크리닝으로 수행하였다. 여기서는 MAF≥0.03 기준을 적용하여 약 2,200,500개의 고품질 SNP를 확보하였고(표현형별 상⋅하위 30% 이분화), 이 단계는 마커 유의성 보고가 아닌 feature 랭킹을 위한 절차로만 사용하였다. 이후 Random Forest 기반 feature importance로 후보를 정제하여 최종 핵심 SNP를 단백질 함량 26개, 등숙률 51개, 초장 19개, 수장 20개로 선정하였다.
2. 보정⋅캘리브레이션(보조 분석, GCTA-MLMA)
집단구조 및 친연도 보정을 확인하기 위해 GCTA v1.93.2 (MLMA)를 별도로 적용하였다. 모델에는 ecotype (Japonica⋅Tongil-type), PC1, PC2를 고정 공변량(covariates)으로, GRM 기반 다유전자 효과(g)를 무작위효과(random effect)로 포함하였다. 이 경로에서는 보다 보수적으로 MAF≥0.05를 적용했으며, case⋅control 재정의에 따라 표현형별 유효 SNP 수는 다음과 같이 달라졌다: 단백질 함량 1,543,341개, 등숙률 1,550,409개, 초장 1,560,408개, 수장 1,565,209개로 유의성 평가는 Bonferroni (α=0.05/m; -log10P≈7.49-7.50)와 suggestive (1/m; -log10P≈6.19-6.20) 기준을 사용하였고, Manhattan과 QQ-plot으로 보정 적정성을 제시하였다.
표현형 조사 및 전장유전체 연관 분석(GWAS)
1. 표현형 조사
연구의 표현형 데이터(단백질 함량, 등숙률, 초장, 수장)는 농촌진흥청 국립식량과학원에서 육성년도 1976~2022년 재배기에 표준 절차에 따라 생산⋅관리된 자료를 측정 및 활용하였으며, 최근에는 표현형의 디지털화를 위한 표준 프로토콜 확립 및 딥러닝을 이용한 표현형 조사가 활발히 이루어지고 있다(
Cho et al. 2024,
Yu & Chung 2024). 단백질 함량은 도정 백미를 대상으로 Kjeldahl 법(
Chen et al. 2023)으로 산출하였고, 등숙률(grain filling rate)은 이삭 전체에서 여문 종자 비율(%)로 계산하였다. 초장과 수장은 성숙기에 품종별 무작위 10개체를 측정하여 평균값을 사용하였다. 생태형(Japonica⋅Tongil-type)별 표본 수와 평균을 측정하였다(
Table 2). 표현형별 데이터는 분포도를 통해 시각화하였다(
Fig. 2).
2. 데이터 분류 기준
각 형질의 관측값을 기준으로 상위 30%와 하위 30%를 각각 case와 control로 정의하였다. 이러한 극단값 기반 이분화(extreme-phenotype design)는 집단 간 대비를 높여 효과 크기를 증폭하고, 비유전적 잔차 변동을 상대적으로 줄여 분류 모델 학습에 적합한 입력을 제공하였다. 집단구조의 영향을 점검하기 위해 PLINK 분석과 병행하여 GCTA-MLMA 보정 분석을 보조적으로 수행하고 QQ-plot으로 보정 적정성을 확인하였다.
3. PLINK 기반 기계학습 입력자료 후보 구성
PLINK v1.9를 이용하여 단변량 case⋅control 분석을 수행하고, 후속 기계학습 단계에 사용할 입력 후보 변이 집합을 구성하였다. 본 절의 분석 목적은 마커 수준의 통계적 유의성 보고가 아니라 예측 모델에 투입할 변이 후보를 합리적으로 집계하는 데 있다. 집단구조 보정의 적정성과 통계적 해석은 다음 절(GCTA-MLMA)에서 상세히 기술한다. 공통 품질 기준으로 bi-allelic, genotyping rate≥0.90, MAF≥0.03을 적용하였다. 표현형별 후보 수집 기준은 단백질 함량 –log
10P≥6, 등숙률 -log
10P≥7, 초장 -log
10P≥5, 수장 -log
10P≥18로 설정하였고, 이에 따른 Manhattan plot을 제시하였다(
Fig. 3).
집계된 후보 변이는 Random Forest 기반 변수 중요도로 우선순위를 재평가하여 핵심 SNP로 축약하여 최종적으로 단백질 26개, 등숙률 51개, 초장 19개, 수장 20개의 SNP를 확보하였다. 최종 핵심 패널은 이후 분류 모델 학습 및 평가에 사용하였다.
4. GCTA-MLMA를 통한 보정 적정성 평가
집단구조와 친연도 보정을 점검하기 위해 GCTA v1.93.2 (MLMA)를 별도로 수행하였다. 모형에는 ecotype (Japonica⋅Tongil-type), PC1, PC2를 고정 공변량으로, GRM 기반 다유전자 효과(g)를 무작위 효과로 포함하였다. 이 경로에서는 MAF≥0.05의 보수적 기준을 적용하였다. 이는 본 분석의 목적이 보정 적정성(calibration)과 해석의 타당성을 확인하는 데 있으며, 표본 규모가 n=85로 제한된 상황에서 희귀 대립유전자(low MAF)는 (i) 군집별(case⋅control) 소표본에서의 희소성(sparsity)로 회귀 추정의 불안정성을 높이고, (ii) 분산 성분 및 검정 통계의 분산을 키워 변동성을 증가시킬 수 있기 때문이다. MAF 0.05 기준은 각 분석 서브셋(상⋅하위 30%)에서의 최소 minor allele count (MAC)를 확보해 검정력과 분산 추정의 안정성을 동시에 도모하기 위한 선택이다. case⋅control 재정의에 따라 유효 SNP 수는 형질별로 상이하였으며, 단백질 1,543,341개, 등숙률 1,550,409개, 초장 1,560,408개, 수장 1,565,209개였다. 유의성 평가는 Bonferroni (α=0.05/m; -log
10P≈7.49-7.50) 및 suggestive (1/m; -log
10P≈6.19-6.20) 기준을 적용하였고, 결과는 Manhattan과 QQ-plot으로 제시하였다(
Fig. 4).
표본 규모가 n=85인 조건에서 genome-wide 유의 변이는 검출되지 않았으며, QQ-plot은 과도한 인플레이션 없이 적절한 보정 상태를 시사하였다. 본 분석은 보정 적정성 및 해석상의 타당성 확인을 목적으로 수행되었으며, 머신러닝 입력 변수 선정에는 사용하지 않았다.
5. 생태형 분포 점검 및 기계학습용 최종 데이터셋 구성
각 형질에 대해 상위 30%와 하위 30%를 기준으로 정의한 case와 control의 생태형 분포(Japonica 비율)을 점검하였다. 대부분의 형질에서 두 집단 간 분포 불균형은 크지 않았으나, 수장에서는 control 집단의 Japonica 비율이 96.9%로 높게 관찰되었다. 이 분포 특성은 이후 결과 해석과 검증 설계에 반영하였다(
Table 3).
후속 예측 모델 학습을 위해 품질 기준을 충족한 변이들 가운데 PLINK 단변량 분석으로 도출된 후보를 바탕으로 입력 특성 집합을 구성하였으며, 여기서 얻은 후보 변이는 Random Forest 기반 변수 중요도로 우선순위를 재평가하여 핵심 SNP 패널로 축약하였다. 최종 패널의 크기는 형질별로 단백질 26개, 등숙률 51개, 초장 19개, 수장 20개였으며, 이들 데이터를 SVM, KNN, RF, C5.0, PLS 등 여러 분류기를 학습⋅평가하였다.
기계학습 기반 표현형 예측 모델 개발
형질별 GWAS에서 선별된 SNP 마커를 이용해 유전형 데이터를 기반으로 한 표현형 예측 모델을 구축하였다. 각 반복에서 훈련 세트 70% 내에서 PLINK 후보를 바탕으로 RF 중요도 산정 및 핵심 SNP 선정을 수행하고, 내부 5-fold 교차검증으로 튜닝을 마쳤다. 최종 평가는 검증 세트 30%에서 산출했으며, 이 과정을 seed (1-5)×반복(≥100)으로 수행해 평균 성능을 보고하였으며, 머신러닝 활용 연구사례가 급증하고 있다(
Lee et al. 2025).
적용된 기계학습 알고리즘으로는 Caret package v6.0 (
Kuhn 2008)에서 지원되는 Support Vector Machine (SVM), K-Nearest Neighbors (KNN), Random Forest (RF), C5.0 Decision Tree Algorithm (C5.0), Partial Least Squares Regression (PLS) 등 총 5가지의 모델이 사용되었다(
Fig. 5). R에서 작동하는 Caret 패키지는 다양한 머신러닝 알고리즘을 적용할 수 있어 유전체 정보를 활용하여 예측 정확도를 비교할 수 있는 장점이 있다(
Cho et al. 2022,
Yu et al. 2021). 각 모델은 내부 튜닝(caret::tuneLength=10) 및 5-fold cross-validation을 수행하였고, 모델의 예측 성능은 TP, FP, TN, FN 값을 기반으로 민감도(sensitivity), 특이도(specificity), 정확도(accuracy)를 validation set 평가하였다. 최종적으로 각 형질별로 가장 높은 예측 정확도를 보이는 모델과 해당 SNP 마커 세트를 선발하였다.
결과 및 고찰
형질별 예측 성능 요약
본 연구는 85개 품종의 전장유전체 변이를 바탕으로, PLINK 단변량(case⋅control) 스캔을 기계학습(ML) 입력용 사전 스크리닝에 활용하고, GCTA v1.93.2 (MLMA)를 별도로 수행해 집단구조⋅친연도 보정의 적정성(QQ-plot)을 확인하였다. 각 형질의 관측값을 상위 30%/하위 30%로 이분화하여(case⋅control) 후보 마커를 수집한 뒤, Random Forest 기반 feature importance로 핵심 피처를 최종 압축하였다. 이 과정에서 선별된 핵심 SNP 수는 형질별로 단백질 26개, 등숙률 51개, 초장 19개, 수장 20개였으며, 이들만을 입력으로 SVM, KNN, RF, C5.0, PLS 등 여러 분류기를 학습⋅평가하였다. Seed 반복을 통해 교차검증 결과, 평균 예측 정확도는 단백질 81.8%, 등숙률 81.0%, 초장 73.1%, 수장 94.0%로 확인되었다. 형질별 대표 모델은 단백질 PLS, 등숙률 KNN, 초장 PLS, 수장 SVM이었다(
Table 4).
이와 같은 결과는 소수의 정보성 SNP 조합만으로도 실용 수준의 조기 분류가 가능함을 시사한다.
1. 단백질 함량 결과
조사된 63개 품종의 평균 단백질 함량은 6.6%였고, 생태형별로 Japonica 6.5%, Tongil-type 6.9%였다(
Table 2). PLINK에서는 단백질 함량 값을 상⋅하위 30%로 이분화한 뒤(case⋅control) PLINK v1.9 단변량 스캔을 수행하여 후속 ML 입력용 후보를 수집했고, 수집 임계값은 -log
10P≥6로 두었다(
Fig. 3A). 동시에 GCTA-MLMA로 ecotype⋅PC1, PC2 (고정), GRM (무작위)를 포함한 혼합모형을 적용해 보정 적정성(QQ-plot)을 확인했다. 표본 규모(n=85)에서는 genome-wide 유의 신호는 검출되지 않았으며, MLMA는 보정 확인용으로만 사용했다(
Fig. 4).
후속 feature selection에서는 PLINK 후보를 대상으로 Random Forest 중요도를 계산해 핵심 SNP 26개로 압축했고, 이를 입력으로 SVM, KNN, RF, C5.0, PLS를 비교했다. 다중 seed 반복 및 학습폴드 내 feature 선정 등 교차검증 결과, 평균 정확도 81.8%를 얻었으며, PLS와 SVM이 가장 안정적 성능을 보였다(
Table 4,
Fig. 6). 이로써 소수(26개)의 정보성 마커 조합만으로도 실용적 수준의 조기 분류가 가능함을 확인했다.
2. 등숙률 결과
등숙률은 생태형 간 평균이 Japonica 83.1%, Tongil-type 76.7%로 차이가 관찰되었다(
Table 2). PLINK에서는 등숙률 값을 상⋅하위 30%로 이분화한 뒤(case⋅control) PLINK v1.9 단변량 스캔을 수행하여 후속 ML 입력용 후보를 수집했고, 수집 임계값은 -log
10P≥7로 두었다(
Fig. 3B). 동시에 GCTA-MLMA를 적용하여 보정 적정성(QQ-plot)을 확인했으며, 표본 규모(n=85)에서는 genome-wide 유의 신호는 관찰되지 않았다(
Fig. 4).
후속 feature selection에서는 PLINK 후보를 대상으로 Random Forest 중요도를 계산해 핵심 SNP 51개로 압축했고, 이를 입력으로 SVM, KNN, RF, C5.0, PLS를 비교했다. 다중 seed 반복과 교차검증(학습 폴드 내 피처 선정)을 통해 평가한 결과, 평균 정확도는 81.0%였으며, KNN과 PLS가 상위 성능을 보였다(
Table 4,
Fig. 7). 이 결과는 등숙률에 대해서도 소수의 정보성 마커 조합만으로 실용 수준의 조기 분류가 가능함을 뒷받침한다.
3. 초장 결과
평균 초장은 Japonica 76.7 cm, Tongil-type 82.9 cm였다(
Table 2). PLINK 에서는 초장 값을 상⋅하위 30%로 이분화한 뒤 PLINK v1.9 단변량 스캔으로 후속 ML 입력용 후보를 수집했고, 수집 임계값은 -log
10P≥5로 두었다(
Fig. 3C). GCTA-MLMA로 ecotype⋅주성분을 고정효과, GRM을 무작위효과로 포함한 혼합모형을 적용해 QQ-plot 정합성을 점검했으며, 본 표본 규모에서는 genome-wide 유의 신호는 관찰되지 않았다(
Fig. 4).
후속 feature selection에서는 PLINK 후보를 대상으로 Random Forest 중요도를 계산해 핵심 SNP 19개로 압축했고, 이를 입력으로 SVM, KNN, RF, C5.0, PLS를 비교했다. 다중 seed 반복 교차검증 결과, 평균 정확도 73.1%였으며 PLS가 가장 안정적이었다(
Table 4,
Fig. 8).
초장은 다유전자성이고 환경 기여도가 상대적으로 크기 때문에 다른 형질 대비 정확도가 낮게 나타난 것으로 해석된다. 문헌에서 보고된 sd1 인근의 신호는 본 패널에서도 suggestive 수준으로 관찰되었으나, 보수적 임계선(7.49)을 넘는 피크는 제한적이었다(
Dong et al. 2024).
4. 수장 결과
수장은 수량성과 밀접한 형질로, 본 패널의 평균은 22.3 cm였다(
Table 2). 수장에서는 control 집단의 생태형 분포가 Japonica 96.9%로 편중되어 있었고(
Table 3), 이를 고려해 GCTA-MLMA로 ecotype (고정)과 주성분/GRM (보정요인)을 포함한 혼합모형을 적용해 QQ-plot 정합성을 점검했다. 본 표본 규모에서는 genome-wide 유의 신호는 관찰되지 않았으며, 혼합모형은 보정 확인용으로만 사용했다(
Figs. 3D,
4D).
PLINK에서는 수장 값을 상⋅하위 30%로 이분화한 뒤(case⋅control) PLINK v1.9 단변량 스캔으로 후속 ML 입력용 후보를 모았고, 수집 임계값은 -log
10P≥6으로 두었다. 이후 Random Forest 중요도로 후보를 압축해 핵심 SNP 20개를 선정했고, 이를 입력으로 SVM, KNN, RF, C5.0, PLS를 비교했다. 다중 seed 반복 교차검증 결과, 평균 정확도 94.0%로 매우 높았으며 SVM/RF/KNN이 최상위 성능을 보였다(
Table 4,
Fig. 9). 이 결과는 본 패널에서 수장 변이가 비교적 크게 남아 있고, 소수의 정보성 마커 조합만으로도 실용 수준의 고정확도 조기 분류가 가능함을 보여준다.
본 연구는 85개 국내 육성 벼 품종의 전장유전체 변이를 바탕으로, PLINK 단변량(case⋅control) 스캔을 주 분석의 사전 스크리닝으로 활용하고, GCTA v1.93.2 (MLMA)를 결과 해석을 뒷받침하는 진단적 분석으로 수행하여 집단구조⋅친연도 보정의 적정성(Manhattan⋅QQ-plot)을 확인하였다. Resequencing으로 확보한 12,718,879개 변이는 공통 QC (bi-allelic, call rate≥0.90, MAF≥0.03) 후 약 2,200,500개의 고품질 SNP로 축약되었고, 각 형질의 관측값을 상⋅하위 30%로 이분화해 PLINK로 후속 ML 입력용 후보를 수집하였다. 이어 Random Forest 기반 중요도로 후보를 재평가하여 핵심 SNP 패널(단백질 함량 26개, 등숙률 51개, 초장 19개, 수장 20개)로 압축하고, SVM, KNN, RF, C5.0, PLS 등 분류기를 학습⋅평가하였다. 다중 seed 교차검증 결과, 단백질 81.8%, 등숙률 81.0%, 초장 73.1%, 수장 94.0%의 평균 정확도를 달성하여, 소수의 정보성 SNP 조합만으로도 실용 수준의 조기 분류가 가능함을 확인하였다.
MLMA는 ecotype, PC1, PC2 (고정효과), GRM (무작위효과)을 포함하여 모델 보정의 적합성을 점검하기 위해 수행되었으며, 표본 규모 n=85에서는 genome-wide 유의 변이가 검출되지 않았다. QQ-plot은 과도한 인플레이션 없이 보정이 양호함을 시사하였다. 최종 ML 입력 변수의 정의는 PLINK 경로(MAF≥0.03, ~2.20M SNP)에 기반하며, MLMA 결과는 해석의 타당성 검증을 위한 보조 근거로 제시되었다.
교차-생태형 일반화에 대하여. 본 연구는 표본 수와 계통 구성의 제약으로 인해 교차-생태형 일반화(Japonica 학습→Tongil-type 검증, 역도)의 정량적 추정이 불안정해질 가능성이 높다. 과대⋅과소 추정의 위험을 최소화하기 위해 본문에는 해당 수치를 포함하지 않았다. 대신 PLINK 스크리닝 이후의 혼합모형 분석에서 ecotype과 주성분, GRM을 반영해 보정의 적합성을 점검하고, QQ-plot으로 과도한 인플레이션이 없음을 확인했으며, case⋅control의 생태형 분포를 표로 제시하고, ecotype 라벨을 ML 입력 변수로 사용하지 않았음을 명확히 했다. 향후 연구에서는 표본을 확대하고 생태형 층화를 균형화한 뒤, 교차-생태형 일반화 성능을 정량 평가할 필요가 있다.
형질별로 보면, 수장은 표현형 변이가 상대적으로 크게 남아 있고 정보성 마커가 집중되어 핵심 20개 SNP만으로 94%의 높은 정확도를 보였다. 초장은 다유전자성과 환경 기여가 큰 특성, 및 단간화에 따른 집단 내 분산 축소로 인해 73.1%로 낮게 나타났다. 단백질 함량과 등숙률은 일반적으로 유전력이 낮은 편임에도, 극단값 이분화(top/bottom 30%) 및 피처 압축 전략으로 각각 81.8%, 81.0%의 실용적 성능을 확보하였다.
본 연구의 표본 수(n=85)가 제한적이어서 GWAS의 검정력과 ML 모델의 일반화 가능성이 제약될 수 있으며, 일부 서브셋의 생태형 편중이 신호 해석을 어렵게 만들 가능성이 있다. 향후 표본 수를 대폭 확대하고 생태형 층화 표집을 병행하여 보정된 혼합모형 분석의 검출력과 모델의 외부 일반화 성능을 개선하는 연구를 수행할 필요가 있다. 그럼에도 본 연구에서 제시한 핵심 SNP 패널은 Fluidigm 설계 기준을 충족하며, 표준화 칩으로 구현 시 유전형 기반 조기 선발에 적용 가능한 컴팩트 패널로 기능할 수 있음을 확인하였다. 향후 표본 수 확대를 통해 통계적 검출력과 예측 안정성을 동시에 끌어올리고, 교차-생태형⋅다환경 검증을 체계화함으로써 디지털 육종의 실용 구현을 한층 앞당길 수 있을 것으로 기대한다.
적요
전장유전체 resequencing으로 얻은 변이를 PLINK 기준(bi-allelic, call rate≥0.90, MAF≥0.03)으로 품질 관리해 약 220만개의 SNP을 확보하였고, 각 형질의 상⋅하위 30%를 이용한 단변량 스캔으로 ML 입력 후보를 확보하였다. Random Forest 기반 피처 랭킹과 압축으로 단백질 함량 26개, 등숙률 51개, 초장 19개, 수장 20개의 핵심 SNP 패널을 확정하고, 다중 알고리즘 비교에서 수장 94.0%, 단백질 81.8%, 등숙률 81.0%, 초장 73.1%의 평균 정확도를 달성했다. GCTA-MLMA는 ecotype⋅PC1, PC2⋅GRM을 포함해 보정의 적정성을 점검하기 위한 보조 분석으로 수행되었으며, 본 연구의 표본 수 85개체 에서의 genome-wide 유의 신호는 검출되지 않았다. PLINK를 통해 확보된 최종 선정 마커는 Fluidigm SNP 칩 설계 요건을 충족해, 유전형 기반 조기 선발에 적용 가능한 컴팩트 패널을 제시한다. 이는 국내 벼 품종의 디지털 육종 실현을 앞당길 실용적 단계로서 의의가 있다.
사사
본 연구는 2024년도 농촌진흥청 연구사업(과제번호: PJ01721001)의 지원으로 수행되었습니다.
Fig. 1Workflow for selecting key SNP markers in rice.
Fig. 2Data distribution for each phenotype.
Fig. 3GWAS Manhattan plots using PLINK (case⋅control design). Black line: -log10P=6 (protein content), 7 (grain filling rate), 5 (height), and 18 (panicle length).
Fig. 4GWAS results for four traits in 85 rice cultivars using GCTA v1.93.2 (MLMA; fixed covariates: ecotype, PC1, PC2; random effect: GRM). Left: Manhattan plots; right: QQ-plots. SNPs tested: (a) 1,543,341; (b) 1,550,409; (c) 1,560,408; (d) 1,565,209. Red line: Bonferroni threshold α=0.05/m (-log10P≈7.49 for a-c, 7.50 for d). Blue line: suggestive threshold 1/m (-log10P≈6.19 for a-c, 6.20 for d).
Fig. 5Schematic representation of machine learning models used in phenotype prediction: (a) Support Vector Machine (SVM), (b) K-Nearest Neighbors (KNN), (c) Random Forest (RF), (d) C5.0 Decision Tree (C5.0), and (e) Partial Least Squares Regression (PLS).
Fig. 6Protein content model performance across five random seeds. Panels (a)-(e) denote independent repeats. Within each panel: (A) ROC curves comparing SVM, KNN, RF, C5.0, and PLS; (B) genotype-based PCA with points colored by true case⋅control and shaped by ecotype; (C) predicted probability (score) distributions for case vs. control.
Fig. 7Grain filling rate model performance across five random seeds. Panels (a)-(e) denote independent repeats. Within each panel: (A) ROC curves comparing SVM, KNN, RF, C5.0, and PLS; (B) genotype-based PCA with points colored by true case⋅control and shaped by ecotype; (C) predicted probability (score) distributions for case vs. control.
Fig. 8Height model performance across five random seeds. Panels (a)-(e) denote independent repeats. Within each panel: (A) ROC curves comparing SVM, KNN, RF, C5.0, and PLS; (B) genotype-based PCA with points colored by true case⋅control and shaped by ecotype; (C) predicted probability (score) distributions for case vs. control.
Fig. 9Panicle length model performance across five random seeds. Panels (a)-(e) denote independent repeats. Within each panel: (A) ROC curves comparing SVM, KNN, RF, C5.0, and PLS; (B) genotype-based PCA with points colored by true case⋅control and shaped by ecotype; (C) predicted probability (score) distributions for case vs. control.
Table 1Rice varieties used in this study.
Table 1
|
No. |
Variety (English) |
Variety (Korean) |
Systematic Name |
Issued year |
Ecotype |
|
01 |
Chamdongjin |
참동진 |
Jeonju623 |
2020 |
japonica |
|
02 |
Odae |
오대 |
Suweon303 |
1980 |
japonica |
|
03 |
Samgwang |
삼광 |
Suweon474 |
2003 |
Tongil-type |
|
04 |
Milyang23 |
밀양23 |
Milyang23 |
1976 |
Tongil-type |
|
05 |
Giho |
기호 |
Suweon306 |
1980 |
japonica |
|
06 |
Yechan |
예찬 |
Iksan583 |
2017 |
japonica |
|
07 |
Cheongcheongbyeo |
청청벼 |
Milyang46 |
1979 |
Tongil-type |
|
08 |
Mimyeon |
미면 |
Milyang260 |
2012 |
Tongil-type |
|
09 |
Sangnambatbyeo |
상남밭벼 |
Milyang93 |
1988 |
japonica |
|
10 |
Yongmunbyeo |
용문벼 |
Suweon332 |
1985 |
Tongil-type |
|
11 |
Hangangchal 1 |
한강찰1호 |
Milyang167 |
2006 |
Tongil-type |
|
12 |
Seopyeong |
서평 |
Gyehwa22 |
2003 |
japonica |
|
13 |
Nokwoo |
녹우 |
Suweon560 |
2014 |
Tongil-type |
|
14 |
Mogwoo |
목우 |
Suweon519 |
2009 |
Tongil-type |
|
15 |
Aromi |
아로미 |
Milyang302 |
2016 |
japonica |
|
16 |
Jangseongbyeo |
장성벼 |
Iri362 |
1985 |
Tongil-type |
|
17 |
Nokyang |
녹양 |
Suweon490 |
2006 |
Tongil-type |
|
18 |
Samgangbyeo |
삼강벼 |
Milyang55 |
1982 |
Tongil-type |
|
19 |
Jinbo |
진보 |
Yeongdeog45 |
2009 |
japonica |
|
20 |
Geonganghongmi |
건강홍미 |
Milyang234 |
2010 |
japonica |
|
21 |
Hyangmibyeo 1 |
향미벼1호 |
Suweon393 |
1993 |
Tongil-type |
|
22 |
Joeunheukmi |
조은흑미 |
Cheolwoen80 |
2011 |
japonica |
|
23 |
DASAN1HO |
다산1호 |
Suweon499 |
2006 |
Tongil-type |
|
24 |
Namil |
남일 |
Suweon472 |
2002 |
japonica |
|
25 |
Joan |
조안 |
Suweon478 |
2003 |
japonica |
|
26 |
Borami |
보라미 |
Milyang211 |
2007 |
japonica |
|
27 |
Heukjinjubyeo |
흑진주 |
Suweon415 |
1997 |
japonica |
|
28 |
Nampungbyeo |
남풍벼 |
Suweon294 |
1981 |
Tongil-type |
|
29 |
Hangangchalbyeo |
한강찰벼 |
Suweon290 |
1979 |
Tongil-type |
|
30 |
Saemimyeon |
새미면 |
Milyang278 |
2014 |
Tongil-type |
|
31 |
Boseogheugchal |
보석흑찰 |
Suweon512 |
2008 |
japonica |
|
32 |
MY299BK |
MY299BK |
Milyang299 |
2016 |
japonica |
|
33 |
Hangaru |
한가루 |
Suweon594 |
2016 |
japonica |
|
34 |
Heugnambyeo |
흑남벼 |
Iri427 |
1997 |
japonica |
|
35 |
Taebaekbyeo |
태백벼 |
Suweon287 |
1979 |
Tongil-type |
|
36 |
Cheongdam |
청담 |
Suweon498 |
2006 |
japonica |
|
37 |
Yeongdeogbyeo |
영덕벼 |
Yeongdeog3 |
1985 |
japonica |
|
38 |
Seolbaek |
설백 |
Cheolweon76 |
2010 |
japonica |
|
39 |
Gangbaek |
강백 |
Iksan478 |
2006 |
japonica |
|
40 |
Mogyang |
목양 |
Suweon525 |
2010 |
Tongil-type |
|
41 |
Hwaseongbyeo |
화성벼 |
Suweon330 |
1985 |
japonica |
|
42 |
Cheongwoo |
청우 |
Suweon585 |
2016 |
Tongil-type |
|
43 |
CW92MR |
CW92MR |
Cheolweon92 |
2016 |
japonica |
|
44 |
Hanareumchal |
한아름찰 |
Milayng288 |
2015 |
Tongil-type |
|
45 |
Jinbuolbyeo |
진부올벼 |
Jinbu11 |
1991 |
japonica |
|
46 |
Miwoo |
미우 |
Suweon597 |
2017 |
Tongil-type |
|
47 |
Jinmibyeo |
진미벼 |
Suweon349 |
1989 |
japonica |
|
48 |
Seomyeong |
서명 |
Gyehaw30 |
2009 |
japonica |
|
49 |
Nonganbyeo |
농안벼 |
Suweon392 |
1993 |
japonica |
|
50 |
Areumbyeo |
아름벼 |
Milyang160 |
1999 |
Tongil-type |
|
51 |
Jeogjinjuchal |
적진주찰 |
Suweon524 |
2010 |
japonica |
|
52 |
Dacheong |
다청 |
Iksan495 |
2008 |
japonica |
|
53 |
Cheongbaekchal |
청백찰 |
Cheolweon77 |
2010 |
japonica |
|
54 |
HONGJINJU |
홍진주 |
Suweon501 |
2006 |
japonica |
|
55 |
Joryeongbyeo |
조령벼 |
Milyang107 |
1992 |
japonica |
|
56 |
Yeongwoo |
영우 |
Suweon573 |
2015 |
Tongil-type |
|
57 |
Baegjinju1ho |
백진주1호 |
Suweon491 |
2005 |
japonica |
|
58 |
Palbangmi |
팔방미 |
Suweon543 |
2012 |
Tongil-type |
|
59 |
Anmi |
안미 |
Suweon523 |
2010 |
japonica |
|
60 |
Geonyangmi |
건양미 |
Suweon533 |
2011 |
japonica |
|
61 |
Baegokchal |
백옥찰 |
Milyang225 |
2008 |
japonica |
|
62 |
Andabyeo |
안다벼 |
Suweon431 |
1998 |
Tongil-type |
|
63 |
Jinseolchal |
진설찰 |
Jinbu50 |
2012 |
japonica |
|
64 |
Dasan2 |
다산2호 |
Suweon518 |
2009 |
Tongil-type |
|
65 |
Nunkeunheugchal1ho |
눈큰흑찰1호 |
Milyang282 |
2014 |
japonica |
|
66 |
MY298BB |
MY298BB |
Milyang298 |
2016 |
japonica |
|
67 |
Daeripbyeo1 |
대립벼1호 |
Suweon391 |
1993 |
japonica |
|
68 |
Seonhyangheukmi |
선향흑미 |
Suweon532 |
2011 |
japonica |
|
69 |
Unilchal |
운일찰 |
Unbong52 |
2014 |
japonica |
|
70 |
Hanareum4 |
한아름4호 |
Milyang295 |
2016 |
Tongil-type |
|
71 |
Ilpumbyeo |
일품벼 |
Suweon355 |
1990 |
japonica |
|
72 |
Jungwonbyeo |
중원벼 |
Suweon325 |
1984 |
Tongil-type |
|
73 |
Manmi |
만미 |
Milyang162 |
2002 |
japonica |
|
74 |
Jogwang |
조광 |
Milyang213 |
2007 |
japonica |
|
75 |
Cheongun |
청운 |
Suweon537 |
2012 |
japonica |
|
76 |
Saeodae1ho |
새오대1호 |
Cheolwon103 |
2022 |
japonica |
|
77 |
Chamjinmi |
참진미 |
Jeonju654 |
2022 |
japonica |
|
78 |
Hangadeuk |
한가득 |
Suweon639 |
2022 |
japonica |
|
79 |
Dangchanjinmi |
당찬진미 |
Jeonju632 |
2022 |
japonica |
|
80 |
Amissal |
아미쌀 |
Jeonju653 |
2022 |
Tongil-type |
|
81 |
Dapyeong |
다평 |
Milyang363 |
2022 |
japonica |
|
82 |
Amimyeon |
아미면 |
Milyang355 |
2022 |
Tongil-type |
|
83 |
Jinokchal |
진옥찰 |
Milyang366 |
2022 |
japonica |
|
84 |
Chamnuri |
참누리 |
Jeonju636 |
2021 |
japonica |
|
85 |
Chohong |
초홍 |
Milyang356 |
2021 |
japonica |
Table 2Number of samples (left) and phenotype measurements (right) for protein content, grain filling rate, height, and panicle length according to rice ecotype.
Table 2
Ecotype
phenotype |
Japonica |
Tongil-type |
Total |
|
Ecotype
phenotype |
Japonica |
Tongil-type |
Average |
|
Protein content |
45 |
18 |
63 |
|
Protein content (%) |
6.5 |
6.9 |
6.6 |
|
Grain Filling Rate |
54 |
26 |
80 |
|
Grain Filling Rate (%) |
83.1 |
76.7 |
81.0 |
|
Height |
56 |
28 |
84 |
|
Height (cm) |
76.7 |
82.9 |
78.7 |
|
Panicle length |
45 |
21 |
66 |
|
Panicle length (cm) |
20.9 |
25.2 |
22.3 |
Table 3Ecotype composition of case⋅control groups defined by the top/bottom 30% per phenotype.
Table 3
|
phenotype |
Bottom 30% |
Top 30% |
Case
n |
Case
Japonica n |
Case
Japonica (%) |
Control n |
Control Japonica n |
Control Japonica (%) |
|
Protein content (%) |
6.1 |
7.1 |
20 |
13 |
65 |
20 |
16 |
80 |
|
Grain filling rate (%) |
79.07 |
86 |
26 |
20 |
76.9 |
24 |
13 |
54.2 |
|
Height (cm) |
74 |
83 |
26 |
15 |
57.7 |
28 |
21 |
75 |
|
Panicle length (cm) |
21 |
23.5 |
20 |
5 |
25 |
32 |
31 |
96.9 |
Table 4Summary of the primary ML model per phenotype.
Table 4
|
Phenotype |
Core SNPs |
Primary model |
Accuracy |
Sensitivity |
Specificity |
|
Protein content |
26 |
PLS |
0.818 |
0.833 |
0.840 |
|
Grain Filling Rate |
51 |
KNN |
0.810 |
0.956 |
0.714 |
|
Height |
19 |
PLS |
0.731 |
0.850 |
0.711 |
|
Panicle length |
20 |
SVM |
0.940 |
0.800 |
0.985 |
References
- 1. Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114-2120.
- 2. Chen S, Qi G, Zhang L, Duan X, Bai M, Hu M, Li P, Zhao W, Sun X, Guo Y, Chen W, Wang Z. 2023. Detection of salmon meat freshness using QCM gas sensor array combined with physicochemical method. Microchem J 194: 109353.
- 3. Cho E, Cho S, Kim M, Ediriweera TK, Seo D, Lee SS, Cho CY, Lee JH. 2022. Single nucleotide polymorphism marker combinations for classifying Yeonsan Ogye chicken using a machine learning approach. J Anim Sci Technol 64: 830
- 4. Cho Y, Kim JY, Ha J. 2024. Application of deep learning technology for phenotyping tissue specific length of sprout vegetables using YOLOv8. Korean Society of Breeding Science 56: 416-424.
- 5. De Summa S, Malerba G, Pinto R, Mori A, Mijatovic V, Tommasi S. 2017. GATK hard filtering: tunable parameters to improve variant calling for next generation sequencing targeted gene panel data. BMC Bioinformatics 18(Suppl 5):119
- 6. Ding H, Wang C, Cai Y, Yu K, Zhao H, Wang F, Li T, Bai Z, Tang L, Cui F. 2024. Characterization of a wheat stable QTL for spike length and its genetic effects on yield-related traits. BMC Plant Biol 24: 292
- 7. Dong J, Ma Y, Hu H, Wang J, Yang W, Fu H, Wang Z, Zhang S. 2024. The function of SD1 on shoot length and its pyramiding effect on shoot length and plant height in rice (Oryza sativa L.). Rice 17: 21
- 8. Kang HM, Sul JH, Service SK, Zaitlen NA, Kong SY, Freimer NB, Sabatti C, Eskin E. 2010. Variance component model to account for sample structure in genome-wide association studies. Nat Genet 42: 348-354.
- 9. Kim EG, Jadamba C, Yoo SC. 2025. Identification of genes conferring nitrogen deficiency tolerance by GWAS. Plant Breed Biotechnol 13: 33-52.
- 10. Koboldt DC. 2020. Best practices for variant calling in clinical sequencing. Genome Medicine 12: 91
- 11. Kuhn M. 2008. Building predictive models in R using the caret package. J Stat Softw 28: 1-26.
- 12. Lee SY, Han JH, Bak HJ, Ha SK, Lee HS, Lee G, Kim K, Mo Y. 2025. Machine learning-based heading date QTL detection in rice. Plant Breed Biotechnol 13: 108-118.
- 13. Li H. 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv: 1303.3997.
- 14. Liu M, Fan F, He S, Guo Y, Chen G, Li N, Wang H, Li S. 2022. Creation of elite rice with high-yield, superior-quality and high resistance to brown planthopper based on molecular design. Rice 15: 17
- 15. Murray MG, Thompson W. 1980. Rapid isolation of high molecular weight plant DNA. Nucl Acids Res 8: 4321-4326.
- 16. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. 2006. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet 38: 904-909.
- 17. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, Sham PC. 2007. PLINK: A tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics 81: 559-575.
- 18. Quille-Mamani J, Ramos-Fernández L, Huanuqueño-Murillo J, Quispe-Tito D, Cruz-Villacorta L, Pino-Vargas E, Silva-López J, Ángel Ruiz L. 2025. Rice yield prediction using spectral and textural indices derived from UAV imagery and machine learning models in Lambayeque, Peru. Remote Sensing 17: 632
- 19. Van der Auwera GA, Carneiro MO, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A, Jordan T, Shakir K, Roazen D, Thibault J, Banks E, Garimella K, Altshuler D, Gabriel S, DePristo MA. 2013. From FastQ data to high-confidence variant calls: The genome analysis toolkit best practices pipeline. Current Protocols in Bioinformatics 43: 11.10.1-11.10.33.
- 20. VanRaden PM. 2008. Efficient methods to compute genomic predictions. J Dairy Sci 91: 4414-4423.
- 21. Yang J, Lee SH, Goddard ME, Visscher PM. 2011. GCTA: A tool for genome-wide complex trait analysis. Am J Hum Genet 88: 76-82.
- 22. Yu GE, Shin Y, Subramaniyam S, Kang SH, Lee SM, Cho C, Kim D, Kim CK. 2021. Machine learning, transcriptome, and genotyping chip analyses provide insights into SNP markers identifying flower color in Platycodon grandiflorus. Sci Rep 11: 8019
- 23. Yu JK, Chung YS. 2024. Image-based digitalization of germplasm stock: Overcoming its limitations. Korean Society of Breeding Science 56: 293-300.
- 24. Zhou X, Stephens M. 2012. Genome-wide efficient mixed-model analysis for association studies. Nat Genet 44: 821-824.
