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

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