ATAC-seq(6) -- 寻找motif及差异peak分析

寻找motif

构造HOMER软件指定使用的peak文件格式

mkdir ../homer
ls *.sorted.narrowPeak | while read id;
do 
file=$(basename $id )
sample=${file%%.*} 
echo $sample 
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > ../homer/${sample}.homer_peaks_bed
findMotifsGenome.pl ../homer/${sample}.homer_peaks.bed  mm10 ../homer/${sample}_motifDir -len 8,10,12
done

差异peak分析

构造计数矩阵

#合并bed文件
cat 2-cell-1_peaks.narrowPeak 2-cell-2_peaks.narrowPeak 2-cell-4_peaks.narrowPeak 2-cell-5_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge | uniq > merge.bed

awk '{print "peak_"NR"\t"$1"\t"$2"\t"$3"\t""."}' merge.bed > merge.bed.saf
featureCounts -T 4 -a  merge.bed.saf -F SAF -o counts_subread.txt ../bowtie2/*.bam

head counts_subread.txt
# Program:featureCounts v2.0.1; Command:"featureCounts" "-T" "4" "-a" "merge.bed.saf" "-F" "SAF" "-o" "counts_subread.txt" "../bowtie2/2-cell-1.bam" "../bowtie2/2-cell-2.bam" "../bowtie2/2-cell-4.bam" "../bowtie2/2-cell-5.bam" 
Geneid  Chr Start   End Strand  Length  ../bowtie2/2-cell-1.bam ../bowtie2/2-cell-2.ba../bowtie2/2-cell-4.bam   ../bowtie2/2-cell-5.bam
peak_1  chr1    3516001 3516332 .   332 6   10  2   30
peak_2  chr1    3661830 3662104 .   275 3   9   1   14
peak_3  chr1    3673035 3673612 .   578 16  36  5   21
peak_4  chr1    3728391 3728824 .   434 4   23  3   14
peak_5  chr1    3930572 3930989 .   418 19  21  0   24
peak_6  chr1    3958616 3959549 .   934 103 155 46  130
peak_7  chr1    4051851 4052293 .   443 7   10  0   26
peak_8  chr1    4091560 4092122 .   563 7   20  0   4

拿到表达矩阵后，后面的流程和RNA-Seq分析流程一样，用DESeq2进行差异分析。