Hi everyone,
In the original paper for ATAC-seq, the authors shifted the reads +4 bp for the +strand and -5 bp for the -strand: http://www.nature.com/nmeth/journal/v10/n12/full/nmeth.2688.html
How can I easily shift paired-end alignments for my ATAC-seq samples as done in the paper?
Thanks!
3 answers
The bedtools command should extract the paired-end alignments as bedpe format, then the awk command should shift the fragments as needed:
bedtools bamtobed -i reads.bam -bedpe | awk -v OFS="\t" '{if($9=="+"){print $1,$2+4,$6+4}else if($9=="-"){print $1,$2-5,$6-5}}' > fragments.bed
Note, the BAM file should be sorted by read name beforehand:
samtools sort -n -T aln.sorted -o aln.sorted.bam aln.bam
You could convert the BED file back to a BAM file using bedtools bedtobam command
Hi, besides pos/start/end, if I also want to extract the sequence inforamtion about the fragment, what should I do?
You can use bedtool's getfasta command:
bedtools bamtobed -i input.bam | bedtools getfasta -fi genome.fasta -bed stdin > fragments.fasta
Hello James,
If I understand your awk command correctly, $2 represents 5' end of R1 (/left-most) read and $6 represents 5' end of the R2 (/right-most) read. If that is correct, then my question is, why do we add 4 to both ends? Shouldn't it be $2+4 and $6-5?
Thanks.
Hi Ravi, I have the same confusion as yours. Have you understood it? On the first ATAC paper published in 2013, they said "For peak-calling and footprinting, we adjusted the read start sites to represent the center of the transposon binding event." And the ENCODE atac papline showed they adjust one end of the fragment.So I'm mostly in favour of this idea, what do you think?
enter code here
cmd = 'zcat -f {} | '
cmd += 'awk \'BEGIN {{OFS = "\\t"}}'
cmd += '{{ if ($6 == "+") {{$2 = $2 + 4}} '
cmd += 'else if ($6 == "-") {{$3 = $3 - 5}} print $0}}\' | '
cmd += 'gzip -nc > {}'
Hi, I used awk in linux to trim the end of bed file data. $ awk 'BEGIN {OFS = "\t"} ; {if ($6 == "+") print $1, $2 + 4, $3 + 4, $4, $5, $6; else print $1, $2 - 5, $3 - 5, $4, $5, $6}' input.bed >output.bed
To prepare macs2 BEDPE, I believe for paired end
bedtools bamtobed -i reads.bam -bedpe | perl -lane 'if ($F[0] ne $F[3]) {print STDERR "Warnings: read pairs mapped to different chroms: $_"; next;} if ($F[7] eq "+" and $F[8] eq "-") {$F[1]+=4;$F[5]-=5; print "$F[0]\t$F[1]\t$F[5]";}elsif($F[7] eq "-" and $F[8] eq "+"){$F[4]+=4;$F[2]-=5; print "$F[0]\t$F[4]\t$F[2]";}else {print STDERR "Warnings: invalid line: $_";}' | sort -k1,1 -k2,2n -k3,3n > reads.sorted.bed
*if pairs mapped to dierent chrom, it will report Warning
*reads(+) would +4
*reads(-) would -5
*if pairs mapped to the same strand, also report Warning
*coordinates -sorted
Please correct me if I am wrong
Log in to answer this question.
-
ATAC-seq +4 -5 shift in single-end setting before motif analysis
written by sarahmanderni
-
Why should ATAC-seq mapped reads be shifted +4 and -5 for +strand and -strand, respectively
written by progistar •
-
ATAC-seq +4 -5 shift
written by Jingyue •
-
ATAC-SEQ primers order
written by clizama •
-
single-cell RNA-seq power analysis
written by niu.shengyong •
-
Forum: Update: CRISPR off-target effects causing unexpected mutations - good science or media hype?
written by BioinfGuru
-
Published scRNA data with bechmarks?
written by Vivian •
-
Tutorial: Rna-Seq Review Papers
written by Istvan Albert
-
How Pca Was Done To Relate Health Status And Abundance Of Bacteria?
written by Tianyang Li •
-
How Many Human Drug Targets Are There?
written by Andrew Su
Just out of curiosity: did it make a difference on a genome-wide scale, e.g. for annotations, if you do that shifting?
Hi,robm, Does single-end ATAC-seq data need to shift the reads?Do you have any idea about this?
Depends really on what you want to do with the data. What is your goal?
Thank you for you reply. I want to do footprint analysis .
Then I would definitely do it, to be sure to have reads positioned as exact as possible. There is a package, chromVAR, for footprinting analysis. Think you should check it out, and also check if the program does the shifting internally (which I assume it does probably).
Thanks a lot! I used pyDNase to find footprints before, and I will check the chromVAR out.Thank you again!
The chromVAR is especially designed for sparse ATAC-seq data. I never used it, so if you find it helpful, it would be great if you could post here a little summary on your experience. Cheers!