今天的分享從conda安裝的bowtie2一直這樣報錯說起。
(chipseq) vip13t16@bw-X11DAi-N:~$ bowtie2 --help
/home/data/vip13t16/biosoft/miniconda/yes/envs/chipseq/bin/bowtie2-align-s: error while loading shared libraries: libtbb.so.2: cannot open shared object file: No such file or directory
(ERR): Description of arguments failed!
Exiting now ...然後我發現群裡這樣報錯的不止我一個人！
bowtie2報錯經過一番探索在簡書上找到這樣的帖子：Linux下bowtie2安裝(非conda)和配置
養成習慣，安裝在固定目錄下比如家目錄下，建biosoft文件夾，下面建各個軟體的文件夾
mkdir biosoft/bowtie2 && cd biosoft/bowtie2
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/bowtie2-2.3.5.1-linux-x86_64.zip/download#
unzip bowtie2-2.3.5.1-linux-x86_64.zip
#下面把bowtie2寫入bashrc,以便以後隨便哪個目錄都可以調用
vim ~/.bashrc
#最後一行加入
export PATH="$PATH:/home/data/vip13t16/biosoft/bowtie2/bowtie2-2.3.5.1-linux-x86_64/"按Esc鍵退出編輯:wq保存並退出，最後source ~/.bashrc
然後再構建索引進行比對
合併bam文件因為一個樣品分成了多個lane進行測序，所以在進行peaks calling的時候，需要把bam進行合併。
如果不用循環：
samtools merge control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam通常我們用批處理：
cd ~/project/epi/
mkdir mergeBam
source activate chipseq
cd ~/project/epi/align
ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ../mergeBam/$id.merge.bam $id*.bam ;done得到全新的bam文件如下：一共是9個。
merge.bam是否做PCR重複（要根據文章）PCR重複如果存在PCR重複，起始點、終止點都一樣，但是有兩條帶。在http://www.bio-info-trainee.com 有兩篇可以學習的帖子。
兩篇推文cd ~/project/epi/mergeBam
source activate chipseq
ls *merge.bam | while read id ;do (nohup samtools markdup -r $id $(basename $id ".bam").rmdup.bam & );done
# 針對新生成的.rmdup.bam文件構建索引、統計
ls *.rmdup.bam |xargs -i samtools index {}
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
使用macs2進行找peaksmacs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用實例
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01一般而言，我們照葫蘆畫瓢，按照這個實例替換對應部分就行了，介紹一下各個參數的意義
-f: -t和-c提供文件的格式，可以是」ELAND」, 「BED」, 「ELANDMULTI」, 「ELANDEXPORT」, 「ELANDMULTIPET」 (for pair-end tags), 「SAM」, 「BAM」, 「BOWTIE」, 「BAMPE」 「BEDPE」 任意一個。如果不提供這項，就是自動檢測選擇。-g: 基因組大小， 默認提供了hs, mm, ce, dm選項， 不在其中的話，比如說擬南芥，就需要自己提供了。-B: 會保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores-q: q值，也就是最小的PDR閾值， 默認是0.05。q值是根據p值利用BH計算，也就是多重試驗矯正後的結果。-p：這個是p值，指定p值後MACS2就不會用q值了。-m: 和MFOLD有關，而MFOLD和MACS預構建模型有關，默認是5：50，MACS會先尋找100多個peak區構建模型，一般不用改，因為你很大概率上不會懂。cd ~/project/epi/mergeBam
source activate chipseq
#如果不存在 ${id}_summits.bed這個文件
ls *merge.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.bam -t $id.merge.bam -f BAM -B -g mm -n $id --outdir ../peaks 2> $id.log &
fi
done
#對去PCR重複的再做一次
mkdir dup
mv *rmdup* dup/
cd dup/
ls *.merge.rmdup.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_rmdup_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.rmdup.bam -t $id.merge.rmdup.bam -f BAM -B -g mm -n ${id}_rmdup --outdir ../peaks 2> $id.log &
fi
done
jimmy大神的我的但是我對比了一下我和大神的結果，我的peak文件夾裡空空如也。是的，小菜狗又有哪裡出問題了(T~T)
打開了我的log於是然後去找了MACS的官方文檔，https://github.com/macs3-project/MACS，都已經MACS3了啊，看了一眼前幾天的推文：ChIP-Seq數據分析上下遊打通，人家也是用的MACS3，那就來下載吧~
MACS (3.0.0a6)UsageExample for regular peak calling on TF ChIP-seq:
macs3 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01Example for broad peak calling on Histone Mark ChIP-seq:
macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1Example for peak calling on ATAC-seq (paired-end mode):
macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01There are currently twelve functions available in MAC3 serving as sub-commands.
SubcommandDescriptioncallpeakMain MACS3 Function to call peaks from alignment results.bdgpeakcallCall peaks from bedGraph output.bdgbroadcallCall broad peaks from bedGraph output.bdgcmpComparing two signal tracks in bedGraph format.bdgoptOperate the score column of bedGraph file.cmbrepsCombine BEDGraphs of scores from replicates.bdgdiffDifferential peak detection based on paired four bedGraph files.filterdupRemove duplicate reads, then save in BED/BEDPE format.predictdPredict d or fragment size from alignment results.pileupPileup aligned reads (single-end) or fragments (paired-end)randsampleRandomly choose a number/percentage of total reads.refinepeakTake raw reads alignment, refine peak summits.callvarCall variants in given peak regions from the alignment BAM files.安裝&使用：（此處省略Python環境配置）
pip install MACS3
#下載成功會有：
#Successfully built MACS3
#Installing collected packages: cykhash, MACS3
#Successfully installed MACS3-3.0.0a4 cykhash-1.0.2
ls *merge.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_summits.bed ];
then
echo $id
nohup macs3 callpeak -c Control.merge.bam -t $id.merge.bam -f BAM -B -g mm -n $id --outdir ../peaks 2> $id.log &
fi
done
現在的文件在~/project/epi/peaks裡面
ls -lh *.bed |cut -d " " -f 5-
#以下是結果
0 8月 8 14:30 Control_summits.bed
69K 8月 8 14:30 H2Aub1_summits.bed
1.1M 8月 8 14:31 H3K36me3_summits.bed
1.2M 8月 8 14:31 Ring1B_summits.bed
1.9M 8月 8 14:30 RNAPII_8WG16_summits.bed
814K 8月 8 14:31 RNAPII_S2P_summits.bed
2.0M 8月 8 14:29 RNAPII_S5PRepeat_summits.bed
2.7M 8月 8 14:30 RNAPII_S5P_summits.bed
3.1M 8月 8 14:30 RNAPII_S7P_summits.bed下一次，將帶來激動人心的 可視化 部分的分享