ATAC-seq测序概述

# better alignment results are frequently achieved with --very-sensitive
# use -X 2000 to allow larger fragment size (default is 500)
bowtie2 --very-sensitive -X 2000 -x $Bowtie2Index -1 ${sample}*_1.fq.gz -2 ${sample}*_2.fq.gz -p 10 2> $${sample}.bowtie2.log |
samtools sort -@ 10 -O bam -o ${sample}.sorted.bam

samtools index -@ 10 $WORKDIR/bowtie2/${sample}.sorted.bam

samtools view -@ 10 -h ${sample}.bam | grep -v chrM | samtools sort -@ 10 -O bam -o ${sample}.rmChrM.bam

# REMOVE_DUPLICATES=false: mark duplicate reads, not remove.
# Change it to true to remove duplicate reads.
java -XX:ParallelGCThreads=10 -Djava.io.tmpdir=/tmp -jar picard MarkDuplicates \
  QUIET=true INPUT=${sample}.bam OUTPUT=${sample}.marked.bam METRICS_FILE=${sample}.sorted.metrics \
  REMOVE_DUPLICATES=false CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENT TMP_DIR=/tmp

## dupliate reads labeled as 1024
samtools view -h -b -F 1024 ${sample}.bam > ${sample}.rmDup.bam

# Remove multi-mapped reads (i.e. those with MAPQ < 30, using -q in SAMtools)
samtools view -h -q 30 ${sample}.bam > ${sample}.rmMulti.bam

samtools merge -@ 10 condition1.merged.bam sample1.bam sample2.bam sample3.bam
samtools index -@ 10 condition1.merged.bam

# -f BAMPE, use paired-end information
# --keep-dup all, keep all duplicate reads.
macs2 callpeak -f BAMPE -g hs --keep-dup all --cutoff-analysis -n sample1 -t sample1.shifted.bam --outdir macs2/sample1 2> macs2.log

# bam to bigwig, normalize using 1x effective genome size
# effective genome size: https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html
bamCoverage --numberOfProcessors 8 --binSize 10 --normalizeUsing RPGC --effectiveGenomeSize $effect_genome_size --bam sample1.shifted.bam -o sample1.shifted.bw

samtools view ATAC_sorted.bam | awk '$9>0' | cut -f 9 | sort | uniq -c | sort -b -k2,2n | sed -e 's/^[ \t]*//' > fragment_length_count.txt

samtools index ATAC_sorted.bam
mtReads=$(samtools idxstats ATAC_sorted.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats ATAC_sorted.bam | awk '{SUM += $3} END {print SUM}')

echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'

# the bam used here is sorted bam after duplicates marking
# sort bam by names
samtools sort -@ 10 -n -O BAM -o tmp.bam ATAC.bam

# calculate PBC metrics
bedtools bamtobed -bedpe -i tmp.bam | awk 'BEGIN{OFS="\t"}{print $1,$2,$4,$6,$9,$10}' \
  | grep -v 'chrM' | sort | uniq -c | awk 'BEGIN{mt=0;m0=0;m1=0;m2=0}($1==1){m1=m1+1} \
  ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} \
  END{printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n", mt,m0,m1,m2,m0/mt,m1/m0,m1/m2}' > ATAC_sorted.pbc.qc
rm tmp.bam

# total reads
total_reads=$(samtools view -c ATAC.bam)

# reads in peaks
reads_in_peaks=$(bedtools sort -i ATAC_peaks.narrowPeak \
  | bedtools merge -i stdin | bedtools intersect -u -nonamecheck \
  -a ATAC.bam -b stdin -ubam | samtools view -c)

# FRiP score
FRiP=$(awk "BEGIN {print "${reads_in_peaks}"/"${total_reads}"}")

## reference: https://bedops.readthedocs.io/en/latest/content/usage-examples/master-list.html

summit_bed=(sample1_summits.bed sample2_summits.bed sample3_summits.bed)

out=fAdrenal.master.merge.bed

tmpd=/tmp/tmp$$
mkdir -p $tmpd

## First, union all the peaks together into a single file.
bedlist=""
for bed in ${beds[*]}
do
    bedlist="$bedlist $summit_bed"
done

# extended summits to 500-bp windows (±250 bp)
bedops --range 250 -u $bedlist > $tmpd/tmp.bed

## The master list is constructed iteratively.  For each pass through
## the loop, elements not yet in the master list are merged into
## non-overlapping intervals that span the union (this is just bedops
## -m).  Then for each merged interval, an original element of highest
## score within the interval is selected to go in the master list.
## Anything that overlaps the selected element is thrown out, and the
## process then repeats.
iters=1
solns=""
stop=0
while [ $stop == 0 ]
do
    echo "merge steps..."

    ## Condense the union into merged intervals. This klugey bit
    ## before and after the merging is because we don't want to merge
    ## regions that are simply adjacent but not overlapping
    bedops -m --range 0:-1 $tmpd/tmp.bed \
        | bedops -u --range 0:1 - \
        > $tmpd/tmpm.bed

    ## Grab the element with the highest score among all elements forming each interval.
    ## If multiple elements tie for the highest score, just grab one of them.
    ## Result is the current master list.  Probably don't need to sort, but do it anyway
    ## to be safe since we're not using --echo with bedmap call.
    bedmap --max-element $tmpd/tmpm.bed $tmpd/tmp.bed \
        | sort-bed - \
        > $tmpd/$iters.bed
    solns="$solns $tmpd/$iters.bed"
    echo "Adding `awk 'END { print NR }' $tmpd/$iters.bed` elements"

    ## Are there any elements that don't overlap the current master
    ## list?  If so, add those in, and repeat.  If not, we're done.
    bedops -n 1 $tmpd/tmp.bed $tmpd/$iters.bed \
       > $tmpd/tmp2.bed

    mv $tmpd/tmp2.bed $tmpd/tmp.bed

    if [ ! -s $tmpd/tmp.bed ]
    then
        stop=1
    fi

    ((iters++))
done

## final solution
bedops -u $solns \
   > $out

## Clean up
rm -r $tmpd

exit 0