featureCounts output

Question

[QX]

Hi all,

I try to use featureCounts for making count table from the macs2 narrowPeaks files. However, I found that some rows of chr, start, end have repeat value:

Do anyone know what is wrong from my output?

Best

[ATpoint]

What is your featureCounts command? The usual way (at least how I go about this) is to make a SAF file first and then couns reads over the intervals of the SAF.

Converting from BED to SAF/GFF

How to create SAF from text file for FeatureCounts

[QX]

Hi, here is my code (I actually take it from one of your post years ago :))

awk 'OFS="\t" {print $1"."$2+1"."$3, $1, $2+1, $3, "+"}' ${peak_all_dir}peak_all_peaks.narrowPeak > ${peak_all_dir}featureCounts_peaks.saf

featureCounts -a ${peak_all_dir}featureCounts_peaks.saf \
    -F SAF \
    --read2pos 5 \
    -T $SLURM_CPUS_PER_TASK \
    -p \
    -o ${peak_all_dir}featureCounts_peaks_countMatrix.txt ${dedup_dir}*.bam

[QX]

Do you have any idea why it happens to my count table?

Answer 1

[rfran010]

Your narrowPeak files has repeated peaks. Usually I think happens due to different summits within peaks or something. There is a MACS2 option to disable this type of output

[rfran010]

Alternatively, you can filter your peak or SAF file for unique lines. (Can be done with basic command line tools, or bedtools with the peak file).

[QX]

thank you, I will try!

[QX]

hi @rfran010, can you elaborate more on "different summits within peaks"; why this is the case?

[rfran010]

Not sure what type of answer you're interested in, but you can have a region that has multiple peaks that are apparently connected.

Take this random example of ATAC-seq signal, normal behavior will call this as one peak, but you can see three sharp "summits". I believe --call-summits is the relevant option where MACS3 will then re-analyze peaks for summits and instead of calling this as peak1 (red bar) it will report peak1a, peak1b, and peak1c (blue arrows), however each of these will have the exact same chr, start, and end of the red peak. [[please note, this is a random shot from a browser I was looking at. It may or may not be a call-able peak, but either way it illustrates the point]]

See also this post: Why are there duplicate peaks in ENCODE ATAC-seq output?