def sample_cleanup(data, sample):
    """ stats, cleanup, and link to samples """

    ## get maxlen and depths array from clusters
    maxlens, depths = get_quick_depths(data, sample)

    try:
        depths.max()
    except ValueError:
        ## If depths is an empty array max() will raise
        print("    no clusters found for {}".format(sample.name))
        return

    ## Test if depths is non-empty, but just full of zeros.
    if depths.max():
        ## store which min was used to calculate hidepth here
        sample.stats_dfs.s3["hidepth_min"] = data.paramsdict["mindepth_majrule"]

        ## If our longest sequence is longer than the current max_fragment_length
        ## then update max_fragment_length. For assurance we require that
        ## max len is 4 greater than maxlen, to allow for pair separators.
        hidepths = depths >= data.paramsdict["mindepth_majrule"]
        maxlens = maxlens[hidepths]

        ## Handle the case where there are no hidepth clusters
        if maxlens.any():
            maxlen = int(maxlens.mean() + (2.*maxlens.std()))
        else:
            maxlen = 0
        if maxlen > data._hackersonly["max_fragment_length"]:
            data._hackersonly["max_fragment_length"] = maxlen + 4

        ## make sense of stats
        keepmj = depths[depths >= data.paramsdict["mindepth_majrule"]]
        keepstat = depths[depths >= data.paramsdict["mindepth_statistical"]]

        ## sample summary stat assignments
        sample.stats["state"] = 3
        sample.stats["clusters_total"] = depths.shape[0]
        sample.stats["clusters_hidepth"] = keepmj.shape[0]

        ## store depths histogram as a dict. Limit to first 25 bins
        bars, bins = np.histogram(depths, bins=range(1, 26))
        sample.depths = {int(i):v for i, v in zip(bins, bars) if v}

        ## sample stat assignments
        ## Trap numpy warnings ("mean of empty slice") printed by samples
        ## with few reads.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=RuntimeWarning)
            sample.stats_dfs.s3["merged_pairs"] = sample.stats.reads_merged
            sample.stats_dfs.s3["clusters_total"] = depths.shape[0]
            try:
                sample.stats_dfs.s3["clusters_hidepth"] = int(sample.stats["clusters_hidepth"])
            except ValueError:
                ## Handle clusters_hidepth == NaN
                sample.stats_dfs.s3["clusters_hidepth"] = 0
            sample.stats_dfs.s3["avg_depth_total"] = depths.mean()
            LOGGER.debug("total depth {}".format(sample.stats_dfs.s3["avg_depth_total"]))
            sample.stats_dfs.s3["avg_depth_mj"] = keepmj.mean()
            LOGGER.debug("mj depth {}".format(sample.stats_dfs.s3["avg_depth_mj"]))
            sample.stats_dfs.s3["avg_depth_stat"] = keepstat.mean()
            sample.stats_dfs.s3["sd_depth_total"] = depths.std()
            sample.stats_dfs.s3["sd_depth_mj"] = keepmj.std()
            sample.stats_dfs.s3["sd_depth_stat"] = keepstat.std()

    else:
        print("    no clusters found for {}".format(sample.name))

    ## Get some stats from the bam files
    ## This is moderately hackish. samtools flagstat returns
    ## the number of reads in the bam file as the first element
    ## of the first line, this call makes this assumption.
    if not data.paramsdict["assembly_method"] == "denovo":
        refmap_stats(data, sample)

    log_level = logging.getLevelName(LOGGER.getEffectiveLevel())

    if not log_level == "DEBUG":
        ## Clean up loose files only if not in DEBUG
        ##- edits/*derep, utemp, *utemp.sort, *htemp, *clust.gz
        derepfile = os.path.join(data.dirs.edits, sample.name+"_derep.fastq")
        mergefile = os.path.join(data.dirs.edits, sample.name+"_merged_.fastq")
        uhandle = os.path.join(data.dirs.clusts, sample.name+".utemp")
        usort = os.path.join(data.dirs.clusts, sample.name+".utemp.sort")
        hhandle = os.path.join(data.dirs.clusts, sample.name+".htemp")
        clusters = os.path.join(data.dirs.clusts, sample.name+".clust.gz")

        for f in [derepfile, mergefile, uhandle, usort, hhandle, clusters]:
            try:
                os.remove(f)
            except:
                pass
def sample_cleanup(data, sample):
    

    
    maxlens, depths = get_quick_depths(data, sample)

    try:
        depths.max()
    except ValueError:
        
        print("    no clusters found for {}".format(sample.name))
        return

    
    if depths.max():
        
        sample.stats_dfs.s3["hidepth_min"] = data.paramsdict["mindepth_majrule"]

        
        
        
        hidepths = depths >= data.paramsdict["mindepth_majrule"]
        maxlens = maxlens[hidepths]

        
        if maxlens.any():
            maxlen = int(maxlens.mean() + (2.*maxlens.std()))
        else:
            maxlen = 0
        if maxlen > data._hackersonly["max_fragment_length"]:
            data._hackersonly["max_fragment_length"] = maxlen + 4

        
        keepmj = depths[depths >= data.paramsdict["mindepth_majrule"]]
        keepstat = depths[depths >= data.paramsdict["mindepth_statistical"]]

        
        sample.stats["state"] = 3
        sample.stats["clusters_total"] = depths.shape[0]
        sample.stats["clusters_hidepth"] = keepmj.shape[0]

        
        bars, bins = np.histogram(depths, bins=range(1, 26))
        sample.depths = {int(i):v for i, v in zip(bins, bars) if v}

        
        
        
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=RuntimeWarning)
            sample.stats_dfs.s3["merged_pairs"] = sample.stats.reads_merged
            sample.stats_dfs.s3["clusters_total"] = depths.shape[0]
            try:
                sample.stats_dfs.s3["clusters_hidepth"] = int(sample.stats["clusters_hidepth"])
            except ValueError:
                
                sample.stats_dfs.s3["clusters_hidepth"] = 0
            sample.stats_dfs.s3["avg_depth_total"] = depths.mean()
            LOGGER.debug("total depth {}".format(sample.stats_dfs.s3["avg_depth_total"]))
            sample.stats_dfs.s3["avg_depth_mj"] = keepmj.mean()
            LOGGER.debug("mj depth {}".format(sample.stats_dfs.s3["avg_depth_mj"]))
            sample.stats_dfs.s3["avg_depth_stat"] = keepstat.mean()
            sample.stats_dfs.s3["sd_depth_total"] = depths.std()
            sample.stats_dfs.s3["sd_depth_mj"] = keepmj.std()
            sample.stats_dfs.s3["sd_depth_stat"] = keepstat.std()

    else:
        print("    no clusters found for {}".format(sample.name))

    
    
    
    
    if not data.paramsdict["assembly_method"] == "denovo":
        refmap_stats(data, sample)

    log_level = logging.getLevelName(LOGGER.getEffectiveLevel())

    if not log_level == "DEBUG":
        
        
        derepfile = os.path.join(data.dirs.edits, sample.name+"_derep.fastq")
        mergefile = os.path.join(data.dirs.edits, sample.name+"_merged_.fastq")
        uhandle = os.path.join(data.dirs.clusts, sample.name+".utemp")
        usort = os.path.join(data.dirs.clusts, sample.name+".utemp.sort")
        hhandle = os.path.join(data.dirs.clusts, sample.name+".htemp")
        clusters = os.path.join(data.dirs.clusts, sample.name+".clust.gz")

        for f in [derepfile, mergefile, uhandle, usort, hhandle, clusters]:
            try:
                os.remove(f)
            except:
                pass
def sample_cleanup ( data , sample ) : ## get maxlen and depths array from clusters maxlens , depths = get_quick_depths ( data , sample ) try : depths . max ( ) except ValueError : ## If depths is an empty array max() will raise print ( "    no clusters found for {}" . format ( sample . name ) ) return ## Test if depths is non-empty, but just full of zeros. if depths . max ( ) : ## store which min was used to calculate hidepth here sample . stats_dfs . s3 [ "hidepth_min" ] = data . paramsdict [ "mindepth_majrule" ] ## If our longest sequence is longer than the current max_fragment_length ## then update max_fragment_length. For assurance we require that ## max len is 4 greater than maxlen, to allow for pair separators. hidepths = depths >= data . paramsdict [ "mindepth_majrule" ] maxlens = maxlens [ hidepths ] ## Handle the case where there are no hidepth clusters if maxlens . any ( ) : maxlen = int ( maxlens . mean ( ) + ( 2. * maxlens . std ( ) ) ) else : maxlen = 0 if maxlen > data . _hackersonly [ "max_fragment_length" ] : data . _hackersonly [ "max_fragment_length" ] = maxlen + 4 ## make sense of stats keepmj = depths [ depths >= data . paramsdict [ "mindepth_majrule" ] ] keepstat = depths [ depths >= data . paramsdict [ "mindepth_statistical" ] ] ## sample summary stat assignments sample . stats [ "state" ] = 3 sample . stats [ "clusters_total" ] = depths . shape [ 0 ] sample . stats [ "clusters_hidepth" ] = keepmj . shape [ 0 ] ## store depths histogram as a dict. Limit to first 25 bins bars , bins = np . histogram ( depths , bins = range ( 1 , 26 ) ) sample . depths = { int ( i ) : v for i , v in zip ( bins , bars ) if v } ## sample stat assignments ## Trap numpy warnings ("mean of empty slice") printed by samples ## with few reads. with warnings . catch_warnings ( ) : warnings . simplefilter ( "ignore" , category = RuntimeWarning ) sample . stats_dfs . s3 [ "merged_pairs" ] = sample . stats . reads_merged sample . stats_dfs . s3 [ "clusters_total" ] = depths . shape [ 0 ] try : sample . stats_dfs . s3 [ "clusters_hidepth" ] = int ( sample . stats [ "clusters_hidepth" ] ) except ValueError : ## Handle clusters_hidepth == NaN sample . stats_dfs . s3 [ "clusters_hidepth" ] = 0 sample . stats_dfs . s3 [ "avg_depth_total" ] = depths . mean ( ) LOGGER . debug ( "total depth {}" . format ( sample . stats_dfs . s3 [ "avg_depth_total" ] ) ) sample . stats_dfs . s3 [ "avg_depth_mj" ] = keepmj . mean ( ) LOGGER . debug ( "mj depth {}" . format ( sample . stats_dfs . s3 [ "avg_depth_mj" ] ) ) sample . stats_dfs . s3 [ "avg_depth_stat" ] = keepstat . mean ( ) sample . stats_dfs . s3 [ "sd_depth_total" ] = depths . std ( ) sample . stats_dfs . s3 [ "sd_depth_mj" ] = keepmj . std ( ) sample . stats_dfs . s3 [ "sd_depth_stat" ] = keepstat . std ( ) else : print ( "    no clusters found for {}" . format ( sample . name ) ) ## Get some stats from the bam files ## This is moderately hackish. samtools flagstat returns ## the number of reads in the bam file as the first element ## of the first line, this call makes this assumption. if not data . paramsdict [ "assembly_method" ] == "denovo" : refmap_stats ( data , sample ) log_level = logging . getLevelName ( LOGGER . getEffectiveLevel ( ) ) if not log_level == "DEBUG" : ## Clean up loose files only if not in DEBUG ##- edits/*derep, utemp, *utemp.sort, *htemp, *clust.gz derepfile = os . path . join ( data . dirs . edits , sample . name + "_derep.fastq" ) mergefile = os . path . join ( data . dirs . edits , sample . name + "_merged_.fastq" ) uhandle = os . path . join ( data . dirs . clusts , sample . name + ".utemp" ) usort = os . path . join ( data . dirs . clusts , sample . name + ".utemp.sort" ) hhandle = os . path . join ( data . dirs . clusts , sample . name + ".htemp" ) clusters = os . path . join ( data . dirs . clusts , sample . name + ".clust.gz" ) for f in [ derepfile , mergefile , uhandle , usort , hhandle , clusters ] : try : os . remove ( f ) except : pass
stats, cleanup, and link to samples
stats cleanup and link to samples