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 포맷이
수월하니깐요.. :)

그럼....