너도 AB1파일 가지고 pdf 만들 수 있어 (2)

너도 AB1 파일 사용해서 pdf 만들 수 있어 첫번째 글에 이은 두번째!!!

chatGPT한테 물어 물어 좀 쉽게 작성을 해보려고 했으나...

아.. 이것 또한 쉽지 않습니다.

예제를 보여줘서 학습을 시키면 좋을텐데... 잘 모르게서

그냥 간단히 matplotlib 코드만 간단히 활용해 보았습니다.

ab1 파일을 분석하는 것은 그냥 노가다로... ㅎㅎ 

우선 그냥 chatGPT한테 

"ab1파일 parsing해서 pdf로 만드는 스크립트 작성해줘" 하면 

from Bio import SeqIO

import matplotlib.pyplot as plt

# ab1 파일 경로 설정

abif_file = "example.ab1"

# ab1 파일 파싱하여 전기영동 결과 얻기

abif_record = SeqIO.read(abif_file, "abi")

electropherogram = abif_record.annotations['abif_raw']['DATA2']

# matplotlib를 이용하여 전기영동 결과 그래프 그리기

plt.plot(electropherogram)

plt.title("Electropherogram")

plt.xlabel("Base Position")

plt.ylabel("Signal Strength")

plt.show()

이렇게 알려줍니다.

이거 그대로 작성하면 우리가 기대했던 이미지하고 영 딴판의 결과가 저장됩니다.

그래서 세부 내용들은 Applied Biosystems에서 제공하는 Data File Format pdf를 좀 들쳐봐야합니다.

그래서 코드를 좀 수정 하면...

from Bio import SeqIO
import matplotlib.pyplot as plt

# ab1 파일 경로 설정
abif_file = "example.ab1"

# ab1 파일 파싱하여 전기영동 결과 얻기
abif_record = SeqIO.read(abif_file, "abi")
poc = record.annotations['abif_raw']['PLOC1']
a = record.annotations['abif_raw']['DATA10']
c = record.annotations['abif_raw']['DATA12']
g = record.annotations['abif_raw']['DATA9']
t = record.annotations['abif_raw']['DATA11']

data = {"A":a_seq, "C":c_seq,"G":g_seq,"T":t_seq}

plt.figure(figsize=(len(poc)/10,5))
for base, color in zip("ACGT",["g","b","k","r"]):
    plt.plot(data[base],color=color)

tmp = [None]*len(a_seq)

i=0
for pnt in poc:
    tmp[pnt]=seq[i]
    i+=1

plt.xticks(range(len(tmp)),tmp, fontsize=6)
plt.savefig('output.pdf')

이 코드를 사용하면 응? 좀 이상하지만 약간 그럴싸한 이미지가 보이실겁니다.

시퀀칭 업체에 Sanger Sequencing에 맡기면 fasta파일과 ab1파일과 함께 오는 pdf파일과 다르긴 하지만 얼추 비슷한...

그럼 다음 기회에는 Sanger Sequencing 맡기면 함께 받아 볼 수 있는 pdf 파일을 만들어보기로 해봐요. 내년 쯤에는 할 수 있지 않을까 합니다. :)

출처: @ye._.vely618

genpept 파일을 parsing하고 싶을때

genpept는 무슨 파일인고?

자세한 설명은 >여기<

간단히 얘기하자면 예전에 사용하셨던 genbank 파일 포맷입니다.
이름만 바꾼건지 모 그렇습니다. ㅎㅎ
biopython에서도 genpept가 아니라 genbank를 사용해서 접근하면 잘 parsing됩니다.

import os,sys
from Bio import SeqIO
try:
    inFile = sys.argv[1]
except:
    print ''
    exit(1) 

for seq_record in SeqIO.parse(inFile, "genbank"):
    print (seq_record.id)
    print (seq_record.name)
    seq_anno = seq_record.annotations
    print (seq_anno['accessions'])
    print (seq_record.seq)

이렇게 코드 짜서 사용하시면 됩니다. annotations안의 정보는 dict형식이라서 dict사용하는방식으로 확인 할 수 있습니다.

간만에 gp파일 다운받아서 정리하는 김에 biopython의 genbank 사용법 정리해보았습니다.  :)

크.. WSL이 있으니 참 편하긴 합니다.

출처: @sana_twice.09

Microbiome Database를 만들어볼까? -NCBI편 4-

지난번 글에서 받기시작했던 nt.gz파일은 잘 받아졌나요?

그럼 이제 이 파일에서 무엇인가 뽑아내야 겠죠?

무엇을 뽑아내느냐?
16S rRNA를 뽑아낼겁니다.
어떻게?
다음 스크립트를 작성해 봅시다.

$vi parser_nt.py

import glob,sys,re,gzip
from Bio import SeqIO

try:
        input_fa = sys.argv[1]
except:
        print "fasta_split.py <in.fa.gz>"
        exit(1)
for rec in SeqIO.parse(gzip.open(input_fa),format='fasta'):
        desc = rec.description
        seq = str(rec.seq)
        name = desc.strip().split('\x01')[0]
        if name.upper().find(' 16S RIBOSOMAL RNA ') != -1 or name.find(' 16S RIBOSOMAL RNA,') != -1:
                if name.upper().find('MITOCHONDRIA')!= -1:
                        pass
                elif name.upper().find('CHLOROPLAST')!= -1:
                        pass
                elif name.upper().find(' 23S ') != -1:
                        pass
                else:
                        if name.upper().find(' SIMILAR ') != -1:
                                pass
                        elif name.upper().find(' INTRON') != -1:
                                pass
                        elif name.upper().find(' PLASTID') != -1:
                                pass
                        else:
                                print '>%s\n%s\n'%(desc,seq)

※ 위의  스크립트를 수정해서 입맛에 맞게 교정하시고 사용하시면되겠습니다.

$python parser_nt.py nt.gz > nt.fasta
위에 스크립트 실행시키면 떡 하니 수 gb 짜리 파일이 하나 나올겁니다.

이 파일안에는 nt서열 중에서 16S rRNA 서열 (대신 mitochondria와 chloroplast의 16S rRNA는 제외하고 이것도 16S rRNA지만 저한테는 일단 필요없어서 뺐습니다. 사용하고 싶으시면 사용하셔도 됩니다. :) )이 500만개 정도 들어 있습니다.

이 파일에는 온갖 종의 16S rRNA 서열이 있지만 문제가 있다는 점!

16S rRNA 서열이 품질이...
어떤 녀석은 full 서열이 있지만 어떤녀석은 서열의 일부만 가지고 있는 경우가 있습니다.

그래서 그런 녀석들을 잘 확인해서 제거를 하던지 merge를 하던지... 

그건 개인 취향으로 남겨 놓도록 하겠습니다. :)

aka 필터링을 하던 안하던 문제 생기면 그 문제의 책임은 오롯이 당신의 것!!

그리고 즐거운 추석보내시기 바랍니다. :)

출처 SM

Microbiome Database를 만들어볼까? -NCBI편 3-

지난 시간에 이곳에서 귀한 자료를 받아 봤을 겁니다.
그게 무엇이냐!!!

지금까지 공개된 bacteria의 서열들이죠 정확하게 말하면
현재까지 수많은 연구자들이 자발적(자의타의)으로
공개해준 것을 NCBI가 아름답게 정리한 RNA서열들만 다운받았습니다.

근데 우리는 RNA 서열이 아니라 rRNA 서열이 필요하죠

이전 글에서 우리는 rna라는 폴더 안에 12만개에 달하는 파일들을 다운받았습니다.
아마 그냥 ls하시면 어쩌구 저쩌구 long 할겁니다.
ls로 불러오기에 item이 너무 많다 이거죠
그럼 어쩌지?? @.@

우리에겐 귀도 훃님의 파이썬이 있지 않겠습니꽈

간단하게 다음과 같은 스크립트를 뚝딱 뚝딱 만들어보죠

import os,glob
from Bio import SeqIO
for files in glob.glob('rna/*gz`):
    for rec in SeqIO.parse(gzip.open(files), format='fasta'):
        name =rec.description
        seq = rec.seq
        if name.find('[product=16S ribosomal RNA]') != -1:
            print '>%s\n%s\n'.format(name,seq)

(python2.7에 Biopython이 설치되어 있어야 하고 *_rna_form_genomic.fna.gz 파일이 rna폴더 밑에 위치하고 있어야 합니다.)

$python script.py > ncbi.fa

라고 해주면 헤더에 "[product=16S ribosomal RNA]"가 포함된 aka 우리가 원하는 바로 그것! 16S rRNA 서열을 각 종에서 샤샤샥 ncbi.fa라는 파일에 저장할 수가 있습니다.

-12만개 파일 읽어 오는거라 순차적으로 하면 3-4시간 걸릴것이고
subprocess로 잘 해주시면 적어도 십수분? 30분이내면 충분히 끝날각 되겠습니다.

출처: SM Town

Fasta 형식의 파일에서 빈서열 제거

가끔씩 Blast를 수행하고자 formatdb를 수행 할 때,
다음과 같은 에러를 접한 적이 있으리라 본다.

[formatdb] WARNING: Cannot add sequence number XXXXX XX.XXX.XXX.
 because it has zero-length.
[formatdb] FATAL ERROR: Fatal error when adding sequence to BLAST database.

formatdb를 수행하려는 fasta 서열에
빈 서열을 가지고 있기 때문에 나오는 에러로
빈 서열을 제거하면 OK!

vi check.py

import glob,sys
from Bio import SeqIO


file = sys.argv[1]


ow = open(file.split('.')[0]+'.check','w')
seqs = SeqIO.parse(open(file), format='fasta')
 for rec in seqs:
        name = rec.description
        seq = rec.seq.tostring()


        if len(seq.strip()) != 0:
                ow.write('>'+name+'\n')
                ow.write(seq+'\n')


ow.close()

[test]$python check.py test.fasta