Python HTSeq.mixed_SAM_alignments_with_buffer方法代码示例

# 需要导入模块: import HTSeq [as 别名]
# 或者: from HTSeq import mixed_SAM_alignments_with_buffer [as 别名]
def count_reads_in_features( sam_filename, gff_filename, samtype, order, stranded,
        overlap_mode, feature_type, id_attribute, quiet, minaqual, out,
        file_regex, sample_col, debug ):


    # Keep track of what features we see and
    # also the counts from those features
    features = HTSeq.GenomicArrayOfSets( "auto", stranded != "no" )
    counts = {}
    lens = {}

    # Try to open samfile to fail early in case it is not there
    if sam_filename != "-":
        open( sam_filename ).close()

    # OPen and read the GFF file (features)
    gff = HTSeq.GFF_Reader( gff_filename )
    i = 0
    try:
        for f in gff:
            if f.type == feature_type:
                try:
                     # Grab out the feature id so we can make our count dict
                    feature_id = f.attr[ id_attribute ]
                except KeyError:
                    raise ValueError, ( "Feature %s does not contain a '%s' attribute" %
                        ( f.name, id_attribute ) )
                if stranded != "no" and f.iv.strand == ".":
                    raise ValueError, ( "Feature %s at %s does not have strand information but you are "
                        "running htseq-count in stranded mode. Use '--stranded=no'." %
                        ( f.name, f.iv ) )
                features[ f.iv ] += feature_id
                counts[ f.attr[ id_attribute ] ] = 0
                lens[f.attr[id_attribute]] = f.iv.length
            i += 1
            # progress report
            if i % 100000 == 0 and not quiet:
                sys.stderr.write( "%d GFF lines processed.\n" % i )
    except:
        sys.stderr.write( "Error occured when processing GFF file (%s):\n" % gff.get_line_number_string() )
        raise

    if not quiet:
        sys.stderr.write( "%d GFF lines processed.\n" % i )

    if len( counts ) == 0:
        sys.stderr.write( "Warning: No features of type '%s' found.\n" % feature_type )

    # now find the alignment file
    if samtype == "sam":
        SAM_or_BAM_Reader = HTSeq.SAM_Reader
    elif samtype == "bam":
        SAM_or_BAM_Reader = HTSeq.BAM_Reader
    else:
        raise ValueError, "Unknown input format %s specified." % samtype

    try:
        if sam_filename != "-":
            read_seq_file = SAM_or_BAM_Reader( sam_filename )
            read_seq = read_seq_file
            first_read = iter(read_seq).next()
        else:
            read_seq_file = SAM_or_BAM_Reader( sys.stdin )
            read_seq_iter = iter( read_seq_file )
            first_read = read_seq_iter.next()
            # oops! have to put that first read back on
            read_seq = itertools.chain( [ first_read ], read_seq_iter )
        #pe_mode = first_read.paired_end
    except:
        sys.stderr.write( "Error occured when reading beginning of SAM/BAM file.\n" )
        raise

    try:
        # Do a mixture of paired and single end reads
        read_seq = HTSeq.mixed_SAM_alignments_with_buffer( read_seq )

        # Successful Counts
        paired_end = 0
        single_end = 0
        collapsed = 0
        collapsed_trimmed = 0

        # Total counts for FPKM
        total_counts = 0

        empty = 0
        ambiguous = 0
        notaligned = 0
        lowqual = 0
        nonunique = 0
        i = 0
        for r in read_seq:
            i += 1 # This handles reporting at read 0
            if i % 100000 == 0 and not quiet:
                sys.stderr.write( "%d SAM alignment records processed.\n" % ( i,) )
                if debug:
                    break
            # Figure out if single or paired end

            # SINGLE END
#.........这里部分代码省略.........