Hi everyone! I'm performing some ChIP seq experiments and I want to read depth coverage over the region -500 to 3000pb of the TSS of all Ensemble protein coding transcript.
I have different approach, but as I am totally naive, I appreciate very much your advises.
a)First of all I've normalized the bam files using bamCoverage (deepTools), the output is a bedgraph, then I will do the anntotaion peak over TSS using ChIPseeker (R)
b)Also, I can do the peakannot over the bed file (after bam file conversion) and I will normalize the peakannot when I will plot it.
Other questions is, Is necessary make the peakcalling to calculate the read depth over TSS?
I'm looking forward to hear from you
Than you!!!
2 answers
a) looks good to me.
b) I don't really understand what you want to achieve here. What is 'doing a peakannot' ? Do you refer to the annotatePeak function of the ChIPseeker R package ?
It is not necessary to make a peakcalling to calculate read depth over TSS. You only need genome coverage and TSS annotation.
There are many many tools to do that. Since you are already familiar with deeptools, you could try their computeMatrix tool.
Could it be that your ChIP-seq data is from a female and that there is no Y ? This is a warning, not an error, and it could be ok to proceed if the non-Y coordinates are considered correct by the script.
More generally, check that the chromosome names in the bed file are the same as in the bedgraph file.
Carlo
a normalization equal to 25 million reads
Do you mean that you want to normalize as if you had 25 million reads mapped ? If so, you need to use a third option : --scaleFactor set to 25000000/number_of_reads_mapped. To calculate the number of reads mapped, you can just use samtools flagstat on the .bam file.
25000000/31567890 = 0.7919 : this is your scaleFactor for this particular bam file :)
last question, when I've tried to run reads <- import.bed(con="NGS974Norm")
NGSNorm is the bedgraph file, 4 cols and size ~ 400MB
cov <- coverage(reads) peakAnno <- annotatePeak(reads, tssRegion=c(-500, 2000), TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene, annoDb="org.Hs.eg.db")
I have the next error Error in rsqlite_send_query(conn@ptr, statement) : out of memory
any ideas?
Log in to answer this question.
-
best way to do a boxplot for abundnace matrix (microbiome data)
written by Lila M
-
How to calculate read depth for a specific gene list?
written by Lila M
-
how can I get the coverage over introns only in a RNA-seq experiment?
written by Lila M
-
What is the definition of "read depth" vs "coverage"? (again...)
written by ariel.balter
-
Select TSS to know if genes are differentially transcribed
written by Lila M
-
out of memory in R after run ChIPSeeker
written by Lila M
-
The best option for ChIP seq normalization.
written by Lila M
-
Is there any option in annotatePEak (chipseeker) to normalize the reads?
written by Lila M
-
MNase-seq / fragment size selection
written by Constantine
-
Deeptools sample scaling
written by Constantine
