Genrich for single-end ATAC-seq data

[amitjavilaventura]

Hello,

I am trying to run Genrich for single-cell ATAC-seq data but doesn't work.

I am aligning the fastq files to the mm9 with Bowtie and them I am sorting the bam files by query name (as required for Genrich) with samtools sort. Then I am using the sorted bam files as input for Genrich but it is returning Error! Experimental sample has no analyzable fragments.

The code I used is that one:

samtools sort -n -m 5G -@ 5 -o sample_sorted.bam sample_.bam

Genrich -j \
-t sample_sorted.bam \
-o sample_.narrowPeak \
-p 0.0005 \
-y

Can anybody give me some advice?

Thank you very much.

[Maarten-vd-Sande]

Make sure to always check other issues (see #6 and #8)

[jsh58]

Thanks for the question. Could you please post the output when you run the command in verbose mode (-v)?

[amitjavilaventura]

Yes, thank you. The commands I am running are this ones:

samtools sort -n -@ 5 -m 5G -o bamfiles/${sample}/${sample}_sorted_name.bam bamfiles/${sample}/${sample}_.bam #genrich needs a bamfile sorted by query name as an input.

Genrich -j \
-t bamfiles/${sample}/${sample}_sorted_name.bam \
-o ${trial}/results/genrich/${sample}/${sample}_.narrowPeak \
-p 0.0005 \
-y \
-r \
-v

It is returning this:

[bam_sort_core] merging from 0 files and 5 in-memory blocks...
Processing experimental file #0: bamfiles/WTplusWNT/WTplusWNT_sorted_name.bam
  BAM records analyzed:   21551493
    Unmapped:             11235897
    Paired alignments:           0
    Unpaired alignments:  10315596
  PCR duplicates --
    Paired aln sets:             0
      duplicates:                0 (0.0%)
    Discordant aln sets:         0
      duplicates:                0 (0.0%)
    Singleton aln sets:          0
      duplicates:                0 (0.0%)
  Fragments analyzed:            0
    Full fragments:              0
    Half fragments:              0
    ATAC-seq cut sites:          0
      (expanded to length 100bp)
Error! Experimental sample has no analyzable fragments

I think that the problem is the sorted bam, because samtools sort is returning this: [bam_sort_core] merging from 0 files and 5 in-memory blocks... When normally samtools sort returns: merging from X files... with X being different than 0.

[jsh58]

Yes, that is strange, and likely a problem. I suggest you eliminate -@ 5 from samtools sort. Also, depending on your version of samtools, you may need to specify -O or -T.

[amitjavilaventura]

Still doesn't work. It is strange, because with PE files samtools sort and Genrich work well.

[jsh58]

Does samtools think it is producing a valid BAM file? Try running this:

$ samtools view -c <BAM>

[amitjavilaventura]

I have run this with my bam file and it returned the following number 21551493. I do not know what does it mean.

[jsh58]

If there were no error/warning messages, it means that samtools does think it is a valid BAM file. And the count matches what Genrich said:

  BAM records analyzed:   21551493

So I don't know why Genrich did not analyze any fragments, since you ran it with -y. Can you post the first few lines of the BAM file? samtools view -F0x4 <BAM> | head -n4

[amitjavilaventura]

Sure, I can.

These are the first lines, which I obtained by running what you told me samtools view -F0x4 bamfiles/WTplusWNT/WTplusWNT_.bam | head -n4:

A00302:59:H53NCDRXX:2:1101:13313:1251   0   chr13   48456896    255 51M *   0   0   ATTGCCATGCCCAAATGTACATGCCCTATCGTCATGCCCAAATGTACATGC FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:0  MD:Z:51 NM:i:0  XM:i:2
A00302:59:H53NCDRXX:2:1101:13476:1251   0   chrM    2266    255 51M *   0   0   GAATAAAAAATCCTCCGAATGATTATAACCTAGACTTACAAGTCAAAGTAA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:0  MD:Z:51 NM:i:0  XM:i:2
A00302:59:H53NCDRXX:2:1101:16984:1251   0   chrM    15058   255 51M *   0   0   CTTTCCTTCATACCTCAAAGCAACGAAGCCTAATATTCCGCCCAATCACAC :FFFFFF:FFFF:FFFFF:FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF XA:i:0  MD:Z:51 NM:i:0  XM:i:2
A00302:59:H53NCDRXX:2:1101:17074:1251   16  chr9    33925314    255 51M *   0   0   CTGCACTGTGTAATAGCACTTTCATAGAAATGAATATAGTCTGTACAAAGG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:0  MD:Z:51 NM:i:0  XM:i:2  

Thank you.

[jsh58]

Thanks, that's very helpful! The problem was a bug when singleton alignments do not have alignment scores (AS) in the BAM records. The bug has now been fixed (hopefully). Please let me know if you have any other problems.

Note that if you want to take advantage of Genrich's ability to analyze multimapping reads, you will need to use an aligner that reports alignment scores with AS (such as Bowtie2).

[amitjavilaventura]

I have redownloaded Genrich and I have ran the same command as the last time and it returned the same:

Processing experimental file #0: bamfiles/WTplusWNT/WTplusWNT_sorted_name.bam
  BAM records analyzed:    21551493
    Unmapped:              11235897
    Paired alignments:            0
    Unpaired alignments:   10315596
  PCR duplicates --
    Paired aln sets:              0
      duplicates:                 0 (0.0%)
    Discordant aln sets:          0
      duplicates:                 0 (0.0%)
    Singleton aln sets:           0
      duplicates:                 0 (0.0%)
  Fragments analyzed:             0
    Full fragments:               0
    Half fragments:               0
    ATAC-seq cut sites:           0
      (expanded to length 100bp)
Error! Experimental sample has no analyzable fragments

I do not know if I have done it the right way or not.

[jsh58]

To get the latest update, you must run this:

$ git clone https://github.com/jsh58/Genrich.git
$ cd Genrich
$ make

And then run the newly compiled Genrich program.

[amitjavilaventura]

Thank you. It has run normally!