filter_blast_results.py

blastoutput = "../input/R_2013_11_30_15_23_49_user_CTL-148-Metagenoma_Mauro_Eduardo2.CTLT_PGM.IonXpress_001_Q20_TrimRL_AT_80pb_cdhit97.all.blast_output"

from subprocess import call

#http://biopython.org/DIST/docs/api/Bio.SearchIO.BlastIO-module.html
from Bio import SearchIO

blast_results = SearchIO.parse(blastoutput, 'blast-text')

results = {
  '16S': [],
  '18S': [],
  '28S': [],
  'COI': [],
  'left out': [],
  # '18S and 28S': []
}

coi_list = [" COI ", "(COI)", "cytocrome oxidase", " COX ", "(COX)", " CO1 ", "(CO1)"]
s28_list = ["28S", "large subunit ribosomal", "large ribosomal subunit"]


for qresult in blast_results:

    found_one = False

    for hit in qresult.hits:        
        if "16S" in hit.description:
            results['16S'].append(qresult.id)
            found_one = True
        if "18S" in hit.description:
            results['18S'].append(qresult.id)
            found_one = True
        # if "28S" in hit.description:
        if any(t in hit.description for t in s28_list):
            results['28S'].append(qresult.id)
            found_one = True
        #COI
        if any(t in hit.description for t in coi_list):
            results['COI'].append(qresult.id)
            found_one = True        
    if not found_one:
        results['left out'].append(qresult.id)


#remove duplicated
for group in results:
    results[group] = set(results[group])

results['18S and 28S'] = results['18S'].intersection(results['28S'])

# records = []

records = {
  '16S': [],
  '18S': [],
  '28S': [],
  'COI': [],
  '18S and 28S': [],
  'left out': []
}

parser = open(blastoutput)

groups = ["16S","18S","28S","COI","18S and 28S", "left out"]

record = ["as"]
found_query = False

for line in parser:    
    if line.startswith("Query="):
        found_query = True
        #This is used to ignore the header of the blast file
        if record[0].startswith("Query="):
            id = record[0].split(" ")[1].strip()
            for group in groups:
                if id in results[group]:
                    records[group].append(record)
        record = []
    if found_query:
        record.append(line)
#treat last record
# if record[0].startswith("Query="):
id = record[0].split(" ")[1].strip()
for group in groups:
    if id in results[group]:
        records[group].append(record)


# print(len(records))
# print(records[1])
for group in groups:
    outputfile = open("%s.blast_output" % (group), "w")
    blastheader = open("../input/blastheader.txt", "r").read()
    outputfile.write(blastheader)
    for query in records[group]:
        # for line in query:
        outputfile.writelines(query)