PePr相关内容:
提前安装的软件和文件:sratoolkit;fastq;bowtie2;samtools;PePr;
参考基因组:hg38
从原始数据到差异Peak
原始数据下载:
prefetch SRRXXXXXX
##数据批量下载
prefetch --option-file SRR_Acc_List.txt
数据处理:
1. SRRXXXXXX.sra--->SRRXXXXXX.fastq.gz+质量控制
fastq-dump --split-3 --gzip --outdir /outpath1 --split-files /path/SRRXXXXXX.sra
fastqc -f fastq -o /outpath2/ /outpath1/SRRXXXXXX.fastq.gz
2.SRRXXXXXX.fastq.gz--->SRRXXXXXX.sam
建立参考基因组索引:
bowtie2 -build /refgenomePath/hg38.fa /refgenomePath/hg38
单端:
bowtie2 -p 16 -x /refgenomePath/hg38 -U /outpath1/SRRXXXXXX.fastq.gz -S /outpath3/SRRXXXXXX.sam
双端:
bowtie2 -p 16 -x /refgenomePath/hg38 -1 /outpath1/SRRXXXXXX_1.fastq.gz -2 /outpath1/SRRXXXXXX_2.fastq.gz -S /outpath3/SRRXXXXXX.sam
3.sam-->bam--->sorted.bam
samtools view -bS /outpath3/SRRXXXXXX.sam >/outpath3/SRRXXXXXX.bam
samtools sort -o /outpath3/SRRXXXXXX.sort.bam /outpath3/SRRXXXXXX.bam
4.生成bam.bai文件
samtools index /outpath3/SRRXXXXXX.sort.bam
5.识别差异Peak
这里需要注意:PePr的输入文件必须是具有生物学重复的,如果没有生物学重复可以把把同一个样本写两次放进去,也会识别处差异peak
例如:
PePr -c SRR6418912.sort.bam,SRR6418912.sort.bam --chip2 SRR6418918.sort.bam,SRR6418918.sort.bam -f bam --diff
bed文件注释
#安装BiocManager::install("ChIPseeker")
#加载library(ChIPseeker)
#安装人的注释包
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
#读取chipseq峰的bed文件
Scr_peak <- readPeakFile("/home/wanghongli/THOR_example_data/NA__PePr_chip2_peaks.bed")
#注释，TSS的范围可自定义
#加载人基因组注释包
require(TxDb.Hsapiens.UCSC.hg38.knownGene)
#对txdb进行指定
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
#进行注释
Scr_peakAnno <- annotatePeak(Scr_peak, tssRegion = c(-3000, 3000), TxDb = txdb)
#输出结果
write.table(Scr_peakAnno, file = "Scr_peak.txt",sep = '\t', quote = FALSE, row.names = FALSE)
合并bed文件
# 两个BED文件，第一步合并成一个文件，但来自两个文件的区域是分开的
cat A.bed B.bed > C.bed
# 合并后的文件PEAK数是合并前两个文件的总和
$ sort -k1,1 -k2,2n C.bed > D.bed
