서론
Pseudomonas aeruginosa (PAO1) 처리에 따른 애기장대 전사체 변화를 분석한 RNA-Seq 프로젝트에서, 본격적인 DEG 분석에 들어가기 전 가장 먼저 챙긴 것은 "이 분석을 다른 사람도 똑같이 재현할 수 있는가"였습니다.
분석 요청서에는 이런 원칙이 명시되어 있었습니다.
분석에 사용되는 소프트웨어의 버전 및 주요 파라미터는 별도의 지시가 없는 경우 기본(Default) 설정을 사용한다.
그리고 분석을 누구나 재현할 수 있게 분석을 워크플로우 형태로 구성할 것.
많은 내용은 아니지만, 이걸 실제로 지키려면 생각보다 많은 걸 기록하고 확인해야 합니다.
어떤 도구를 어떤 버전으로 썼는지, 어떤 옵션을 왜 그 값으로 정했는지, 그리고 다른 사람이 똑같은 순서로 다시 실행할 수 있는 구조인지까지 신경 써야 하기 때문입니다.
이 작업이 필요했던 이유는 단순히 규칙이나 가이드로 제시되었기 때문만은 아니었습니다.
이번 프로젝트는 신뢰할 수 있는 분석 체계를 세우고, 그 체계로 데이터를 검증하고, 핵심 발견을 하고, 그 발견을 다시 두 방향(공간적 확장, 방법론적 견고성)으로 재검증하는 구조로 설계되어 있었습니다.
재현 가능한 파이프라인을 구축하는 것은 이 분석을 진행하는 데 있어서 전체 구조를 구성하는데 첫 번째 전제였습니다.
첫 시작이 불안하면 뒤에 나오는 모든 검증도 의미가 없어진다고 생각하기 때문입니다.
워크플로우: 두 경로로 나눠진 이유
전체 파이프라인은 Tuxedo Protocol의 Quantification Mode를 기반으로, fastp(전처리) → TopHat(정렬) → Cuffquant(정량화)까지는 하나의 경로로 진행되다가 이후 두 갈래로 나뉩니다.
- 발현량 비교용: Cuffnorm(16개 샘플 전체 발현량) / Cuffdiff(Root 8개 샘플 DEG)
- 방법론 비교용: htseq-count → TCC(Root 8개 샘플, 통계 모델을 바꿔도 같은 결론이 나오는지 확인)
경로를 나눈 이유는 명확합니다.
Cuffdiff와 TCC는 완전히 다른 통계 모델을 쓰는 도구인데, 같은 데이터로 같은 결론이 나온다면 그 결론은 특정 도구의 우연한 산물이 아니라는 근거가 됩니다.
(이 비교는 시리즈 뒷부분에서 자세히 다룹니다.)
여기서 한 가지 설계 포인트가 있습니다.
DEG 분석(Cuffdiff/TCC)에는 16개 전체가 아니라 Root 8개 샘플만 사용했습니다.
Root와 Shoot을 한꺼번에 넣고 비교하면 조직 간 차이라는 훨씬 큰 변동 요인이 통계 모델에 섞여 들어가, PAO1 처리 효과만을 순수하게 보기 어려워지기 때문입니다.
TopHat 결과물인 accepted_hits.bam이 두 경로 모두의 공통 출발점이 되고, 이후 각 단계는 *.cxb, *.fpkm_table, *.diff, *.count 같은 실제 파일로 연결됩니다.
이렇게 파일 단위로 흐름을 명확히 해두면, 나중에 문제가 생겼을 때 어디까지는 정상이고 어디부터 문제인지 좁혀가기 쉬워집니다.
(다음 글에서 이 내용이 실제로 어떻게 쓰였는지 또한 설명할 예정니다.)
조건을 기록한다는 것의 구체적인 의미
버전과 파라미터를 기록한다는 것은 단순히 이 프로그램을 사용했다고 기록하는 게 아닙니다.
bowtie2, fastp, TopHat, Cufflinks, HTSeq, TCC, samtools까지 각 단계별로 사용한 패키지와 버전을 표로 남겼고, Reference Genome(TAIR10)과 GTF 파일도 정확히 명시했습니다.
파라미터 원칙은 "요청서에 명시된 옵션은 그대로 반영하고, 명시되지 않은 옵션은 소프트웨어 기본값을 쓴다"였습니다.
이 중 몇 가지는 기본값이 아니라 직접 계산해서 넣은 값이라 짚어볼 만합니다.
TopHat의 -r(mate-inner-dist) 250은 라이브러리의 insert size(450bp)에서 read length의 두 배(100bp × 2 = 200bp)를 뺀 값입니다.
Paired-end 시퀀싱에서 두 read 사이의 실제 간격을 TopHat에 명시해주는 옵션인데, 시퀀싱 조건에 맞춰 직접 계산해서 확인한 후에 적용했습니다.
mask.gtf도 별도로 만들었습니다.
TAIR10 GTF에서 gene_biotype이 protein_coding이 아닌 유전자(rRNA, tRNA 등)를 제외한 파일입니다.
rRNA나 tRNA는 발현량이 극단적으로 높아서, 이런 유전자가 정규화 계산에 그대로 섞이면 전체 발현량 정규화가 왜곡될 수 있습니다.
이 파일을 Cuffquant와 Cuffdiff 단계에 동일하게 적용해 정규화 왜곡을 막았습니다.
이렇게 정리한 실행 명령어는 rnaseq.cmd.sh(메인 파이프라인)와 add_analysis.cmd.sh(보조 분석)에 순서대로 기록해서 소스 코드로 함께 제출했습니다.
사람이 아니라 스크립트가 반복하게 만들기
기록을 남기는 것뿐만 아니라 기록이 실제로 실행 가능한 구조로 남아있어야 재현성이 의미가 있습니다.
가장 핵심적인 설계 결정은 16개 샘플을 하나의 반복문으로 묶지 않고 샘플별 독립 명령어로 구성한 것입니다.
반복문이 더 깔끔해 보일 수 있지만, 독립적으로 구성하면 특정 샘플에서 오류가 나도 다른 샘플 처리에 영향을 주지 않습니다.(반복문 또한 스크립트에 작성은 해둠.)
이 구조 덕분에 실제로 문제가 생긴 샘플 몇 개만 골라서 재실행할 수 있었는데, 이 이야기는 다음 글에서 이어집니다.
모든 파라미터는 명령어에 직접 명시적으로 기재해서 별도 설정 파일 없이도 스크립트 자체가 곧 파라미터 기록이 되도록 했습니다.
그리고 단계별로 필요한 부분만 골라 재실행할 수 있도록 sed로 특정 구간을 지정해서 실행했습니다.
mkdir -p logs
nohup bash -c "sed -n '{시작},{끝}p' rnaseq.cmd.sh | /bin/bash" >> logs/run.log 2>&1 &
전체를 처음부터 다시 돌리고 싶을 때는 /bin/bash rnaseq.cmd.sh 한 줄이면 됩니다.
결과적으로 다른 분석자가 이 파이프라인을 재현하려면 필요한 건 세 가지뿐입니다. conda 환경(edurnaseq), 스크립트 두 개, 그리고 원본 Raw data·Reference·GeneModel 데이터가 존재하는 디렉토리입니다.
다음 편 예고
이렇게 구조를 짜놓은 게 실제로 쓸모가 있었던 순간이 있었습니다.
TopHat을 실행하고 결과물 개수를 세어봤는데, 16개가 나와야 할 BAM 파일이 12개밖에 나오지 않았습니다.
다음 편에서는 이 문제를 로그만으로 어떻게 추적해서 원인을 찾고 해결했는지 다룹니다.
🔗 전체 분석 코드와 실행 로그는 GitHub 저장소에서 확인할 수 있습니다. github.com/ourkofe/rna_analysis_report
'한국바이오협회 유전체 분석 과정 > RNA-seq 분석 프로젝트' 카테고리의 다른 글
| [RNA-Seq 분석 프로젝트 로그 #6] RNA-Seq 분석, 결론을 두 방향으로 재검증하기 (Root→Shoot 확장 & Cuffdiff vs TCC) (0) | 2026.08.30 |
|---|---|
| [RNA-Seq 분석 프로젝트 로그 #5] RNA-Seq DEG는 어떻게 골라내고 해석하는가 (0) | 2026.08.30 |
| [RNA-Seq 분석 프로젝트 로그 #4] RNA-Seq 이상치 샘플, 포함할까 뺄까 : 판단 기준과 실제 사례 (0) | 2026.08.30 |
| [RNA-Seq 분석 프로젝트 로그 #3] RNA-Seq 파이프라인, 산출물 개수가 안 맞을 때 로그로 원인을 찾는 법 (0) | 2026.08.30 |
| [RNA-Seq 분석 프로젝트 로그 #1] PAO1과 애기장대, 무엇을 분석하려 했을까 (0) | 2026.08.18 |