1）自己实验获得的原始测序文件（.fastq.gz），测序文件格式为.fastq格式，而.gz格式是压缩格式，为了节省存储空间。
可以将测序文件修改为自己样品的名字
2）从NCBI的GEO数据库中下载sra文件，需要转换为.fastq文件。
$ wget -c ‘下载地址’ #wget下载数据，-c参数为断点续传
从GEO数据库中下载的.sra文件
#如果是单端测序文件就只需要fastq-dump命令将sra文件转换为fastq测序文件。
$ fastq-dump SRR10485695.sra
#如果是双端测序文件转换，需要加入–split-e参数，将.fastq测序文件转换为*_1.fastq和*_2.fastq。
$ fastq-dump --split-e SRR10485695.sra
#转换为fastq文件后需要压缩
$ gzip *.fastq
2.数据处理
1）trim_galore质量过滤，这一个非常流行的用于「去接头序列」的软件，用于处理高通量测序得到的原始数据。通常我们从测序公司拿到数据后，第一步就是评估数据的质量以及对raw data去接头处理。公司拿来的数据通常附带了clean data以及去接头的说明文件。
具体用法可见：
FelixKrueger/TrimGalore
?github.com/FelixKrueger/TrimGalore
参数说明：-o 输出到trim文件夹中；–cores 8线程数为8，默认是20 ；–fastqc对去完接头的数据进行fastqc质量控制。具体参数可以使用$trim_galore -help查看。
fastqc具体用法见：
FastQC A Quality Control tool for High Throughput Sequence Data
?www.bioinformatics.babraham.ac.uk/projects/fastqc/
$ trim_galore .fastq.gz -gzip -o trim --cores 8 --fastqc
2）bowtie2序列比对 (比对软件还有bwa和bowtie<和bowtie2不同>，这里使用的是bowtie2)。Bowtie 2是一个超快且内存高效的工具，用于将序列读取对齐到参考基因组上。参数说明：-x是使用bowtie2索引，索引位置会根据不同物种发生改变，索引文件在/home/reference/下。-p为线程数(需根据服务器的总线程数考虑)；-U为需要比对的文件；|为通道符，将前面的值赋予后面的命令；samtools sort是对bam文件进行排序，一般情况按照序列在fasta文件中的顺序和序列从左到右的位点排序；-@为设置线程数，同上-p；-O输出文件夹，-o输出文件名。
$ bowtie2 -x /home/reference/Dmelanogaster/UCSC/dm6/Sequence/Bowtie2Index/dm6 -p 30 -U _trimmed.fq.gz | samtools sort -@ 30 -O bam -o bam/_sort.bam
3）gatk去除重复。在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates。主要的重复来源于PCR的扩增，还有部分的测序产生的重复。去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。参数说明：-I为输入需要去除重复的样本。–ADD_PG_TAG_TO_READS为是否添加PG tag； --REMOVE_SEQUENCING_DUPLICATES为是否去除重复序列；–CREATE_INDEX为建造索引。-O为输出设定；-M为创造矩阵文本。
$ gatk MarkDuplicates -I _sort.bam --ADD_PG_TAG_TO_READS false --REMOVE_SEQUENCING_DUPLICATES true --CREATE_INDEX true -O bam/.rmdup.bam -M bam/.rmdup.matrix.txt
4）bamCompare格式转换。为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。
$ bamCompare -p 10 -b1 .rmdup.bam -b2 ._IN.rmdup.bam --skipNAs --scaleFactorsMethod readCount --operation subtract --outFileFormat bedgraph -o …/bw_data/.compareInput.bedgraph --extendReads 200
参数说明：-p为线程数；-b1为ip；-b2为input；–skipNA是跳过bin中没有覆盖的部分； --scaleFactorsMethod是选择缩放样本的方法，可选参数为：readCount/SES/ None； --operation默认情况下输出两个示例的log2比率，此脚本选用subtract为做差值；–outFileFormat设置输出文件格式，可选bigwig及bedgraph；-o为输出设置； --extendReads此参数允许将读取扩展为片段大小。通常不建议对拼接读取的数据(比如RNA-seq)使用此特性，因为它会在跳过的区域上扩展读取。默认参数为200。
5）compareinput to move0 to rpm.bedgraph：上一步做完差值后，可能会存在负值，所以这一步需要将其矫正为0，为之后的统计做准备。
$ awk ‘{if($4<0){$4=0};print}’ .compareInput.bedgraph > .move0.bedgraph
$ totalreads=awk '{sum+=$4}END{print sum}' <sample>.move0.bedgraph |echo 'reads left after remove duplication: '$totalreads
$ awk -v totalreads=$totalreads ‘{$4=$4/totalreads*1000000;print}’ .move0.bedgraph > .rpm.bedgraph
参数说明：awk是一种处理文本文件的语言，是一个强大的文本分析工具。行匹配语句 awk ’ ’ 只能用单引号。首先将compareInput.bedgraph中第四列小于0的值转换为0；sum+=$4表示sum=sum+$4，再输出sum=totalreads；-v是赋值一个用户定义变量。将第二步计算出的值赋予此命令中，将原来第四列的值比上totalreads再×1000000，最后输出到第四列(改变第四列的值)。
6）sort排序：sort将文件的每一行作为一个单位，相互比较，比较原则是从首字符向后，依次按ASCII码值进行比较，最后将他们按升序输出。参数说明：-k是指定列数；-k1,1 -k2,2n是先按照第一列排序，如果相同就按照第二列，n表示数字类型。n后面跟着r的，表示reverse。
$ sort -k1,1 -k2,2n .rpm.bedgraph > .sort.bedgraph
7）将bedgraph转换为bw文件。参数说明：根据dm6的染色体具体位置，将bedgraph转为bw。Bedgraph文件中的染色体位置是相对位置，需要实际染色体的文件，转换为绝对位置，即真实位置。
$ bedGraphToBigWig .sort.bedgraph /public/software/ucsc-tools/dm6.chrom.sizes .rpm.bw
8）bamCoverage标准化处理。将bam文件转为bw文件，同时进行标准化。参数说明：-p设置线程数；-b输入bam文件；–skipNAs同上；–normalizeUsing标准化处理设置为CPM，-o输出。
$ bamCoverage -p n T h r e a d ? b b a m / nThread -b bam/ nThread?bbam/{sample}.rmdup.bam --skipNAs --normalizeUsing CPM -o bw/${sample}.CPM.bw
9）call peak：支持bam/sam/bed等文件格式。
$ macs2 callpeak -t ._sort.bam -c ._sort.bam
参数说明：-t实验组，IP的数据文件；-c对照组，input。脚本会将输出的文件放在当前文件夹的peak文件夹中。
MACS2输出文件解读：
① NAME_peaks.xls——包含peak信息的tab分割的文件，前几行会显示callpeak时的命令。输出信息包含：染色体号、peak起始位点、peak结束位点、peak区域长度、peak的峰值位点(summit position)、peak 峰值的高度(pileup height at peak summit, -log10(pvalue) for the peak summit)、peak的富集倍数(相对于random Poisson distribution with local lambda)。<XLS里的坐标和bed格式的坐标还不一样，起始坐标需要减1才与narrowPeak的起始坐标一样>
② NAME_peaks.narrowPeak——1：染色体号；2：peak起始位点；3：结束位点；4：peak name；5：int(-10*log10qvalue)；6：正负链；7：fold change；8：-log10pvalue；9：-log10qvalue；10：relative summit position to peak start。
3.结果可视化
1）利用bedtools intersect取overlap （可选步骤！！ 这个是求两个peak的交集）
$ bedtools intersect -a .bed -b .bed -wa -wb > overlap_.txt
$ bedtools window -w n #利用设置一定的大小的window进行overlap，-w为设置window大小为n
参数说明：bedtools intersect为两个取交集；bedtools window是取一定大小的窗口，进行overlap，即可以存在最多窗口大小的差距再取overlap。-a -b是指定两个文件；-wa -wb是将两个文件中overlap的列出来；>输出到指定文件中。
2）利用bedtools intersect做peak注释。也可以用R进行peak注释，可以参考：
使用ChIPseeker进行peak注释_庐州月光的博客-CSDN博客
?blog.csdn.net/weixin_43569478/article/details/108079450
$ bedtools intersect -a stage6_peak.bed -b UCSC.D.melanogaster.dm6.repeats.gtf -wa -wb > overlap_repeat_anno.txt
参数说明：-a为peak文件；-b为注释文件。输出txt文件为注释结果文件。(-a -b的表格尽量不要有tittle，不然容易出bug)
3）peak数据制作矩阵(scale-regions/reference-point)
$ computeMatrix scale-regions -S .bw input..bw -R .bed -a 1000 -b 1000 -o ip_input.computeMatrix
$ computeMatrix reference-point --referencePoint center -b .bed -S .bw .bw --skipZeros -o .gz --outFileNameMatrix
参数说明：scale-regiuons mode简单来说会将给定bed file范围内的结合信号做一个统计（指的是一段长度），并将基因长度统一scale到设定regionBdoyLength的长度，加上统计基因上游和下游X bp的信号。reference-point mode则是给定一个bed file，以某个点为中心开始统计信号（TSS/TES/center）。但实际上我在尝试的时候regionBdoyLength参数也还是可以用的，所以估计和scale-regions区别也不是太大，主要是作图的一点区别。
4）作图plotProfile & plotHeatmap
$ plotHeatmap -m .gz -o heatmap_.pdf -zMin 0 --zMax 8 --colorMap coolwarm --missingDataColor 1
参数说明：-m输入矩阵；-o输出pdf文件；-zMin 0 --zMax 8是热图的最大值和最小值；–colorMap是选择热图颜色类型。做出的pdf热图可以直接在Ai里编辑。
ps. 温馨提示，具体操作都需要自己实践，实践的路上总会出现bug。以上内容如果有什么建议和意见，可以私信，欢迎讨论。