from 生信技能樹
paper: Brookes, E. et al. Polycomb associates genome-wide with a specific RNA polymerase II variant, and regulates metabolic genes in ESCs. Cell Stem Cell 10, 157–170 (2012).

GSE accession: GSE34518侍郭，總共是9個樣本，但是很多樣本都分開在多個lane測序的，所以每個樣本其實是有多個sra文件麸恍，需要進(jìn)行合并校辩。

Use prefetch to download them all, then transform those SRA files to fastq files by sra-toolkits, then align them to mm10, and call peaks.

測序用的是： Illumina Genome Analyzer II 測序儀，測序策略是 SE50。
data="/home/zzz/data/test"

#查找sra accesion list窍霞，用vim命令寫入sra.list文件
SRR391032
SRR391033
SRR391034
SRR391035
SRR391036
SRR391037
SRR391038
SRR391039
SRR391040
SRR391041
SRR391042
SRR391043
SRR391044
SRR391045
SRR391046
SRR391047
SRR391048
SRR391049
SRR391050

下載sra并且轉(zhuǎn)換為fastq

cd $data
conda install -c bioconda sra-tools #安裝sra-tools
cd $data/sra
cat sra.list |while read i; do prefetch $i; done #download sra files
## 默認(rèn)下載目錄：~/ncbi/public/sra/ ,但是我的是下載在當(dāng)前文件夾下的罩旋，而且每個sra文件都是一個文件夾啊央。
ls -lh ~/ncbi/public/sra/

第一步需要制作配置文件，需要先從ncbi下載sra的information文件sra.table

## 直接用excel制作config文件涨醋，或者寫代碼, 我當(dāng)然是用的excel瓜饥，這個代碼太麻煩了啦。浴骂。乓土。
cut -f 4,7 sra.table |cut -d":" -f 2 |sed 's/ChIPSeq//g' | sed 's/MockIP//g'|sed  's/^ //' |tr ' ' '_' |perl -alne '{$h{$F[0]}++ if exists $h{$F[0]}; $h{$F[0]}=1 unless exists $h{$F[0]};print "$F[0]$h{$F[0]}\t$F[1]"}' > config

得到內(nèi)容如下，命名為config.txt

RNAPII_S5P_1    SRR391032
RNAPII_S5P_2    SRR391033
RNAPII_S2P_1    SRR391034
RNAPII_S7P_1    SRR391035
RNAPII_8WG16_1    SRR391036
RNAPII_8WG16_2    SRR391037
RNAPII_S2P_2    SRR391038
RNAPII_S2P_3    SRR391039
RNAPII_S7P_2    SRR391040
H2Aub1_1    SRR391041
H2Aub1_2    SRR391042
H3K36me3_1    SRR391043
H3K36me3_2    SRR391044
Control_1    SRR391045
Control_2    SRR391046
Ring1B_1    SRR391047
Ring1B_2    SRR391048
Ring1B_3    SRR391049
RNAPII_S5PRepeat_1    SRR391050

批量sra轉(zhuǎn)fq文件

## 下面用到的 config 文件溯警，就是上面自行制作的趣苏。
cat config.txt |while read i;
do echo $i
arr=($i)
srr=${arr[1]}
sample=${arr[0]}
# 單端測序數(shù)據(jù)的sra轉(zhuǎn)fasq
fastq-dump -A  $sample -O $data/raw  --gzip --split-3  $data/sra; done
# -A: Replaces accession derived from <path> in filename(s) and deflines (only for single table dump)
#不知道為什么我的這個命令就是跑不通，所以我用了自己的代碼：
fastq-dump --gzip --split-3 $data/sra/*sra
#然后把fastq.gz文件手動改了名稱愧膀，其實可以下載的時候直接用-O參數(shù)重新命名拦键。

qc

cd $data/raw/
ls *gz | xargs fastqc -t 10 -o $data/raw/qc
multiqc $data/raw/qc -o $data/raw/qc

用trim_galore進(jìn)行數(shù)據(jù)過濾，此處是單端測序

cd $data/raw/
ls *gz | while read i;
do 
nohup trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $data/clean  $i & 
done 

trimmed_data進(jìn)行QC

mkdir $data/clean/qc
ls $data/clean/*gz | xargs fastqc -t 10 -o $data/clean/qc
multiqc $data/clean/qc -o $data/raw/qc

使用bowtie2進(jìn)行比對

直接用bowtie2進(jìn)行比對和統(tǒng)計比對率, 需要提前下載參考基因組然后使用命令構(gòu)建索引檩淋，或者直接就下載索引文件芬为，這里用常用的mm10

# 索引大小為3.2GB萄金， 不建議自己下載基因組構(gòu)建，可以直接下載索引文件媚朦，代碼如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip

單端測序數(shù)據(jù)進(jìn)行比對

ls ../clean/*gz |while read i;
do 
file=$(basename $i)
sample=${file%%.*}
# echo $file $sample
bowtie2 -p 5 -x $bowtie2_index -U  $i | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done #最后只輸出bam文件氧敢，沒有sam文件

注意到這里沒有用samtools 的view命令將sam轉(zhuǎn)換成bam文件，然后再sort询张，我就比較了一下這兩個命令的區(qū)別孙乖。下面是結(jié)果：

Activated-24h-1.bam.count.txt:ENSMUSG00000000028    3092
Activated-24h-2.bam.count.txt:ENSMUSG00000000028    3483
Activated-72h-1.bam.count.txt:ENSMUSG00000000028    2260
Activated-72h-2.bam.count.txt:ENSMUSG00000000028    1989
Resting-1.bam.count.txt:ENSMUSG00000000028  157
Resting-2.bam.count.txt:ENSMUSG00000000028  185

對過濾后mapping得到的bam文件進(jìn)行QC

cd $data/bam
ls *.bam  |xargs -i samtools index {}
ls *.bam  |while read i ;do (nohup samtools flagstat $i > $(basename $i ".bam").stat & );done 

grep '%' *stat #得到所有mapping率
Control_1_trimmed.stat:7415799 + 0 mapped (87.76% : N/A)
Control_2_trimmed.stat:7202987 + 0 mapped (86.18% : N/A)
H2Aub1_1_trimmed.stat:8949983 + 0 mapped (97.19% : N/A)
H2Aub1_2_trimmed.stat:13195346 + 0 mapped (97.27% : N/A)
H3K36me3_1_trimmed.stat:11732645 + 0 mapped (98.85% : N/A)
H3K36me3_2_trimmed.stat:13404798 + 0 mapped (98.38% : N/A)
Ring1B_1_trimmed.stat:4620324 + 0 mapped (93.31% : N/A)
Ring1B_2_trimmed.stat:4633085 + 0 mapped (93.57% : N/A)
Ring1B_3_trimmed.stat:22884650 + 0 mapped (95.12% : N/A)
RNAPII_8WG16_1_trimmed.stat:7465754 + 0 mapped (96.07% : N/A)
RNAPII_8WG16_2_trimmed.stat:20652996 + 0 mapped (95.30% : N/A)
RNAPII_S2P_1_trimmed.stat:24966425 + 0 mapped (97.05% : N/A)
RNAPII_S2P_2_trimmed.stat:6095629 + 0 mapped (94.73% : N/A)
RNAPII_S2P_3_trimmed.stat:8659690 + 0 mapped (96.81% : N/A)
RNAPII_S5P_1_trimmed.stat:11791542 + 0 mapped (97.45% : N/A)
RNAPII_S5P_2_trimmed.stat:12171421 + 0 mapped (98.09% : N/A)
RNAPII_S5PRepeat_1_trimmed.stat:4158664 + 0 mapped (82.71% : N/A)
RNAPII_S7P_1_trimmed.stat:6378735 + 0 mapped (80.81% : N/A)
RNAPII_S7P_2_trimmed.stat:5962539 + 0 mapped (82.54% : N/A)

有幾個文件建index的時候報錯，查了網(wǎng)站之后份氧，用了下面的命令唯袄，發(fā)現(xiàn)正常的bam文件返回值是VALID，而報錯的文件不是蜗帜。因此我重新跑了QC命令之后就好了恋拷。

gunzip -t H2Aub1_1_trimmed.bam && echo "VALID"
gzip: H2Aub1_1_trimmed.bam: decompression OK, trailing garbage ignored
gunzip -t Control_1_trimmed.bam && echo "VALID"
VALID

合并bam文件
samtools 的merge命令只能用來merge sort過的bam文件
當(dāng)有多個bam文件時，一般思路就是對每一個bam進(jìn)行sort厅缺、index后蔬顾，再merge成一個整體merged.bam，然后對merged.bam再進(jìn)行sort湘捎、index诀豁，才算能用了，得到最終結(jié)果應(yīng)該是是sorted.merge.bam

## 如果不用循環(huán)：
## samtools merge  control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam
## 循環(huán)命令
mkdir  ~/project/epi/ mergeBam
cd ~/project/epi/align
ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge $data/merge.bam/$i.merge.bam $i*.bam ;done
#sort -u 的意思是輸出的時候去除重復(fù)

mapping過的bam文件窥妇，進(jìn)行篩選舷胜，qc

#去除duplicate
ls  $data/merge.bam/*merge.bam  | while read i ;do (nohup samtools markdup -r $i $(basename $id ".bam").rmdup.bam & );done #這樣寫nohup命令雖然可以幾個同時跑，比較快秩伞，但是我不知道nohup命令的log文件被寫到哪里了逞带。
#去除低mapping的reads，以及mapping多次的reads
ls $data/merge.bam/*rmdup.bam | while read i; 
do
file=$(basename $i)
sample=${file%%.*}
echo $file $sample
samtools view $i -bF 4 -q 10 > ${sample}.uniq.bam
done
#過濾后的reads進(jìn)行qc
ls $data/merge.bam/*.uniq.bam  |xargs -i samtools index {}
ls $data/merge.bam/*.uniq.bam  |while read i ;do
samtools flagstat $i > $(basename $i ".uniq.bam").stat;done 
grep '%' *stat
control.stat:10120184 + 0 mapped (100.00% : N/A)
H2Aub1.stat:14453515 + 0 mapped (100.00% : N/A)
H3K36me3.stat:20416595 + 0 mapped (100.00% : N/A)
Ring1B.stat:20587109 + 0 mapped (100.00% : N/A)
RNAPII_8WG16.stat:19619979 + 0 mapped (100.00% : N/A)
RNAPII_S2P.stat:21815365 + 0 mapped (100.00% : N/A)
RNAPII_S5PRepeat.stat:3324476 + 0 mapped (100.00% : N/A)
RNAPII_S5P.stat:7593988 + 0 mapped (100.00% : N/A)
RNAPII_S7P.stat:8086221 + 0 mapped (100.00% : N/A)

使用macs2進(jìn)行call peak

cd  $data/merge.bam 
ls  *uniq.bam |cut -d"." -f 1 |while read i;
do 
macs2 callpeak -c  control.uniq.bam -t $i.uniq.bam -f BAM -B -g mm -n $i --outdir $data/peaks  2 > $id.log
done  

-t: 實驗組的輸出結(jié)果
-c: 對照組的輸出結(jié)果
-f: -t和-c提供文件的格式纱新，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一個。如果不提供這項穆趴，就是自動檢測選擇脸爱。
-g: 基因組大小， 默認(rèn)提供了hs, mm, ce, dm選項未妹， 不在其中的話簿废，比如說擬南芥，就需要自己提供了络它。
-n: 輸出文件的前綴名
-B: 會保存更多的信息在bedGraph文件中族檬，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores, 但這個參數(shù)意義不大，得到的bedgraph文件沒啥用化戳。
-q: q值单料，也就是最小的PDR閾值埋凯， 默認(rèn)是0.05。q值是根據(jù)p值利用BH計算扫尖，也就是多重試驗矯正后的結(jié)果白对。
-p： 這個是p值，指定p值后MACS2就不會用q值了换怖。
-m: 和MFOLD有關(guān)甩恼，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)，默認(rèn)是5：50沉颂，MACS會先尋找100多個peak區(qū)構(gòu)建模型条摸，一般不用改，因為你很大概率上不會懂铸屉。

對比了一下只去掉adaptor屈溉，去掉adaptor和duplicate，去除adaptor以及沒有比對上抬探、多重比對子巾、duplicate的peaks（raw_peaks，rmdup_peaks小压，peaks）

wc -l *bed
#raw_peaks
       0 ../raw_peaks/control_summits.bed
    1182 ../raw_peaks/H2Aub1_summits.bed
   40034 ../raw_peaks/H3K36me3_summits.bed
   26029 ../raw_peaks/Ring1B_summits.bed
   41628 ../raw_peaks/RNAPII_8WG16_summits.bed
   20029 ../raw_peaks/RNAPII_S2P_summits.bed
   38659 ../raw_peaks/RNAPII_S5PRepeat_summits.bed
   56672 ../raw_peaks/RNAPII_S5P_summits.bed
   72203 ../raw_peaks/RNAPII_S7P_summits.bed
  296436 total
#rmdup_peaks
       0 control_summits.bed
    1182 H2Aub1_summits.bed
   39815 H3K36me3_summits.bed
   26029 Ring1B_summits.bed
   41628 RNAPII_8WG16_summits.bed
   20029 RNAPII_S2P_summits.bed
   38659 RNAPII_S5PRepeat_summits.bed
   56750 RNAPII_S5P_summits.bed
   72203 RNAPII_S7P_summits.bed
  296295 total
#peaks
       0 ../peaks/control_summits.bed
    2203 ../peaks/H2Aub1_summits.bed
   23031 ../peaks/H3K36me3_summits.bed
   25062 ../peaks/Ring1B_summits.bed
   34926 ../peaks/RNAPII_8WG16_summits.bed
   20118 ../peaks/RNAPII_S2P_summits.bed
   33668 ../peaks/RNAPII_S5PRepeat_summits.bed
   54452 ../peaks/RNAPII_S5P_summits.bed
   62553 ../peaks/RNAPII_S7P_summits.bed
  256013 total

發(fā)現(xiàn)去不去PCR duplicate對于最后call peaks的結(jié)果不會影響太大线梗，但是去不去多重比對和沒有比對對結(jié)果影響比較大。

根據(jù)生信技能樹的經(jīng)驗：
前幾個月處理這個數(shù)據(jù)集的時候使用的過濾低質(zhì)量reads參數(shù)是短于 35bp的全部丟棄怠益，現(xiàn)在是短于25bp的全部拋棄仪搔，導(dǎo)致了得到的peaks從數(shù)量上千差別不小。

MACS2結(jié)果

_peaks.xls

This file is generated by MACS version 2.2.6
Command line: callpeak -c Activated_input.bam -t Activated_H2AK9ac.bam -f BAM -B -g mm -n Activated_H2AK9ac
ARGUMENTS LIST:
name = Activated_H2AK9ac
format = BAM
ChIP-seq file = ['Activated_H2AK9ac.bam']
control file = ['Activated_input.bam']
effective genome size = 1.87e+09
band width = 300
model fold = [5, 50]
qvalue cutoff = 5.00e-02
The maximum gap between significant sites is assigned as the read length/tag size.
The minimum length of peaks is assigned as the predicted fragment length "d".
Larger dataset will be scaled towards smaller dataset.
Range for calculating regional lambda is: 1000 bps and 10000 bps
Broad region calling is off
Paired-End mode is off

tag size is determined as 50 bps
total tags in treatment: 25164589
tags after filtering in treatment: 25164589
maximum duplicate tags at the same position in treatment = 1
Redundant rate in treatment: 0.00
total tags in control: 18816543
tags after filtering in control: 18816543
maximum duplicate tags at the same position in control = 1
Redundant rate in control: 0.00
d = 214
alternative fragment length(s) may be 214 bps
chr     start   end     length  abs_summit      pileup  -log10(pvalue)  fold_enrichment -log10(qvalue)  name
chr1   4784036 4784310    275    4784120         17.20     8.44714        4.54950 6.37316      Activated_H2AK9ac_peak_1

#length: peak區(qū)域長度
#abs_summit: peak的峰值位點（summit position）
#pileup: peak 峰值的高度（pileup height at peak summit, -log10(pvalue) for the peak summit）
#fold_enrichment: peak的富集倍數(shù)（相對于random Poisson distribution with local lambda）

Coordinates in xls is 1-based which is different with BED format
xls里的坐標(biāo)和bed格式的坐標(biāo)還不一樣蜻牢，起始坐標(biāo)需要減1才與narrowPeak的起始坐標(biāo)一樣烤咧。

_peaks.narrowPeak

narrowPeak文件是BED6+4格式，可以上傳到UCSC瀏覽抢呆。輸出文件每列信息分別包含：
1煮嫌；染色體號
2：peak起始位點
3：結(jié)束位點
4：peak name
5：int(-10*log10qvalue)
6 ：正負(fù)鏈
7：fold change
8：-log10pvalue
9：-log10qvalue
10：relative summit position to peak start（？）

chr1    4784035 4784310 Activated_H2AK9ac_peak_1        63      .       4.54950 8.44714 6.37316 84
chr1    4784829 4785734 Activated_H2AK9ac_peak_2        586     .       12.83651        61.82663        58.671

_summits.bed

BED格式的文件抱虐，包含peak的summits位置昌阿，第5列是-log10pvalue。如果想找motif恳邀，推薦使用此文件懦冰。

chr1    3670687 3670688 Activated_H3K27me3_peak_1       7.21477
chr1    3671910 3671911 Activated_H3K27me3_peak_2       16.00848

補充：

bed格式簡介
BED格式能夠非常簡潔的表示基因組特征和注釋，盡管BED格式描述中定義了12列谣沸，但是僅僅只有3列必須刷钢，因此BED格式按照列數(shù)繼續(xù)細(xì)分為BED3,BED4,BED5,BED6,BED12。
BED12定義的12列分別為：chrom, start, end, name(BED代表的特征名),score(范圍為0~1000乳附，可以是pvalue, 或者是字符串,如"up"), strand(正負(fù)鏈), thickstart, thickednd(額外著色位置, 比如說表示外顯子), itemRgb(RGB顏色,如255,0,0), blockCount(區(qū)塊數(shù)量, 如外顯子), blockSizes(由逗號隔開的區(qū)塊大小), blockStarts(由逗號隔開的區(qū)塊起始位點)内地。
對BED文件的細(xì)分格式進(jìn)行舉例說明
BED3：chr1 11873 14409
BED4: chr1 11873 14409 uc001aaa.3
BED5: chr1 11873 14409 uc001aaa.3 0
BED6: chr1 11873 14409 uc001aaa.3 0 +
BED12: chr1 11873 14409 uc001aaa.3 0 + 11873 12000 123,123,123 3 354,109,1189, 0,739,1347,

1-based coordinate system:序列的第一個堿基設(shè)為數(shù)字1伴澄，如SAM, VCF, GFF, wiggle格式
0-based coordinate system :序列的第一個堿基設(shè)為數(shù)字0，如BAM, BCFv2, BED, PSL格式

.bdg

bedGraph格式瓤鼻，可以導(dǎo)入UCSC或者轉(zhuǎn)換為bigwig格式秉版。兩種bfg文件：treat_pileup, and control_lambda.
NAME_peaks.broadPeak
BED6+3格式與narrowPeak類似，只是沒有第10列茬祷。

image.png

使用deeptools進(jìn)行可視化

還有很多其他用法清焕，參考：https://vip.biotrainee.com/d/226-如何使用deeptools處理bam數(shù)據(jù)

deeptools提供bamCoverage和bamCompare進(jìn)行格式轉(zhuǎn)換，為了能夠比較不同的樣本祭犯，需要先將基因組分成等寬分箱(bin)秸妥，統(tǒng)計每個分箱的read數(shù)，最后得到描述性統(tǒng)計值沃粗。對于兩個樣本粥惧，描述性統(tǒng)計值可以是兩個樣本的比率，或是比率的log2值最盅，或者是差值突雪。如果是單個樣本，可以用SES方法進(jìn)行標(biāo)準(zhǔn)化涡贱。

bamCoverage的基本用法

bamCoverage -e 170 -bs 10 -b chip_sorted.bam -o chip.bw
# chip_sorted.bam是前期比對得到的BAM文件

得到的bw文件就可以送去IGV/Jbrowse進(jìn)行可視化咏删。 這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來的read長度，且設(shè)置分箱的大小问词。其他參數(shù)還有

--filterRNAstrand {forward, reverse}: 僅統(tǒng)計指定正鏈或負(fù)鏈
--region/-r CHR:START:END:選取某個區(qū)域統(tǒng)計
--smoothLength:通過使用分箱附近的read對分箱進(jìn)行平滑化
如果為了其他結(jié)果進(jìn)行比較督函，還需要進(jìn)行標(biāo)準(zhǔn)化，deeptools提供了如下參數(shù)：

--scaleFactor:縮放系數(shù)
--normalizeUsingRPKMReads: Per Kilobase per Million mapped reads (RPKM)標(biāo)準(zhǔn)化
--normalizeTo1x: 按照1x測序深度(reads per genome coverage, RPGC)進(jìn)行標(biāo)準(zhǔn)化
--ignoreForNormalization： 指定那些染色體不需要經(jīng)過標(biāo)準(zhǔn)化
如果需要以100為分箱激挪，并且標(biāo)準(zhǔn)化到1x辰狡，且僅統(tǒng)計某一條染色體區(qū)域的正鏈，輸出格式為bedgraph,那么命令行可以這樣寫

mouse測序數(shù)據(jù)垄分，2號染色體中間會浪費一堆reads
bedgraph轉(zhuǎn)換成tdf格式宛篇，igvtools可以完成轉(zhuǎn)換，tdf格式文件可以快速被igv打開

查看TSS附近的信號強度

#先得到bw文件锋喜，之前已經(jīng)得到些己，這里就不再贅述
bed="/home/zzz/data/reference/annotation/mm10_refseq.bed"
for i in  $data/bw/*rmdup.bw;
do 
echo $i
file=$(basename $i )
sample=${file%%.*} 
echo $sample  
computeMatrix reference-point  --referencePoint TSS   \
-b 2000 -a 2000    \
-R  $bed \
-S $i  \
--skipZeros  -o matrix1_${sample}_TSS_2K.gz  \
--outFileSortedRegions ${sample}_TSS_2K.bed
# 輸出的gz為文件用于plotHeatmap, plotProfile

## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.png
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720 

done 

image.png

畫genebody的圖

bed="/home/zzz/data/reference/annotation/mm10_refseq.bed"
for i in  $data/bw/*.rmdup.bw ;
do 
echo $i
file=$(basename $i )
sample=${file%%.*} 
echo $sample 

computeMatrix scale-regions  -p 5  \
-b 3000 -a 3000    \
-R  $bed \
-S $i  \
--regionBodyLength 15000 --skipZeros \
-o matrix1_${sample}_genebody.gz  \
--outFileSortedRegions regions1_${sample}_genebody.bed

## both plotHeatmap and plotProfile will use the output from   computeMatrix
plotHeatmap -m matrix1_${sample}_genebody.gz  -out ${sample}_Heatmap_genebody.png
plotHeatmap -m matrix1_${sample}_genebody.gz  -out ${sample}_Heatmap_genebody.pdf --plotFileFormat pdf  --dpi 720  #輸出文件有兩種格式，一個是png嘿般，一個是pdf，實操的時候可以只保留一個涯冠。
plotProfile -m matrix1_${sample}_genebody.gz  -out ${sample}_Profile_genebody.png
plotProfile -m matrix1_${sample}_genebody.gz  -out ${sample}_Profile_genebody.pdf --plotFileFormat pdf --perGroup --dpi 720 
done 

上面的批量代碼其實就是為了統(tǒng)計全基因組范圍的peak在基因特征的分布情況炉奴，也就是需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進(jìn)行可視化蛇更，用plotProfile以折線圖的方式展示覆蓋情況瞻赶。

computeMatrix具有兩個模式:scale-region和reference-point赛糟。前者用來信號在一個區(qū)域內(nèi)分布，后者查看信號相對于某一個點的分布情況砸逊。無論是那個模式璧南，都有兩個參數(shù)是必須的，-S是提供bigwig文件师逸，-R是提供基因的注釋信息司倚。還有更多個性化的可視化選項。

使用R包對找到的peaks文件進(jìn)行注釋

bedPeaksFile ='8WG16_summits.bed'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile )  
keepChr= !grepl('_',seqlevels(peak)) #將chr_去掉篓像，只留下常染色體和性染色體
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Mm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)

更詳細(xì)的教程參見另一篇分享:http://www.reibang.com/p/a4f497608485

homer

homer安裝

conda install homer

homer 軟件配置

mkdir ~software/homer && cd homer
wget http://homer.ucsd.edu/homer/configureHomer.pl

# conda安裝的Homer實際沒有包含參考序列或注釋數(shù)據(jù) 动知；但是可以使用 configureHomer.pl下載數(shù)據(jù)
perl /home/zzz/data/software/homer/configureHomer.pl -install mm10

Usage: configureHomer.pl [options]

Options:
        -list 列出可用的包
        -install 安裝Homer或homer需要用到的數(shù)據(jù)包
        -version 安裝homer包時，可以指定包版本
        -remove 移除包
        -update 更新所有包到最新版本
        -reinstall 強制重新安裝所有已經(jīng)安裝過的包
        -all 安裝所有包
        -getFacts (add humor to HOMER - to remove delete contents of data/misc/)
        -check 檢查第三方軟件：samtools, DESeq2, edgeR
        -make 重新配置和編譯可執(zhí)行文件
        -sun SunOS系統(tǒng)员辩，使用gmake 和 gtar代替make 和 tar
        -keepScript 不更新configureHomer.pl
        -url 安裝時盒粮，使用的資源地址，默認(rèn)：http://homer.ucsd.edu/homer/
        Hubs & BigWig settings (with read existing settings from config.txt if upgrading):
            -bigWigUrl <base urls for bigWigs> (Setting for makeBigWigs.pl)
            -bigWigDir <base directory for bigWigs> (Setting for makeBigWigs.pl)
            -hubsUrl <base urls for hubs> (Setting for makeMultiWigHub.pl)
            -hubsDir <base directory for hubs> (Setting for makeMultiWigHub.pl)
    Configuration files: 下載 update.txt奠滑，更新config.txt

homer 使用

cd  $data/peaks
for i in *.bed;
do
echo $i
file=$(basename $i )
sample=${file%%.*} 
echo $sample  
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $i >${sample}.bed  ## 將bed文件的第4列丹皱，第1列，第2列宋税，第3列打印出來摊崭，是為了讓bed文件符合homer軟件的輸入
findMotifsGenome.pl ${sample}.bed mm10 ${sample}_homer_motif -len 8,10,12
annotatePeaks.pl ${sample}.bed mm10  1>${sample}.peakAnn.xls 2>${sample}.annLog.txt 
done 

findMotifs常用參數(shù)：

-bg：自定義背景序列；
-size: 用于motif尋找得片段大小弃甥，默認(rèn)200bp爽室；-size given 設(shè)置片段大小為目標(biāo)序列長度；越大需要得計算資源越多淆攻；
-len：motif大小設(shè)置阔墩，默認(rèn)8,10,12；越大需要得計算資源越多瓶珊；
-S：結(jié)果輸出多少motifs, 默認(rèn)25啸箫；
-mis：motif錯配堿基數(shù)，默認(rèn)2bp伞芹；
-norevopp：不進(jìn)行反義鏈搜索motif忘苛；
-nomotif：關(guān)閉重投預(yù)測motif；
-rna: 輸出RNA motif唱较，使用RNA motif數(shù)據(jù)庫扎唾；
-h：使用超幾何檢驗代替二項式分布；
-N：用于motif尋找得背景序列數(shù)目南缓，default=max(50k, 2x input)胸遇；耗內(nèi)存參數(shù)
參考鏈接：http://www.reibang.com/p/93f45acff1f3

不僅僅找了motif，還順便把peaks注釋了一下汉形。得到的后綴為peakAnn.xls 的文件就可以看到和使用R包注釋的結(jié)果是差不多的（homer official web: http://homer.ucsd.edu/homer/ngs/peakMotifs.html）纸镊。

image.png

MEME尋找motif

需要通過bed格式的peaks的坐標(biāo)來獲取fasta序列倍阐。這個是在線的，我應(yīng)該不會去用逗威。MEME峰搪，鏈接：http://meme-suite.org/

參考鏈接：http://www.reibang.com/p/1384173c353b