Reproduce genome coverage of ENCODE-CSHL Long RNA-seq data, but failed?

[josephyan]

I've try to reproduce the graphic view (genome-wide coverage) of the ENCODE-CSHL long RNA-seq http://hgdownload.cse.ucsc.edu/golde...shlLongRnaSeq/

Let's say, K562-Chromatin-Rep4 as an example. http://hgdownload.cse.ucsc.edu/golde...talAlnRep4.bam

Luckily, it provides the bigwig file as well.Thus I've an idea to test whether I can reproduce the same data visualization.

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeCshlLongRnaSeq/wgEncodeCshlLongRnaSeqK562ChromatinTotalPlusRawSigRep4.bigWig

1. bigWig provided by paper
I upload the Plus strand bigwig into UCSC, with the very default options:

track name='K562-Chromatin4-Pos' type="bigWig" color=0,0,255 bigDataUrl=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeCshlLongRnaSeq/wgEncodeCshlLongRnaSeqK562ChromatinTotalPlusRawSigRep4.bigWig

2. genomeCoverage generate bigWig
genomeCoverage tool is one of the BEDtools toolbox.

samtools sort wgEncodeCshlLongRnaSeqK562ChromatinTotalAlnRep4.bam wgEncodeCshlLongRnaSeqK562ChromatinTotalAlnRep4.bam.sort

genomeCoverage -bg -ibam wgEncodeCshlLongRnaSeqK562ChromatinTotalAlnRep4.bam.sort -g hg19.chromInfo -strand + >K562chromatin.POS.bedgraph

genomeCoverageBed -bg -ibam wgEncodeCshlLongRnaSeqK562ChromatinTotalAlnRep4.bam.sorted.bam -g ../hg19.chromInfo -strand -|awk '{$4 = - $4;print $0}' >K562ChromatinTotalAlnRep4.REV.bedgraph

Then using bedGraphToBigwig (UCSC exe) to convert bedgraph into bigwig:

bedGraphToBigwig K562chromatin.POS.bedgraph hg19.chromInfo K562chromatin.POS.bigwig

Finally upload them into Galaxy and then send to UCSC genome browser.


Thus, in the image, from top to bottom:
Black: Forward_strand_my_version
Blue: Forward_strand_ENCODE_default
Red: Reverse_strand_my_version
Grey: Reverse_strand_ENCODE_version

The result is strange: General expression wave is similar, they share same peak. But my-version-bigwig give higher basal level (notice the y-axis).

And please note the gene I give an example here, Xist. It seems that if I my self calculate the coverage, the Xist intron is NOT spliced off the primary mRNA (indicated by green box).

Any suggestions?

[schlesin]

CoverageBed probably does not consider the split reads (RNA-seq reads mapping across splice junctions) correctly and instead 'fills-in' in corresponding alignment positions. You will see them when you look in the CIGAR string in the BAM files.

Some bedtools have a -split option to deal with this.

Felix