Skip to content

bamCoverage with DNAse-seq shift #1415

@JavierMenRev

Description

@JavierMenRev

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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions