FASTX Tool Kit에서 phred33 사용하기

FAST{A/Q} 파일을 조작(그 조작말고요 고갱님~ ㅎㅎ)을
보다 수월하게 하는 명령어들을 모아놓은 아주 좋은 패키지 입니다.
(사실 전 몇개 안쓰지만요;;; 잇힝~ )

그러나 지금은 이 툴보다 좋은 프로그램들이 많이나와 있는 상태인데..

그래도 끄적 끄적... ㅎㅎ
혹시 다른 페이지들이 폐쇄되는걸 대비해서
백업용으로.....

fasta tool에서 fastx_quality_stats를 사용할 때 quality score관련한 error 메세지를 볼 수 있는데 이것은 fasta tool kit이 phred기준 +64를 좋아해서 그런다는...

그래서 fastx_quality_stats -i input.fastq -N -Q 33 -o quality.txt

이렇게 하면 된다고 이곳에 나와 있었습니다.

출처: shengliblog

NGSQCToolkit 사용기

아.... 이제서야
지난번에 언급했던
NGS QC Toolkit을 사용해 봤습니다. :)

라이브러리와 perl 모듈을 잘 설치해주면
큰 문제없이 잘 돌아가는것을 확인했고
multi-thread로 실행하는 경우 음.. 빠르더라구요 ㅎㅎ
시간 체크는 못해봤는데...
지금 시간체크 하면서 돌리는게 있으니
정리해서 올리도록 하겠습니다. :)

1. 들어가기전
일단 시스템에 gd관련 라이브러리가 있는지 확인하시고
gd와 libgd-graph 등등 관련 라이브러리를 설치해주시기 바랍니다.
그리고 perl 모듈들이 모두 설치되어 있는지 확인해서 안되어 있다면
설치해주시면 되겠습니다.
gd라이브러리가 없으면 펄의 GD::Graph 설치할때 설치가 안되더군요;;
에러가 나서 몬가 하고 있었는데.. ㅎㅎ 여하튼...
모 이런저런 라이브러리와 모듈을 확인하시고 잘 설치하면
사용하는데 문제 없습니다. :)

2. 사용하기
NGQQCToolkit에는 크게 4가지의 서브 카테고리로 구분되어져 있더군요
1) 포맷 변경
2) QC
3) 자료 통계
4) Trimming

2.1 Format Convert
Fastq -> {454 | Fasta}: Fastq를 454(Fastq,Qual), Fasta 포맷으로 변환

{SangerFastq | SolexaFastq} -> IlluFastq: Sanger와 Solexa의 qual를 Illumina의 통일된 qual score range로 변환 (다만, 1.5+ 로 하는지 1.8+ 로 하는지는 확인 못했습니다.)


2.2 QC
454{QC|QC_PE|QC_PRLL}: 454 데이터를 input으로 하는 QC tools
Illumina와 다르게 QC_PE가 있는건 454의 경우 paired-end로 sequencing 하는 경우는 좀 특별해서 구분해둔듯.. :) (단, input은 SFF 포맷이 아닌 서열 파일과 Quality score파일로 구분해서 입력해야 사용가능하다.)

Ill{QC|QC_PRLL}: 일루미나 read를 처리하는 tools, 454와는 다르게 single-end와 paired-end를 따로 구분하지 않고 -se, -pe 옵션으로 처리하도록 만들어놨다는 점~ :)

PRLL 접미사는 병렬처리를 지원하는 스크립트입니다.
PRLL tools에서는 -c를 사용해서 multi-core를 사용하는데에 반해
일반 tools는 -p 옵션을  사용해도 multi-core를 사용하지 않는 점이 있었습니다.


2.3 Statistics
AvgQuality.pl: quality score 파일을 입력받아 점수를 계산하는 tool
N50Stat.pl: fasta파일을 input으로 받아 N50을 계산하는 tool

2.4 Trimming
AmbiguityFiltering.pl:
HomopolymerTrimming.pl:
TrimmingReads.pl:

결과로 제공되는 figure도 나름 괜찮습니다. :)
속도도 multi-core를 사용하던 안하던 만족할만한 수준이었습니다.
(제가 in-House로 제작한 script가 느린것도 있겠지만요.. ㅎㅎ )


자세한 사용법은 저보다 영어 못하시는 분은 없을테니 메뉴얼 보세요~ ㅎㅎ

>>메뉴얼보러가기<<


추가정보
paired-end fastq raw파일로 3-4g정도의 파일을 single cpu로 처리하는데
2시간에서 2시간 반내외정도로 확인되었습니다. :)

SFF 파일은 어떻게 작업을 해야하나...

아....

딱히 인연이 없던 454 파일을 작업할 기회가 생겨서..
다음과 같이 스크립트를 좀 작성했습니다.

454에서 제공하는 Data analysis를 이용하지 않아서 좀 거시기합니다.
(Homopolymer trimming은 제공하지 못하고 있습니다. ㅎㅎ)

convertsff.py 

SFF파일을 fastq 혹은 fasta, qual 파일로 변환하는 스크립트
fastq로 변환하는 경우 illumina 1.3/1.5+ score로 변환 됩니다.
biopython이 설치되어 있어야 합니다.


SFF_Filter.py

illumina 데이터와는 일단 길이가 차이가 나니.... ㅎㅎ
filtering 스크립트를 간단히 만들었습니다.
cutoff base quality와 cutoff read length는 사용자가 설정 할 수 있게 하였습니다. :)
그리고 N의 포함 정도와 cutoff base quality 포함 정도는 고갱님의 의견을 반영하여
박하지 않게 설정해서 fixed시켜 놨습니다.
맘대로 수정하셔도 무방합니다. :)
다만 더 좋은 옵션이나 방법으로 업데이트 하셨다면 공유를 해주시면 더더욱 감사드리겠습니다.

데이터 변환 및 필터링이 끝나고 나면
QC를 해봐야 겠죠?

좀 간단히 결과를 이쁘게 그려주는게 어디 없을까 하고 있었는데
prinseq라는 프로그램이 있어 잠시 사용해봤습니다.
사용방법은 어렵지 않아요~ :)


Base Quality가 Phred +33인 경우

perl prinseq-lite.pl -verbose -fastq <input.fq> -graph_data <output.gd> -out_good null -out_bad null

Base Quality가 Phred +64인 경우

perl prinseq-lite.pl -verbose -fastq <input.fq> -graph_data <output.gd> -phred64 -out_good null -out_bad null

Quality Check 결과물을 Html로 확인하는 경우
perl prinseq-graphs.pl -i <output.gd> -html_all -o <output_name>

Quality Check 결과물을 png로 확인하는 경우
perl prinseq-graphs.pl -i <output.gd> -png_all -o <output_name>





추후에 시간이되면
Data analysis 프로그램을 설치해서 작업하는 단계나 방법에 대해서 설명하도록 하겠습니다. 그리고 추가적으로 NGS QC toolkit을 사용해서 QC하는 것도..
S대 L모군이 찾아논건데 괜찮아보여서 테스트 해볼까 합니다.
Illumina 외에 454도 지원하고 제일 매력적인건 multi-thread를 지원한다는것!!

모 여하튼...

다음기회에~ :)

파일의 포맷을 변환하는데 필요한 것들

내가 아니란 말이닷!!! ㅋㅋ

python에서 Biopython을 이용하여
간단하게 convert하는 샘플 코드를 제공하고 있으니
여러분들도 쉽게 만들수 있어요~ :)
Biopython에서 제공하는 Tutorial 


오늘 문의가 들어온 파일은 sff파일
Roche의 454 GS FLX? sequencing 결과파일로....
ABI와 함께 illumina한테 밀려서 뒷방으로 들어앉은 파일 포맷입니다.
그러나 아직도 쓰는 이유는 read 길이가 길기때문 :)

그렇습니다. PacBio도 Nanopore다 디립다 길게 sequencing해준다는
애들이 있습니다. 그런데 왜 옛날꺼 쓰냐?? PacBio는 base quality가 안습이고,
Nanopore는.... 언제 출시일지 전 잘 모르겠습니다. 업자가 아닌관계로 ㅎㅎ

그래서 위의 길게 sequencing 해준다는 시퀀서를 제외하고는 Roche의 454가 read 길이가 가장 길다고 할 수 있겠습니다. NGS중에선 말이죠

그런데 sff파일을 보려고 하면 문제가 생깁니다.
권모씨께서 문의를 한것이 그것때문인지는 모르겠지만 걍 일반인이
sff파일을 걍 직접 볼수가 없습니다. 왜냐 binary파일이니깐요(sff파일이 binary라고
알고 있는데  직접 다뤄본적이 없어서... ㅎㅎ )

그래서 사람이 볼수 있게 파일을 변환시켜줘야 한다는 겁니다.

convertSff.py

#!/usr/bin/python

import os, sys
from Bio import SeqIO

try:
inputSFF = sys.argv[1]
outputPREFIX = sys.argv[2]

except:
print "Usage: python convertSFF <input.sff> <output_name>"
print ""
exit(1)


SeqIO.convert(inputSFF,"sff","%s.fasta"%(outputPREFIX), "fasta")
SeqIO.convert(inputSFF,"sff","%s.quality"%(outputPREFIX), "qual")
SeqIO.convert(inputSFF,"sff","%s.fastq"%(outputPREFIX), "fastq")

권모씨의 요청으로 급조한 날림 convert python 코드 ㅋㅋ
이 스크립트를 수행하면 세개의 파일이 나오게 될것으로 예상됩니다. ㅎㅎ
안나오면 어쩔수없고... ㅎㅎ


아.. 그리고 사족으로 LT사의 SOLiD의 경우 우리가 알고 있는 서열과 달리
첫 염기 서열만 서열이고 그 다음부터는 A/G/T/C 알파벳이 아닌 숫자로 되어있는데..
이걸 굳이 변환해서 reference geneome에 mapped 작업하지 말라고 합니다.
Re-sequencing하는 경우라면 변환해서 mapping하지 말고 원래 원본 파일 그대로를
input으로 하는 align 프로그램을 사용해서 mapped한 다음에 그 다음 작업을
일반적으로 사용하는 samtools나 GATK같은 프로그램을 사용하라고 합니다.
(다들 알고있는거 한번더 상기 시켜드렸습니다. 혹시 아나요 SOLiD 포맷을 분석하게 될지.. ㅎㅎ)

분석시 raw 파일을 사용해야 하는 이유는 SOLiD만의 월등한 quality 효과를 볼수 있어서
그러지 않겠나하는....  믿거나 말거나 저 혼자만의 생각입니다.. ㅎㅎ
다만, 타사 제품과 다르게 복잡하게 숫자로 표현한건 아니겠죠...
나름의 숨은 뜻이.... 쿨럭.. (설마... 간지용;;;;; )

Re-sequencing이 아닌 denovo일 경우 모 어쩔수 없이 fastq파일로 변환을 해야 하지 않을까 합니다. assembly 프로그램을 작동시키려면 아무래도 SOLiD format보다는 fastq 포맷이
수월하니깐요.. :)

그럼....

Sequence Quality

근래 많은 차세대 서열 분석기의 등장으로
서열정보가 너무 넘쳐나고 있는 가운데
각 시퀀서에서 나오고 있는 서열이 과연 믿고 쓸만한가에 대한
관심이 없을수 없게 되었다.
다행히 서열 결과에는 각 서열이 쓸만한 것인지 이건 연구/진단에
사용하기 부적합한 것인지 판단할 수 있는 품질 정보가 ASCII 코드 값으로
들어가 있다.

참고: Fastq format

시퀀서마다 약간씩 다르지만....
일단 점수 높은 놈은 좋다는건 진리.. ^^