小王的ChIP-seq傻瓜学习教程(纠错更新中)

2022-10-07 本文已影响0人 pudding815

一套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

#################################################################

其实我觉得到这里就结束了，前面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 ：使用的线程数