SOAPaligner 결과에 Header넣고 SAM/BAM 변환

NGS alignment tool 중 BGI의 SOAPaligner가 있는데, 그 중 SAOP2 (SOAP3는 GPU로 인해 제한적으로 사용가능하기 때문에 SOAP2를 주력으로 사용하고 있다.) 작업 후 sam->bam으로 만드는 작업 중 문제에 대해서 얘기하고자 한다.

SOAPaligner 작업을 하면 sam/bam이 아닌 text파일이 생성된다.

그러나 이 text파일은 우리가 routine하게 사용하는 프로그램의 input으로는
사용할수가 없다는 것이 함정!!!

그래서 일련의 변환 작업을 필요로 한다.

오늘은 SOAP의 sam/bam 변환 작업에 대해서 알아보도록 하겠습니다.

이제와서 갑자기 SOAP을 얘기하냐고 물으신다면
요즘 SOAP을 가까이 하는 관계로 ㅎㅎㅎㅎ
팔은 안으로 굽지만 포스팅은 많이 사용하고 있는 tools 중심이니깐~ :)


첫번째 SOAPaligner는 알아서 잘 alignment 하도록 합니다.

이제 본격적인 변환 작업입니다.
다행히 BGI에서 soap output 형식을 sam 형식으로 변환시켜주는 perl 스크립트를 제공하고 있습니다.
그이름도 아름답도다 soap2sam.pl

근데 함정은 soap으로 alignment된 결과만 sam 형식으로 변환해준다는 점.
sam에 있어야할 header같은 정보는 무시한다는 것!!!


soap에는 없는 sam의 header파일 만드는 방법은
그냥 간단히 reference파일을 사용하여 header로 변경하는 script코드 작성하여도 되고

import os,sys

from Bio import SeqIO

try:

    input_fa = sys.argv[1]

except:

    print "python convert_ref_sam_header.py <ref> > <outpu.header>"

    exit(1)

print "@HD\tVN:1.0\tSO:unsorted"

for rec in SeqIO.parse(open(input_fa), format='fasta'):

    name = rec.description

    seq = rec.seq.tostring()

    name = (name.split()[0]).strip()

    print "@SQ\tSN:%s\tLN:%s" % (name,len(seq.strip()))

print "@PG\tID:soap2\tPN:soap2\tVN:2.2.1"

아니면

samtools faidx reference.fa

작업 후

cut -f 1,2 reference.fa.fai | awk '{print "@SQ\tSN:"$1"\tLN:"$2}' > sam.header

필요한 컬럼과 문자열만 조합하여 추출하여도 된다.
(이 경우 상하단에 @HD와 @PG 헤더를 추가해주어야 한다는 불편함이 있다. )


이런 작업을 거치면
soap을 sam을 거쳐 bam파일을 만들어서 다음 작업을 이어나갈수 있게 된다. :)


혹 fai파일을 만들기 위해 samtools faidx 작업을 수행하다가 에러가 발생하는 경우
fai파일을 만드려는 reference파일을 확인해보기 바란다.
fasta형식의 서열간 빈 line이 있는 경우 에러를 발생시킨다.
이런 경우 간단히 vi로 들어가서":g/^$/d"로 공란을 제거해주면 정상적으로 fai파일을 만들게 된다.

그리고 가장 쉬운 방법은 나중에 알려주는 법
(역시 메뉴얼을 잘 읽어야 한다는... 아놔....)

samtools view -bT reference.fa input.sam > output.bam

위에와 같이 뻘짓 필요없다는... ㅋㅋㅋㅋㅋㅋ ㅠ.ㅜ 아...;;; 제길..

Velvet으로 Assembly하기

결론은 de novo는 흠좀무라는... ㅋㅋ

species가 무엇이 되던간에 말입니다. ㅎㅎ


요즘 몇몇 곰팡이를 조립중에 있는데
SOAPdenovo, Velvet, SPAdes, SSPACE를 중점적으로 사용하고 있습니다.
ALLPATHS-LG는 read 준비하는 것까지는 작동을 하는데
그다음에 본 작업에서 error를 뱉고 죽어버리는 어처구니 없는 상황을 연출해주는 바람에
일단 pass했습니다.

여러 strategy으로 접근 중에 있는데...

아... 본 글은 그게 중요한게 아니라... (다음 기회에 언급 하도록 하겠습니다.)

Velvet assembly할 때 long insert library (Mate-Paired Library)를 한개만 지원하는 것 같아그 문제를 해결하가 위한 tip입니다.

해결 방법은 다음 링크에 나와 있습니다.


단 EMBOSS의 "revseq"를 사용하면 안습이 된다는 것 빼고는 말이죠~ :)

그래서 준비 했습니다.

#!/usr/bin/python
import os,sys
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

def readread(s):
        return [s.readline(),s.readline(),s.readline(),s.readline()]


try:
        fq = sys.argv[1]

except:
        print "python command.py <input_fastq>"
        print ""
        exit(1)

fq_open = open(fq,"r")

fq_write = open("%s.rev" % (fq),"w")

s1read = readread(fq_open)

while(s1read[0]):
        h1, h2 = s1read[0].strip(), s1read[2].strip()
        s1, q1 = s1read[1].strip(), s1read[3].strip()

        rev = Seq(s1,IUPAC.unambiguous_dna).reverse_complement()
        qual = q1[::-1]

        fq_write.write("%s\n%s\n%s\n%s\n" % (h1,rev,h2,qual))

        s1read = readread(fq_open)

fq_write.close()

-위 소스는 python 2.6 이상에서 작동합니다. :)

velvet으로 조립 시 mate-paired Library가 여러개라면 서열을 변경 한 이후 사용해 보시기 바랍니다. :)

The Great Good Place는 폐쇄중;;

다름이 아니라 무료로 호스팅하고 있는 Hostinger 측에서
제가 사용하고 있는 계정에서 다수의 이메일이 발송되어
자신들의 IP가 블랙리스트에 등록될까 염려되어 제 계정을 막았다는
소식을 접했습니다.

그래서 문의 결과 어떻게 할 수 없다는 답변을 받고 일단
제 도메인들을 모두 이 gwlee.blogspot.kr 포워딩 시켜 놓은 상황입니다.

어차피 아직 공식 개장 하지 않아서 적어놓은 글은 많지 않았지만..
WP셋팅하고 블로깅을 하려고 했는데 이런 시련이 닥치네요 ㅎㅎ

조만간 다른 안식처를 마련하여 공개하도록 하겠습니다. :)