Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # ~/NCBI_record_parser.py
- #
- #FILE_INPUT = open("blast_output.xml")
- #MINIMUM_INT = 20
- #from Bio.Blast import NCBIXML
- #PARSE_FILE = NCBIXML.parse(FILE_INPUT)
- #READ_LINE = PARSE_FILE.next()
- #
- #
- #
- #for READ_LINE in PARSE_FILE:
- # for alignment in READ_LINE.alignments:
- # for hsp in alignment.hsps:
- # if hsp.identities > MINIMUM_INT:
- # print \n***Seperator*** \n\n ALIGNMENT_TITLE: , alignment.title
- # print length: , alignment.length
- # print HSP_SCORE: , hsp.score
- # print HSP_BITS: , hsp.bits
- # print HSP_EXPECT: , hsp.expect
- # print HSP_NUM_ALIGNMENTS: , hsp.num_alignments
- # print HSP_IDENTITIES: , hsp.identities
- # print HSP_POSITIVES: , hsp.positives
- # print HSP_GAPS: , hsp.gaps
- # print HSP_STRAND: , hsp.strand
- # print HSP_FRAME: , hsp.frame
- # print HSP_QUERY: , hsp.query[0:75] + ...
- # print HSP_: QUERY_START, hsp.query_start
- # print HSP_MATCH: , hsp.match[0:75] + ...
- # print HSP_SBJCT: , hsp.sbjct[0:75] + ...
- # print HSP_SBJCT_START: , hsp.sbjct_start
- #in the format below i wrote the variable names in all upper case thought it might be a little easier to distinguish at a glance.
- #
- #The following variables are declared below:
- #
- # --- FILE_INPUT
- # --- MINIMUM_INT
- # --- VARIABLE_ID
- # --- PARSE_FILE
- # --- READ_LINE
- #
- #
- #
- # set the variable FILE_INPUT equal to the data contained in the file blast_output.xml.
- FILE_INPUT = open("blast_output.xml")
- #
- # established the minimum threshold integer to 20.
- MINIMUM_INT = 20
- #
- # from BioPython import the module NCBIXML
- from Bio.Blast import NCBIXML
- #
- # in the string below the parse() function contained in the module NCBIXML was called to run against the data contained in the variable FILE_INPUT.
- PARSE_FILE = NCBIXML.parse(FILE_INPUT)
- #
- # now we take the string we just created above and we are going to add the next() function. Thie next function tells the interpreter that when parsing the input, move just one line at a time.
- READ_LINE = PARSE_FILE.next()
- #
- #
- # so a horrible visual for this might look like:
- #
- #
- # --- NCBIXML (Module) > parse() > next()
- #
- #
- # its like we are gradually narrowing down the amount of data that will be look at by the interpreter at any given time.
- #
- #
- # --- parse() == read entire file at a time
- #
- # --- next() == read a single line from a file at a time
- #
- #
- #
- # now we begin a for loop to go through our data line by line and process the data. So since above we essentially made READ_LINE a parsing function, when it gets called below, it is being told to run against the data contained in variable PARSE_FILE.
- for READ_LINE in PARSE_FILE:
- #
- #
- # so since:
- #
- # --- FILE_INPUT = open("blast_output.xml")
- #
- # --- A.K.A FILE_INPUT is equal to the data contained in blast_output.xml
- #
- # --- PARSE_FILE = NCBIXML.parse(FILE_INPUT)
- #
- # --- A.K.A PARSE_FILE is equal to the data contained in blast_output.xml and the parse() function is like it stating it is going to stream edit the data in some way.
- #
- # --- READ_LINE = PARSE_FILE.next()
- #
- # --- A.K.A READ_LINE is equal to the data from blast_output.xml, the intention to stream edit, and finally the next() function was chosen as the parsing method. next() is a line-by-line parser
- #
- # now this is where things can get confusing.
- #
- # the command begins what is called commonly refered to in programming as nesting. this would be considered a for loop nest. what happens in for loop nesting can sometimes be difficult to follow.
- #
- # i am going try to break this loop down, by explaining each piece. In the command above:
- #
- #
- # --- for READ_LINE in PARSE_FILE: ; // 1st loop
- #
- # we know its parsing line by line. however when the line below is nested inside:
- #
- # --- for alignment in READ_LINE.alignments: ; // 2nd loop
- #
- # it says to the interpreter, "parse the whole file, one line at a time, but at the end of each line, i want you to pause and run the command below."
- #
- # so then the second loop begins. this one was a little tricky for me too figure out until i did some research. While reading the biopython documentation I ran across this pdf document of the biopython cookbook.
- #
- # --- http://biopython.org/DIST/docs/tutorial/Tutorial.pdf ; // Biopython Tutorial and Cookbook Link
- #
- # there is a section that begins on page 87 of that document called Parsing BLAST output. This chapter really helped put things into focus for me. the most important aspect of this chapter to me was the illustrated database on pages 91-92.
- #
- # what had me the most confused was that i did not realize that the alignment variables used below are set globally when we import the NCBIXML module. additonionally, even though the variables look almost identical, the data they contain references two different tiers of there record database.
- #
- # --- alignment | alignments ; // only one letter difference try not to confuse the two.
- #
- # so the below loop could be read something like this:
- #
- # --- read the document blast_output.xml, at the end of each line pause, then start to parse the alignment entries in the database. this part of the loop parses out individual records at a time.
- #
- # remember nesting is all about narrowing your search results
- #
- #
- for alignment in READ_LINE.alignments:
- #
- #
- # quick recap of wher this program has brought us:
- #
- #FILE_INPUT = open("blast_output.xml") ; // Open file blast_output.xml
- #MINIMUM_INT = 20 ; // Set minimum Threshold value
- #from Bio.Blast import NCBIXML ; // import the biopython module NCBIXML
- #PARSE_FILE = NCBIXML.parse(FILE_INPUT) ; // define PARSE_FILE as the parsed output of blast_output.xml
- #READ_LINE = PARSE_FILE.next() ; // define line-by-line parsing of the xml file by using the builtiin parser tools next() function. read file one line at a time.
- #for alignment in READ_LINE.alignments: ; // at the end of each line, loop through and parse the alignment entries into their own individual records.
- #for hsp in alignment.hsps: ; // each time an individual record is read above, take that record and parse it. this final for loop parses the records into their smallest form, individual strings of data.
- #if hsp.identities > MINIMUM_INT: ; // finally parse each string holds an identities value, and compare that number to the threshold. if the value of the identities string is greater than the threshold integer, continue executing this program. otherwise ignore the string and continue the loop until completion.
- #
- #
- for hsp in alignment.hsps:
- if hsp.identities > MINIMUM_INT:
- print '\n***Seperator*** \n\n ALIGNMENT_TITLE: ', alignment.title
- print 'length: ', alignment.length
- print 'HSP_SCORE: ', hsp.score
- print 'HSP_BITS: ', hsp.bits
- print 'HSP_EXPECT: ', hsp.expect
- print 'HSP_NUM_ALIGNMENTS: ', hsp.num_alignments
- print 'HSP_IDENTITIES: ', hsp.identities
- print 'HSP_POSITIVES: ', hsp.positives
- print 'HSP_GAPS: ', hsp.gaps
- print 'HSP_STRAND: ', hsp.strand
- print 'HSP_FRAME: ', hsp.frame
- print 'HSP_QUERY: ', hsp.query[0:75] + '...'
- print 'HSP_: QUERY_START', hsp.query_start
- print 'HSP_MATCH: ', hsp.match[0:75] + '...'
- print 'HSP_SBJCT: ', hsp.sbjct[0:75] + '...'
- print 'HSP_SBJCT_START: ', hsp.sbjct_start
Advertisement
Add Comment
Please, Sign In to add comment