<s> import os <EOL> import sys <EOL> import argparse <EOL> import itertools <EOL> import math <EOL> import time <EOL> import shlex <EOL> from fastqp import FastqReader , padbases , percentile , mean , bam_read_count , gc , window <EOL> from fastqp . plots import qualplot , qualdist , qualmap , nucplot , depthplot , gcplot , gcdist , kmerplot , mismatchplot , adaptermerplot <EOL> from fastqp . adapters import all_adapter_sequences <EOL> from collections import defaultdict <EOL> from simplesam import Reader , Sam <EOL> from subprocess import Popen , PIPE <EOL> from scipy import stats <EOL> from operator import mul <EOL> from six . moves import reduce <EOL> def run ( args ) : <EOL> """<STR_LIT>""" <EOL> time_start = time . time ( ) <EOL> if args . input . name != '<STR_LIT>' : <EOL> bsize = os . path . getsize ( args . input . name ) <EOL> est_counter = int ( ) <EOL> sample_lengths = list ( ) <EOL> sample_binsizes = list ( ) <EOL> act_nlines = int ( ) <EOL> name , ext = os . path . splitext ( args . input . name ) <EOL> if ( args . leftlimit > <NUM_LIT:0> ) and ( args . rightlimit > <NUM_LIT:0> ) : <EOL> if args . rightlimit < args . leftlimit : <EOL> sys . exit ( "<STR_LIT>" ) <EOL> if args . type : <EOL> ext = '<STR_LIT:.>' + args . type <EOL> if ext not in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] and args . input . name != '<STR_LIT>' : <EOL> sys . exit ( "<STR_LIT>" ) <EOL> if args . name : <EOL> sample_name = args . name <EOL> else : <EOL> sample_name = args . input . name <EOL> if ext in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> with FastqReader ( open ( args . input . name ) ) as fh : <EOL> for read in fh : <EOL> sample_lengths . append ( len ( read ) ) <EOL> sample_binsizes . append ( len ( str ( read ) ) ) <EOL> est_counter += <NUM_LIT:1> <EOL> if est_counter == <NUM_LIT> : <EOL> break <EOL> mean_bentry = mean ( sample_binsizes ) <EOL> mean_len = mean ( sample_lengths ) <EOL> est_nlines = int ( bsize / mean_bentry ) <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" <EOL> "<STR_LIT>" . format ( bytes = mean_bentry , <EOL> len = mean_len , <EOL> est = est_nlines ) ) <EOL> elif ext == '<STR_LIT>' : <EOL> with Reader ( open ( args . input . name ) ) as fh : <EOL> for read in fh : <EOL> sample_lengths . append ( len ( read ) ) <EOL> sample_binsizes . append ( len ( str ( read ) ) ) <EOL> est_counter += <NUM_LIT:1> <EOL> if est_counter == <NUM_LIT> : <EOL> break <EOL> mean_bentry = mean ( sample_binsizes ) <EOL> mean_len = mean ( sample_lengths ) <EOL> est_nlines = int ( bsize / mean_bentry ) <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" <EOL> "<STR_LIT>" . format ( bytes = mean_bentry , <EOL> len = mean_len , <EOL> est = est_nlines ) ) <EOL> elif ext == '<STR_LIT>' : <EOL> est_nlines = sum ( bam_read_count ( args . input . name ) ) <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( est = est_nlines ) ) <EOL> elif ext == '<STR_LIT>' : <EOL> if args . binsize : <EOL> n = args . binsize <EOL> est_nlines = None <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( binsize = n ) ) <EOL> else : <EOL> sys . stderr . write ( "<STR_LIT>" ) <EOL> p1 = Popen ( shlex . split ( '<STR_LIT>' % args . input . name ) , stdout = PIPE ) <EOL> p2 = Popen ( shlex . split ( '<STR_LIT>' ) , stdin = p1 . stdout , stdout = PIPE ) <EOL> est_nlines , _ = p2 . communicate ( ) <EOL> est_nlines = int ( est_nlines ) // <NUM_LIT:4> <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( est = est_nlines ) ) <EOL> elif name == '<STR_LIT>' : <EOL> if args . binsize : <EOL> n = args . binsize <EOL> else : <EOL> n = <NUM_LIT:1> <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( binsize = n ) ) <EOL> est_nlines = None <EOL> if est_nlines is not None : <EOL> if args . binsize : <EOL> n = args . binsize <EOL> else : <EOL> nf = math . floor ( est_nlines / args . nreads ) <EOL> if nf >= <NUM_LIT:1> : <EOL> n = int ( nf ) <EOL> else : <EOL> n = <NUM_LIT:1> <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( binsize = n ) ) <EOL> if ext in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> infile = Reader ( args . input ) <EOL> else : <EOL> infile = FastqReader ( args . input , ext = ext ) <EOL> read_len = defaultdict ( int ) <EOL> cycle_nuc = defaultdict ( lambda : defaultdict ( int ) ) <EOL> cycle_qual = defaultdict ( lambda : defaultdict ( int ) ) <EOL> cycle_gc = defaultdict ( int ) <EOL> cycle_kmers = defaultdict ( lambda : defaultdict ( int ) ) <EOL> cycle_mismatch = { '<STR_LIT:C>' : defaultdict ( lambda : defaultdict ( int ) ) , <EOL> '<STR_LIT>' : defaultdict ( lambda : defaultdict ( int ) ) , <EOL> '<STR_LIT:A>' : defaultdict ( lambda : defaultdict ( int ) ) , <EOL> '<STR_LIT:T>' : defaultdict ( lambda : defaultdict ( int ) ) } <EOL> if args . count_duplicates : <EOL> try : <EOL> from pybloom import ScalableBloomFilter <EOL> bloom_filter = ScalableBloomFilter ( mode = ScalableBloomFilter . SMALL_SET_GROWTH ) <EOL> except ImportError : <EOL> sys . exit ( "<STR_LIT>" ) <EOL> duplicates = <NUM_LIT:0> <EOL> percent_complete = <NUM_LIT:10> <EOL> reads = infile . subsample ( n ) <EOL> for read in reads : <EOL> if isinstance ( read , Sam ) : <EOL> if args . aligned_only and not read . mapped : <EOL> continue <EOL> elif args . unaligned_only and read . mapped : <EOL> continue <EOL> if read . reverse : <EOL> seq = read . seq [ : : - <NUM_LIT:1> ] <EOL> qual = read . qual [ : : - <NUM_LIT:1> ] <EOL> else : <EOL> seq = read . seq <EOL> qual = read . qual <EOL> else : <EOL> seq = read . seq <EOL> qual = read . qual <EOL> if ( args . leftlimit == <NUM_LIT:1> ) and ( args . rightlimit < <NUM_LIT:0> ) : <EOL> pass <EOL> elif ( args . leftlimit >= <NUM_LIT:1> ) and ( args . rightlimit > <NUM_LIT:0> ) : <EOL> try : <EOL> seq = seq [ args . leftlimit - <NUM_LIT:1> : args . rightlimit ] <EOL> qual = qual [ args . leftlimit - <NUM_LIT:1> : args . rightlimit ] <EOL> except IndexError : <EOL> act_nlines += n <EOL> continue <EOL> elif ( args . leftlimit > <NUM_LIT:1> ) and ( args . rightlimit < <NUM_LIT:0> ) : <EOL> try : <EOL> seq = seq [ args . leftlimit - <NUM_LIT:1> : ] <EOL> qual = qual [ args . leftlimit - <NUM_LIT:1> : ] <EOL> except IndexError : <EOL> act_nlines += n <EOL> continue <EOL> if len ( seq ) == <NUM_LIT:0> : <EOL> act_nlines += n <EOL> continue <EOL> cycle_gc [ gc ( seq ) ] += <NUM_LIT:1> <EOL> if args . count_duplicates : <EOL> if seq in bloom_filter : <EOL> duplicates += <NUM_LIT:1> <EOL> else : <EOL> bloom_filter . add ( seq ) <EOL> for i , ( s , q ) in enumerate ( zip ( seq , qual ) ) : <EOL> cycle_nuc [ args . leftlimit + i ] [ s ] += <NUM_LIT:1> <EOL> cycle_qual [ args . leftlimit + i ] [ q ] += <NUM_LIT:1> <EOL> read_len [ len ( qual ) ] += <NUM_LIT:1> <EOL> for i , kmer in enumerate ( window ( seq , n = args . kmer ) ) : <EOL> cycle_kmers [ args . leftlimit + i ] [ kmer ] += <NUM_LIT:1> <EOL> if isinstance ( read , Sam ) and read . mapped : <EOL> try : <EOL> ref = read . parse_md ( ) <EOL> for i , ( s , r ) in enumerate ( zip ( seq , ref ) ) : <EOL> if s != r : <EOL> try : <EOL> cycle_mismatch [ r ] [ args . leftlimit + i ] [ s ] += <NUM_LIT:1> <EOL> except KeyError : <EOL> pass <EOL> except KeyError : <EOL> pass <EOL> if est_nlines is not None : <EOL> if ( act_nlines / est_nlines ) * <NUM_LIT:100> >= percent_complete : <EOL> sys . stderr . write ( "<STR_LIT>" <EOL> "<STR_LIT>" . format ( percent_complete , <EOL> act_nlines , <EOL> time . strftime ( '<STR_LIT>' , <EOL> time . gmtime ( time . time ( ) - time_start ) ) ) ) <EOL> percent_complete += <NUM_LIT:10> <EOL> act_nlines += n <EOL> positions = [ k for k in sorted ( cycle_qual . keys ( ) ) ] <EOL> depths = [ read_len [ k ] for k in sorted ( read_len . keys ( ) ) ] <EOL> basecalls = [ cycle_nuc [ k ] . keys ( ) for k in sorted ( cycle_nuc . keys ( ) ) ] <EOL> bases = set ( list ( itertools . chain . from_iterable ( basecalls ) ) ) <EOL> map ( padbases ( bases ) , cycle_nuc . values ( ) ) <EOL> quantile_values = [ <NUM_LIT> , <NUM_LIT> , <NUM_LIT:0.5> , <NUM_LIT> , <NUM_LIT> ] <EOL> quantiles = [ ] <EOL> for _ , v in sorted ( cycle_qual . items ( ) ) : <EOL> for q in tuple ( v . keys ( ) ) : <EOL> v [ ord ( str ( q ) ) - <NUM_LIT> ] = v . pop ( q ) <EOL> line = [ percentile ( v , p ) for p in quantile_values ] <EOL> quantiles . append ( line ) <EOL> adapter_kmers = set ( ) <EOL> for adapter in all_adapter_sequences : <EOL> for kmer in window ( adapter , n = args . kmer ) : <EOL> adapter_kmers . add ( kmer ) <EOL> observed_expected = dict ( ) <EOL> all_kmers = [ cycle_kmers [ k ] . keys ( ) for k in sorted ( cycle_kmers . keys ( ) ) ] <EOL> kmers = set ( list ( itertools . chain . from_iterable ( all_kmers ) ) ) <EOL> bad_kmers = [ ] <EOL> sequenced_bases = sum ( ( l * n for l , n in read_len . items ( ) ) ) <EOL> priors = tuple ( map ( float , args . base_probs . split ( '<STR_LIT:U+002C>' ) ) ) <EOL> for kmer in kmers : <EOL> kmer_counts = [ ( i , cycle_kmers [ i ] [ kmer ] ) for i in sorted ( cycle_kmers . keys ( ) ) ] <EOL> expected_fraction = reduce ( mul , ( p ** kmer . count ( b ) for b , p in zip ( ( '<STR_LIT:A>' , '<STR_LIT:T>' , '<STR_LIT:C>' , '<STR_LIT>' , '<STR_LIT:N>' ) , priors ) ) , <NUM_LIT:1> ) <EOL> expected = expected_fraction * sequenced_bases <EOL> observed_expected [ kmer ] = sum ( ( n for _ , n in kmer_counts ) ) / expected <EOL> slope , _ , _ , p_value , _ = stats . linregress ( * zip ( * kmer_counts ) ) <EOL> if abs ( slope ) > <NUM_LIT:2> and p_value < <NUM_LIT> : <EOL> bad_kmers . append ( ( kmer , slope , p_value ) ) <EOL> bad_kmers = sorted ( bad_kmers , key = lambda x : x [ <NUM_LIT:2> ] ) [ : <NUM_LIT:10> ] <EOL> pos_gc = [ sum ( [ cycle_nuc [ i ] [ '<STR_LIT:C>' ] , cycle_nuc [ i ] [ '<STR_LIT>' ] ] ) / sum ( [ cycle_nuc [ i ] [ '<STR_LIT:C>' ] , <EOL> cycle_nuc [ i ] [ '<STR_LIT>' ] , <EOL> cycle_nuc [ i ] [ '<STR_LIT:A>' ] , <EOL> cycle_nuc [ i ] [ '<STR_LIT:T>' ] ] ) * <NUM_LIT:100> for i in positions ] <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , column = '<STR_LIT>' , pos = '<STR_LIT:None>' , value = act_nlines ) ) <EOL> for cycle , count in read_len . items ( ) : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , column = '<STR_LIT>' , pos = cycle , <EOL> value = count ) ) <EOL> for i , position in enumerate ( positions ) : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = quantiles [ i ] [ <NUM_LIT:0> ] ) ) <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = quantiles [ i ] [ <NUM_LIT:1> ] ) ) <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = quantiles [ i ] [ <NUM_LIT:2> ] ) ) <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = quantiles [ i ] [ <NUM_LIT:3> ] ) ) <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = quantiles [ i ] [ <NUM_LIT:4> ] ) ) <EOL> for base in bases : <EOL> for position in positions : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = base , pos = position , <EOL> value = cycle_nuc [ position ] [ base ] ) ) <EOL> for position in positions : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = position , <EOL> value = cycle_gc [ position ] ) ) <EOL> for i in range ( <NUM_LIT> ) : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = '<STR_LIT>' , pos = i , <EOL> value = cycle_gc [ i ] ) ) <EOL> for kmer , obs_exp in sorted ( observed_expected . items ( ) , key = lambda x : x [ <NUM_LIT:1> ] ) : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , <EOL> column = kmer , pos = '<STR_LIT:None>' , <EOL> value = obs_exp ) ) <EOL> if args . count_duplicates : <EOL> sys . stdout . write ( "<STR_LIT>" . format ( row = sample_name , column = '<STR_LIT>' , pos = '<STR_LIT:None>' , value = duplicates / act_nlines ) ) <EOL> from zipfile import ZipFile <EOL> with ZipFile ( args . output + '<STR_LIT>' , mode = '<STR_LIT:w>' ) as zip_archive : <EOL> fig_kw = { '<STR_LIT>' : ( <NUM_LIT:8> , <NUM_LIT:6> ) } <EOL> qualplot ( positions , quantiles , zip_archive , fig_kw ) <EOL> median_qual = qualdist ( cycle_qual . values ( ) , zip_archive , fig_kw ) <EOL> qualmap ( cycle_qual , zip_archive , fig_kw ) <EOL> depthplot ( read_len , zip_archive , fig_kw ) <EOL> gcplot ( positions , pos_gc , zip_archive , fig_kw ) <EOL> gcdist ( cycle_gc , zip_archive , fig_kw ) <EOL> nucplot ( positions , bases , cycle_nuc , zip_archive , fig_kw ) <EOL> kmerplot ( positions , cycle_kmers , zip_archive , [ fields [ <NUM_LIT:0> ] for fields in bad_kmers ] , fig_kw ) <EOL> adaptermerplot ( positions , cycle_kmers , adapter_kmers , zip_archive , fig_kw ) <EOL> if isinstance ( infile , Reader ) : <EOL> mismatchplot ( positions <EOL> , cycle_mismatch , zip_archive , fig_kw ) <EOL> time_finish = time . time ( ) <EOL> elapsed = time_finish - time_start <EOL> if not args . quiet : <EOL> sys . stderr . write ( "<STR_LIT>" . format ( counts = act_nlines , <EOL> sec = time . strftime ( '<STR_LIT>' , <EOL> time . gmtime ( elapsed ) ) <EOL> ) ) <EOL> if len ( bad_kmers ) > <NUM_LIT:0> : <EOL> for kmer in bad_kmers : <EOL> sys . stderr . write ( "<STR_LIT>" % ( kmer ) ) <EOL> if median_qual < args . median_qual : <EOL> sys . stderr . write ( "<STR_LIT>" % median_qual ) <EOL> def main ( ) : <EOL> parser = argparse . ArgumentParser ( prog = '<STR_LIT>' , description = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT:input>' , type = argparse . FileType ( '<STR_LIT:r>' ) , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , action = "<STR_LIT:store_true>" , default = False , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , help = '<STR_LIT>' ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = str , help = '<STR_LIT>' ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , default = <NUM_LIT> , help = '<STR_LIT>' ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = str , default = '<STR_LIT>' , help = '<STR_LIT>' ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , default = <NUM_LIT:5> , choices = range ( <NUM_LIT:2> , <NUM_LIT:8> ) , help = '<STR_LIT>' ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = str , default = '<STR_LIT>' , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = str , default = None , choices = [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , default = <NUM_LIT:1> , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , default = - <NUM_LIT:1> , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , type = int , default = <NUM_LIT:30> , help = "<STR_LIT>" ) <EOL> align_group = parser . add_mutually_exclusive_group ( ) <EOL> align_group . add_argument ( '<STR_LIT>' , action = "<STR_LIT:store_true>" , default = False , help = "<STR_LIT>" ) <EOL> align_group . add_argument ( '<STR_LIT>' , action = "<STR_LIT:store_true>" , default = False , help = "<STR_LIT>" ) <EOL> parser . add_argument ( '<STR_LIT>' , '<STR_LIT>' , action = "<STR_LIT:store_true>" , default = False , help = "<STR_LIT>" ) <EOL> args = parser . parse_args ( ) <EOL> run ( args ) <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> main ( ) </s>
