RNA-seq 전처리 #3 Read counting

RNA-seq 전처리 #3 Read counting #

#2026-01-25

#1 read counting

alignment output file인 accepted_hits.bam 파일을 사용해서 gene-level count 테이블을 생성하기.

HTSeq-count 도구를 이용했고 run_htseq_count.sh 스크립트로 작업.

#!/bin/bash
set -euo pipefail

# ===== 설정 =====
TOPHAT_DIR="TopHat"
GTF="./hg38.refGene.gtf"
OUT_DIR="Count"

mkdir -p "${OUT_DIR}"

# 대상 샘플 (요청한 5개)
SAMPLES=(
  "150-2"
  "150-3"
  "33-1"
  "33-2"
  "33-3"
)

echo "[INFO] HTSeq-count started."

for s in "${SAMPLES[@]}"; do
  bam="${TOPHAT_DIR}/tophat_out_${s}/accepted_hits.bam"
  out="${OUT_DIR}/count_${s}.txt"

  if [[ ! -f "${bam}" ]]; then
    echo "[WARN] BAM not found: ${bam} (skip)"
    continue
  fi

  echo "[COUNT] ${bam} -> ${out}"

  # header 포함 SAM 스트림 생성 후 HTSeq-count로 카운팅
  samtools view -H "${bam}" > header.sam
  samtools view "${bam}" | cat header.sam - | python -m HTSeq.scripts.count - "${GTF}" > "${out}"

done

rm -f header.sam
echo "[DONE] HTSeq-count completed."

chmod +x run_htseq_count.sh
./run_htseq_count.sh > htseq_count.log 2>&1

실행하면 샘플에 대해 Count/count_150-2.txt와 같은 gene-level count 파일이 생성된다.

#

#2 HTSeq-count 오류 트러블슈팅

accepted_hits.bam 파일을 HTSeq-count에 입력했는데, 아래와 같은 경고/에러가 발생하며 count 생성이 중단됨

$ samtools view -H ../2306_tophat/data/Bowtie2Index/tophat_out_150-1/accepted_hits.bam > header.sam
samtools view ../2306_tophat/data/Bowtie2Index/tophat_out_150-1/accepted_hits.bam | \
cat header.sam - | \
python -m HTSeq.scripts.count - ./hg38.refGene.gtf > Count/count_150-1.txt

Warning: Read A01056:397:H35NJDSX7:3:1239:2718:33223 claims to have an aligned mate which could not be found in an adjacent line.
Error occured when processing input (record #981 in file <stdin>):
  'pair_alignments' needs a sequence of paired-end alignments
  [Exception type: ValueError]

HTSeq-count는 paired-end 모드에서 read1/read2가 입력 스트림에서 연속된 줄로 들어오는 것을 전제로 mate를 매칭함. accepted_hits.bam은 일반적으로 좌표 기준 정렬(coordinate-sorted) 상태이기 때문에, mate 쌍이 파일 내에서 떨어져 존재할 수 있고(서로 다른 위치로 정렬되면 더 흔함), 그 결과 HTSeq가 “인접 라인에서 mate를 못 찾았다”고 판단하며 오류가 발생한다.

입력 BAM을 read name 기준으로 정렬(name-sorted) 한 뒤, HTSeq-count에 전달하도록 변경함 즉, mate끼리 항상 연속된 상태로 들어가도록 보장해서 pair_alignments 에러를 방지함.

# 1) name-sort BAM 생성
samtools sort -n -o accepted_hits.name_sorted.bam \
  ../2306_tophat/data/Bowtie2Index/tophat_out_150-1/accepted_hits.bam

# 2) name-sorted BAM을 header 포함 SAM 스트림으로 변환 후 HTSeq-count 수행
samtools view -H accepted_hits.name_sorted.bam > header.sam
samtools view accepted_hits.name_sorted.bam | \
cat header.sam - | \
python -m HTSeq.scripts.count - ./hg38.refGene.gtf > Count/count_150-1.txt

#

#3 정리

  • Alignment 결과(accepted_hits.bam)를 입력으로 gene-level read counting 수행
  • hg38.refGene.gtf annotation 기준으로 HTSeq-count 적용
  • header 포함을 위해 samtools view -H + cat 파이프라인 사용
  • BAM을 SAM으로 스트리밍 변환하면서 paired-end 처리 중 mate가 인접하지 않은(coordinate-sorted) BAM 입력으로 인해 pair_alignments 오류가 발생함
  • samtools sort -n을 이용해 name-sorted BAM으로 변환 후 HTSeq-count를 수행하여 mate adjacency를 보장하고 오류를 해결.