I'm wondering if anyone out there can help me. I am trying to calculate mean coverage over exons and introns in my RNAseq dataset using bedtools coverage -mean. I have very large bam files that are actually combined files from two sequencing runs- one single end and one paired-end. I am having a heck of a time making bed files from these bam files. I don't want to loose the second read in a pair (as I understand I would if I just made a regular bed from bamtobed), since the second read is legitimate coverage in my regions. However, if I try to use the -bedpe option the program freaks out because not all the reads in my file are paired (maybe less than half). For this analysis I basically just need a split (on the N cigar) bedfile with intervals for all reads that contribute to coverage in the dataset with the correct strand in the 6th column (which gets complicated when R2 is actually backwards).
Thanks in advance for any suggestions!