<s> """<STR_LIT>""" <EOL> from ruffus import * <EOL> import sys <EOL> import glob <EOL> import gzip <EOL> import os <EOL> import itertools <EOL> import CGAT . CSV as CSV <EOL> import re <EOL> import math <EOL> import types <EOL> import optparse <EOL> import shutil <EOL> import sqlite3 <EOL> import CGAT . GFF as GFF <EOL> import CGAT . GTF as GTF <EOL> import CGAT . Experiment as E <EOL> import CGAT . Pipeline as P <EOL> import CGAT . IOTools as IOTools <EOL> import CGAT . Genomics as Genomics <EOL> import CGAT . Database as Database <EOL> import CGAT . FastaIterator as FastaIterator <EOL> import PipelineGeneset as PGeneset <EOL> import PipelineGO as PGO <EOL> import scipy . stats <EOL> import CGAT . Stats as Stats <EOL> import alignlib <EOL> import CGAT . Mali as Mali <EOL> PARAMS = P . getParameters ( ) <EOL> @ files ( ( ( "<STR_LIT>" , "<STR_LIT>" ) , ) ) <EOL> def createProfile ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> outfile = open ( outfile , "<STR_LIT:w>" ) <EOL> for line in open ( infile ) : <EOL> if line . startswith ( "<STR_LIT:#>" ) : continue <EOL> data = re . split ( "<STR_LIT>" , line [ : - <NUM_LIT:1> ] , <NUM_LIT:1> ) <EOL> print data <EOL> pid , sequence = data <EOL> outfile . write ( "<STR_LIT>" % ( pid , sequence ) ) <EOL> outfile . close ( ) <EOL> def getParts ( src ) : <EOL> '''<STR_LIT>''' <EOL> result = None <EOL> r = [ ] <EOL> last_s = src . getColTo ( ) <EOL> for p in range ( src . getRowFrom ( ) , <EOL> src . getRowTo ( ) ) : <EOL> s = src . mapRowToCol ( p ) <EOL> if s < <NUM_LIT:0> : continue <EOL> if last_s >= s : <EOL> if result : <EOL> r . append ( result ) <EOL> result = alignlib . makeAlignmentVector ( ) <EOL> last_s = s <EOL> result . addPair ( s , p , <NUM_LIT:0> ) <EOL> if result : <EOL> r . append ( result ) <EOL> return r <EOL> def alignIndels ( all_alleles , colcounts , extend_by = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> aa = alignlib . makeAlignatorDPFull ( alignlib . ALIGNMENT_LOCAL , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> alignator = alignlib . makeMultipleAlignatorSimple ( aa ) <EOL> ids = all_alleles . keys ( ) <EOL> for x , c in enumerate ( colcounts ) : <EOL> if c <= <NUM_LIT:1> : continue <EOL> sequences = alignlib . StringVector ( ) <EOL> for sid in ids : <EOL> for allele in all_alleles [ sid ] : <EOL> sequences . append ( allele [ x ] ) <EOL> mali = alignlib . makeMultAlignment ( ) <EOL> alignator . align ( mali , sequences ) <EOL> realigned = [ ] <EOL> for line in str ( alignlib . MultAlignmentFormatPlain ( mali , sequences ) ) . split ( "<STR_LIT:\n>" ) [ : - <NUM_LIT:1> ] : <EOL> data = line [ : - <NUM_LIT:1> ] . split ( "<STR_LIT:\t>" ) <EOL> realigned . append ( data [ <NUM_LIT:1> ] ) <EOL> assert len ( realigned ) == len ( sequences ) <EOL> l = max ( [ len ( r ) for r in realigned ] ) <EOL> i = <NUM_LIT:0> <EOL> for sid in ids : <EOL> for allele in all_alleles [ sid ] : <EOL> if realigned [ i ] : allele [ x ] = realigned [ i ] <EOL> else : allele [ x ] = "<STR_LIT:->" * l <EOL> i += <NUM_LIT:1> <EOL> colcounts [ x ] = l <EOL> def _alignToProfile ( infile , outfile , <EOL> min_score = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> mali = Mali . Mali ( ) <EOL> mali . readFromFile ( open ( "<STR_LIT>" ) ) <EOL> src_mali = Mali . convertMali2Alignlib ( mali ) <EOL> E . debug ( "<STR_LIT>" % ( mali . getNumSequences ( ) , mali . getNumColumns ( ) ) ) <EOL> profile_mali = mali . getClone ( ) <EOL> n = profile_mali . getNumColumns ( ) <EOL> for x in "<STR_LIT>" : <EOL> for y in range ( <NUM_LIT:0> , <NUM_LIT:2> ) : <EOL> profile_mali . addSequence ( "<STR_LIT>" % ( x , y ) , <NUM_LIT:0> , n , x * n ) <EOL> profile_mali = Mali . convertMali2Alignlib ( profile_mali ) <EOL> alignlib . setDefaultEncoder ( alignlib . getEncoder ( alignlib . DNA4 ) ) <EOL> alignlib . setDefaultLogOddor ( alignlib . makeLogOddorUniform ( ) ) <EOL> profile = alignlib . makeProfile ( profile_mali ) <EOL> alignment_mode = alignlib . ALIGNMENT_WRAP <EOL> alignator = alignlib . makeAlignatorDPFull ( alignment_mode , <EOL> - <NUM_LIT> , <EOL> - <NUM_LIT:0.5> ) <EOL> map_seq2profile = alignlib . makeAlignmentVector ( ) <EOL> map_rseq2profile = alignlib . makeAlignmentVector ( ) <EOL> profile . prepare ( ) <EOL> build_mali = alignlib . makeMultAlignment ( ) <EOL> m = alignlib . makeAlignmentVector ( ) <EOL> m . addDiagonal ( <NUM_LIT:0> , n , <NUM_LIT:0> ) <EOL> build_mali . add ( src_mali , m ) <EOL> outf = open ( outfile , "<STR_LIT:w>" ) <EOL> outf_log = open ( outfile + "<STR_LIT>" , "<STR_LIT:w>" ) <EOL> outf_log . write ( "<STR_LIT>" ) <EOL> sequences , aa = alignlib . StringVector ( ) , alignlib . AlignandumVector ( ) <EOL> ids = [ ] <EOL> for pid in mali . getIdentifiers ( ) : <EOL> sequences . append ( re . sub ( "<STR_LIT:->" , "<STR_LIT>" , mali [ pid ] ) ) <EOL> ids . append ( pid ) <EOL> c = E . Counter ( ) <EOL> for s in FastaIterator . FastaIterator ( open ( infile ) ) : <EOL> E . debug ( "<STR_LIT>" % s . title ) <EOL> c . input += <NUM_LIT:1> <EOL> rsequence = Genomics . complement ( s . sequence ) <EOL> seq = alignlib . makeSequence ( s . sequence ) <EOL> rseq = alignlib . makeSequence ( rsequence ) <EOL> alignator . align ( map_seq2profile , seq , profile ) <EOL> alignator . align ( map_rseq2profile , rseq , profile ) <EOL> if map_seq2profile . getScore ( ) > map_rseq2profile . getScore ( ) : <EOL> m , seq , sequence = map_seq2profile , seq , s . sequence <EOL> else : <EOL> m , seq , sequence = map_rseq2profile , rseq , rsequence <EOL> if m . getLength ( ) == <NUM_LIT:0> : <EOL> c . skipped += <NUM_LIT:1> <EOL> continue <EOL> if m . getScore ( ) < min_score : <EOL> c . skipped += <NUM_LIT:1> <EOL> continue <EOL> r = getParts ( m ) <EOL> covered = <NUM_LIT:0> <EOL> for mm in r : <EOL> build_mali . add ( mm ) <EOL> sequences . append ( sequence ) <EOL> ids . append ( s . title ) <EOL> covered += mm . getLength ( ) - mm . getNumGaps ( ) <EOL> mali_covered = m . getColTo ( ) - m . getColFrom ( ) <EOL> outf_log . write ( "<STR_LIT:\t>" . join ( map ( str , ( <EOL> s . title , <EOL> len ( s . sequence ) , <EOL> m . getRowFrom ( ) , <EOL> m . getRowTo ( ) , <EOL> len ( r ) , <EOL> covered , <EOL> "<STR_LIT>" % ( <NUM_LIT> * covered / len ( s . sequence ) ) , <EOL> m . getScore ( ) , <EOL> m . getColFrom ( ) , <EOL> m . getColTo ( ) , <EOL> mali_covered , <EOL> "<STR_LIT>" % ( ( <NUM_LIT> * mali_covered ) / mali . getNumColumns ( ) ) <EOL> ) ) ) + "<STR_LIT:\n>" ) <EOL> c . output += <NUM_LIT:1> <EOL> result = str ( alignlib . MultAlignmentFormatPlain ( build_mali , <EOL> sequences , <EOL> alignlib . UnalignedStacked ) ) <EOL> for pid , data in zip ( ids , result . split ( "<STR_LIT:\n>" ) ) : <EOL> start , sequence , end = data . split ( "<STR_LIT:\t>" ) <EOL> outf . write ( "<STR_LIT>" % ( pid , int ( start ) + <NUM_LIT:1> , int ( end ) , sequence ) ) <EOL> outf . close ( ) <EOL> outf_log . close ( ) <EOL> E . info ( "<STR_LIT>" % str ( c ) ) <EOL> @ follows ( createProfile ) <EOL> @ files ( [ ( x , "<STR_LIT>" % ( x [ : - <NUM_LIT:3> ] , f ) , f ) <EOL> for x , f in itertools . product ( glob . glob ( "<STR_LIT>" ) , ( <NUM_LIT:0> , <NUM_LIT> ) ) ] ) <EOL> def alignToProfile ( infile , outfile , threshold ) : <EOL> _alignToProfile ( infile , outfile , min_score = threshold ) <EOL> @ transform ( alignToProfile <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def importMaliStats ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> table = P . toTable ( outfile ) + "<STR_LIT>" <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ transform ( alignToProfile <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def buildMalisCodon ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ transform ( alignToProfile <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def buildMalisAA ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ transform ( alignToProfile <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def computeColumnStatsNA ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ transform ( buildMalisAA <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def computeColumnStatsAA ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ transform ( ( computeColumnStatsNA , computeColumnStatsAA ) <EOL> , suffix ( "<STR_LIT>" ) <EOL> , "<STR_LIT>" ) <EOL> def importColumnStats ( infile , outfile ) : <EOL> '''<STR_LIT>''' <EOL> table = P . toTable ( outfile ) <EOL> statement = '''<STR_LIT>''' <EOL> P . run ( ) <EOL> @ follows ( alignToProfile <EOL> , importMaliStats <EOL> , computeColumnStatsNA <EOL> , buildMalisAA <EOL> , computeColumnStatsAA <EOL> , importColumnStats <EOL> ) <EOL> def full ( ) : pass <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> sys . exit ( P . main ( sys . argv ) ) </s>
