Set up environment and copy over files
export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/chip_seq_V2/bin:$PATH"Call peaks using MACS3. Call peaks for the other samples: mut_rep2.bam,rescue_rep1.bam,rescue_rep2.bam
macs3 callpeak -g hs -f BAMPE --extsize 200 --shift -100 --nomodel -B -p 1e-3 --keep-dup all -t /home/hgen_share/lec11/mut_rep1.bam -n mut_rep1 --outdir .Run IDR for both conditions
idr --samples mut_rep1_peaks.narrowPeak mut_rep2_peaks.narrowPeak \
--input-file-type narrowPeak \
--rank p.value \
--output-file mut-idr \
--output-file-type bed \
--idr-threshold 0.05 \
--plot \
--log-output-file mut.idr.logPool together peaks. Make a .saf file to use with featureCounts
module load StdEnv/2020 bedtools/2.30.0
cat *-idr | sort -k1,1 -k2,2n | bedtools merge | tee merged.bed | wc -l
bedtools intersect -a merged.bed -b blacklist.hg19.bed -v \
> merged.filtered.bed
awk 'OFS="\t" {print $1"."$2"."$3, $1, $2, $3, "."}' merged.filtered.bed > merged.safRun featureCounts
module load StdEnv/2020 gcc/9.3.0 subread/2.0.1
featureCounts -a merged.saf -F SAF -o xie2022.fc.txt -p /home/hgen_share/lec11/mut_rep1.bam /home/hgen_share/lec11/mut_rep2.bam /home/hgen_share/lec11/rescue_rep1.bam /home/hgen_share/lec11/rescue_rep2.bamAggregate signal centered on ATAC-seq peaks using
computeMatrix reference-point -R open_peaks.bed \
-S mut_rep1.cpm.bw mut_rep2.cpm.bw rescue_rep1.cpm.bw rescue_rep2.cpm.bw \
-bl blacklist.hg19.bed \
--referencePoint center \
--binSize 50 -a 2500 -b 2500 \
--missingDataAsZero --skipZeros \
--samplesLabel mut_1 mut_2 rescue_1 rescue_2 \
-o ATAC.mat.gzDue to memory/speed limitations, you can cancel it with Ctrl + C Instead, the compressed matrix from this analysis has been provided for you in /home/hgen_share/lec11
Choosing the appropriate reference regions influences how the signal behaves. Plot the output based on different reference regions: union_peaks.mat.gz, open_peaks.mat.gz, closed_peaks.mat.gz, mut_peaks.mat.gz, rescue_peaks.mat.gz. Visualize the output through plotHeatmap and/or plotProfile. How do the plots differ using the different reference regions?
plotHeatmap -m ATAC.mat.gz -o heatmap.png --dpi 600 --colorMap Reds --verbose -x "ATAC" -y "" --startLabel "" --endLabel "" --regionsLabel ""
plotProfile -m ATAC.mat.gz -out aggregate_plot.png --perGroup --refPointLabel "ATAC" -T "ATACseq" --dpi 600 --regionsLabel ""Perform motif enrichment analysis using findMotifsGenome.pl
findMotifsGenome.pl open_peaks.bed hg19 HOMER.open_peaks -size given -mask -bg union_peaks.bedYou can cancel it with Ctrl + C