Hi guys,
For visualization purpose (like in a genome browser), I would like to get the per-base genome-wide coverage of paired-end RNA-seq data. I thought this would be easy, turns out it is not. Here is a picture showing what I want to do :
And here is where I came so far :
accepted_hits.bam)And it is where I'm stuck. I thought of using bedtools bamtobed with the -split option along with -bedpe, but those two options seems not compatible. A possible solution would be to extract the halfs of the reads spanning introns that are not linked to the mate pair, save them as singleton and add them later to the coverage, but I have no idea how to do that. Any suggestion, even for a completely different method? Thank you!
samtools sort -no <(samtools view -b -f 2 fwd.bam) temp | bamToBed -i stdin -bedpe | cut -f 1,2,6 | sort -k 1,1 > fwd_merged_mates.bed
bedtools genomecov -bga -i fwd_merged_mates.bed -g chrom_summary > WT1_fwd.bg
•
written
by
Carlo Yague
QUESTION: How to get per base coverage in bigWig format from mapped paired-end reads with splice junctions (bam files such as tophat2 output).
ISSUES: The bedtools suite can do the trick with either paired end reads or single end with splice junctions, but not with both paired end and splice junctions.
SOLUTION: First split your reads according to their strand (see this tutorial). Then simply use samtools depth to compute coverage:
samtools depth fwd.bam > fwd_temp
samtools depth rev.bam > rev_temp
The output will be a three column <chrom> <position> <value> text file. However to be able to load it easily into a genome browser, I want my coverage in bigWig format. To do that I first manually* convert the output into a 4 column <chrom> <start> <end> <value> bedgraph file with [start]=[position -1] and [end] = [position]. Then:
bedGraphToBigWig fwd_temp chrom_summary fwd.bw
bedGraphToBigWig rev_temp chrom_summary rev.bw
That's it, you got your bw file!
Thanks for all the help to answer this one and see ya!
*Edit: When I say manually, I mean with R, python or whatever.
•
written
by
Carlo Yague
Exactly! I'll look into STAR RNA. Will update.
•
written
by
Carlo Yague
Log in to answer this question.
Your goal is get genome coverage from only exons ( and stranded ) ?
No, I want genome coverage for the whole genome. But if a read if split by an intron (mapping with tophat allows it), I need to take that into account when I compute the coverage, i.e, I can't just count +1 between the start and the end of a fragment if one of the read forming the fragment is split by an intron... Hope this clarify my question :)
I guess you simply needs number of reads per gene/Exon. You could use HTSeq-count for counting the number of reads per gene, which takes care of split reads without double counting them. If you need exon level counts, DEXSeq provides some scripts that will flatten your GTF file and can be used with HTSeq-count to get exon level information.
No, I thank you for your interest but this is not what I want. Instead I want to generate coverage tracks that cover the whole genome for visualization in a genome browser. We could call it "stranded per base count" if u want, in contrast with exon/gene count. I edited the question to make this point clear.
For me is still unclear what you are trying to do and frankly this:
Sounds like a bad idea. There is a rule that says: if it is hard to explain it is probably a bad idea.
IMHO If all you want is coverage then first compute the whole genome coverage ignoring genes for each strand then slice and dice that file with bedtools.
Ok there seems to be some misunderstanding. When I say "intron", I'm not talking about actual intron, but about reads that are mapped discontinuously by tophat (although the two concepts are linked). And indeed, all I want is the whole coverage ignoring genes. But because of those discontinuity, I find it hard to compute it :
There isn't really a good solution for the case you mention, since you can't know with certainty that there's no splicing going on between your mates (yes, you can just assume that the GTF defines all legitimate splicing events, but that is an unjustifiable assumption in many situations).
Very nice point.
Currently I'm trying samtools depth. I need little more testing, but it looks like it does what I want.