PePr相關(guān)內(nèi)容:
提前安裝的軟件和文件:sratoolkit;fastq;bowtie2;samtools;PePr;
參考基因組:hg38?
從原始數(shù)據(jù)到差異Peak
原始數(shù)據(jù)下載:
prefetch SRRXXXXXX
##數(shù)據(jù)批量下載
prefetch --option-file SRR_Acc_List.txt
數(shù)據(jù)處理:
1. SRRXXXXXX.sra--->SRRXXXXXX.fastq.gz+質(zhì)量控制
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.識(shí)別差異Peak
這里需要注意:PePr的輸入文件必須是具有生物學(xué)重復(fù)的,如果沒有生物學(xué)重復(fù)可以把把同一個(gè)樣本寫兩次放進(jìn)去,也會(huì)識(shí)別處差異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)
#對(duì)txdb進(jìn)行指定
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene?
#進(jìn)行注釋
Scr_peakAnno <- annotatePeak(Scr_peak, tssRegion = c(-3000, 3000), TxDb = txdb)
#輸出結(jié)果
write.table(Scr_peakAnno, file = "Scr_peak.txt",sep = '\t', quote = FALSE, row.names = FALSE)
合并bed文件
# 兩個(gè)BED文件，第一步合并成一個(gè)文件，但來自兩個(gè)文件的區(qū)域是分開的
?cat A.bed B.bed > C.bed
# 合并后的文件PEAK數(shù)是合并前兩個(gè)文件的總和
$ sort -k1,1 -k2,2n C.bed > D.bed
