add scripts for preparing nt and wgs databases.

update blastdb location file, wgs is available now.

scripts/metagenomics/convert_title.py

@@ -0,0 +1,39 @@
+#! /usr/bin/python
+
+"""
+convert nt and wgs data (fasta format) to giNumber_seqLen
+run formatdb in the command line: gunzip -c nt.gz |formatdb -i stdin -p F -n "nt.chunk" -v 2000
+"""
+
+import os, sys, math
+
+if __name__ == '__main__':
+    
+    seq = [] 
+    len_seq = 0
+    invalid_lines = 0
+    
+    for i, line in enumerate(sys.stdin):
+        line = line.rstrip('\r\n')
+        if line.startswith('>'):
+            if len_seq > 0:
+                print ">%s_%d" %(gi, len_seq)
+                print "\n".join(seq)
+            title = line
+            fields = title.split('|')
+            if len(fields) >= 2 and fields[0] == '>gi':
+                gi = fields[1]
+            else:
+                gi = 'giunknown'
+                invalid_lines += 1
+            len_seq = 0
+            seq = []
+        else:
+            seq.append(line)
+            len_seq += len(line)
+    if len_seq > 0:
+        print ">%s_%d" %(gi, len_seq)
+        print "\n".join(seq)
+        
+    print >> sys.stderr, "Unable to find gi number for %d sequences, the title is replaced as giunknown" %(invalid_lines)
+     

scripts/metagenomics/process_BLAST_db.sh

@@ -0,0 +1,18 @@
+echo "This will take several hours to finish due to the size of the databases (about 30GB)..." 
+echo "Getting nt database from NCBI..."
+wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz
+echo "Changing fasta title lines to >ginumber_seqlength..."
+echo "Formatting nt database to chunks of 2GB each..."
+gunzip -c nt.gz | python convert_title.py | formatdb -i stdin -p F -n "nt.chunk" -v 2000
+echo "Remove the zip file, keep the formatted files."
+rm nt.gz
+
+echo "Getting wgs database from NCBI..."
+wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/wgs.gz
+echo "Changing fasta title lines to >ginumber_seqlength..."
+echo "Formatting wgs database to chunks of 2GB each..."
+gunzip -c wgs.gz | python convert_title.py | formatdb -i stdin -p F -n "wgs.chunk" -v 2000
+echo "Remove the zip file, keep the formatted files."
+rm wgs.gz
+
+echo "Job finished"