2015년 8월 25일 화요일

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/

2014년 5월 15일 목요일

OpenMPI 를 설치하자.

ABySS 를 사용하기 위해 OpenMPI 를 먼저 설치해야 된다. OpenMPI 는 병렬연결된 클러스터 서버를 이용하는 프로그램이다.
ABySS 는 여러 대의 병렬 연결된 서버를 하나로 묶어서 사용할 수 있다. RAM이 모자르던 예전에 유용하게 쓰이던 방법. 요즘은 RAM을 1.5 TB 까지 올릴 수도 있고 1개의 CPU 에 12개의 core 를 넣는 경우도 있으니 굳이 필요한가 싶지만 이런 시스템을 구성할 수 없는 가난한 학교 연구소에서는 여전히 쓰이고 있다.
물론 슈퍼컴퓨터 중에서도 잘만 쓰는 듯...

OpenMPI 는 http://www.open-mpi.org/ 여기서 다운로드 할 수 있다.

압축을 해제한 뒤

./configure --prefix=/target_directory --enable-mpi-thread-multiple

mpi 의 가장 기본적인 기능만 살렸다. 기타 기능은 저도 모르는 관계로 넘어감.

make all
make install

모든 일이 끝나면 target_directory 가 생겨 있을 거다. 여기에 있는 bin/ 디렉토리를 PATH 에 걸고 lib 역시 걸어야 한다.

vim ~/.bashrc

export PATH=/target_directory/bin:$PATH:
export LD_LIBRARY_PATH=/target_directory/lib/:$LD_LIBRARY_PATH:
export OMP_NUM_THREADS=32
export OMP_THREAD_LIMIT=32

ABYSS assembler 를 설치하자.

NGS 로 나온 염기서열의 에러를 수정하는 quake 를 설치하였다.
이제는 수정된 data 를 이용해 genome sequence 를 assemble 하는 프로그램인 abyss 를 설치해보자.

다운로드 받는 곳
http://www.bcgsc.ca/platform/bioinfo/software/abyss

간단한 메뉴얼과 설치방법의 소개
https://github.com/bcgsc/abyss#abyss

설치방법을 보면 먼저
Boost, sparsehash, OpenMPI 가 필요하다고 나온다.

예전에 설치해 둔 적이 있긴 하지만 신버전이 나왔을 테니 다른 글에 설치방법을 설명해두겠다.

sparsehash 는 구글에서 찾을 수 있다.
https://code.google.com/p/sparsehash/

2.0.2 버전 이후 업데이트가 없어서 기존에 설치된 것을 사용하기로 한다. 딱히 설치 방법이 있는게 아니고 그냥 압축을 풀기만 하면 되는지라... 디렉토리 위치만 잘 기억해두도록 하자.

Boost 는 quake 를 설치하면서 1.55.0 version 을 설치하였으니 그 글을 참조하도록..
abyss 매뉴얼에는 1.51.0 과 1.52.0 version 에서 에러가 난다고 한다.
사실 Boost 는 quake 를 위해 설치하였을 뿐 ABySS 에서는 굳이 필요로 하진 않는다. 그냥 압축만 풀어도 된다.

OpenMPI 는 요즘 나오는 CPU 는 모두 multi-core 이므로 이것을 지원하기 위해 나온 툴의 일종이다.
이곳에서 다운 받을 수 있다.
http://www.open-mpi.org/

2014년 5월 현재 1.8 버전이 stable 로 나와 있는데... 오래 전 기억이라 가물가물 하지만 신버전이 잘 안 맞아서 1.4.5 버전으로 설치한 것 같다. 혹시 모르니 신버전으로 다시 재설치 해보기로 했다. 다른 글에 작성 중.


2014년 5월 14일 수요일

quake 를 설치하자.

이제 마지막으로 목표로 했던 quake error correction 을 설치하자.

여기서 프로그램을 받고..
http://www.cbcb.umd.edu/software/quake/

적당한 위치에 압축을 풀고 안으로 들어가 잠시 README 를 감상하자.
쉬워보이는데.. 아무것도 모르는 나로서는 정말 힘들었다.

먼저 Quake/bin 디렉토리 안으로 들어가서 jellyfish/bin 안에 있는 jellyfish 프로그램의 링크 파일을 만들자.

ln -s target_directory/jellyfish/bin/jellyfish jellyfish

그리고 Quake/src 로 들어가서 Makefile 을 조금 고쳐야 한다.

vim Makefile

CFLAGS=-O3 -fopenmp -I/opt/local/include -I --> 요 부분을 찾아라.. 아마 세번째 줄일거다.

이걸 아까 설치한 boost 디렉토리의 위치를 알려줘야 한다. 마지막 boost/ 는 생략해라.
CFLAGS=-O3 -fopenmp -I/target_directory/include -I

저장하며 vim 을 빠져나오고, make 를 실행하면 된다.
여기까지 고생하고 나니 이 부분은 참 쉽다.

마지막으로 PATH 에 Quake/bin 을 추가해주면 된다.

R package 를 설치해보자.

quake 를 사용하기 위해 필요한 또다른 package 인 R 이다.
quake README 에는 쉬우니까 알아서 설치하라고 되어 있는데 이거 설치하면서 정말 그 인간 누군지 때려주고 싶었다. 이거에만 사흘 매달린 듯..

R 은 통계 프로그램으로 오픈소스라서 전세계 사람들이 자유롭게 통계처리를 쉽게 하도록 만들어졌다고 한다.

이곳에서 버전에 맞는 것을 찾아 다운 받자.
http://cran.r-project.org/

그리고 quake 에는 VGAM 이라는 추가적인 package 가 필요하다.
https://www.stat.auckland.ac.nz/~yee/VGAM/
이건데 source 파일로 설치하는 방법을 모르겠다. 그냥 그런가보다 하자.

마찬가지로 소스 파일을 받아 적당한 위치에 압축을 해제하고..

1) ./configure --prefix=/data/ympahk/software/R --with-x=no

--with-x=no 서버에 설치할 때 x11 로 display export 가 안된다고 에러 뜨는 것을 방지하기 위한 용도이다. 아마, 데스크탑에서는 필요없지만 관리자 권한이 없고, x11 을 지원하지 않는 서버의 경우 필요하다. 바로 내 경우 였다.

2) make prefix=/data/ympahk/software/R

뭔가 한참 걸려서 설치를 한다.

저 1번을 몰라서 엄청 헤메었다. 쉬운 건데.. ㅜ_ㅜ

그리고 R 을 path 에 추가시키고 R 을 실행하여 VGAM 을 추가 설치하자.

R 을 실행하면 '>' 만 뜨는데 여기에 install.packages("VGAM") 하면..
미러 서버를 고르라고 하며 나라 이름이 잔뜩 뜬다. 왠지 한가해보이는 서버를 지정하면 설치가 된다. 설치가 안되면 R 에서 빠져나가 다시 들어와서 다시 하자. 반복하다보니 되더라.. -_-..

앗! 아침이 되니 머리가 맑아졌나보다. 다시 구글링을 하다보니 VGAM 을 command line 에서 설치하는 방법을 알게 되었다.

R CMD INSTALL -l ./R-3.1.0/library VGAM_0.9-3.tar.gz

그리고 하면서 느낀 건데.. 관리자 권한 없으면 prefix 쓰지 말자.. 괜히 복잡하다.

jellyfish 설치하기

quake 를 설치하기 위한 2단계 jellyfish 를 설치해보자.

jellyfish 는 DNA sequence 를 kmer 로 분석하는 프로그램이다. quake 는 이 jellyfish 를 사용하여 메모리 사용량을 줄이면서도 더욱 빠르게 error correction 을 수행할 수 있다.

jellyfish 는 http://www.genome.umd.edu/jellyfish.html 이곳에서 받을 수 있고 간단하게 설치법도 나와있다.

설치에 필요한 것은 gcc 4.4 이상과 g++ 4.4 이상 버전이다. 이거면 된다. 따로 설정할 것도 없다.

압축을 풀고 디렉토리 안으로 들어가서...
./configure --prefix=/target_directory
make
make install

하면 target_directory 에 jellyfish/bin 디렉토리에 프로그램이 있다.

직접 quake 를 실행하다가 알게 된 사실 하나를 추가한다.

http://www.cbcb.umd.edu/software/jellyfish/

여기에 보면 jellyfish 2.x 버전과 1.x 버전이 따로 업데이트 되고 있는게 보인다.
2.x 버전은 jellyfish 본연의 기능만을 넣고 있지만 1.x 버전에는 quake 용 명령어가 있다.
즉, quake 용으로 jellyfish 를 사용하려면 1.x 버전을 사용해야 된다는 거다.