def create_otu_table(output_fp, deblurred_list,
                     outputfasta_fp=None, minreads=0):
    """Create a biom table out of all files in a directory

    Parameters
    ----------
    output_fp : string
        filepath to output BIOM table
    deblurred_list : list of (str, str)
        list of file names (including path), sampleid of all deblurred
        fasta files to add to the table
    outputfasta_fp : str, optional
        name of output fasta file (of all sequences in the table) or None
        to not write
    minreads : int, optional
        minimal number of reads per bacterial sequence in order to write
        it to the biom table and fasta file or 0 to write all
    """
    logger = logging.getLogger(__name__)
    logger.info('create_otu_table for %d samples, '
                'into output table %s' % (len(deblurred_list), output_fp))

    # the regexp for finding the number of reads of a sequence
    sizeregexp = re.compile('(?<=size=)\w+')
    seqdict = {}
    seqlist = []
    sampset = set()
    samplist = []
    # arbitrary size for the sparse results matrix so we won't run out of space
    obs = scipy.sparse.dok_matrix((int(1E9), len(deblurred_list)), dtype=np.double)

    # load the sequences from all samples into a sprase matrix
    sneaking_extensions = {'fasta', 'fastq', 'fna', 'fq', 'fa'}
    for (cfilename, csampleid) in deblurred_list:
        if csampleid.rsplit('.', 1)[-1] in sneaking_extensions:
            csampleid = csampleid.rsplit('.', 1)[0]

        # test if sample has already been processed
        if csampleid in sampset:
            warnings.warn('sample %s already in table!', UserWarning)
            logger.error('sample %s already in table!' % csampleid)
            continue
        sampset.add(csampleid)
        samplist.append(csampleid)
        csampidx = len(sampset)-1
        # read the fasta file and add to the matrix
        for chead, cseq in sequence_generator(cfilename):
            cseq = cseq.upper()
            if cseq not in seqdict:
                seqdict[cseq] = len(seqlist)
                seqlist.append(cseq)
            cseqidx = seqdict[cseq]
            cfreq = float(sizeregexp.search(chead).group(0))
            try:
                obs[cseqidx, csampidx] += cfreq
            except IndexError:
                # exception means we ran out of space - add more OTUs
                shape = obs.shape
                obs.resize((shape[0]*2,  shape[1]))
                obs[cseqidx, csampidx] = cfreq

    logger.info('for output biom table loaded %d samples, %d unique sequences'
                % (len(samplist), len(seqlist)))

    # and now make the sparse matrix the real size
    obs.resize((len(seqlist), len(samplist)))

    # do the minimal reads per otu filtering
    if minreads > 0:
        readsperotu = obs.sum(axis=1)
        keep = np.where(readsperotu >= minreads)[0]
        logger.info('keeping %d (out of %d sequences) with >=%d reads' %
                    (len(keep), len(seqlist), minreads))
        obs = obs[keep, :]
        seqlist = list(np.array(seqlist)[keep])
        logger.debug('filtering completed')

    # convert the matrix to a biom table
    table = Table(obs.tocsr(), seqlist, samplist,
                  observation_metadata=None,
                  sample_metadata=None, table_id=None,
                  generated_by="deblur",
                  create_date=datetime.now().isoformat())
    logger.debug('converted to biom table')

    # remove samples with 0 reads
    filter_minreads_samples_from_table(table)

    # save the merged otu table
    write_biom_table(table, output_fp)
    logger.info('saved to biom file %s' % output_fp)

    # and save the fasta file
    if outputfasta_fp is not None:
        logger.debug('saving fasta file')
        with open(outputfasta_fp, 'w') as f:
            for cseq in seqlist:
                f.write('>%s\n%s\n' % (cseq, cseq))
        logger.info('saved sequence fasta file to %s' % outputfasta_fp)
def create_otu_table ( output_fp , deblurred_list , outputfasta_fp = None , minreads = 0 ) : logger = logging . getLogger ( __name__ ) logger . info ( 'create_otu_table for %d samples, ' 'into output table %s' % ( len ( deblurred_list ) , output_fp ) ) sizeregexp = re . compile ( '(?<=size=)\w+' ) seqdict = { } seqlist = [ ] sampset = set ( ) samplist = [ ] obs = scipy . sparse . dok_matrix ( ( int ( 1E9 ) , len ( deblurred_list ) ) , dtype = np . double ) sneaking_extensions = { 'fasta' , 'fastq' , 'fna' , 'fq' , 'fa' } for ( cfilename , csampleid ) in deblurred_list : if csampleid . rsplit ( '.' , 1 ) [ - 1 ] in sneaking_extensions : csampleid = csampleid . rsplit ( '.' , 1 ) [ 0 ] if csampleid in sampset : warnings . warn ( 'sample %s already in table!' , UserWarning ) logger . error ( 'sample %s already in table!' % csampleid ) continue sampset . add ( csampleid ) samplist . append ( csampleid ) csampidx = len ( sampset ) - 1 for chead , cseq in sequence_generator ( cfilename ) : cseq = cseq . upper ( ) if cseq not in seqdict : seqdict [ cseq ] = len ( seqlist ) seqlist . append ( cseq ) cseqidx = seqdict [ cseq ] cfreq = float ( sizeregexp . search ( chead ) . group ( 0 ) ) try : obs [ cseqidx , csampidx ] += cfreq except IndexError : shape = obs . shape obs . resize ( ( shape [ 0 ] * 2 , shape [ 1 ] ) ) obs [ cseqidx , csampidx ] = cfreq logger . info ( 'for output biom table loaded %d samples, %d unique sequences' % ( len ( samplist ) , len ( seqlist ) ) ) obs . resize ( ( len ( seqlist ) , len ( samplist ) ) ) if minreads > 0 : readsperotu = obs . sum ( axis = 1 ) keep = np . where ( readsperotu >= minreads ) [ 0 ] logger . info ( 'keeping %d (out of %d sequences) with >=%d reads' % ( len ( keep ) , len ( seqlist ) , minreads ) ) obs = obs [ keep , : ] seqlist = list ( np . array ( seqlist ) [ keep ] ) logger . debug ( 'filtering completed' ) table = Table ( obs . tocsr ( ) , seqlist , samplist , observation_metadata = None , sample_metadata = None , table_id = None , generated_by = "deblur" , create_date = datetime . now ( ) . isoformat ( ) ) logger . debug ( 'converted to biom table' ) filter_minreads_samples_from_table ( table ) write_biom_table ( table , output_fp ) logger . info ( 'saved to biom file %s' % output_fp ) if outputfasta_fp is not None : logger . debug ( 'saving fasta file' ) with open ( outputfasta_fp , 'w' ) as f : for cseq in seqlist : f . write ( '>%s\n%s\n' % ( cseq , cseq ) ) logger . info ( 'saved sequence fasta file to %s' % outputfasta_fp )
Create a biom table out of all files in a directory

    Parameters
    ----------
    output_fp : string
        filepath to output BIOM table
    deblurred_list : list of (str, str)
        list of file names (including path), sampleid of all deblurred
        fasta files to add to the table
    outputfasta_fp : str, optional
        name of output fasta file (of all sequences in the table) or None
        to not write
    minreads : int, optional
        minimal number of reads per bacterial sequence in order to write
        it to the biom table and fasta file or 0 to write all