forked from TScottLUC/rmsystemsproject
-
Notifications
You must be signed in to change notification settings - Fork 0
/
RebaseSetup.py
40 lines (34 loc) · 1.44 KB
/
RebaseSetup.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
import os
import sys
from Bio import SeqIO
from ftplib import FTP
print('Retrieving sequences from REBASE')
ftp = FTP('ftp.neb.com')
ftp.login(passwd='[email protected]')
ftp.cwd('pub/rebase')
# ftp.retrlines('LIST')
with open('dna_seqs.txt', 'wb') as seqs_file:
ftp.retrbinary('RETR dna_seqs.txt', seqs_file.write)
seqs_file.close()
# clean up lines from file received from REBASE
print('Formatting sequences received from REBASE')
with open('dna_seqs.txt') as infile, open('rebase_seqs.fasta', 'w') as outfile:
for line in infile:
if not line.strip(): continue # skip the empty line
#if 'ThisgenbanknumberdoesnothaveaDNAsequenceassociatedwithit.' in line: continue # skip line w/ no seq
line = line.replace("<>", "") # remove <> from lines
outfile.write(line) # non-empty line. Write it to output
infile.close(); outfile.close()
# writes to file in fasta format, makeblast command fails without this
with open('rebase_seqs.fasta', "r") as input_handle, open("formattedSeqs.fasta", "w") as output_handle:
sequences = SeqIO.parse(input_handle, "fasta")
SeqIO.write(sequences, output_handle, "fasta")
output_handle.close()
input_handle.close()
# Delete temp. seqs file
os.system("rm -f rebase_seqs.fasta")
# Make local blast DB
print('Making local BLAST db')
input_file = 'formattedSeqs.fasta'
makeblast_command = "makeblastdb -in "+input_file+" -out LocalRebaseDB -title LocalRebaseDB -dbtype nucl"
os.system(makeblast_command)