本文目录一览：
- 1、deeptools使用过程中出现问题集
- 2、ATAC-seq专题---生信分析流程
- 3、deepTools 的使用（二）
- 4、day15ChIP-seq开始学起来
- 5、我的ChIP-Seq(4):MAnorm差异分析
deeptools使用过程中出现问题集
1.想要画一个chip-seq的热图，首先要通过bamcompare获得针对IP和input样品的处理后的bigwig文件，
通过运行 命令
bamcompare -b1 IP. bam -b2 input.bam -o log2ratio_IP1.bw，
结果报错
“IP. bam does not appear to have an index. You must index the file first!”
在网上找了下解决方法是运行命令：
samtools index IP. bam
后，重新运行上述bamcompare命令
结果可行。
注：input. bam 也需要进行同样的处理。其实这个index文件也就是测序公司给的*.bam.bai文件。由于我这个分析是在服务器上运行的，上传的时候没有把bai文件也传上去，所以出现了这个问题。-wig文件chipseq
2.前面没有考虑到生物学重复的问题，如果ChIP有两个生物学重复，是不是要合并两个生物学重复后作图？这里就涉及了将两个生物学重复的BAM文件进行merge，需要使用samtools的merge命令：-wig文件chipseq
Usage: samtools merge [-nr] [-h inh.sam] out.bam in1.bam in2.bam[...]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to out.bam [in1.bam]
ATAC-seq专题---生信分析流程
ATAC-seq信息分析流程主要分为以下几个部分：数据质控、序列比对、峰检测、motif分析、峰注释、富集分析，下面将对各部分内容进行展开讲解。
下机数据经过过滤去除接头含量过高或低质量的reads，得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件，使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列，也可以不填这两个参数，软件会自动识别接头并进行剪切。如：-wig文件chipseq
fastp \
--in1 A1_1.fq.gz \ # read1原始fq文件
--out1 A1_clean_1.fq.gz \ # read1过滤后输出的fq文件
--in2 A1_2.fq.gz \ # read2原始fq文件
--out2 A1_clean_2.fq.gz \ # read2过滤后输出的fq文件
--cut_tail \ #从3’端向5’端滑窗，如果窗口内碱基的平均质量值小于设定阈值，则剪切
--cut_tail_window_size=1 \ #窗口大小
--cut_tail_mean_quality=30 \ #cut_tail参数对应的平均质量阈值
--average_qual=30 \ #如果一条read的碱基平均质量值小于该值即会被舍弃
--length_required=20 \ #经过剪切后的reads长度如果小于该值会被舍弃
fastp软件的详细使用方法可参考：。fastp软件对于trim结果会生成网页版的报告，可参考官网示例和，也可以用FastQC软件对trim前后的数据质量进行评估，FastQC软件会对单端的数据给出结果，如果是PE测序需要分别运行两次来评估read1和read2的数据质量。-wig文件chipseq
如：
fastqc A1_1.fq.gz
fastqc A1_2.fq.gz
FastQC会对reads从碱基质量、接头含量、N含量、高重复序列等多个方面对reads质量进行评估，生成详细的网页版报告，可参考官网示例：
经过trim得到的reads可以使用BWA、bowtie2等软件进行比对。首先需要确定参考基因组fa文件，对fa文件建立索引。不同的软件有各自建立索引的命令，BWA软件可以参考如下方式建立索引：
bwa index genome.fa
建立好索引后即可开始比对，ATAC-seq推荐使用mem算法，输出文件经samtools排序输出bam：
bwa mem genome.fa A1_clean_1.fq.gz A1_clean_2.fq.gz
| samtools sort -O bam -T A1 A1.bam
值得注意的是，在实验过程中质体并不能完全去除，因此会有部分reads比对到质体序列上，需要去除比对到质体上的序列，去除质体序列可以通过samtools提取，具体方法如下：首先将不含质体的染色体名称写到一个chrlist文件中，一条染色体的名称写成一行，然后执行如下命令即可得到去除质体的bam-wig文件chipseq
samtools view -b A1.bam $chrlist A1.del_MT_PT.bam
用于后续分析的reads需要时唯一比对且去重复的，bwa比对结果可以通过MAPQ值来提取唯一比对reads，可以用picard、sambamba等软件去除dup，最终得到唯一比对且去重复的bam文件。-wig文件chipseq
比对后得到的bam文件可以转化为bigWig（bw）格式，通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件，可以用以下命令实现bam向bw格式的转换：-wig文件chipseq
bamCoverage -b A1.bam -o A1.bw
此外，可以使用deeptools软件展示reads在特定区域的分布，如：
computeMatrix reference-point \ # reference-pioint表示计算一个参照点附近的reads分布，与之相对的是scale-regions，计算一个区域附近的reads分布-wig文件chipseq
--referencePoint TSS \#以输入的bed文件的起始位置作为参照点
-S A1.bw \ #可以是一个或多个bw文件
-R gene.bed \ #基因组位置文件
-b 3000 \ #计算边界为参考点上游3000bp
-a 3000 \ #计算边界为参考点下游3000bp，与-b合起来就是绘制参考点上下游3000bp以内的reads分布
-o A1.matrix.mat.gz \ #输出作图数据名称
#图形绘制
plotHeatmap \
-m new_A1.matrix.mat.gz \ #上一步生成的作图数据
-out A1.pdf \ # 输出图片名称
绘图结果展示：
MACS2能够检测DNA片断的富集区域，是ATAC-seq数据call peak的主流软件。峰检出的原理如下：首先将所有的reads都向3'方向延伸插入片段长度，然后将基因组进行滑窗，计算该窗口的dynamic λ，λ的计算公式为：λlocal = λBG（λBG是指背景区域上的reads数目），然后利用泊松分布模型的公式计算该窗口的显著性P值，最后对每一个窗口的显著性P值进行FDR校正。默认校正后的P值（即qvalue）小于或者等于0.05的区域为peak区域。需要现在linux环境中安装macs2软件，然后执行以下命令：-wig文件chipseq
macs2 callpeak \
-t A1.uni.dedup.bam \ #bam文件
-n A1 \ # 输出文件前缀名
--shift -100 \ #extsize的一半乘以-1
--extsize 200 \ #一般是核小体大小
--call-summits #检测峰顶信息
注：以上参数参考文献（Jie Wang，et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications）-wig文件chipseq
ATAC分析得到的peak是染色质上的开放区域，这些染色质开放区域常常预示着转录因子的结合，因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例，首先在linux环境中安装homer，然后用以下命令进行motif分析：-wig文件chipseq
findMotifsGenome.pl \
A1_peaks.bed \ #用于进行motif分析的bed文件
genome.fa \ #参考基因组fa文件
A1 \ #输出文件前缀
-size given \ #使用给定的bed区域位置进行分析，如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析
homer分析motif的原理及结果参见：
根据motif与已知转录因子的富集情况可以绘制气泡图，从而可以看到样本与已知转录因子的富集显著性。
差异peak代表着比较组合染色质开放性有差异的位点，ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount，可以选择edgeR、DESeq2等模型进行差异分析。-wig文件chipseq
在科研分析中我们往往需要将peak区域与基因联系起来，也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例，需要在R中安装ChIPseeker包和GenomicFeatures包，然后就可以进行分析了。-wig文件chipseq
library(ChIPseeker)
library(GenomicFeatures)
txdb- makeTxDbFromGFF(‘gene.gtf’)#生成txdb对象，如果研究物种没有已知的TxDb,可以用GenomicFeatures中的函数生成
peakfile -readPeakFile(‘A1_peaks.narrowPeak’)#导入需要注释的peak文件
peakAnno - annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)
# 用peak文件和txdb进行peak注释，这里可以通过tssRegion定义TSS区域的区间
对于peak注释的结果，也可以进行可视化展示，如：
p - plotAnnoPie(peakAnno)
通过注释得到的peak相关基因可以使用goseq、topGO等R包进行GO富集分析，用kobas进行kegg富集分析，也可以使用DAVID在线工具来完成富集分析。可以通过挑选感兴趣的GO term或pathway进一步筛选候选基因。-wig文件chipseq
deepTools 的使用（二）
此工具基于multiBamSummary或multiBigwigSummary的输出文件。主要用于同一批样本之间的相关性比较，查看是否差异很大，也可做聚类图，来观察哪些样本相似。需求不大，所以先不讲。-wig文件chipseq
此工具主要用于Chip-seq数据，通过比较input和实验组观察实验质量。因为Chip实验有高达90%的假阳性。
----numberOfSamples 参数主要是随机从样本中抽取的箱数来计算相对的覆盖reads数，-bs 可以设置每个箱的base数。
做出的图如上图，大概有三种情况。横坐标是按照bin里read数排序，纵坐标代表bin里reads数的总和百分比。第一个图中，在排序前97%的bin中所有reads数只占了55%，所以有3%的bin中包含了45%的reads。也就是说H3K4me3在某些地方强烈富集，正是我们想看到的。第二个图也很美，而第三个图很难区分input和实验组有什么区别，是个失败的实验。-wig文件chipseq
该工具可用于评估给定样品的测序深度。它采样100万bp，计算重叠读数的数量，并可以报告一个直方图，告诉你有多少碱基被覆盖多少次。接受多个BAM文件，但它们都应该对应于相同的基因组组装。
此工具可以看两个样本的覆盖深度相关性。
此工具主要估计片段长度和频率。
计算GC偏差。
day15ChIP-seq开始学起来
linux学的半半拉拉，觉得需要用实战来激励一下，之后再回过头还要再学linux的哦。所以先来ChIP-seq吧，这个挺急需的，去年的数据估计也需要用这一套流程来跑。感谢Jimmy大神无私地录教程。 -wig文件chipseq
raw data转换，quality control。
call peak-通常是用MACS2来找到peaks。有的可以IGV可视化。
包括结合峰相关基因的注释、GO分析KEGG通路富集分析和motif分析，主要是 motif和peaks的注释。
包括差异结合峰相关基因的注释、GO分析、KEGG通路富集分析和motif分析。即多次重复的相关性，或者处理前后peaks差异的对比。
还能对目的蛋白的调控机制做进一步的深入分析。
从知乎大神-生信宝典这里学习到
这些之前都知道啦。bam和sam文件可以帮助我们探索reads在参考基因组中的比对情况。
只包含区域和区域的覆盖度信息，文件更小，用于可视化更方便，可以导入基因组浏览器（Genome Browser）中进行可视化，以查看reads在参考基因组各个区域的覆盖度并检测测序深度。这几个文件在[ChIP-seq数据分析]、Call Peak阶段会生成，可以利用[IGV]等工具进行可视化，方便查看组蛋白修饰的程度。-wig文件chipseq
wig文件包括染色体长度，步长多少，span多少
bw文件，bigWig是 wig文件的二进制压缩格式，可通过wigToBigWig工具转换。bw文件可以直接导入IGV。
bedGraph: bed与wig的结合，更省空间和灵活，展示信息与wig类似。
分析过程中的bed文件一般代表区域信息，如表示Peak位置的bed文件，表示基因注释的bed12文件。
下文可以详读。
Pedro Rosmaninho et al, Zeb1 potentiates genome‐wide gene transcription with Lef1 to promote glioblastoma cell invasion, The EMBO Journal (2018). DOI: 10.15252/embj.201797115-wig文件chipseq
我的ChIP-Seq(4):MAnorm差异分析
哈哈，搜了一圈没发现网上有关于MAnorm的中文教程或者是说明，本文将是第一篇~撒花✿✿ヽ(°▽°)ノ✿那就要用心写了，感到鸭梨.jpg==
首先，MAnorm是什么，可以做什么呢？
简单地说，这是一款寻找两个ChIP-Seq样本之间差异peak的软件。一般ChIP的流程中，若是单一处理的细胞系，那么callpeak之后可能会做binding motif的分析或是peak相关gene的功能分析等；但若是两种处理的细胞系（比如饥饿组和对照组），我们肯定想要知道两种处理下，组蛋白修饰的差异，类似于RNA-Seq中差异表达基因的分析，所以这时就需要进行差异分析。MAnorm就可以实现这样的分析需要。-wig文件chipseq
一般来说，上述差异分析不一定要在peaks水平进行，完全可以在reads水平，这个就叫做“一步法”；而通过先分别callpeak再比较peaks的density或者depth等，就是所谓的“两步法”。不同方法有不同类型的软件可供选择，这就是ChIP分析成熟的地方，不过技术流大可根据自己的目的写脚本进行个性化处理，这个暂且不表。-wig文件chipseq
那么差异分析软件如何选择呢？根据组蛋白修饰类型、样品是否有重复、是否需要callpeak（即predefined region set），下图一目了然：
我的样品有宽峰窄峰两种修饰、无重复，项目时间紧张尽量想用一个软件实现，所以选择了MAnorm。
话不多说，直接看图：
1.1.4版本 ：
conda/PyPi
需要注意的是，此版本只支持bed格式且不支持paired-end模式，会把所有reads当成single-end处理。若reads文件想用支持更多的格式（sam/bam/bedpe等），请用v1.2.0。-wig文件chipseq
1.2.0版本 ：
暂时只能从Github复制源码进行安装。方法：
建议首先阅读使用说明，最好从linux中 manorm --help ，或者在Github中找到相应版本的 附带说明 ，这一点很重要，因为有时网上搜到的说明和你实际用的版本不一致，会走弯路，不要问我咋知道的。-wig文件chipseq
所以要准备的文件有4个：
将如上文件移动至新文件夹下待用。***tips：这里不再需要对照组In的文件了
基本命令（--p1 --p2 --r1 --r2 -o是5个必需参数，注意是两个-）：
建议试运行一组数据先，根据报错文件调整格式。软件还不太成熟，需要多调整格式。
运行约10min，产生4个结果文件：
sample1peaks_vs_sample2peaks_all_MAvalues.xls ：这个是主要的结果文件，Excel格式，里面的peak_group有标注是common/1unique/2unique的。-wig文件chipseq
output_figures 文件夹 ：4个图，计算的Mvalue Avalue(MA)及校正之后的MA，大概就是这个意思，还需要读文献琢磨
output_filters 文件夹 ：3个peaks.bed文件，可能就是条件严格了点之后的结果，两个biased包括的peaks很少，一个unbiased包括的peaks很多跟all那个文件差不了多少。-wig文件chipseq
output_tracks 文件夹 ：3个wig文件，是M A values的，UCSC可视的文件类型。
综上，决定用main output file即第一个结果，进行后面的分析。
