# HTSeq raw count 읽기

> htseq-count가 BAM과 GTF를 이용해 read를 유전자에 배정하는 방법, raw count 행렬과 특수 카운터를 읽는 법.

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

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

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

```text
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**라고 부릅니다.

:::note[raw count의 raw는 FASTQ라는 뜻이 아니다]
FASTQ는 시퀀서가 만든 원시 read이고, HTSeq raw count는 FASTQ를 정렬하고 유전자에 배정한 뒤 만든 가공 데이터입니다. 여기서 `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에는 다음과 같은 줄이 있다고 해봅시다.

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

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

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

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

## 출력은 유전자 × 샘플 행렬

```text
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 결과 파일의 끝부분에서는 다음과 같은 행을 볼 수 있습니다.

```text
__no_feature
__ambiguous
__too_low_aQual
__not_aligned
__alignment_not_unique
```

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

```python
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를 비교 가능한 값으로 바꾸기](/lessons/rna-seq-quantification/)에서 이어집니다.

## 받은 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에서 제외합니다.

### 공식 문서

- [HTSeq 공식 문서: htseq-count로 feature별 read 세기](https://htseq.readthedocs.io/en/latest/htseqcount.html)
- [HTSeq 공식 튜토리얼: read를 feature에 배정하는 내부 과정](https://htseq.readthedocs.io/en/latest/tutorials/exon_example.html)