Ⅰ chapter1.高通量序列实验简介：设计与生物信息学分析
2021/4/16
1、设计高通量测序实验的步骤
2、介绍了最广泛使用的应用，并描述了基本的测序概念。
3、可用于生物信息学分析的各种软件程序，以理解测序数据。
1、Insert ：用于测序的DNA片段
2、Read ： Insert 被测序到的部分
3、Single Read（SR） ：一种只从 Insert 序列一端测序的测序程序
4、Pair Read（PR） ：一种从 Insert 序列两端测序的测序程序
5、Flowcell ：连接DNA芯片并进行测序的一种小玻璃芯片。 Flowcell 被探针覆盖，允许与DNA片段连接的接头杂交。
6、Lane ： Flowcell 由8个物理分离的通道组成，称为 Lane 。在所有 Lane 上并行进行测序。
7、Multiplexing/Demultiplexing ：在同一 Lane 上对几个茄燃样本进行测序称为多路复用 Multiplexing ，在一条 Lane 上测序的 Reads 的分离称为分路复用 Demultiplexing ，通过一个识别每个Reads*的索引将其与已知样本的索引进行比较。
8、Pipeline ：一系列的计算过程
（一）reading
1、Resequencing ：在一个给定的样本中找到相对于参考基因组的变体
实验细节 ：从相关细胞中提取DNA，进行由DNA碎片化和测序组成的样品制备
基本分析总结 ：将序列片段映射到参考基因组，并通过总结片段与其基因组位点的差
异来识别相对于参考基因组的变异对应的“地图”
2、Target-enriched sequencing ：靶点富集测序是一种特定的 Resequencing 形式，只
关注特定的基因组基因座。
实验细节 ：在从细胞中提取DNA并进行样品制备后，进行一个富集过程来捕获相关的
位点，靶富集可以使用“定制的”靶富集探针在基因组的特定区域进行，或
使用可用的试剂盒，如exome-enrichment kits。
基本分析总结 ：与 Resequencing 相同
3、De novo assembly ：识别一个基因组序列，而无需任何额外的参考
实验细节 ：与 Resequencing 相同
基本分析总结 ：组装过程依赖于DNA片段的重叠。这些重叠被合并成一致序列，称为
contigs 和 scaffolds 。
（二）counting
1、ChIP-Seq/RIP-Seq ：找到RNA或DNA结合蛋白的结合位置
实验细节 ：（1）首先，进行了ChIP/RIP实验：蛋白质与DNA/RNA结合，并与之交
联。然后DNA/RNA被分裂。
（2）蛋白质 pull down 经历免疫沉淀过程，交联被逆转
（3）对富集于蛋白结陆纳虚合位点中的DNA/RNA片段进行测序
基本分析总结 ：被序列排列的片段被映射到基因组中。基因组中丰富的位置是通过检
测基因组的映射片段的“ peaks (峰)”发现的,这些峰值应该明显高于在
周围的位点中已映射的片段，并且与对照样本相比要高得多----通常
是ChIP实验的输入DNA或其他由非特异性抗体进行的免疫沉淀样
本。
2、RNA-Seq ：检测和比较基因表达水平
实验细节 ：从细胞中提取总RNA，在样品制备过程中，mRNA被 pull down 并破碎。
然后，mRNA片段被逆转录成cDNA，cDNA片段测序。
基本分析总结 ：cDNA片段被映射到参考基因组中。映射到每个基因的片段被计数和
标准化，以便比较不同的基因和不同的样本。在一个RNA-Seq实验
中，通过检测映射到一个未注释区域的基因组上的片段束，可以找到
未标记的基因和转录本。
（三）reading/counting
microRNA-Seq ：检测和计数 microRNAs
实验细节 ：从细胞中提取总RNA，通过识别大多数已知的microRNA分子共同的自然
早燃 结构来分离microRNA，然后对microRNA片段进行逆转录和测序。
基本分析总结 ：被测序的片段被映射到基因组中，然后，微RNA可以被检测和计数。
1、在 reading 中，覆盖范围对应于 平均覆盖基因组中每个碱基 的 reads 数量。
一般来说， 30X覆盖率 被认为是识别基因组变异的最小值，而 de novo 通常需要一个更高的覆盖范围。
2、在 counting 中，覆盖的概念并不简单，因为the number of reads along the genome is not expected to be uniform.
可帮助评估是否有足够的reads序列的分析是“ *saturation report* (饱和度报告)”，使用所有的reads确定表达水平，表达水平与 取一部分reads重新计算的表达水平 比较。
1、基因组的重复性
要唯一地对重复区域的read映射进行评分，它必须 比 重复区域 或 边界相邻的非重复序列 更长。更长的reads或PE reads允许“拯救”非唯一端，也映射到基因组中的非唯一区域。
2、差异剪接变异
同一基因表达的转录本不同时：
3、测序样本与参考基因组的遗传距离
如果被测序的样本与参考基因组有遗传距离，可能需要更长的reads来确定基因组中每个read的来源。
4、寻找结构变异
基因组的结构变化，如长的插入或缺失，倒位和易位可以通过Paired-End信息找到。
5、De Novo 装配
挑战：测序错误、低复杂度区域和重复区域等
更长的PE reads会导致更好的装配，使用一些具有不同insert length的序列库可以改进组装过程。
1、 Resequencing：有遗传距离。。。
2、RNA-Seq：使用来自不同重复的数据，并将其合并为一个具有更高统计显著性的值。
3、ChIP-Seq：+控制样本
1、Raw Data 处理
此步骤的可用软件：Illumina’s CASAVA software，Illumina运行会生成“base-calling”文件(*.bcl)，它们只有在转换为通用fastq格式时才会非常有用,在此文件转换过程中，还执行解复用过程，即从同一lane上排序的不同样本分离读取。
2、质量控制和read操作
此步骤的可用软件：CASAVA和FastQC
测序运行完成后，在开始分析之前，应检查运行的质量是否以下参数，这些参数可能说明样品和运行的质量。
3、为De Novo Assembly组装 Contigs 和 Scaffolds
此步骤的可用软件：SOAPdenovo，ABySS，Velvet，ALL-PATHS
4、mapping
此步骤的可用软件：BWA ，Bowtie，TopHat
5、 Variant Calling and Filtering
此步骤的可用软件： SAMtools，GATK，MAQ
帮助检测变异的两个基本参数如下：
（1）Coverage at the loci
（2）被测序的等位基因的频率
6、Assembling Transcripts
7、 Gene Expression Analysis
此步骤的可用软件：Cufflinks，Myrna
一种常见的归一化方法FPKM，计算如下：
8、 Peak Detection
此步骤的可用软件：MACS，SICER
Ⅱ VarBen：NGS变异数据模拟软件
2021年3月3号，来自 国家卫生健康委临床检验中心、中国科学院计算技术研究所等单位的研究人员开发了肿瘤突变数据模拟软件VarBen，包括SNV、InDel、CNV、SV等多种变异形式的模拟。同时该软件支持目前临床上常用的测序平台，包括Illumina、华大智造MGISEQ以及Ion torrent测序平台。
据介绍，VarBen采用的是比对到参考基因组特定位点的测序reads进行编辑的方式来进行突变模拟，该方法可保留测序过程“湿实验”部分核酸提取、靶向捕获、文库制备以及测序过程中产生的错误分布模式，从而保证模拟数据更加的接近真实。研究人员根据不同高通量测序平台的原理以及不同种类基因变异的特点，使得VarBen能够对几乎所有的变异类型进行模拟，包括SNV、Indel、Complex insertion-deletion、CNV以及SV。此外，VarBen同时支持全基因组、全外显子组以及靶向panel测序数据的模拟，并且适用于多个测序平台。
文献名称： VarBen: Generating in Silico Reference Data Sets for Clinical Next-Generation Sequencing Bioinformatics Pipeline Evaluation
期刊名称：The Journal of Molecular Diagnostics （国际分子诊断领域权威期刊） ( IF 5.553 )(JMD)
发表时间：2020-12-18
文章介绍1： https://www.x-mol.com/paper/1340047108494508032
文章介绍2： https://www.sohu.com/a/453971998_120054289
github： https://github.com/nccl-jmli/VarBen
数据路径： SRP171108
软件亮点（和其他软件相比）：
测试软件：
注意：bwa、samtools、bedtools 要加入环境变量中，流程里面直接使用的时bwa，samtools、bedtools，没有给传软件路径的参数
最终生成的bam文件存储在 outdir/edit.sorted.bam
i)从BAM文件中根据用户指定的变异等位基因分数(VAF)的指定位置或区域随机选择Read；
ii)根据指定的变体类穗旁型(包括)修改、插入或删除重叠位置或区域的Reads中的SNV和InDel；
iii)编辑后的Reads被重新映射到参考基因组以获得适当的比对信息并与从原始输入BAM中剩下的的Read进行合并。
1） 获得单倍型 "step1: deal with mutfile and get haplotypes"：varben.deal_mut.checkMutInput.get_haplotypes
① 判悔轿断突变类型 varben.deal_mut.checkMutInput.check_haplotype -> varben.deal_mut.checkMutInput.check_mut_info
2）获得序列： "step2: deal haplotypes and get total_chosen_reads, total_chosen_reads_muts"
varben.deal_mut.dealHaplotype.deal_haplotype_multi
1）读碧族肆取结构变化(SV)模拟的类型分类：
类型1：Read2与断点A重合或者Read1与断点B重合
类型2：Read1与断点A重合或者Read2与断点B重合
类型3：断点落在配对Reads中间
类型4: 配对Read完全落在两个断点之间
类型5：断点同时与Read1和Read2重合
2）SV变异模拟：
i）根据BAM文件中指定的VAF，选择与指定断点有Overlap的Read；
ii）选择的Read根据删除区域分为五种类型，这些Read被编辑以产生 soft-clipper Reads和另一端Read；
iii）编辑后的Reads重新进行比对；
iv）合并重新比对的Read和的原始BAM中的剩余Read；
0）检测输入文件格式 "step0: prepare sv list" ：varben.deal_sv.checkSVInput.check_sv_file
1）获得插入序列 "step1: get insert size of paired reads"：varben.common.methods.get_insertSize_range
都是 [100, 1000]
2）获得编辑序列："step2: deal with sv and get total edited reads"
①varben.deal_sv.dealSVType.deal_sv -->
②varben.deal_sv.dealSVType.deal_cnv -->
③varben.deal_sv.getReadsType.pos_type_classify --> ###里面设置了process=20,
④varben.deal_sv.getReadsType.posType_sub_paired -->
varben.deal_sv.dealSVType.get_write_reads --> ###输出需要增加或者删除或者修改的序列ID
3）将需要处理的序列写入Bam ："step3: get reads by region bed and write bam file"
①varben.deal_sv.getReadsByRegion.get_reads_by_region -->
②varben.common.bamconvert.getRegionReads --> ###生成bam
samtools view raw.bam -b -h
-o cnv_out/tempDir/used_tmp.bam
-U cnv_out/tempDir/exclude_tmp.bam
-L cnv_out/tempDir/consider_region.bed
③varben.common.bamconvert.bamIndex ---> ###创建bam索引
④varben.deal_sv.writeBamByChr.write_sub_bam
4）合并编辑后的read到原始文件中，并重新比对生成新的bam文件："step4: merge edited reads and remap to new bam, consider about the tag, RG, life reads"
varben.deal_sv.mergeEditBam.merge_edit_bam
5）合并bam文件："step5: remapped edit reads and merge"
varben.common.bamconvert.remap --> varben.common.bamconvert.bamToFastq
(edit.bam -> fq1,fq2 -> edit.rmap.bam -> edit.remap.sort.bam -> edit.sort.bam)
Ⅲ DNA/RNA序列比对软件整理
文章仅是记录自己的学习使用，有错误请指出，我立刻改正
在对比对工具进行比较时，通常将其分为DNA比对工具（DNA-seq）和RNA比对工具（RNA-seq）。它们的区别在于是否会考虑跨外显子的比对，即：是否会将没有比对上的reads劈开，对劈开后的两部分再次比对）。
随着现在各种seq测序的出现，我们已经不能简单的根据是比对DNA还是RNA来判断。比对工具的选择主要依据reads的比对是否需跨外显子。（PRO-seq/GRO-seq，它们虽然在建库时捕获的RNA，但是它们的比对并不需要考虑跨外显子。）
常用工具：
DNA-seq：BWA；bowtie&bowtie2
RNA-seq：STAR；HISAT2；Tophat&Tophat2
BWA主要应用二代测序后的大量短小片段与参考基因组之间的定位比对。需要先嫌燃对参考序列建建立索引，BWA也是基于 BWT和 FM-Index 理论来对参考基因组做索引。根据测序方法的不同，有单末端序列(Single-end，SE)比对和双末端序列(Pair-end，PE)比对。
bowtie出现在测序行业还不成熟的时候，序列长度普遍在50bp以下，bowtie的只满足长度在50bp以下的reads的比对。官方称其可以把短的DNA序列（35bp）快速的比对到人类基因组上。
Bowtie2 是一款经典的短读长序列（ 50-100 bp，最多可到1000 bp ) 比对软件，节约内存且灵活与成熟的短序列比对软件，比较适合下一代测序技术。支持单端测序（unpaired) 和双端测序的比对。支持全局比对（end-to-end align ) 和 局部比对( local align )。其通常使用全文分索引(FM-index)以及Burrows-Wheeler 变换(BWT)索引基因组使得比对非常快速且内存高效，但是这种方法不适合于找到较长的、带缺口的序列比对
结论：bowtie和bowtie2，是两个不同类型的比对工具，bowtie2并非是bowtie的升级。尺有所长寸有所短，bowtie适合长度在50b长度以内的reads比对，而bowtie2适合50-100b，甚至更长的reads比对。但是这两个都属DNA-seq比对工具
RNA-Seq测序的特性，天然的会有一部分数据延伸到内含子区，这部分跨越外显子和内含子的reads就称为『junction reads』，所以RNA-Seq比对软件需要针对此进行优化。
（ junction：转录组reads比对不同于基因组reads比对(如ChIP-seq、WES等)的地方在于，比对的reads可能来源于2个被内含子隔开的外显子区域，导致reads一端比对在第一个外显子的后面部分，另一端比对在第二个外显子的前面部分，芹衡虚即跨剪切位点，从而形成exon-exon junction (剪接点)。这些reads又称为junction reads，对转录本的拼接、鉴定和差异分析具有重要的意义。）
（soft-clip事件: 即reads末端存在低质量碱基或接头导致比对不上的， STAR会自动尝试截去未比对部分，只保留比对上的部分。）
STAR是ENCODE皇家御用的RNA-seq比对工具，ENCODE计划(ENCyclopedia Of DNA Elements)又称人类基因组DNA元件网络全书计划，是2003年在人类基因组计划完成之后紧接着的又一个大型国际科研项目。
Tophat2的原作者们也不知道是出于什么考虑，不再更新Tophat2，转而开发了一个新的比对工具HISAT2，更是推荐人们使用HISAT2，声称其速度更快，内存占用率更小，准确率更高。
此外，HISAT2不仅支持RNA-seq的比对还支持DNA-seq比对，唯一需要做的就是加上一个参数--no-spliced-alignment。但是就目前拦散来看，大部分人都是使用HISAT2做RNA-seq，没人使用它做DNA-seq
Tophat/Tophat2工具本身不能进行比对，它是通过调用bowtie/bowtie2进行比对的。划重点，bowtie2不是bowtie的升级版，但是Tophat2是Tophat2的升级版。因此Tophat只可以调用bowtie，而Tophat2不仅可以调用bowtie2（默认）还可以更改设置调用bowtie。
Tophat/Tophat2调用bowtie/bowtie2后，会首先使用bowtie/bowtie2对序列进行比对，对于那些没有比对上的，会考虑其跨外显子的可能性，将reads劈开重新比对。
全长转录组（Full-length transcriptome）是基于PacBio和Nanopore三代测序平台，无需打断拼接，直接获得包含5’UTR、3’UTR、polyA尾的mRNA全长序列及完整结构信息，从而准确分析有参考基因组物种可变剪接及融合基因等结构信息，克服无参考基因组物种转录本拼接较短、信息不完整的难题。同时还可以借助二代测序数据，进行转录本特异性表达分析，获得更加全面的注释信息。
传统的使用比较多的长读长比对软件是GMAP，05年发表公布，最开始是用来比对低通量的est序列的，后来也有进一步升级为GSNAP支持高通量的二代测序。PacBio测序技术出现后，常用于Iso-seq转录本的鉴定，目前仍是相关研究引用量最高的比对软件，该软件也一直在持续更新升级。其可以将转录本序列与参考基因组序列比对，输出gff文件，比对速度稍慢。
Minimap2是生信大牛李恒18年用C语言开发的可以用于三代数据（subreads、iso-seq）比对的长序列比对软件，与传统的三代比对工具GMAP相比，其速度有非常显著的提升，当然同时消耗的内存也比较大。使用方法也比较简单，近几年引用次数增长的也很迅速，所以大家可以试试用minimap2进行Iso-seq的比对。
Ⅳ 如何利用python进行高通量测序数据的分析 知乎
建议读两本书：
1、集体智慧编程 (豆瓣)
因为Python是一门不需要花太多精力（甚至可以说很少），就可以基本掌握的一门语言，所以推荐这本书。题主提到以后想学机器学习，这是一本非常好的入门书，书中的例子源码都是Python实现的，并且能帮你迅速熟悉Python相关的各种计算库。
2、统计学习方法 (豆瓣)
考虑到题主要学得踏实，这本书深入浅出地讲了和机器学习有关的一切数学基础知识，一整本的干货，没有废话，非常值得一读。题主数学专业的话，读起来应该会比我更顺畅。
Ⅳ 以下哪些工具可以获得序列标识图
序列标识图是指将蛋白质序列上的结构和功能信息映射到一个二维图上，以便更好地理解和分析蛋白质序列的结构和功能。以下是一些可以获得序列标识图的工具:
1. InterProScan: InterProScan是一种基于序列的功能注释工具，可以预测蛋白质序誉察列的结构和功能，并输出序列标识图。它使用多种算法和数据库，包括Pfam、PRINTS、ProSite、SUPERFAMILY等，以提高预测准确性。
2. Protter: Protter是一种在线的序列标识图生成工具，它可以将FASTA格式的蛋白质序列映射到一个可视化的二维图上，显示出蛋白喊饥质的各种结构和功能信息，包括域、保守区域、跨膜区域、信号肽等。
3. ESPript: ESPript是一种在线的序列标识图生成工具，它可以将蛋白质序列映射到一个二维图上，以显示出蛋白质的二级结构、保守区域、氨基酸残基间的相互作用庆渗茄等信息。它还支持多种输出格式和样式，用户可以根据需要进行调整。
4. WebLogo: WebLogo是一种在线的序列标识图生成工具，它可以将多个蛋白质序列的保守区域映射到一个二维图上，以显示出每个位置上的氨基酸分布情况，并计算每种氨基酸出现的频率。它还支持多种输出格式和样式，用户可以根据需要进行调整。
总之，以上工具都可以用于生成蛋白质序列标识图，用户可以根据自己的需要选择合适的工具，并根据输出结果进行分析和解释。
Ⅵ ATAC-seq专题---生信分析流程
ATAC-seq信息分析流程主要分为以下几个部分：数据质控、序列比对、峰检测、motif分析、峰注释、富集分析，下面将对各部分内容进行展开讲解。
下机数据经过过滤去除接头含量过高或低质量的reads，得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件，使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列，也可以不填这两个参数，软件会自动识别接头并进行剪切。如：
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软件的详细使用方法可参考：https://github.com/OpenGene/fastp。fastp软件对于trim结果会生成网页版的报告，可参考官网示例http://opengene.org/fastp/fastp.html和http://opengene.org/fastp/fastp.json，也可以用FastQC软件对trim前后的数据质量进行评估，FastQC软件会对单端的数据给出结果，如果是PE测序需要分别运行两次来评估read1和read2的数据质量。
如：
fastqc A1_1.fq.gz
fastqc A1_2.fq.gz
FastQC会对reads从碱基质量、接头含量、N含量、高重复序列等多个方面对reads质量进行评估，生成详细的网页版报告，可参考官网示例：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html
经过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
samtools view -b A1.bam $chrlist > A1.del_MT_PT.bam
用于后续分析的reads需要时唯一比对且去重复的，bwa比对结果可以通过MAPQ值来提取唯一比对reads，可以用picard、sambamba等软件去除p，最终得到唯一比对且去重复的bam文件。
比对后得到的bam文件可以转化为bigWig（bw）格式，通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件，可以用以下命令实现bam向bw格式的转换：
bamCoverage -b A1.bam -o A1.bw
此外，可以使用deeptools软件展示reads在特定区域的分布，如：
computeMatrix reference-point \ # reference-pioint表示计算一个参照点附近的reads分布，与之相对的是scale-regions，计算一个区域附近的reads分布
--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软件，然后执行以下命令：
macs2 callpeak \
-t A1.uni.dep.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）
ATAC分析得到的peak是染色质上的开放区域，这些染色质开放区域常常预示着转录因子的结合，因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例，首先在linux环境中安装homer，然后用以下命令进行motif分析：
findMotifsGenome.pl \
A1_peaks.bed \ #用于进行motif分析的bed文件
genome.fa \ #参考基因组fa文件
A1 \ #输出文件前缀
-size given \ #使用给定的bed区域位置进行分析，如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析
homer分析motif的原理及结果参见：http://homer.ucsd.e/homer/motif/index.html
根据motif与已知转录因子的富集情况可以绘制气泡图，从而可以看到样本与已知转录因子的富集显著性。
差异peak代表着比较组合染色质开放性有差异的位点，ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount，可以选择edgeR、DESeq2等模型进行差异分析。
在科研分析中我们往往需要将peak区域与基因联系起来，也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例，需要在R中安装ChIPseeker包和GenomicFeatures包，然后就可以进行分析了。
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进一步筛选候选基因。