peaks差异分析 |
您所在的位置：网站首页 › deseq2标准化方法 › peaks差异分析 |
|
1. 获得peak集
image.png
这里的逻辑是：把四个样品合并，call peaks，然后获得peaks文件与前面idr 处理后的peaks进行overlap，都overlap的peaks，作为最终的peaks。 gsize=199000000 ## callpeak and idr analysis of sample A s1_r1=SRR6322534 s1_r2=SRR6322535 s1=SRR6322534_SRR6322535 input1=SRR6322538 # call peaks for replicte 1 macs2 callpeak -t ./${s1_r1}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1_r1} -g $gsize --keep-dup all -p 0.01 # call peaks for replicte 2 macs2 callpeak -t ./${s1_r2}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1_r2} -g $gsize --keep-dup all -p 0.01 # call peaks for combined dataset macs2 callpeak -t ./${s1_r1}.ff.bam ./${s1_r2}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1} -g $gsize --keep-dup all -p 0.01 # run idr idr --samples ${s1_r1}_peaks.narrowPeak ${s1_r2}_peaks.narrowPeak --peak-list ${s1}_peaks.narrowPeak --input-file-type narrowPeak --output-file ./${s1}.idr --rank p.value --soft-idr-threshold 0.05 --plot # get idr produced narrowPeak file cut -f 1-10 ${s1}.idr | sort -k1,1 -k2,2n -k3,3n >${s1}.idr.narrowPeak ## callpeak and idr analysis of sample B s2_r1=SRR8423051 s2_r2=SRR8423052 s2=SRR8423051_SRR8423052 input2=SRR8423055 # call peaks for replicte 1 macs2 callpeak -t ./${s2_r1}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2_r1} -g $gsize --keep-dup all -p 0.01 # call peaks for replicte 2 macs2 callpeak -t ./${s2_r2}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2_r2} -g $gsize --keep-dup all -p 0.01 # call peaks for combined dataset macs2 callpeak -t ./${s2_r1}.ff.bam ./${s2_r2}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2} -g $gsize --keep-dup all -p 0.01 # run idr idr --samples ${s2_r1}_peaks.narrowPeak ${s2_r2}_peaks.narrowPeak --peak-list ${s2}_peaks.narrowPeak --input-file-type narrowPeak --output-file ./${s2}.idr --rank p.value --soft-idr-threshold 0.05 --plot # get idr produced narrowPeak file cut -f 1-10 ${s2}.idr | sort -k1,1 -k2,2n -k3,3n >${s2}.idr.narrowPeak ## Peaks in combined sample bams macs2 callpeak -t ${s1_r1}.ff.bam ${s1_r2}.ff.bam ${s2_r1}.ff.bam ${$s2_r2}.ff.bam -c $input1.ff.bam $input2.ff.bam -f BAM -n ${s1}_${s2} -g $gsize --keep-dup all -p 0.01 获得最终的peaks文件 Cat ${s1_r1}_${s1_r2}.idr.narrowPeak ${s2_r1}_${s2_r2}.idr.narrowPeak | sort -k1,1 -k2,2n -k3,3n | bedtools intersect -a ${s1}_${s2}_peaks.narrowPeak -b - -F 0.5 -u >${s1}_${s2}_peaks.f.narrowPeak ## filter peaks against blacklist bedtools intersect -a ${s1}_${s2}_peaks.f.narrowPeak -b ../ref/blacklist.bed -f 0.25 -v >temp && mv temp ${s1}_${s2}_peaks.f.narrowPeak image.png 2. 统计peaks区域的counts，主要利用deeptools的 multiBamSummary s1=SRR6322534 s2=SRR6322535 s3=SRR8423051 s4=SRR8423052 peak=./peak.narrowPeak ## for comparison between samples without replicate cut -f 1-3 $peak >peak.1.bed #把peaks转换成 bed 文件 multiBamSummary BED-file --bamfiles ../callpeak/downsample/$s2.ff.bam ../callpeak/downsample/$s3.ff.bam --BED peak.1.bed --smartLabels -p 4 --outRawCounts peak_counts.1.txt --extendReads 134 ## for comparison between samples with replicates cut -f 1-3 $peak >peak.2.bed multiBamSummary BED-file --bamfiles ../callpeak/downsample/$s1.ff.bam ../callpeak/downsample/$s2.ff.bam ../callpeak/downsample/$s3.ff.bam ../callpeak/downsample/$s4.ff.bam --BED peak.2.bed --smartLabels -p 4 --outRawCounts peak_counts.2.txt --extendReads 134The multiBamSummary command is part of the deepTools package and is used to generate summary statistics for multiple BAM files. Here's an explanation of the options in the example command you provided: BED-file: The name of the BED file containing the genomic regions of interest. --bamfiles: A list of BAM files to analyze. --smartLabels: Automatically generate labels for the BAM files based on their file names. -p: The number of threads to use for parallel processing. --outRawCounts: The name of the output file containing the raw counts for each genomic region. --extendReads: The number of base pairs to extend the reads in each direction. 3. R里面进行差异分析 1. read in sample information library(DESeq2) library(tidyverse) library(pheatmap) cat meta.2.txt SRR6322534 rabbit, anti-IR_beta, sc-711 HepG2_IRb 19272489 SRR6322535 rabbit, anti-IR_beta, sc-711 HepG2_IRb 37475005 SRR8423051 rabbit, anti-IR_beta, sc-711 HepG2_IRb_starve 28559621 SRR8423052 rabbit, anti-IR_beta, sc-711 HepG2_IRb_starve 45616747 meta |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |