HTSeq raw count 읽기

HTSeq는 Python 생물정보학 소프트웨어 패키지이고, htseq-count 는 그 안에서 read를 유전자별로 세는 명령줄 프로그램입니다. HTSeq raw count는 htseq-count가 만든 정규화 전 데이터 표입니다.

이 표에는 각 RNA-seq 샘플에서 유전자별로 몇 개의 서열 조각(read 또는 read pair)이 배정됐는지가 적힙니다. read를 기준 유전체에 정렬한 다음, 그 위치를 유전자 위치 정보인 주석(annotation)과 대조해 만듭니다.

여러 샘플을 합친 count matrix에서는 보통 행이 유전자이고 열이 샘플입니다. 예를 들어 sample_A 열의 한 부분이 다음과 같다고 해봅시다.

gene_id             sample_A
ENSG00000000003     9896

sample_A: 실험자가 붙인 샘플 이름
ENSG00000000003: 유전자를 구분하는 Ensembl ID
9896: 이 유전자에서 왔다고 판단된 read, 즉 시퀀서가 읽어 낸 짧은 서열 조각의 수

paired-end 데이터라면 양쪽 read를 한 쌍으로 묶어 세므로, 9896은 read pair 9,896개가 이 유전자에 배정됐다는 뜻입니다. RNA 분자가 정확히 9,896개 있었다는 뜻은 아닙니다. 아직 샘플별 시퀀싱 깊이나 유전자 길이를 보정하지 않은 숫자라서 raw count라고 부릅니다.

HTSeq와 GTF가 맡는 역할

정렬이 끝난 read에는 chr1의 11870번 위치에 붙었다는 정보만 있습니다. 그 위치가 어느 유전자인지는 정렬 결과만 보고 알 수 없습니다. 별도의 유전자 위치표와 대조해야 합니다.

HTSeq는 대량 시퀀싱(high-throughput sequencing) 데이터를 처리하는 Python 기반 생물정보학 소프트웨어 패키지입니다. 그 안의 htseq-count 명령이 정렬된 read를 유전자별로 세는 일을 합니다.
GTF(Gene Transfer Format)는 게놈 위에 유전자, 전사체, 엑손이 어디 있는지 적어 둔 주석 파일입니다. 여기서 주석(annotation)은 서열을 바꾸는 메모가 아니라, 게놈 좌표에 생물학적 의미를 붙인 정보입니다.

예를 들어 정렬 결과에는 read가 chr1:1000-1075에 붙었다고 기록되어 있고, GTF에는 다음과 같은 줄이 있다고 해봅시다.

chr1  source  exon  900  1100  .  +  .  gene_id "GENE1";

이 줄은 chr1의 900번부터 1100번까지가 GENE1에 속한 엑손이라는 뜻입니다. read의 위치가 이 구간과 겹치므로 htseq-count는 GENE1의 count를 올립니다. 즉, 하는 일은 두 좌표표를 대조하는 것입니다.

단계	파일 또는 프로그램	담긴 정보와 역할
입력	SAM/BAM/CRAM	각 read가 게놈 어디에 정렬됐는가
입력	GTF	각 유전자와 엑손이 게놈 어디에 있는가
처리	`htseq-count`	두 파일의 좌표를 겹쳐 보고 read를 유전자에 배정
출력	raw count	유전자별로 배정된 read 또는 fragment의 수

RNA-seq에서는 보통 GTF에서 exon이라고 표시된 줄을 읽습니다. 여러 엑손이 같은 gene_id를 가지면 한 유전자에 속한 것으로 묶고, 그 엑손들에 배정된 관측값을 모두 더합니다. GTF의 9개 컬럼과 gene_id를 자세히 읽는 법은 GTF: 유전자 주석 파일에서 이어집니다.

출력은 유전자 × 샘플 행렬

gene_id              sample_1  sample_2
ENSG00000000003           9896      6279
ENSG00000000005             32        49
ENSG00000000419           5129      2406

일반적인 bulk RNA-seq 결과를 합치면 행은 유전자, 열은 샘플, 값은 count인 행렬이 됩니다. 분석 도구에 따라 입력 방향이 다를 수 있으므로, 행과 열을 추측하지 말고 파일의 첫 행과 분석 도구의 요구 형식을 확인해야 합니다.

숫자 하나가 올라가는 조건

기본 union 모드에서는 read가 겹치는 유전자 집합을 확인합니다.

read가 겹치는 유전자	처리
정확히 하나	그 유전자의 count를 1 올림
없음	`__no_feature`에 기록
둘 이상	기본 설정에서는 `__ambiguous`에 기록

paired-end 데이터에서는 두 mate를 별개의 read 두 개가 아니라 한 read pair, 즉 한 cDNA fragment의 증거로 세므로 count가 한 번 올라갑니다. --nonunique=fraction처럼 일부 옵션은 소수 count를 만들 수 있지만, 기본적인 gene-level count는 정수입니다.

`__`로 시작하는 행은 유전자가 아니다

HTSeq는 유전자에 들어가지 않은 read도 이유별로 집계합니다. 이 특수 카운터는 모두 __로 시작하므로 유전자 행과 구분할 수 있습니다.

특수 카운터	의미	많이 나올 때 확인할 것
`__no_feature`	어떤 annotation feature와도 겹치지 않음	genome build, 염색체 이름, strandedness, GTF
`__ambiguous`	둘 이상의 feature와 겹쳐 하나를 고를 수 없음	겹치는 유전자와 annotation, overlap mode
`__too_low_aQual`	최소 alignment quality보다 낮아 제외됨	aligner의 MAPQ와 HTSeq의 `-a` 설정
`__not_aligned`	정렬되지 않은 레코드	upstream alignment 결과
`__alignment_not_unique`	게놈의 여러 위치에 정렬됨	반복서열, aligner의 multi-mapping 설정

HTSeq 결과 파일의 끝부분에서는 다음과 같은 행을 볼 수 있습니다.

__no_feature
__ambiguous
__too_low_aQual
__not_aligned
__alignment_not_unique

이 행들은 유전자가 아니므로 DESeq2 입력에서는 제외합니다. 하지만 파일을 읽자마자 버리기 전에 샘플별 값을 확인하는 편이 좋습니다. 특정 샘플에서 __no_feature나 __alignment_not_unique의 비율이 유난히 크다면 정렬, annotation, 샘플 품질 차이를 먼저 점검해야 합니다.

special = genes_by_samples.loc[
    genes_by_samples.index.str.startswith("__")
]
gene_counts = genes_by_samples.loc[
    ~genes_by_samples.index.str.startswith("__")
]

print(special)

특수 카운터가 크다는 사실만으로 샘플을 실패로 판정하지는 않습니다. protocol, annotation 범위, multi-mapping 처리 방식에 따라 기대값이 달라지므로 다른 샘플과의 상대적 차이와 실행 옵션을 함께 봅니다.

raw count를 그대로 비교하면 안 되는 이유

같은 조직(tissue)도 더 깊게 시퀀싱하면 대부분의 유전자 count가 함께 커집니다. 유전자가 길수록 read가 걸릴 기회도 많습니다. 따라서 서로 다른 유전자의 숫자나 서로 다른 샘플의 총량을 raw count만 보고 직접 비교하면 안 됩니다.

값	보정 상태	주된 용도
HTSeq raw count	보정 전	DESeq2·edgeR 같은 count 모델의 입력
DESeq2 normalized count	샘플별 깊이 보정	시각화와 샘플 내 패턴 확인
TPM	유전자 길이와 샘플 깊이 보정	한 샘플의 상대적 발현 구성 확인

DESeq2는 raw count를 받은 뒤 샘플별 시퀀싱 깊이와 유전자별 변동을 모델 안에서 추정합니다. 같은 유전자를 조건 사이에서 비교할 때 유전자 길이는 모든 샘플에서 같으므로, DEG 검정 전에 TPM으로 바꿀 필요가 없습니다. 정규화 값의 용도 차이는 raw count를 비교 가능한 값으로 바꾸기에서 이어집니다.

받은 count matrix에서 확인할 것

값이 음수가 아닌 정수인가
행 ID가 gene_id, gene symbol, transcript ID 중 무엇인가
열 하나가 샘플 하나인지, 샘플 이름이 metadata와 맞는가
__ 특수 카운터가 있는가, 샘플 사이에서 유난히 큰 값이 있는가
사용한 genome build와 GTF annotation 버전은 무엇인가
single-end인지 paired-end인지, strandedness와 overlap mode는 무엇인가
정규화된 값이 아니라 실제 raw count가 맞는가

공개 저장소의 processed data라는 분류만으로는 이 내용을 알 수 없습니다. 파일의 행과 값을 직접 열고, 데이터 페이지의 처리 방법과 논문 Methods를 함께 확인해야 합니다.

사례: GSE251845 파일 읽기

GSE251845의 GSE251845_htseq_raw_counts.csv.gz도 유전자가 행이고 조직 샘플이 열인 HTSeq count matrix입니다. 예를 들어 24C는 24번 환자의 종양 조직에 붙인 샘플 이름입니다. 한 유전자 행과 24C 열이 만나는 숫자는 그 조직의 정렬 결과에서 해당 유전자에 배정된 fragment 수를 뜻합니다.

파일 끝의 __no_feature, __ambiguous, __alignment_not_unique 같은 행은 유전자가 아닙니다. 따라서 유전자별 차등 발현 분석에서는 이 행들을 별도로 확인한 뒤 count matrix에서 제외합니다.