Metaseq: Average ChIP-seq signal over promoters error

Question

[Sudhir Jadhao]

Hello everyone,

I am new to python,

I am trying to plot chipseq signal over promoter using metaseq . The script available at metaseq site, It uses gtf file in the from of database. I don't have gtf file in the from of database, its giving error for that . Anyone have metaseq script for plotting chipseq signal over promoter which uses bam file data around 1.5GB (use the user defined data, not example data) and .gtf file

[Fidel]

Hi, I recommend you to use deepTools. You just need three commands and a bed file that you can download from UCSC:

$ bamCoverage -b file.bam -o file.bigwig
$ computeMatrix reference-point --scoreFile file.bigwig --regionsFile genes.bed -- beforeRegionStartLength 2000 --afterRegionStartLength 500 --referencePoint TSS -out matrix.tab.gz
$ profiler -m matrix.tab.gz -o profile.png

However, is not recommended to directly plot chip-seq signals as biases are frequently found. E.g. human and mouse promoters tend to have an enrichment of reads just because they are GC rich. Instead the log2ratio of ChIP-seq vs. input is preferable.

[endrebak]

deeptools seems great. Note: in the code above profiler has become plotProfile

[hlsz.laszlo]

Dear all,

How to create a meta-gene profile using metaseq? Unfortunately, I didn't find any help for this.

Answer 1

[Ryan Dale]

edit: For more standard cases like what you're asking about, try deepTools, ngsplot, or HOMER for a much simpler approach. Use metaseq when you need a lot of customization or are integrating with other tools.

metaseq's example 1 uses a database to show how to get TSSs. You can see this question and others for how to get TSSs or promoters using other methods, and just substitute that file for tsses.gtf in the example.

For a simpler example, see Plotting Reads Around Tss, which refers to the following script that you can modify for your data. Make sure to change the "processes" argument if you have multiple cores on your computer:

	# From http://www.biostars.org/p/83800/:
	# "What I want to do is to plot reads of my histone marks (in bam file)
	# around TSS with CpG and TSS without CpG (Essentially a coverage profile)."
	#
	#
	# To install metaseq and dependencies, see:
	#
	#
	# https://pythonhosted.org/metaseq/install.html
	#
	#
	# To download the example data used here, make sure you're in the directory
	# this script is saved in, and then use:
	#
	# git clone https://gist.github.com/a2e63a2fb93d05341de5.git demo_data
	#
	# (Or see https://gist.github.com/daler/a2e63a2fb93d05341de5 and download the
	# files individually)
	#
	import metaseq
	import pybedtools
	import numpy as np
	from matplotlib import pyplot as plt
	bam = metaseq.genomic_signal('demo_data/h3k4me3-chr21.bam', 'bam')
	cpg = pybedtools.BedTool('demo_data/cpg-chr21.bed.gz')
	tss = pybedtools.BedTool('demo_data/tss-chr21.bed.gz')
	# extend by 5 kb up/downstream
	tss = tss.slop(b=5000, genome='hg19')
	tss_with_cpg = tss.intersect(cpg, u=True)
	tss_without_cpg = tss.intersect(cpg, v=True)
	# change this to as many CPUs as you have in order to run in parallel
	processes = 1
	# each read will be extended 3' to a total size of this many bp
	fragment_size = 200
	# the region +/-5kb around each TSS will be split into a total of 100 bins,
	# change as needed
	bins = 100
	x = np.linspace(-5000, 5000, bins)
	# most of the work happens here
	y1 = bam.array(tss_with_cpg, bins=bins, processes=processes, fragment_size=fragment_size)
	y2 = bam.array(tss_without_cpg, bins=bins, processes=processes, fragment_size=fragment_size)
	plt.rcParams['font.size'] = 11
	plt.rcParams['font.family'] = 'Arial'
	plt.plot(x, y1.mean(axis=0), label='with cpg', color='k')
	plt.plot(x, y2.mean(axis=0), label='without cpg', color='r', linestyle='--')
	plt.legend(loc='best')
	plt.xlabel('Distance from TSS (bp)')
	plt.ylabel('Mean H3K4me3 read density')
	plt.show()

view raw metaseq_demo.py hosted with ❤ by GitHub

Answer 2

[Ming Tommy Tang]

See my post here.

[Sudhir Jadhao]

Thanks for reply..

I tried with HTSeq but its giving below error

Traceback (most recent call last):
  File "streaming_throughRead.py", line 17, in <module>
    tsspos[window]+=p
  File "_HTSeq.pyx", line 526, in HTSeq._HTSeq.GenomicArray.__getitem__ (src/_HTSeq.c:9489)
  File "_HTSeq.pyx", line 372, in HTSeq._HTSeq.ChromVector.__getitem__ (src/_HTSeq.c:6649)
IndexError: start too small

Answer 3

[chemcehn]

Just try another program, e.g. HOMER, it is much simpler to use