RNA-seq 실습: FASTQ에서 발현량 행렬까지

1. 무엇을 실습하나

이 실습에서 답할 질문은 하나입니다.

paired-end RNA-seq FASTQ가 믿을 만한 발현량 표로 바뀌려면, 각 단계에서 무엇을 만들고 무엇을 확인해야 하는가?

최종 산출물은 유전자별 TPM 발현량 표와 raw count 표입니다. 그러나 파일 두 개가 생겼다는 사실만으로 분석이 끝나지는 않습니다. 트리밍 전후의 품질, 레퍼런스 정합성, 정렬률과 count 배정률을 차례로 확인해야 뒤의 발현 분석을 신뢰할 수 있습니다.

전체 흐름은 다음과 같습니다.

FastQC → Trim Galore! → FastQC → STAR → RSEM + featureCounts

원리를 먼저 이해하고 싶다면 FASTQ 품질 확인과 트리밍, RNA read 정렬, 발현량 행렬 만들기를 읽은 뒤 시작하세요.

paired-end FASTQ는 무엇을 담고 있나

RNA-seq에서는 조직(tissue)이나 세포에서 꺼낸 RNA를 짧은 fragment로 만들고, sequencer가 각 fragment의 염기서열을 읽습니다. paired-end 방식은 fragment의 한쪽 끝을 R1, 반대쪽 끝을 R2에서 읽습니다. 두 파일의 같은 번째 레코드는 같은 fragment에서 나온 한 쌍입니다.

FASTQ 레코드 하나는 네 줄입니다. 아래 값은 형식을 보여 주기 위한 짧은 가상 예시입니다.

@read0001/1
ACGTTGCA
+
IIIIIHGF

줄	가상 값	의미
1	`@read0001/1`	read 식별자와 pair 방향
2	`ACGTTGCA`	sequencer가 읽은 염기서열
3	`+`	염기서열과 quality 구분선
4	`IIIIIHGF`	각 염기의 품질을 문자로 표현한 Phred quality

두 번째 줄과 네 번째 줄의 문자 수는 같아야 합니다. R1의 첫 레코드와 R2의 첫 레코드는 같은 fragment를 가리켜야 하며, 한쪽 파일만 일부 잘렸다면 짝이 어긋날 수 있습니다. 첫 프롬프트가 파일명뿐 아니라 레코드 구조와 pair를 확인하는 이유입니다.

FASTQ에는 유전자 이름이 없습니다. ACGTTGCA... 같은 짧은 서열만 있으므로, STAR가 이 서열을 레퍼런스 게놈의 위치에 놓아야 어느 유전자에서 왔는지 판단할 수 있습니다. GTF는 그 위치에 어떤 exon과 유전자가 있는지 알려 주는 좌표표입니다.

입력과 출력

입력	역할
`sample_R1.fastq.gz`, `sample_R2.fastq.gz`	같은 fragment의 양쪽 끝을 읽은 paired-end raw read
레퍼런스 게놈 FASTA	STAR가 read의 위치를 찾는 기준
같은 릴리스의 GTF	유전자와 exon의 좌표 주석
STAR 인덱스	게놈 정렬에 사용하는 사전 계산 파일
RSEM 레퍼런스	transcript와 gene 발현량 계산에 사용하는 사전 계산 파일

단계	주요 출력	다음에 하는 일
FastQC	HTML, ZIP 리포트	원본 read의 문제를 찾음
Trim Galore!	`_val_1.fq.gz`, `_val_2.fq.gz`	adapter와 낮은 품질 말단을 정리함
STAR	게놈 BAM, transcriptome BAM, 정렬 로그	위치와 정렬 품질을 확인함
RSEM	`genes.results`	TPM과 expected count를 얻음
featureCounts	gene별 count 파일과 summary	차등 발현 분석용 raw count를 얻음

raw read에서 숫자 하나까지

최종 표의 ENSG... = 1250 같은 숫자는 FASTQ에 직접 들어 있던 값이 아닙니다. 다음 판단을 통과한 read 또는 fragment를 모아 만든 값입니다.

quality가 너무 낮거나 adapter가 남은 부분을 정리합니다.
read pair가 레퍼런스의 어느 위치에서 왔는지 찾습니다.
그 위치가 GTF의 어느 유전자 exon과 겹치는지 확인합니다.
한 유전자에 배정된 fragment를 세거나, 여러 transcript에 걸친 가능성을 모델로 나눕니다.

featureCounts의 raw count는 3번과 4번에서 비교적 직접 배정된 fragment 수입니다. RSEM의 expected_count는 여러 transcript에 걸쳐 애매한 read를 확률적으로 나눈 기대값이라 소수가 될 수 있습니다. TPM은 유전자 길이와 sample 안의 sequencing depth를 함께 보정해, 한 sample 안에서 유전자들의 상대적 발현 비중을 비교하기 쉽게 만든 값입니다.

이 실습에서는 한 sample을 끝까지 처리합니다. 여러 조건의 차등 발현을 검정하려면 같은 절차로 모든 sample을 처리한 뒤 raw count 열을 하나의 행렬로 합쳐야 합니다. 한 sample의 TPM이 높다는 사실만으로 조건 간 증가나 감소를 말할 수는 없습니다.

실습에서 확인할 순서

에이전트에게 한 단계씩 요청하면서 다음 질문을 확인합니다.

R1과 R2는 정말 같은 sample의 완전한 한 쌍인가?
낮은 품질과 adapter 문제는 어느 위치에 있으며 트리밍 뒤 줄었는가?
read 대부분이 기대한 레퍼런스에 정렬되는가?
게놈 BAM과 transcriptome BAM이 각각 올바른 정량 도구로 들어갔는가?
RSEM과 featureCounts 결과의 행과 값은 무엇을 뜻하는가?
다른 컴퓨터에서 같은 분석을 다시 실행할 정보가 남아 있는가?

2. 에이전트에게 작업 환경 점검시키기

Codex CLI나 Claude Code를 입력 FASTQ가 있는 작업 폴더에서 실행한 뒤 아래 프롬프트를 전달합니다. 경로를 알고 있다면 프롬프트의 대괄호 부분을 실제 값으로 바꾸세요. 모른다면 에이전트가 현재 폴더에서 후보를 찾고, 선택이 필요한 지점에서 멈춥니다.

실습 환경 점검 프롬프트

paired-end RNA-seq FASTQ를 발현량 표로 만드는 실습 환경을 점검해줘.

입력 후보:
- FASTQ: 현재 폴더 아래의 *.fastq.gz 또는 *.fq.gz
- STAR index: [경로를 알면 입력, 아니면 후보 탐색]
- RSEM reference prefix: [경로를 알면 입력, 아니면 후보 탐색]
- GTF: [경로를 알면 입력, 아니면 후보 탐색]

요구사항:
1. 아직 분석 명령은 실행하지 마.
2. paired-end FASTQ 후보를 찾고 파일명, 크기, R1/R2 짝을 표로 정리해. sample 이름을 어떻게 정했는지도 설명해.
3. FastQC, Trim Galore!, STAR, RSEM, featureCounts, samtools의 설치 여부와 버전을 확인해. 설치되지 않은 도구를 임의로 설치하지 마.
4. STAR index의 genomeParameters.txt, RSEM reference prefix의 관련 파일, GTF가 실제로 존재하는지 확인해.
5. 확인 가능한 범위에서 게놈 빌드와 annotation 릴리스가 서로 맞는지 검사해. 파일명만으로 추정했다면 추정이라고 표시해.
6. 원본 FASTQ와 레퍼런스는 수정하지 마. 작업 산출물은 qc/raw, trimmed, qc/trimmed, align, quant, counts 아래에 저장하는 계획을 세워.
7. 사용할 sample 이름, 절대 경로, 발견한 도구 버전, 누락된 입력과 도구, 다음 단계에서 실행할 명령의 개요를 한국어로 보고하고 멈춰.

입력이나 R1/R2 짝이 모호하면 임의로 고르지 말고 가능한 후보를 보여준 뒤 멈춰.

환경 점검 결과에서 가장 먼저 볼 것은 R1과 R2의 짝입니다. 서로 다른 sample을 한 쌍으로 넣으면 이후 모든 파일은 정상적으로 생성될 수 있지만 생물학적으로 의미 없는 결과가 됩니다.

도구가 없다면 에이전트에게 사용 중인 운영체제와 패키지 관리 방식에 맞는 설치 계획을 따로 요청하세요. 이 문서는 도구를 전역 환경에 자동 설치하게 하지 않습니다. 기존 연구 환경을 덮어쓰거나 서로 다른 버전이 섞이는 일을 피하기 위해서입니다.

준비 뒤 사용할 폴더

환경 점검이 끝나면 에이전트가 다음 구조를 제안합니다. 아직 실행하지 않았으므로 입력 FASTQ와 레퍼런스를 제외한 출력 폴더는 비어 있어도 됩니다.

qc/raw/: 트리밍 전 FastQC HTML과 ZIP
trimmed/: 원본을 보존한 채 새로 만든 paired FASTQ
qc/trimmed/: 트리밍 후 FastQC HTML과 ZIP
align/: STAR BAM, 정렬 로그와 실행 기록
quant/: RSEM의 유전자·아이소폼 결과와 실행 기록
counts/: featureCounts 결과표, summary와 실행 기록

이 구조에서 raw FASTQ → trimmed FASTQ → BAM → 발현량 표의 입력과 출력이 섞이지 않습니다. sample이라는 이름도 모든 단계에서 같아야 합니다. 예를 들어 입력이 tumor01_R1.fastq.gz와 tumor01_R2.fastq.gz라면 STAR prefix와 RSEM 결과도 tumor01로 유지하는 편이 여러 sample을 합칠 때 안전합니다.

3. 프롬프트로 한 단계씩 처리하기

아래 프롬프트는 위에서부터 하나씩 전달합니다. 각 단계가 끝나면 에이전트의 요약만 읽지 말고 실제 리포트와 로그를 함께 확인한 뒤 다음 단계로 넘어가세요.

3-1. raw FASTQ 구조와 품질 확인

첫 단계에서는 read를 바꾸지 않습니다. gzip 무결성과 FASTQ 레코드 구조를 검사하고, FastQC로 원본 상태를 기록합니다.

1단계 프롬프트: raw FASTQ와 FastQC

앞에서 확인한 paired-end FASTQ에 대해 원본 품질 검사만 진행해줘.

1. 선택한 R1과 R2에 gzip 무결성 검사를 실행해. 원본 파일을 풀거나 수정하지 마.
2. 각 파일의 첫 FASTQ 레코드를 안전하게 읽어 헤더 형식, 염기서열 길이와 quality 문자열 길이가 맞는지 확인해. 출력에는 전체 염기서열이나 개인 식별 가능성이 있는 헤더를 그대로 노출하지 말고 구조만 요약해.
3. qc/raw 디렉터리를 만들고 FastQC를 실행해. 사용한 FastQC 버전과 명령을 기록해.
4. 생성된 HTML과 zip 결과에서 total sequences, sequence length, per base sequence quality, adapter content, duplication 수준을 R1과 R2로 나눠 표로 정리해.
5. 각 항목을 PASS/WARN/FAIL 문자만으로 판정하지 말고 실제 수치와 그래프 패턴을 설명해.
6. 결과 파일의 경로와 크기를 확인해. 오류가 있으면 원인을 조사하되 트리밍이나 정렬로 넘어가지 마.
7. 원본에서 관찰된 문제, 트리밍이 필요한 근거, 다음 단계에서 유지해야 할 최소 read 길이와 quality 기준의 제안을 한국어로 보고하고 멈춰.

아직 Trim Galore!, STAR, RSEM, featureCounts는 실행하지 마.

FastQC에서는 전체 read 수, 염기별 품질, read 길이 분포와 adapter 신호를 함께 봅니다. RNA-seq는 일부 유전자가 매우 많이 발현돼 duplication이 높을 수 있으므로, duplication의 빨간 표시 하나만으로 실패라고 판정하지 않습니다. 각 모듈의 의미는 QC(FastQC)에서 확인할 수 있습니다.

FastQC HTML을 열면 그래프가 많지만 처음에는 다음 네 항목부터 읽습니다.

항목	그래프의 축이나 값	먼저 확인할 것
Basic Statistics	read 수, 길이, GC 비율	R1과 R2의 read 수가 같고 예상한 길이인가
Per base sequence quality	x축은 read 위치, y축은 Phred score	3′ 말단으로 갈수록 품질이 얼마나 떨어지는가
Sequence Length Distribution	x축은 read 길이	raw read 길이가 하나인지 여러 길이가 섞였는가
Adapter Content	x축은 read 위치, y축은 adapter 검출 비율	read 끝에서 특정 adapter 신호가 올라가는가

Phred score 20은 염기 하나가 잘못 읽혔을 확률을 약 1%로 표현한 값이고, 30은 약 0.1%입니다. 그래프의 노란 상자가 특정 위치에서 아래로 내려간다면 그 위치의 많은 read에서 품질이 나빠졌다는 뜻입니다. read 몇 개의 최솟값이 낮다는 사실과 전체 분포가 낮아진 현상을 구분해야 합니다.

R1과 R2의 Total Sequences가 다르면 트리밍 전에 원인을 확인합니다. 정상적인 paired FASTQ라면 두 파일에는 같은 수의 레코드가 있어야 합니다. 파일 전송이 중간에 끊겼거나 서로 다른 처리 단계의 파일을 짝지었을 수 있습니다.

3-2. 트리밍 전후 비교

트리밍의 목적은 점수를 예쁘게 만드는 것이 아니라, 정렬을 방해하는 adapter와 신뢰하기 어려운 말단을 필요한 만큼만 제거하는 것입니다.

2단계 프롬프트: Trim Galore!와 전후 비교

raw FastQC 결과를 유지하고 트리밍과 트리밍 후 품질 비교만 진행해줘.

1. 앞 단계에서 선택한 R1/R2에 Trim Galore! paired-end 모드를 사용해.
2. 별도 근거가 없다면 --quality 20과 --length 20을 사용하고, 다른 값을 쓴다면 raw FastQC의 어떤 관찰 때문에 바꿨는지 먼저 설명해.
3. 출력은 trimmed 디렉터리에 저장하고 원본 FASTQ는 수정하지 마.
4. 생성된 *_val_1.fq.gz와 *_val_2.fq.gz의 gzip 무결성과 paired 파일 존재 여부를 확인해.
5. qc/trimmed에 FastQC를 다시 실행해.
6. raw와 trimmed의 total sequences, 길이 분포, per base quality, adapter content를 같은 표에서 비교해. 제거된 read 수와 비율도 계산해.
7. adapter와 낮은 품질 말단이 줄었는지, read 손실이 과도하지 않은지 실제 수치로 판정해.
8. 사용한 명령, 도구 버전, 결과 경로, 다음 단계에 사용할 정확한 FASTQ 경로를 보고하고 멈춰.

아직 STAR 정렬이나 정량은 실행하지 마.

전후 비교에서는 다음 네 질문에 답할 수 있어야 합니다.

3′ 말단의 낮은 품질이 개선됐는가
adapter 신호가 줄었는가
짧아진 read가 지나치게 많지 않은가
제거된 read의 수와 비율이 예상보다 크지 않은가

트리밍 결과가 나쁘면 같은 옵션으로 다음 단계에 밀어 넣지 않습니다. 원본 데이터의 품질 문제인지, adapter 자동 인식이 맞지 않았는지부터 확인합니다.

트리밍 전후의 숫자는 같은 단위를 놓고 비교해야 합니다. 예를 들어 raw R1이 2,000만 read이고 trimmed R1이 1,940만 read라면 60만 read, 즉 3%가 pair 보존 조건을 만족하지 못해 제거된 셈입니다. 이 숫자는 계산 방법을 보여 주는 가상 예시이며 실제 허용 범위는 실험과 원본 품질에 따라 달라집니다.

paired mode에서는 한쪽 read가 최소 길이보다 짧아질 때 그 fragment의 두 read를 함께 제외할 수 있습니다. 그래서 R1만 보고 손실을 해석하지 않고 Trim Galore! 리포트의 pair 처리 결과와 R1·R2 파일의 레코드 수를 함께 봅니다.

--quality 20은 모든 염기를 Q20 이상으로 만든다는 뜻이 아닙니다. Cutadapt가 말단에서 quality 기준을 적용해 잘라내는 규칙입니다. read 내부에 낮은 quality 염기가 하나 있다는 이유만으로 그 위치 앞뒤를 전부 버리지는 않습니다. --length 20은 트리밍 뒤 20 nt보다 짧은 read를 제외하는 기준입니다.

전후 그래프가 개선됐더라도 read 대부분을 잃었다면 좋은 결과가 아닙니다. 반대로 원본 adapter 신호가 거의 없고 품질이 충분했다면 트리밍 변화가 작아도 실패가 아닙니다. 목표는 최대한 많이 자르는 것이 아니라, 정렬에 쓸 수 있는 정보를 보존하면서 명확한 문제만 제거하는 것입니다.

3-3. STAR로 게놈과 전사체에 정렬

STAR는 같은 입력 read에서 용도가 다른 두 BAM을 만듭니다. 게놈 좌표 BAM은 featureCounts에, transcriptome 좌표 BAM은 RSEM에 들어갑니다.

3단계 프롬프트: STAR 정렬과 로그 판독

트리밍 후 FASTQ를 STAR로 정렬하는 단계만 진행해줘.

1. 앞에서 검증한 STAR index와 trimmed R1/R2의 절대 경로를 다시 확인해.
2. 사용 가능한 CPU와 메모리를 확인하고 다른 작업을 방해하지 않는 runThreadN 값을 정해. 선택 이유를 보고해.
3. align 디렉터리에 sample 이름을 prefix로 사용해 STAR를 실행해.
4. paired gzipped FASTQ이므로 --readFilesCommand zcat을 사용해.
5. --twopassMode Basic, --quantMode TranscriptomeSAM GeneCounts, --outSAMtype BAM Unsorted를 포함해. 종양 fusion 분석용 --chim* 옵션은 이번 기본 실습에 임의로 추가하지 마.
6. 실행한 전체 명령, STAR 버전, index의 메타데이터를 텍스트 파일로 남겨.
7. sample.Aligned.out.bam, sample.Aligned.toTranscriptome.out.bam, sample.Log.final.out이 생성됐는지 확인해. samtools quickcheck는 게놈 BAM에 실행해.
8. Log.final.out에서 input reads, uniquely mapped reads %, multi-mapping %, unmapped 사유와 splice 수치를 표로 정리해.
9. 정렬률이 낮거나 두 BAM 중 하나가 없으면 원인을 조사해. 레퍼런스 불일치, read 길이, 오염 가능성을 구분하고 다음 단계로 넘어가지 마.
10. 결과가 다음 단계에 사용할 수 있다고 판단한 근거와 정확한 BAM 경로를 한국어로 보고하고 멈춰.

아직 RSEM과 featureCounts는 실행하지 마.

Log.final.out의 수치는 데이터 종류와 실험 설계에 따라 달라집니다. uniquely mapped reads가 80% 이상이면 대체로 양호한 출발점이고 RNA-seq에서 일부 multi-mapping은 흔하지만, 이 숫자를 보편적인 통과선으로 사용하지는 않습니다. 값이 예상보다 나쁘면 게놈과 annotation 버전, read 길이, 트리밍 강도와 오염을 함께 점검합니다.

STAR 로그의 비율은 먼저 Number of input reads를 분모로 읽습니다. paired-end에서 input read 하나는 R1 파일의 한 줄이 아니라 R1과 R2로 이루어진 한 fragment pair를 뜻합니다. FASTQ 레코드 수를 셀 때와 STAR 로그를 비교할 때 이 단위가 같은지 확인해야 합니다.

STAR 로그 항목	의미	값이 예상보다 나쁠 때 먼저 볼 것
`Uniquely mapped reads %`	한 위치에 가장 잘 맞은 fragment 비율	종, genome build, read 품질과 오염
`% of reads mapped to multiple loci`	여러 위치에 비슷하게 맞은 비율	반복서열, 유사 유전자와 짧은 read
`% of reads unmapped: too many mismatches`	허용하기 어려울 만큼 불일치한 비율	잘못된 종, 낮은 품질과 오염
`% of reads unmapped: too short`	정렬 근거가 충분하지 않은 비율	과도한 트리밍과 짧은 원본 read
`Number of splices`	exon 사이를 건너 정렬된 junction 수	RNA-seq library와 annotation 맥락

예를 들어 input fragment가 2,000만 개이고 unique 85%, multiple loci 8%, unmapped 7%라면 대략 1,700만 개가 한 위치, 160만 개가 여러 위치에 정렬되고 140만 개가 정렬되지 않은 장면입니다. 가상 예시의 세 비율을 합치면 전체에 가까워지지만, 실제 로그에는 세분된 범주와 반올림이 있으므로 표시된 숫자를 직접 확인해야 합니다.

두 BAM은 같은 정렬 결과를 서로 다른 좌표 체계로 표현합니다.

파일	정렬 좌표	이 실습의 사용처
`sample.Aligned.out.bam`	chromosome의 게놈 좌표	GTF exon과 겹쳐 featureCounts로 집계
`sample.Aligned.toTranscriptome.out.bam`	transcript 서열 좌표	RSEM이 transcript별 가능성을 계산

transcriptome BAM을 featureCounts에 넣거나 genome BAM을 현재 RSEM 명령에 넣으면 파일 형식은 BAM이라도 좌표의 의미가 맞지 않습니다. 확장자만 보고 연결하지 않고 STAR 명령의 산출물 이름을 확인합니다.

samtools quickcheck가 성공했다는 것은 BAM header와 끝부분이 심하게 손상되지 않았다는 뜻입니다. 정렬률이 좋다거나 올바른 sample이라는 뜻은 아닙니다. 로그의 수치, 파일 크기와 reference 정보를 별도로 봐야 합니다.

3-4. RSEM으로 TPM 계산

RSEM은 STAR가 만든 transcriptome 좌표 BAM과 RSEM 레퍼런스를 받아 유전자와 아이소폼 수준의 발현량을 계산합니다.

4단계 프롬프트: RSEM 유전자 발현량

검증한 transcriptome BAM으로 RSEM 정량 단계만 진행해줘.

1. sample.Aligned.toTranscriptome.out.bam과 RSEM reference prefix의 관련 파일이 존재하는지 다시 확인해.
2. STAR index, RSEM reference와 GTF가 같은 게놈 빌드와 annotation 릴리스를 사용했는지 확인 가능한 근거를 기록해.
3. quant 디렉터리에 결과가 생기도록 rsem-calculate-expression을 --alignments --paired-end --estimate-rspd --append-names 옵션으로 실행해. thread 수는 현재 자원에 맞게 정해.
4. 실행한 전체 명령과 RSEM 버전을 텍스트 파일로 남겨.
5. sample.genes.results와 sample.isoforms.results의 존재, 행 수, 열 이름, 결측값과 음수 값 여부를 검사해.
6. genes.results에서 gene_id, expected_count, TPM 앞부분을 보여주고 TPM 합계를 계산해. --append-names로 붙은 이름을 gene ID와 혼동하지 않게 설명해.
7. TPM이 높은 상위 10개 유전자를 출력하되, 이것만으로 샘플 품질이나 생물학적 결론을 확정하지 마.
8. 결과 경로, 핵심 검증값, 경고를 한국어로 보고하고 멈춰.

아직 featureCounts는 실행하지 마.

sample.genes.results는 유전자 수준의 expected_count, TPM과 FPKM을 포함합니다. 기본 발현 탐색은 이 파일에서 시작하고, 특정 전사체의 차이가 중요할 때 sample.isoforms.results를 봅니다.

결과표의 행 하나는 유전자 하나입니다. RSEM 버전과 레퍼런스 설정에 따라 세부 열은 달라질 수 있지만 주로 다음 열을 읽습니다.

열	먼저 읽을 뜻
`gene_id`	annotation에서 가져온 안정적인 유전자 식별자
`transcript_id(s)`	그 유전자에 연결된 transcript 식별자
`length`	annotation에 있는 transcript 길이를 반영한 값
`effective_length`	fragment가 실제로 놓일 수 있는 위치 수를 반영한 유효 길이
`expected_count`	애매한 정렬을 확률적으로 나눈 기대 fragment 수
`TPM`	길이와 sample 내 총량을 보정한 상대 발현 비중

expected_count는 1250.37처럼 소수가 될 수 있습니다. 하나의 fragment가 여러 transcript에 맞을 때 RSEM이 가능성을 나누기 때문입니다. 파일이 잘못된 것이 아닙니다. 반면 뒤의 featureCounts raw count는 직접 배정 규칙으로 센 정수입니다.

TPM은 한 sample 안에서 합이 대체로 100만이 되도록 만든 비율입니다. 어떤 유전자의 TPM이 200이면 그 유전자가 sample의 전체 transcript 신호 가운데 약 200/1,000,000을 차지한다는 뜻입니다. 유전자 길이 효과를 보정했기 때문에 같은 sample 안의 유전자 발현 비중을 살펴보는 데 유용합니다.

하지만 서로 다른 sample의 TPM을 DESeq2 입력으로 사용하지 않습니다. 조건 간 차등 발현 검정은 count가 가진 평균과 분산 관계를 모델링하므로 정수 raw count와 sample별 metadata가 필요합니다. TPM 상위 10개도 그 sample에서 많이 관찰된 유전자일 뿐, 실험 조건의 특징이라고 말하려면 다른 sample과 비교해야 합니다.

--append-names를 사용하면 ENSG..._GENE처럼 ID와 gene symbol이 함께 보일 수 있습니다. symbol은 사람이 읽기 쉽지만 바뀌거나 중복될 수 있으므로, 여러 sample 결과를 합칠 때는 같은 annotation 릴리스의 gene ID를 기준 키로 사용하는 편이 안전합니다.

3-5. featureCounts로 raw count 만들기

featureCounts는 게놈 좌표 BAM에서 GTF의 exon과 겹치는 fragment를 gene ID별로 셉니다. RSEM과 입력 BAM이 다르다는 점이 핵심입니다.

Aligned.toTranscriptome.out.bam → RSEM → TPM
Aligned.out.bam + GTF → featureCounts → raw count

5단계 프롬프트: featureCounts raw count

검증한 genome BAM과 GTF로 featureCounts 단계만 진행해줘.

1. sample.Aligned.out.bam과 GTF의 절대 경로, BAM의 paired-end 여부를 확인해.
2. 설치된 featureCounts 버전의 도움말에서 paired-end fragment 집계 옵션을 확인해. 지원된다면 -p --countReadPairs를 사용하고, 버전에 따른 차이가 있으면 적용한 옵션을 설명해.
3. -t exon -g gene_id --extraAttributes gene_name을 사용하고 counts/sample.featureCounts.txt에 저장해. thread 수는 현재 자원에 맞게 정해.
4. 실행한 전체 명령, Subread/featureCounts 버전과 GTF 식별 정보를 텍스트 파일로 남겨.
5. 결과 표의 주석 열과 sample count 열을 구분하고, gene ID 수, count가 0인 gene 수, 총 assigned fragment 수를 계산해.
6. sample.featureCounts.txt.summary에서 Assigned와 각 미배정 사유의 수와 비율을 표로 정리해.
7. 배정률이 낮다면 GTF와 BAM의 contig 이름, annotation 릴리스, strandedness 설정을 우선 점검해. 임의로 옵션을 바꿔 재실행하지 마.
8. 출력 파일 경로와 이후 여러 sample을 count matrix로 합칠 때 사용할 열을 한국어로 보고하고 멈춰.

featureCounts 표의 raw count는 여러 sample의 차등 발현 분석에 사용합니다. TPM과 raw count는 모두 발현량을 나타내지만, 서로 바꿔 쓸 수 있는 같은 값은 아닙니다. 용도 차이는 read를 발현량 행렬로 바꾸기에서 설명합니다.

결과 파일을 열면 왼쪽에는 유전자 주석, 오른쪽에는 sample count가 있습니다. 아래는 구조를 보여 주기 위한 가상 예시입니다.

Geneid	Chr	Start	End	Strand	Length	`sample.Aligned.out.bam`
`ENSG000001`	`chr1`	1001	1800	`+`	800	1,250
`ENSG000002`	`chr1`	5001	6200	`-`	1,200	47

첫 행의 1,250은 sample의 게놈 BAM에서 ENSG000001의 exon에 featureCounts 규칙대로 배정된 fragment 수입니다. Length = 800은 count를 길이로 나눈 값이 아니라 GTF feature에서 계산한 주석 길이입니다. 마지막 count 열만 여러 sample의 raw count matrix에 들어가고, Chr, Start, End, Strand, Length는 gene별 주석으로 한 번만 유지할 수 있습니다.

summary 파일은 count에 들어오지 못한 fragment도 이유별로 보여 줍니다.

summary 범주	뜻
`Assigned`	지정한 feature와 gene ID에 최종 배정됨
`Unassigned_NoFeatures`	정렬 위치가 지정한 exon과 겹치지 않음
`Unassigned_Ambiguity`	둘 이상의 gene에 걸려 하나로 정하지 못함
`Unassigned_MultiMapping`	여러 게놈 위치에 정렬되어 현재 규칙에서 제외됨
`Unassigned_Unmapped`	레퍼런스에 정렬되지 않음

Unassigned_NoFeatures가 크면 GTF가 BAM과 같은 genome build인지, chromosome 이름이 chr1과 1처럼 다르지 않은지 먼저 확인합니다. library가 strand-specific인데 기본 unstranded 설정을 썼어도 배정률이 크게 달라질 수 있습니다. strandedness는 파일명만으로 추측하지 않고 library preparation 기록이나 RSeQC 같은 별도 추론 결과로 정해야 합니다.

featureCounts와 RSEM의 count가 정확히 같을 필요는 없습니다. RSEM은 transcriptome에서 애매한 read를 확률적으로 배분하고, featureCounts는 지정한 겹침 규칙으로 fragment를 배정합니다. 두 도구가 서로 다른 질문과 규칙을 사용하므로 차이가 생깁니다. 중요한 것은 각 값이 어떤 입력 BAM, annotation과 옵션에서 나왔는지 기록하는 것입니다.

3-6. 전체 산출물 감사하기

마지막 단계는 새 분석을 추가하는 단계가 아닙니다. 지금까지의 입력, 명령, 버전과 결과가 서로 연결되는지 검사합니다.

6단계 프롬프트: 산출물과 재현성 최종 점검

지금까지 만든 RNA-seq 실습 산출물을 읽기 전용으로 최종 점검해줘.

1. raw FASTQ부터 FastQC, trimmed FASTQ, STAR BAM과 로그, RSEM 결과, featureCounts 결과까지 파일 경로와 크기를 단계 순서대로 표로 정리해.
2. 각 단계의 입력 파일이 바로 앞 단계의 검증된 출력인지 확인해. 특히 RSEM은 transcriptome BAM, featureCounts는 genome BAM을 사용했는지 검사해.
3. 실행 기록에서 도구 버전, 전체 명령, reference genome과 annotation 식별 정보가 남아 있는지 확인해.
4. raw 대비 trimmed read 수, STAR input reads, RSEM 입력과 featureCounts 집계 단위 사이에 설명되지 않는 불일치가 없는지 점검해.
5. FastQC 전후 핵심 변화, STAR unique/multi/unmapped 비율, featureCounts Assigned 비율을 한 표로 요약해.
6. 재현에 필요한데 기록되지 않은 정보와 추가 확인이 필요한 경고를 별도 목록으로 작성해.
7. 파일을 수정하거나 분석을 다시 실행하지 말고, 완료된 것과 미완료된 것을 한국어로 보고하고 멈춰.

최종 감사표는 파일 존재 여부만 나열해서는 안 됩니다. 다음 연결 관계가 한 줄씩 추적되어야 합니다.

확인할 연결	맞아야 하는 관계
raw와 trimmed read 수	손실량이 Trim Galore! 리포트로 설명됨
trimmed read와 STAR input	STAR가 실제로 트리밍된 pair를 읽음
STAR와 RSEM	transcriptome BAM이 RSEM 입력으로 사용됨
STAR와 featureCounts	genome BAM과 같은 빌드의 GTF가 사용됨
출력표와 sample 이름	모든 결과가 처음 선택한 같은 sample을 가리킴

예를 들어 FastQC는 tumor01을 검사했는데 STAR prefix가 normal01이라면 파일은 모두 존재해도 추적이 끊깁니다. 파일명, 절대 경로와 명령 기록을 함께 보는 이유입니다.

이 실습의 해석 범위

이 실습이 확인하는 것은 한 paired-end RNA-seq sample의 read를 정렬하고 유전자별 발현량으로 정리하는 기술적 과정입니다. 결과가 만들어졌다고 다음 내용을 자동으로 알 수 있는 것은 아닙니다.

sample이 원래 붙인 생물학적 조건과 일치하는가
조직 안에 어떤 세포 종류가 얼마나 섞였는가
다른 조건보다 특정 유전자가 증가하거나 감소했는가
관찰한 차이가 통계적으로 재현되는가
발현 변화가 질병의 원인인가 결과인가

이 질문에는 여러 biological replicate, 정확한 metadata, sample 간 QC와 차등 발현 모델이 더 필요합니다. 이 페이지의 raw count는 그 다음 분석의 입력이지 최종 결론이 아닙니다.

4. 완료 기준

raw와 trimmed FastQC HTML을 실제로 비교했다.
FASTA, GTF, STAR index와 RSEM reference의 버전 정합성을 기록했다.
STAR Log.final.out의 unique, multi-mapping과 unmapped 사유를 확인했다.
게놈 BAM과 transcriptome BAM의 용도를 구분했다.
RSEM genes.results에서 TPM과 expected count를 확인했다.
featureCounts 표와 summary에서 gene별 raw count와 배정률을 확인했다.
도구 버전, 실행 명령, 입력과 출력 경로가 재현 가능하게 남아 있다.

한 sample에서 여기까지 통과했다면 같은 절차를 다른 sample에 적용하고, sample별 featureCounts 열을 gene ID 기준으로 합쳐 count matrix를 만들 수 있습니다. 여러 sample을 자동화할 때도 먼저 한 쌍을 끝까지 통과시킨 뒤 확장하는 편이 오류를 찾기 쉽습니다.

모두 실행한 뒤 답해 볼 질문

R1과 R2의 같은 번째 레코드가 같은 fragment여야 하는 이유는 무엇인가?
FastQC에서 duplication 경고만으로 RNA-seq sample을 실패라고 할 수 없는 이유는 무엇인가?
트리밍 뒤 adapter 신호는 줄었지만 read의 절반을 잃었다면 무엇을 다시 확인해야 하는가?
STAR가 만든 genome BAM과 transcriptome BAM은 각각 어느 도구에 들어가는가?
RSEM expected_count가 소수이고 featureCounts raw count가 정수일 수 있는 이유는 무엇인가?
TPM을 DESeq2 입력으로 사용하지 않는 이유는 무엇인가?
featureCounts의 Unassigned_NoFeatures가 클 때 레퍼런스와 GTF에서 무엇을 확인해야 하는가?
같은 결과를 재현하려면 FASTQ와 명령 외에 어떤 버전 정보를 기록해야 하는가?