I am trying to make a UCSC Genome Browser hub using bigWig files. My starting file format are BAM files. Briefly, my script does as follows:
- Scale the BAM file according to a specificic pre-determined scaling factor (using samtools).
- From scaled BAM, make a bedgraph representing coverage (using genomeCoverageBed)
- From beGgraph file, make bigwig file (using bedGraphToBigwig)
My script looks as follows:
samtools view -s 0.75 -b file.bam > file.scaled.bam; genomeCoverageBed -bg -split -ibam file.scaled.bam; -g hg38.chrom.sizes > file.bedgraph; sort -k1,1 -k2,2n file.bedgraph > file.sorted.bedgraph; bedGraphToBigWig file.sorted.bedgraph hg38.chrom.sizes file.bw
I am having an issue displaying my bigwig file and noticed this is most likely due to my bedGraph file not having a "chr" in front of each chromosome name (so I have "1" instead of "chr1"). I am confused as how this happened.
Initially, when first testing genomeCoverageBed, I kept getting the following warning:
*****WARNING: Genome (-g) files are ignored when BAM input is provided.
I therefore removed the -g argument in the genomeCoverageBed. However, I just noticed (I should have noticed before) that the "chr" is missing in my bedgraph file. Presumably, this is what is causing my bigwig file to then be ill-formatted.
I am confused, am I missing something? Why is genomeCoverageBed not outputting a bedGraph file with "chr" in it, like it should be doing?