Bash로 시퀀싱 데이터 다루기

이 토픽을 마치면

터미널에서 FASTA/FASTQ 파일을 열어보고, 서열 개수를 세고, 특정 서열을 검색할 수 있습니다. 마우스 없이 키보드만으로 데이터를 다루는 첫 경험입니다.

터미널이란?

연구실에서 피펫으로 시약을 정밀하게 다루듯, 터미널은 컴퓨터의 데이터를 정밀하게 다루는 도구입니다. 클릭 대신 명령어를 입력합니다.

Google Colab에서 셀 앞에 !를 붙이면 Bash 명령어를 실행할 수 있습니다.

파일 만들고 살펴보기

먼저 실습용 FASTA 파일을 만들어봅니다.

bash

# 실습용 FASTA 파일 생성
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

# 파일 내용 확인
cat genes.fasta

bash

# 파일의 줄 수 세기
wc -l genes.fasta

wc는 Word Count의 약자입니다. -l 옵션은 줄(line) 수만 셉니다.

grep: 서열 검색하기

grep은 파일에서 특정 패턴을 찾아주는 명령어입니다. 논문에서 키워드를 검색하는 것과 같습니다.

bash

# FASTA 헤더(>) 줄만 출력
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta

bash

# 서열 수 세기 (헤더 줄 = 서열 수)
grep -c ">" genes.fasta

-c 옵션은 매칭된 줄의 개수(count)만 출력합니다.

bash

# 특정 유전자 찾기
grep "BRCA1" genes.fasta

FASTQ 형식 이해하기

실제 시퀀싱 장비에서 나오는 데이터는 FASTQ 형식입니다. FASTA와 달리 품질 점수(quality score)가 포함되어 있습니다.

bash

# FASTQ 파일 생성 (4줄이 한 세트)
echo "@read_001
ATGCGATCGATCGATCGATCG
+
IIIIIIIIIIIIIIIIIIIII
@read_002
GCTAGCTAGCTAGCTAGCTAG
+
IIIIIIIIIIIIIIIIIIIII
@read_003
NNNNNATCGATCGATCGATCG
+
!!!!!IIIIIIIIIIIIIIIII" > reads.fastq

cat reads.fastq

FASTQ의 4줄 구조:

@로 시작하는 리드 ID
염기서열
+ (구분선)
품질 점수 (높을수록 좋음, I=최고, !=최저)

awk: 데이터 가공하기

awk는 텍스트 데이터를 열(column) 단위로 처리하는 도구입니다. 엑셀에서 특정 열만 선택하는 것과 비슷합니다.

bash

# 탭 구분 데이터 만들기
echo "BRCA1	chr17	43044295	43170245
TP53	chr17	7661779	7687538
EGFR	chr7	55019017	55211628" > gene_locations.tsv

# 유전자 이름(1열)만 출력
awk '{print $1}' gene_locations.tsv

bash

# 유전자 길이 계산 (4열 - 3열)
awk '{print $1, $4 - $3, "bp"}' gene_locations.tsv

파이프(|): 명령어 연결하기

파이프(|)는 한 명령어의 출력을 다음 명령어의 입력으로 연결합니다. 실험에서 시약 A의 결과물을 시약 B에 넣는 것과 같습니다.

bash

# FASTA에서 헤더만 뽑아서 개수 세기
echo ">BRCA1_human
ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA
>EGFR_human
ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG" > genes.fasta

grep ">" genes.fasta | wc -l

직접 해보기 (Faded Example)

아래 빈칸을 채워 FASTQ 파일에서 리드 수를 세는 명령어를 완성하세요.

FASTQ에서 리드 ID는 @로 시작합니다.

빈칸 채우기bash

# @로 시작하는 줄의 개수를 세려면?
grep  "" test.fastq

흔한 에러 & 해결법

Q: command not found가 뜹니다

Colab에서는 ! 접두어를 붙여야 합니다: !grep ">" genes.fasta

Q: 한글이 깨져서 나옵니다

파일 인코딩 문제입니다. file genes.fasta로 인코딩을 확인하고, UTF-8이 아니면 iconv -f EUC-KR -t UTF-8 input.txt > output.txt로 변환합니다.

Q: Permission denied가 뜹니다

파일에 실행 권한이 없을 때 발생합니다. chmod +x script.sh로 권한을 부여하세요.

다음 글에서는 Python으로 이 데이터를 더 정밀하게 분석하는 법을 배웁니다.