Extracting counts from RNA-Seq: trouble with bedtools + multicov

Question

[Kristin Muench]

Hello,

I just aligned my first piece of RNA-Seq data and have a .bam file. I'd like to make a file of counts corresponding to genes that I can use in DESeq2. I downloaded a list of all genes from the UCSC gene browser and saved it as a .bam file. I'm trying to do the following:

bedtools multicov -bams UCSC_knownGene.bam -bed myData_out/accepted_hits.bam

But then I get this error:

Could not open input BAM files.

I've only read about this error on forums where people are using multiple bam files, so I'm not sure why I'm getting it if I'm only submitting one file that was directly downloaded from UCSC Genome Browser. Has anyone had a similar issue? Is there a better way to get .bam files to use for this?

Thank you!

[Varun Gupta]

Your list of genes should be in bed format.

https://genome.ucsc.edu/FAQ/FAQformat.html#format1

Also is your bam file indexed? Use samtools index to index your bam file

[Kristin Muench]

Not sure it was indexed - thank you!

[h.mon]

Try:

file gene_list.bam

it will probably say "ASCII text". If it a text file, then do:

head gene_list.bam

If unsure, post here the result of head. Bam files should be gzipped. I suspect the file you downloaded gene list is either bed (as Varun Gupta pointed out) or gtf format.

Answer 1

[Devon Ryan]

You immediate problem is likely solved by:

bedtools multicov -bams myData_out/accepted_hits.bam -bed UCSC_knownGene.bam

I'm assuming you just misnamed the knownGene file and that it's really a bed12 file.

Your mediate problem is that you shouldn't do any of this to begin with. Bedtools multicov should not be used to get counts for RNAseq. The resulting counts violate the assumptions made by tools like DESeq2 (and edgeR and limma and many more). Use featureCounts or htseq-count and download a GTF file from UCSC.