o Rŀgz@sdZddlZddlmZddlmZddlmZddlmZddlm Z ddlm Z d Z d Z e e Ze e Zdd d Zd dZddZddZddZddZGdddZGdddeZGdddZedkrxddlmZedSdS) a*Bio.SearchIO parser for BLAT output formats. This module adds support for parsing BLAT outputs. BLAT (BLAST-Like Alignment Tool) is a sequence similarity search program initially built for annotating the human genome. Bio.SearchIO.BlastIO was tested using standalone BLAT version 34, psLayout version 3. It should be able to parse psLayout version 4 without problems. More information on BLAT is available from these sites: - Publication: http://genome.cshlp.org/content/12/4/656 - User guide: http://genome.ucsc.edu/goldenPath/help/blatSpec.html - Source download: http://www.soe.ucsc.edu/~kent/src - Executable download: http://hgdownload.cse.ucsc.edu/admin/exe/ - Blat score calculation: http://genome.ucsc.edu/FAQ/FAQblat.html#blat4 Supported Formats ================= BlatIO supports parsing, indexing, and writing for both PSL and PSLX output formats, with or without header. To parse, index, or write PSLX files, use the 'pslx' keyword argument and set it to True. >>> # blat-psl defaults to PSL files >>> from Bio import SearchIO >>> psl = 'Blat/psl_34_004.psl' >>> qresult = SearchIO.read(psl, 'blat-psl') >>> qresult QueryResult(id='hg19_dna', 10 hits) >>> # set the pslx flag to parse PSLX files >>> pslx = 'Blat/pslx_34_004.pslx' >>> qresult = SearchIO.read(pslx, 'blat-psl', pslx=True) >>> qresult QueryResult(id='hg19_dna', 10 hits) For parsing and indexing, you do not need to specify whether the file has a header or not. For writing, if you want to write a header, you can set the 'header' keyword argument to True. This will write a 'psLayout version 3' header to your output file. >>> from Bio import SearchIO >>> qresult = SearchIO.read(psl, "blat-psl") >>> SearchIO.write(qresult, "example.psl", "blat-psl", header=True) (1, 10, 19, 23) >>> import os >>> os.remove("example.psl") Note that the number of HSPFragments written may exceed the number of HSP objects. This is because in PSL files, it is possible to have single matches consisting of noncontiguous sequence fragments. This is where the HSPFragment object comes into play. These fragments are grouped into a single HSP because they share the same statistics (e.g. match numbers, BLAT score, etc.). However, they do not share the same sequence attributes, such as the start and end coordinates, making them distinct objects. In addition to parsing PSL(X) files, BlatIO also computes the percent identities and scores of your search results. This is done using the calculation formula posted here: http://genome.ucsc.edu/FAQ/FAQblat.html#blat4. It mimics the score and percent identity calculation done by UCSC's web BLAT service. Since BlatIO parses the file in a single pass, it expects all results from the same query to be in consecutive rows. If the results from one query are spread in nonconsecutive rows, BlatIO will consider them to be separate QueryResult objects. In most cases, the PSL(X) format uses the same coordinate system as Python (zero-based, half open). These coordinates are anchored on the plus strand. However, if the query aligns on the minus strand, BLAT will anchor the qStarts coordinates on the minus strand instead. BlatIO is aware of this, and will re-anchor the qStarts coordinates to the plus strand whenever it sees a minus strand query match. Conversely, when you write out to a PSL(X) file, BlatIO will reanchor qStarts to the minus strand again. BlatIO provides the following attribute-column mapping: +----------------+-------------------------+-----------------------------------+ | Object | Attribute | Column Name, Value | +================+=========================+===================================+ | QueryResult | id | Q name, query sequence ID | | +-------------------------+-----------------------------------+ | | seq_len | Q size, query sequence full | | | | length | +----------------+-------------------------+-----------------------------------+ | Hit | id | T name, hit sequence ID | | +-------------------------+-----------------------------------+ | | seq_len | T size, hit sequence full length | +----------------+-------------------------+-----------------------------------+ | HSP | hit_end | T end, end coordinate of the last | | | | hit fragment | | +-------------------------+-----------------------------------+ | | hit_gap_num | T gap bases, number of bases | | | | inserted in hit | | +-------------------------+-----------------------------------+ | | hit_gapopen_num | T gap count, number of hit gap | | | | inserts | | +-------------------------+-----------------------------------+ | | hit_span_all | blockSizes, sizes of each | | | | fragment | | +-------------------------+-----------------------------------+ | | hit_start | T start, start coordinate of the | | | | first hit fragment | | +-------------------------+-----------------------------------+ | | hit_start_all | tStarts, start coordinate of each | | | | hit fragment | | +-------------------------+-----------------------------------+ | | match_num | match, number of non-repeat | | | | matches | | +-------------------------+-----------------------------------+ | | mismatch_num | mismatch, number of mismatches | | +-------------------------+-----------------------------------+ | | match_rep_num | rep. match, number of matches | | | | that are part of repeats | | +-------------------------+-----------------------------------+ | | n_num | N's, number of N bases | | +-------------------------+-----------------------------------+ | | query_end | Q end, end coordinate of the last | | +-------------------------+-----------------------------------+ | | | query fragment | | | query_gap_num | Q gap bases, number of bases | | | | inserted in query | | +-------------------------+-----------------------------------+ | | query_gapopen_num | Q gap count, number of query gap | | | | inserts | | +-------------------------+-----------------------------------+ | | query_span_all | blockSizes, sizes of each | | | | fragment | | +-------------------------+-----------------------------------+ | | query_start | Q start, start coordinate of the | | | | first query block | | +-------------------------+-----------------------------------+ | | query_start_all | qStarts, start coordinate of each | | | | query fragment | | +-------------------------+-----------------------------------+ | | len [#]_ | block count, the number of blocks | | | | in the alignment | +----------------+-------------------------+-----------------------------------+ | HSPFragment | hit | hit sequence, if present | | +-------------------------+-----------------------------------+ | | hit_strand | strand, hit sequence strand | | +-------------------------+-----------------------------------+ | | query | query sequence, if present | | +-------------------------+-----------------------------------+ | | query_strand | strand, query sequence strand | +----------------+-------------------------+-----------------------------------+ In addition to the column mappings above, BlatIO also provides the following object attributes: +----------------+-------------------------+-----------------------------------+ | Object | Attribute | Value | +================+=========================+===================================+ | HSP | gapopen_num | Q gap count + T gap count, total | | | | number of gap openings | | +-------------------------+-----------------------------------+ | | ident_num | matches + repmatches, total | | | | number of identical residues | | +-------------------------+-----------------------------------+ | | ident_pct | percent identity, calculated | | | | using UCSC's formula | | +-------------------------+-----------------------------------+ | | query_is_protein | boolean, whether the query | | | | sequence is a protein | | +-------------------------+-----------------------------------+ | | score | HSP score, calculated using | | | | UCSC's formula | +----------------+-------------------------+-----------------------------------+ Finally, the default HSP and HSPFragment properties are also provided. See the HSP and HSPFragment documentation for more details on these properties. .. [#] You can obtain the number of blocks / fragments in the HSP by invoking ``len`` on the HSP N)log) SearchIndexer)Hit)HSP) HSPFragment) QueryResult) BlatPslParserBlatPslIndexer BlatPslWriterz^\d+\s+\d+\s+\d+\s+\d+cs4durdd|dDSfdd|dDS)aCTransform the given comma-separated string into a list (PRIVATE). :param csv_string: comma-separated input string :type csv_string: string :param caster: function used to cast each item in the input string to its intended type :type caster: callable, accepts string, returns object NcSg|]}|r|qSr .0xr r G/var/www/html/myenv/lib/python3.10/site-packages/Bio/SearchIO/BlatIO.py z"_list_from_csv..,csg|]}|r|qSr r r casterr rr)split) csv_stringrr rr_list_from_csvs rcsLt|t|krtdt|t|f|dkr|Sfddt||DS)aZReorients block starts into the opposite strand's coordinates (PRIVATE). :param starts: start coordinates :type starts: list [int] :param blksizes: block sizes :type blksizes: list [int] :param seqlen: sequence length :type seqlen: int :param strand: sequence strand :type strand: int, choice of -1, 0, or 1 z9Unequal start coordinates and block sizes list (%r vs %r)rcsg|] \}}||qSr r )rstartblksizeseqlenr rrz$_reorient_starts..)len RuntimeErrorzip)startsblksizesrstrandr rr_reorient_startss r%cCst|ddkr@|dddkr"|d|ddd|d dkS|ddd kr@|d |d |ddd|d dkSd S)z%Validate if psl is protein (PRIVATE).r$+tendtstarts blocksizes-tstarttsizeFr)pslr r r _is_proteins$ r3c Cs|rdnd}d}||d|d}|d|d}t||}|dkr%dS||}|dkr/dn|}||d|d |d }|dkr[d |d ||d tdtd||}|S) zCalculate millibad (PRIVATE).r,r'rqendqstartr)r/matches repmatches mismatchesi qnuminsert)minroundr) r2 is_proteinsize_mulmillibad qali_size tali_sizeali_sizesize_diftotalr r r_calc_millibads*    rDcCs@|rdnd}||d|dd?||d|d|dS)zCalculate score (PRIVATE).r,r'r6r7r8r9 tnuminsertr )r2r<r=r r r _calc_score#s  rFcst|}|r d}n |dddkrdnd}z|dddkr dnd}Wn ty.d}Ynw|r3dndt|d|d|d |}t|dd kr_t|d fd d |dD|d|}n|d }t|t|krvt|dksyJJtt|dd t||dD}tt|fdd t||dD} d|vrd|vrt|dt|dkrt|krt| ksJJn t|t| ksJg} t|D]L\} } |d} | sdn| | }|d}|sdn|| }t||||d}d|_ | d|_ | d|_ | | d|_ | | d|_ ||_||_| |qt| }|j |dks0J|j |dks:J|j |dksDJ|j |dksNJfdd |jD}||jkrh|dkskJJ|d|_|d|_|d|_|d|_|d|_|d |_|d!|_|d"|_|d|d|_|d|d!|_|d |d"|_||_d#t ||d$|_!t"|||_#t|dd k|_$|S)%z*Create high scoring pair object (PRIVATE).rr$r(r'r+r,qstartsr-qsizer&r*cg|]}|qSr r )riblocksize_multiplierr rrFrz_create_hsp..r0cSsg|]\}}||qSr r rryr r rrPrcsg|] \}}||qSr r rMrKr rrUrtseqsqseqs)hitqueryDNAr5r4r/r)csg|]}|qSr r )rspanrKr rrrr6r8r7ncountr9 qbaseinsertrE tbaseinsertgY@g?)%r3 IndexErrorr%rlistr! enumerategetr molecule_type query_start query_end hit_starthit_end query_strand hit_strandappendr hit_span_allquery_span_all match_num mismatch_num match_rep_numn_numquery_gapopen_num query_gap_numhit_gapopen_num hit_gap_num ident_num gapopen_numgap_numquery_is_proteinrD ident_pctrFscore_has_hit_strand)hidqidr2r<qstrandhstrandrGhstartsquery_range_all hit_range_allfragsidxqcoordshseqlisthseqqseqlistqseqfraghsp hit_spansr rKr _create_hsp/s  ,       &         rc@s:eZdZdZdddZddZddZd d Zd d Zd S)rzParser for the BLAT PSL format.FcCs||_|j|_||_dSInitialize the class.N)handlereadlinelinepslx)selfrrr r r__init__s  zBlatPslParser.__init__ccsd|jsdStt|js#|j|_|jsdStt|jr|D]}d|_|Vq'dS)z1Iterate over BlatPslParser, yields query results.Nblat) rresearch _RE_ROW_CHECKstriprr_parse_qresultprogram)rqresultr r r__iter__s  zBlatPslParser.__iter__cCs|jsJdd|jdD}||i}|d|d<t|d|d<|d|d <t|d |d <t|d |d <t|d|d<t|d|d<t|d|d<t|d|d<t|d|d<t|d|d<t|d|d<|d|d<t|d|d<t|d |d!<t|d"|d#<t|d$|d%<t|d&|d'<t|d(t|d)<t|d*t|d+<t|d,t|d-<|jrt|d.|d/<t|d0|d1<|S)2z6Return a dictionary of parsed column values (PRIVATE).cSr r r r r r rrrz,BlatPslParser._parse_row.. qname rH tnamer0rr6r'r8r&r7r,rVr9rWrErXr$ r5 r4r/r) blockcountr-rGr*rPrO)rrr_validate_colsintrr)rcolsr2r r r _parse_rows:     zBlatPslParser._parse_rowcCsR|jst|dkrtd|jt|fdSt|dkr'td|jt|fdS)z2Validate column's length of PSL or PSLX (PRIVATE).rzAInvalid PSL line: %r. Expected 21 tab-separated columns, found %izBInvalid PSLX line: %r. Expected 23 tab-separated columns, found %iN)rr ValueErrorr)rrr r rrs    zBlatPslParser._validate_colsccsPd}d}d}d}d}d}d}d\}} d\} } d\} } gg}} | dur+| } |} | } |jr;|} | d }| d } n|}d\}} | |krH|}n|}| | ksR||krU|}n|}| durt| | | }||||krzt|}| d |_||g}||ks||krt| d }|D]}||q| d |_|V||krdSg}|j |_q!)z$Yield QueryResult objects (PRIVATE).rr'r,r&rN)NNTrrr0)idrH) rrrrdrseq_lenrabsorbrr)r state_EOFstate_QRES_NEWstate_QRES_SAME state_HIT_NEWstate_HIT_SAME qres_state file_statecur_qidcur_hidprev_qidprev_hidcurprevhit_listhsp_list hit_staterrRrr r rrs`          zBlatPslParser._parse_qresultNF) __name__ __module__ __qualname____doc__rrrrrr r r rrs " rc@s.eZdZdZeZd ddZddZddZd S) r z"Indexer class for BLAT PSL output.FcCstj|||ddS)r)rN)rr)rfilenamerr r rr.szBlatPslIndexer.__init__c cs|j}|dd}d}d}|}|}tt|s3|}|}|s+dStt|r |}dd||D}|durM||}n||} | |krg| |||fV| }|t |}|}|sy| |||fVdSq4)zCIterate over the file handle; yields key, start offset, and length.rrN TcSr r r r r r rrIrz+BlatPslIndexer.__iter__..) _handleseektellrrr_RE_ROW_CHECK_IDXrrdecoder) rr query_id_idx qresult_keytab_char start_offsetr end_offsetrcurr_keyr r rr2s:   zBlatPslIndexer.__iter__c Cs~|j}||d}d}d}d} |}|s |Sdd||D}|dur/||}n ||} | |kr: |S||7}q)zFReturn raw bytes string of a QueryResult object from the given offset.rNrTcSr r r r r r rrfrz*BlatPslIndexer.get_raw..)rrrrr) roffsetrrr qresult_rawrrrrr r rget_rawYs(   zBlatPslIndexer.get_rawNr) rrrrr_parserrrrr r r rr )s   'r c@s2eZdZdZd ddZddZddZd d Zd S) r zWriter for the blat-psl format.FcCs||_||_||_dSr)rheaderr)rrrrr r rrus zBlatPslWriter.__init__cCs|j}d\}}}}|jr|||D],}|rA||||d7}|t|7}|tdd|D7}|tdd|D7}q||||fS)zWrite query results to file.)rrrrr'cs|]}t|VqdSNr1rrRr r r z+BlatPslWriter.write_file..cs|]}t|jVqdSr)r fragmentsrr r rr)rrwrite _build_header _build_rowrsum)rqresultsrqresult_counter hit_counter hsp_counter frag_counterrr r r write_file|s   zBlatPslWriter.write_filecCsd}|dd7}|S)z-Build header, tab-separated string (PRIVATE).zpsLayout version 3 z match mis- rep. N's Q gap Q gap T gap T gap strand Q Q Q Q T T T T block blockSizes qStarts tStarts match match count bases count bases name size start end name size start end count %s z---------------------------------------------------------------------------------------------------------------------------------------------------------------r )rrr r rrs zBlatPslWriter._build_headerc sg}|D];}|jD]4}t|dd}|rdndg}||j||j||j||j||j||j||j ||j fdd|j D}|j |kr^t d|j }|djdkrkd } nd } td d|jD|j |j|dj} |djdkrd} |jr| d 7} nd } | d 7} td d|jD|j |j| } || ||j||j||j||j||j||j||j||j|t||ddd|Dd|ddd| Dd|ddd| Dd|jr2|ddd|jDd|ddd|jDd|ddd|Dq qd|dS)zGReturn a string or one row or more of the QueryResult object (PRIVATE).rrFr,r'crIr r )rsrKr rrrz,BlatPslWriter._build_row..z0HSP hit span and query span values do not match.rr(r.cSg|]}|dqSrr r r r rrrr+cSrrr r r r rrrrcsrrstrr r r rrrz+BlatPslWriter._build_row..csrrrr r r rrrcsrrrr r r rrrcsrrrseqr r r rrrcsrrrr r r rrrrcsrrrr r r rrr )hspsgetattrrdrgrhrirjrkrlrmrnrfrerrbr%r{rrcrur|rr^r_r`rarjoinr query_allhit_all) rr qresult_linesrRrrrreff_query_spans block_sizesr$rGryrzr rKrrsv                        CzBlatPslWriter._build_rowN)FF)rrrrrrrrr r r rr rs   r __main__) run_doctestr)rrmathrBio.SearchIO._indexrBio.SearchIO._modelrrrr__all___PTR_ROW_CHECKcompilerencoderrr%r3rDrFrrr r r Bio._utilsrr r r rs64         hI~