bamCoverage with --extendReads

[Aarthi Mohan]

Hi All,

I am using bamCoverage to get the fragment pileup of my ChIP-seq data (RNAPol2). I made sure my BAM file has only paired-end reads, and used them to generate a bedgraph file using bamCoverage (--extendReads)

bamCoverage -b T2_Ser25_L001_trimmed_CP_pos_sorted.bam --binSize=1 -o T2_Ser25_L001_trimmed_CP_pos_sorted_frag.bed --extendReads --outFileFormat=bedgraph:

chr1    0       596     0.00

chr1    596     610     1.00

chr1    610     745     2.00

chr1    745     758     3.00

chr1    758     788     4.00

chr1    788     790     5.00

chr1    790     795     6.00

chr1    795     799     7.00

chr1    799     877     8.00

chr1    877     967     7.00

I also used bamToBed (with -bedpe option ) and extracted the columns 1,2,6 from the bed file to create a MACS2 BEDPE input. Following is the bedgraph file from 

bedtools genomecov:

chr1    0       745     0

chr1    745     758     1

chr1    758     788     2

chr1    788     790     3

chr1    790     795     4

chr1    795     799     5

chr1    799     863     6

chr1    863     877     5

chr1    877     1028    4

chr1    1028    1035    5

MACS2 treat pileup without SPMR option:

chr1    0       745     0.00000

chr1    745     758     1.00000

chr1    758     788     2.00000

chr1    788     790     3.00000

chr1    790     795     4.00000

chr1    795     799     5.00000

chr1    799     863     6.00000

chr1    863     877     5.00000

chr1    877     1028    4.00000

chr1    1028    1035    5.00000

I understand that 4th column in bamCoverage is different because it extends each mate individually and counts them. But, why is the pileup start site from bamCoverage different from bedtools and MACS2? Is there any way to actually obtain the extended mate pair information as bedfile from bamCoverage?

Thanks for your help!

Regards,

Aarthi

[Devon Ryan]

Hi Aarthi,

As to why there's a slight difference in start position early on we'd
have to see a reproducible example, though it's likely that this is
due deepTools applying the expected fragment length to discordant
alignments. There's no way to obtain extended mate information from
deepTools.

Devon
--
Devon Ryan, Ph.D.
Email: dpr...@dpryan.com
Data Manager/Bioinformatician
Max Planck Institute of Immunobiology and Epigenetics
Stübeweg 51
79108 Freiburg
Germany

> --
> You received this message because you are subscribed to the Google Groups
> "deepTools" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to deeptools+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Aarthi Mohan]

Hi Devon,

Thanks for the reply. I can provide you the BAM file for 'chr1' or all of them, whichever you may find useful. 

I just extracted the 'chr1' reads and used bamCoverage on it. Now, the pileup start sites is altered again. For a paired-end BAM file, isn't the extension done using the right mates' end location?

bamCoverage -b T2_Ser25_L001_trimmed_chr1_CP_pos_sorted.bam -o T2_Ser25_L001_trimmed_chr1_CP_pos_sorted_frag.bed --binSize=1 --outFileFormat=bedgraph --extendReads

chr1    0       606     0.00

chr1    606     620     1.00

chr1    620     745     2.00

chr1    745     758     3.00

chr1    758     788     4.00

chr1    788     790     5.00

chr1    790     795     6.00

chr1    795     799     7.00

chr1    799     863     8.00

chr1    863     873     7.00

Thanks for your help!

Aarthi

[Aarthi Mohan]

Sorry to repost on the same issue, but I thought this might be useful. 

When I create BEDPE file (for 'chr1'), all my fragment length stays below 1kb (i used -X 1000 in bowtie2 to mark discordant pairs). Also, I have filtered out all the discordant pairs from my BAM file (retaining only YT:Z:CP). But, the messages from bamCoverage indicate, it is seeing a maximum fragment length greater than 1kb. 

bamCoverage -b T2_Ser25_L001_trimmed_chr1_CP_pos_sorted.bam -o T2_Ser25_L001_trimmed_chr1_CP_pos_sorted_frag.bed --binSize=1 --outFileFormat=bedgraph --

extendReads
verbose: False
out_file_for_raw_data: None
numberOfSamples: None
bedFile: None
bamFilesList: ['T2_Ser25_L001_trimmed_chr1_CP_pos_sorted.bam']
ignoreDuplicates: False
numberOfProcessors: 32
samFlag_exclude: None
save_data: False
blackList: None
stepSize: 1
smoothLength: None
center_read: False
defaultFragmentLength: 257.0
chrsToSkip: []
region: None
maxPairedFragmentLength: 1028.0
samFlag_include: None
binLength: 1
blackListFileName: None
minMappingQuality: None
zerosToNans: False

[Devon Ryan]

Go ahead and just send me chr1. 

Sent from my iPhone

[Aarthi Mohan]

Hi Devon,

Here it is. Thanks for your help!

Regards,

Aarthi

[Devon Ryan]

Hi Aarthi,

This is looking increasingly light a slight bug that I'm surprised that we never caught (it's been around since December 4th). I'll finish confirming that this is indeed the case and release a bugfix (version 2.3.3) tomorrow if needed.

Devon

-- 
Devon Ryan, PhD
Bioinformatician / Data manager
Bioinformatics Core Facility
Max Planck Institute for Immunobiology and Epigenetics
Email: dpry...@gmail.com

[Aarthi Mohan]

I see, thanks for all the help Devon!

Cheers,

Aarthi

--
You received this message because you are subscribed to a topic in the Google Groups "deepTools" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/deeptools/gPo_qDvj7Ww/unsubscribe.
To unsubscribe from this group and all its topics, send an email to deeptools+...@googlegroups.com.

[Devon Ryan]

2.3.3 is now available via pypi and bioconda, hopefully this solves the issue.

Devon
--
Devon Ryan, Ph.D.
Email: dpr...@dpryan.com
Data Manager/Bioinformatician
Max Planck Institute of Immunobiology and Epigenetics
Stübeweg 51
79108 Freiburg
Germany