一套H3K4me3的数据，GSE167940，单端single，无input!!有了！一翻邮件交流发给我了
一套H3K27ac的数据，GSE165809，单端single，有input
H3K4me3
H3K27ac
1、使用SRA-Toolkit转换SRA为fastq并压缩
module load SRA-Toolkit
fastq-dump SRR*(通配符表示对该文件夹下所有SRR开头的文件操作)
#解压这步的高级用法，fastq-dump默认最多6线程。#fasterq-dump -e 30 SRR*
gzip *.fq
#压缩也有高剂用法，多线程软件pigz，#pigz -p 30 *fastq
用的是SRATools
2、Trim_galore质控及去接头
首先判断片段长度
zcat SRR.fastq.gz |head -n 10
大差不差都50
判断首先判断pread类型
less $1 | head -n 1000 | awk '{if(NR%4==0) printf("%s",$0);}' \
| od -A n -t u1 -v \
| awk 'BEGIN{min=100;max=0;} \
{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END \
{if(max<=126 && min<59) print "Phred33"; \
else if(max>73 && min>=64) print "Phred64"; \
else if(min>=59 && min<64 && max>73) print "Solexa64"; \
else print "Unknown score encoding"; \
print "( " min ", " max, ")";}'
均为33
使用TrimGalore质控
trim_galore -q 25 -j 5 --phred33 --length 25 -e 0.1 --stringency 4 -o ./ *gz
#-j 是线程，最大8就可以了，软件用不了,但必须在py3环境下，否则会自动切换单线程
3、Bowtie2比对
需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件index：
用bowtie2官网提供的index
！！！这块还不会得问师兄
服务器上有index直接用
比对的时候发现一个问题，index存放再bowtie2文件夹下，含有六个文件。-x index的路径不但精确到文件夹，还要精确到6个mm10文件，前缀必须是mm10不然会报错（如下图）
或者可以提前定义一下路径 index=/storage2/anlei/reference/index/bowtie2/mm10/mm10
bowtie2 -p 50 -x /storage2/anlei/reference/index/bowtie2/mm10/mm10 --very-sensitive -X 2000 -U SRR13811453_trimmed.fq.gz -S /storage2/anlei/wangwenjing/ChIP-seq/bowtie2/SRR13811453.sam
附上4个文件的比对结果。
4、Samtools对文件sam转bam并排序
#先转格式
samtools sort -O bam -@50 SRR13579713.sam -o SRR13579713.raw.bam
#再排序
samtools sort -@ 50 -n SRR13579713.raw.bam -o SRR13579713.sort.bam
5、sort.bam文件去PCR重复
https://wenku.baidu.com/view/97dfc4227a563c1ec5da50e2524de518964bd3dd.html
#fixmate 以名称排序的定位alignment填入配对的座标，ISIZE（inferred insert size猜测的插入序列大小）对相应的标签（flag)
samtools fixmate -@ 50 -m SRR13579713.sort.bam SRR13579713.sort.fixed.bam
#再排序
samtools sort -@ 50 SRR13579713.sort.fixed.bam -o SRR13579713.position.fixed.bam
#
samtools markdup -@ 50 -r SRR13579713.position.fixed.bam SRR13579713.rmdup.bam
#
samtools flagstat SRR13579713.rmdup.bam > SRR13579713.rmdup.alignment.flagstat
#
samtools stats SRR13579713.rmdup.bam > SRR13579713.rmdup.alignment.stat
#multiqc软件质控
multiqc *stat
结果会生成一个文件夹和一个html文件，就是质控的结果
#对排序后的序列索引，用于快速的随机处理，此处将产生.bai文件
samtools index SRR13579713.rmdup.bam
#bedtools转为bed格式
bedtools bamtobed -i SRR13579713.rmdup.bam > SRR13579713.rmdup.bed
5、使用MACS2进行Call peaks
##必须有input，不然会假阳性的
https://www.jianshu.com/p/d62b42ab0d51
对于不同的对象（TF、组蛋白、ATAC参数都不一样）
https://www.likecs.com/show-204987447.html
#narrow、broad、gapped peaks format的区别与联系
macs2 callpeak -f BED -n "H3K27ac_hCG_4h" --keep-dup all -t SRR13579715.rmdup.bed -c SRR13579716.rmdup.bed -g mm --broad-cutoff 0.1
macs2 callpeak -f BED -n "H3K27ac_hCG_4h" --keep-dup all -t SRR13579715.rmdup.bed -c SRR13579716.rmdup.bed -g mm --broad-cutoff 0.1
#################################################################
其实我觉得到这里就结束了，前面call peak之后产出的就是相对于input矫正后的富集peak了。以下什么bam啊bw啊之类的，都是没有input矫正的，也就在IGV看看。
deeptools将bam文件转bw标准化
用去除PCR重复的rmdup.bam，需要有bai文件在 参考：https://www.jianshu.com/p/491557586118
bamCoverage --bam SRR13579713.rmdup.bam -o SRR13579713.bw --normalizeUsing CPM --effectiveGenomeSize 2913022398 --binSize 10 -p50
bamCoverage --bam SRR13579715.rmdup.bam -o /storage2/anlei/wangwenjing/ChIP-seq/deeptools/bw/SRR13579715.bw --normalizeUsing CPM --effectiveGenomeSize 2913022398 --binSize 10 -p50
--normalizeUsing师兄说他经常用CPM
-binSize应该是默认10
--effectiveGenomeSize：可比对基因组区域的大小（bp），可以从这里查http://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html
--scaleFactor：数值，可以是由spike-in推算的校正因子
--numberOfProcessors, -p ：使用的线程数
网友评论