FASTQ 파일 내의 base count 하기

fasta 파일의 경우에는 fastasize.py 라는 걸 선구자 분들이 만들어 놔서 간단하게 카운트가 되었다. 그래서 매번 fastq 를 fasta 로 변환하고 count 를 했는데 찾다보니 직접 count 할 수 있다는 걸 알게 되었다.

https://www.biostars.org/p/78043/

위의 웹페이지는 일종의 토론을 한 내용이고.. 거의 마지막에 간단한 명령어가 있어서 이것을 애용한다.

grep "^[ACGTN]" test.fastq | tr -d "\n" | wc -m

awk '{s++}END{print s/4}' test.fastq

위는 base number 를 count 하고 밑은 read number 를 카운트 한다.

----------------------
2018년 3월 21일
위의 내용을 수정함.

위에서 base number 를 count 할 때 안 맞는다. 이상해서 찾아보니 quality 에도 ATGCN 이 존재할 수 있기 때문에 그것들까지 count 되는 듯 싶다. 그래서 찾아보니...

awk 'NR%4==2{c++; l+=length($0)}END{print "Number of reads: "c;print "Number of bases in reads: "l}' fastq.file

이라는 방법이 있었고, 다시 해서 맞춰보니 fastasize.py 와 일치하는 수자가 나온다.

FASTQ 포멧의 전환...

https://en.wikipedia.org/wiki/FASTQ_format

위의 위키피디아에는 fastq 포멧에 대한 설명이 되어 있다.
프로그래머 라면 포멧에 대한 설명도 매우 중요하겠지만 나는 그냥 이용자일 뿐이라서 그다지 중요하지 않다.

보다 중요한 것은 맨 밑에...
FASTQ 를 FASTA 포멧으로 전환하기..

zcat input_file.fastq.gz | awk 'NR%4==1{printf ">%s\n", substr($0,2)}NR%4==2{print}' > output_file.fa

정확하게는 gzip 으로 압축되어 있는 fastq 파일의 압축을 해제하면서 fasta 파일로 전환하는 것이다. 압축이 되어 있지 않다면 zcat 대신 cat 을 쓰면 된다.

fastq 포멧은 Phred+64 와 Phred+33 이 있고 게다가 범위 값도 Sanger, Solexa, Illumina (v1.3, v1.5, v1.8) 버전마다 달라서 그야말로 혼돈의 카오스 이다.
자세한 것은 위키피디아 항목을 참조하자.

그냥 간단하게 전부 무시하고 Phred+64 와 Phred+33 사이의 전환법만 알고 넘어가자.

Illumina CASAVA 1.8 버전을 1.3 버전으로 바꾸기 (Phred+64 를 Phred+33 으로 전환)

sed -e '4~4y/!"#$%&'\''()*+,-.\/0123456789:;<=>?@ABCDEFGHIJ/@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghi/' myfile.fastq


Illumina CASAVA 1.3 버전을 1.8 버전으로 바꾸기 (Phred+33 를 Phred+64 으로 전환)

sed -e '4~4y/@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghi/!"#$%&'\''()*+,-.\/0123456789:;<=>?@ABCDEFGHIJ/' myfile.fastq

대략 2013년 이후로는 Illumina Phred+33 으로 강제 통일된 상태나 마찬가지..
만약 보유하고 있는 파일의 fastq 포멧을 알고 싶다면 fastqc 프로그램을 이용해서 파일을 분석해보자.

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/