This is a test version of Biostars. For the public version, visit https://www.biostars.org.
ATAC-seq peaks centered between TSS and TES

I am analyzing ATAC-seq data and right now trying to confirm my reads enriched around TSS sites however I am obtaining enrichment between TSS and TES for all my samples in two cell lines.

The code I run:

computeMatrix reference-point \
  -S m7_rep1.bw \
  -R hg38.TSS.bed \
  --beforeRegionStartLength 1000 \
  --regionBodyLength 2000 \
  --afterRegionStartLength 1000 \
  --binSize 100 \
  -o m7_rep1_matrix.gz
plotHeatmap \
  -m m7_rep1_matrix.gz \
  -out m7_rep1_TSS_enrichment_heatmap.png

Enrichment

I generated the hg38.TSS.bed

awk '$3 == "transcript"
                {
                    if ($7 == "+") print $1 "\t" $4-1000 "\t" $4+1000; 
                    else if ($7 == "-") print $1 "\t" $5-1000 "\t" $5+1000;}
                ' hg38.refGene.gtf > hg38.TSS.bed

I am also not sure about the black genes.

What could be the thing I am doing or is going wrong?

I appreciate any help, thank you

compute-matrix atac atac-seq tss-enrichment enrichment-graph
• 3,019 views • link
written by buffealo • buffealo

To troubleshoot, maybe try separating the transcripts by strand, instead of mixing the two. In your awk statement, for instance, write one or the other strand if-else case and visualize that. This may help highlight the problem.

• 0 views • link
written by Alex Reynolds Alex Reynolds

I concluded that I initially made the wrong configuration of TSSs initially; thank you so much for your directive reply.

• 0 views • link
written by buffealo • buffealo

2 answers

I think this is a simple labelling issue. I do not know DeepTools well, but based on your awk you simply give it a chr-start-end BED file, so it doesn't know about TSS and TES. All it sees is a 2000bp interval, and you seem to tell it to extend this additionally by 1kb to each direction. The plot shows exactly this. Your actuall TSS (based on the BED file) is centered and perfectly enriched for signal. The borders in your BED file are labelled as TSS/TES (maybe a DeepTools default?) and the additional 1kb stretches up/downstream are called -1.0 and 1.0. I assume the only issue is that you need to change the labelling.

• 0 views • link
written by ATpoint ATpoint

I got your point. This time, I established my configuration like this:

TSSs: To target the exact TSSs, only considered one location and its 1 bp distance;

 awk '$3 == "transcript" { if ($7 == "+") { print $1"\t"$4-1"\t"$4 } else { print $1"\t"$5-1"\t"$5 } }' hg38.refGene.gtf > exact.TSS

Heatmap generation: I pretty much used the computeMatrix in default settings, except that I added --referencePoint TSS and a piece of calmer color --colorMap 'Greens'

computeMatrix reference-point -S m7_rep1.bw -R exact.TSS --referencePoint TSS  -o m7_rep1_matrix_TSS.gz; plotHeatmap -m m7_rep1_matrix_TSS.gz --colorMap 'Greens' -out m7_rep1_TSS_enrichment_heatmap_exact.png

where I think right now I am only considering the exact transcription start sites (correct positioning) and enrichment around TSS (correct labeling)

Now the result is:

heatmap

• 0 views • link
written by buffealo • buffealo

Looks reasonable. The missing values probably come from a bigwig that has not included missing signal as zeros.

• 0 views • link
written by ATpoint ATpoint

If –missingDataAsZero was not set, such cases will be colored in black by default.

see plotHeatmap docs

• 0 views • link
written by gglim • gglim

Log in to answer this question.

More posts like this