bamCoverage with DNAse-seq shift

Question

[JMR]

Hi,

I am trying to create a bigwig file using the bamCoverage function for DNAse-seq data and I wanted to check if my implementation is correct. For ATAC-seq I am doing this:

      alignmentSieve \
        --ATACshift \
        --bam file.bam \
        -o file.shifted.bam

and then

      bamCoverage \
        -b file.shifted.bam \
        -o file.bw \
        --binSize 1 \
        --normalizeUsing CPM \
        --outFileFormat bigwig

I think for DNAse-seq I will first have to do the shifts on a BED file and then call bamCoverage. I was thinking doing something like:

bedtools bamtobed -i file.bam \
  | awk 'BEGIN{OFS="\t"}{
      if ($6=="+") print $1,$2,$2+1,".",1,"+";
      else         print $1,$3,$3+1,".",1,"-";
    }' \
  | bedtools bedtobam -g file.chrom.sizes -i - > file.shifted.bam

then

samtools sort -o file.shifted.sorted.bam file.shifted.bam
samtools index file.shifted.sorted.bam

and lastly the same bamCoverage call

      bamCoverage \
        -b file.file.shifted.sorted.bam \
        -o file.bw \
        --binSize 1 \
        --normalizeUsing CPM \
        --outFileFormat bigwig

Am I missing something? I haven't found an awk command for DNAse-seq but I think the above should work? As an additional note all my data is paired end.

Also, this is the only issue I found with a similar(ish) question. Is using that offset flag equivalent to what I am doing?

Thanks in advance.

Answer 1

[vivekbhr]

If you do not need the shifted BAM file for any other analysis, you can simply use the --offset option to create the desired bigwig file.

[JMR]

I haven't used deeptools too much but based on the description of the function I think if I add --Offset 1 will that be equivalent? Thanks for the help!

[vivekbhr]

Yes indeed if you just need the start of your DNAse protected fragment, use --offset 1