I am new to bioinformatics and python,Now I have downloaded nucleosome-seq data and TSS-seq(which also names CAGE) data,so I try to draw a picture like this:
The difference was I have the TSS data,so I don't need to subtract the TSS feature from Human GTF annotation. But now I face the problems---the question I want to ask is:
I don't know how to write python script to convent my bed file to instead of GTF.
my TSS data are bed format like this:chrom,chromStart,strand.PS:the read length is 27bp;
here is my script，could anybody tell me why I am fault?
import HTSeq import pysam from matplotlib import pyplot bamfile = HTSeq.BAM_Reader("DRR000367.bam" ) bedfile = HTSeq.BED_Reader("Poly.bed") coverage = HTSeq.GenomicArray( "auto", stranded=False, typecode="i" ) for almnt in bamfile: if almnt.aligned: almnt.iv.length = fragmentsize coverage[ almnt.iv ] += 1 list(bedfile) profile = numpy.zeros( 2*halfwinwidth, dtype='i' ) for p in tsspos: window = HTSeq.GenomicInterval( p.chrom, p.pos - halfwinwidth, p.pos + halfwinwidth, "." ) wincvg = numpy.fromiter( coverage[window], dtype='i', count=2*halfwinwidth ) if p.strand == "+": profile += wincvg else: profile += wincvg[::-1] pyplot.plot( numpy.arange( -halfwinwidth, halfwinwidth ), profile ) pyplot.show()
Thanks very much.Any who give me any information would be pretty appreciate.