美文网首页chip-seq生命科学-简书专题科研信息学
Chip-Seq数据处理— 上游篇

Chip-Seq数据处理— 上游篇

作者: 黄晶_id | 来源:发表于2019-03-23 14:48 被阅读278次

Chip-Seq数据是这篇文章里的：Super-Enhancers Promote Transcriptional Dysregulation in Nasopharyngeal Carcinoma
数据下载地址

准备工作——创建环境/安装软件

1.创建 chipseq环境

conda create -n chipseq python=2   #创建 chipseq环境
conda info --envs  #查看当前环境命令，看看是否添加上了 chipseq环境
如图，chipseq环境添加成功

2.启动 chipseq环境

source activate chipseq #启动当前环境
如图，chipseq环境启动成功

3.在chipseq环境下安装软件

还是老规矩这些软件其实可以几个几个一起安装，但是新手上路的我还一个一个来吧，这样保险一点，为什们说这样更保险呢？
答：如果我们放在一起安装 conda install -y sra-tools trim-galore 当前面那个安装出错，后面那个就自动不安装了。

conda install -y sra-tools
conda install -y trim-galore 
conda install -y samtools
conda install -y deeptools  
conda install -y homer
conda install -y meme
conda install -y macs2 
conda install -y bowtie
conda install -y bowtie2

4. 从NCBI上下载SRA数据

过程就是：先下载下来这些SRR号然后导入到Linux里面，最后用一行代码就可以把每一个SRR号所对应的数据下载好。

进入NCBI找到箭头所指，就可以下载该文章的SRR号清单，是一个txt文档
打开后长这样，然后导入到Linux里面
这就导入成功了
之后就要在Linux里用代码把数据库里的数据下载到我们的服务器了
运行下载数据代码之前，如果不做任何前处理就这样直接下载，那会相当相当的慢。那怎么办呢？有一款软件 ——aspera，下载加速器，专治这种从数据库上下载不下来数据的情况。
这个软件不能在conda里下载，下面介绍下载途径
加速下载先装个：asper（这是个下载这个软件的网站） 
#在网站找到最新版本的下载链接
wget -c https://download.asperasoft.com/download/sw/connect/3.8.1/ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
#解压
tar zxvf ibm-aspera-connect-3.7.4.147727-linux-64.tar.gz
#bash之后就安装到自己的环境下了
bash ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.sh

加速下载的软件我们安装好了，下面我们来下载数据：

nohup cat SRR_List.txt|while read id; do (prefetch ${id});done &

默认情况下，这行代码会自动生成 ~/ncbi/public/sra 文档，将下载好的sra数据放进去。

sra数据下载完成

5. sra 数据转成 fastq

第一步一定要将下载好的、默认放到 ncbi 文件夹里的 sra 数据，移动到自己的工作目录，否则是转不成功的。

$ mv ~/ncbi/public/sra ~/project/chipseq

我就直接在自己建立的工作目录~/project/chipseq/sra里运行的代码没有成功！后来我跑去~/ncbi/public/sra目录下运行了一个样本试了一下，可以的，看来果然是文件位置的事情。

在文件所在位置就可以成功转成 faq 压缩文件
把sra数据移动到我们工作目录后，我们开始sra转faq。
正式运行代码之前，必须先拿一个样品测试下代码能否运行成功，这点很关键，因为这步就算成功运行也特别慢，要是代码再出错了就更浪费时间了。 
# 拿第一个样品做测试
ls SRR5315196.sra |fastq-dump -gzip --split-3 -O ./ SRR5315196.sra

sra-tools 里的 fastq-dump工具可以将SRR文件转换为FASTQ格式，--split-3参数表示如果是双端测序就自动拆分，如果是单端不受影响。也就是说，--split-3参数可以将PE的sra文件解压后的fastq文件拆分成_1.fastq和_2.fastq，如果示例数据集是SE测序，不会进行拆分。--gzip转换fastq为压缩文件，节省空间。
单个测试成功，那我们就写循环进行批量转换格式

cat >sra.sh  #写脚本
ls SRR* |while read id; do (fastq-dump -gzip --split-3 -O ./ ${id}) 1>./${id}.sra.log 2>&1;done  #多个循环
nohup bash sra.sh & #挂后台运行

这步一定一定记得挂后台运行（nohup cmd &）因为特别慢，不挂后台你一掉线就功亏一篑了。
我们解读下这个写入脚本里面的循环命令：
ls SRR* 是为了能够把当前文件夹下所有从NCBI上下载的SRR数据列出来

所有的SRR数据
fastq-dump -gzip --split-3 -O ./ ${id} 这条命令是告诉系统，我要转换成fastq格式并且要压缩的；而且表示如果是双端测序就自动拆分，如果是单端不受影响；-O ./ ${id} 表示输出（output）到当前文件夹下，文件名前缀不变（这样的话输出张这样SRR5315196.fastq.gz）
1>./${id}.sra.log 2>&1 （0.标准输入；1.标准输出；2.标准错误）这个命令是说，重定向标准输出到当前文件下SRR5315197.sra.sra.log文件，且标准输出、标准错误到一个文件中（2>&1）
如图所示，SRA转fastq成功

6.制作config文件

制作config文件的作用: 从NCBI上下载SRA数据，之后再转成 fastq.gz格式，在此过程中把原本的文件名（SRR号）改成在跑流程时可以区分各样本的文件名，生信分析中 文件命名很重要，生成fastq文件时，如果不进行更改操作就会直接生成 SRR****.fastq.gz ，我们不能这样，只有这些SRR号我们根本不知道这些样品是些什么，生信分析很重要的一点就是从文件名上我们就要知道这是个什么数据。
首先还去NCBI下载SRA数据的那个界面下载个 txt 文件

image.png
下载到自己的电脑上后，用Xftp导入到服务器上
Xftp
导入成功后，使用命令：
$ head -1 SraRunTable.txt | tr '\t' '\n'|cat -n
image.png
查看该文件的第一行（表头）把它转化成列，并加上行号。
使用命令 $ less -S SraRunTable.txt 查看该表的全部信息，如下图，所有的行都在了，列还有很多，屏幕没有放下，按想有翻页键（End）查看。
$ less -S SraRunTable.txt 命令之后
假装循环我会了，把名字变好之后是这样的
这就是我们所说的 raw data 了
查看SraRunTable.txt文件知道文章是单端测序
查看SraRunTable.txt文件知道文章是单端测序

查看质量报告（QC）

一般过滤（trim_galore）之前都要用fastqc看一下，质量报告决定需不需要下一步的过滤:
ls ../raw/*.gz|xargs fastqc -t 10 -o ./

生成了质量报告
这个过程不是很费时，完成之后会生成.html的文件，使用 Xftp 传到电脑上查看即可。
FastQC Report
其实这篇文章的数据挺好的，不需要过滤，现在我们假装它需要过滤。
使用命令multiqc ./ 把这些质量报告整合到一起，根据整合报告结果选择我们过滤的参数。
整合完成后会生成multiqc_data和multiqc_report.html两个文件。 
drwxrwxr-x 2 jhuang jhuang 4.0K Mar 19 17:23 multiqc_data
-rw-rw-r-- 1 jhuang jhuang 1.3M Mar 19 17:23 multiqc_report.html

我们将multiqc_report.html导入到本地查看，根据整合质量结果，决定过滤 trim_galore 的各种参数。（决定参数需要经验，现在用别人的就好）

本地查看multiqc_report.html

文章里关于chipseq的介绍介绍：
Chromatin immunoprecipitation, sequencing, and analysis
Chromatin immunoprecipitation (ChIP) with antibodies H3K27ac, H3K4me1, H3K4me3, and RNA Pol IIwas performed as described previously. Briefly, NPC cells were crosslinked and sonicated to obtain chromatin DNA fragments between 200 and 250 bp with a Bioruptor Sonicator (Diagenode). Sonicated lysates were cleared and incubated with magnetic beads bound with antibodies. ChIP DNA was eluted, reverse crosslinked, precipitated, and sequenced on Illumina HiSeq 2000.
Single-end 50-bp reads were aligned to hg19 reference genome using Bowtie aligner with parameter –m 2 –k 2 –best –strata –l 49. MACS2 was used for peak calling with parameter –B –SPMR –no-model –ext-size 200 -q 0.01. MACS argument –broad was used for all histone ChIP sequencing (ChIP-seq). Duplicate reads were not considered for peak calling. ChIP-seq heatmaps were drawn using deepTools plotHeatmap function. Pileup signal (bedgraph files) generated by MACS2 were converted to big-wig files using bedGraphToBigWig and visualized in IGV.
SEs were identified using ROSE by stitching peaks within a 12-kb interval. Stitched enhancers were classified as SEs or TEs and assigned to the nearest ensemble genes

从文章里我们知道：1.是单端测序；2.用Bowtie2比对的 3. 参考基因组是 hg19

知道这些我们开始进行下面的流程。

使用trim_galore软件进行质控/过滤

单端过滤的代码如下：

##为了使代码的使用性更加灵活，所以先赋值
analysis_dir=/home/jhuang/project/chipseq
bin_trim_galore="trim_galore"
ls ../raw/*gz | while read fq1;
do 
nohup $bin_trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 5 -o $analysis_dir/clean  $fq1 & 
done

过滤后的文件如下：

trim_galore后生成的文件
过滤完成后理论上再进行QC，以查看过滤前后质量的变化以及我们的过滤效果，但是我们这里数据本来就挺好的，而且上面已经介绍过了如何进行QC，所以这里我们就不做了。

使用bowtie2进行比对

构建索引

bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件。根据文章我们知道人家用的参考基因组是hg19不建议自己下载基因组构建，可以直接下载索引文件，hg19 索引大小为3.5GB，代码如下：（索引比较大，下载了一夜，一定要记得挂后台哦~）

mkdir reference && cd reference
nohup wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/hg19.zip &
unzip hg19.zip

如果物种是小鼠，你就下载小鼠参考基因组的索引和注释文件, 常用的是mm10（看文章人家用的参考基因组是什么,就下载哪个） mm10 索引大小为3.2GB。

wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
索引建立好了以后开始比对

单端测序数据的比对代码如下：

cd /home/jhuang/project/chipseq/align
cat > bidui.sh
#这步要注意，首先一定要试试在该文件夹下能否成功调用 bowtie2 ，如果不行就用软件所安装位置的绝对路径进行赋值。
bin_bowtie2=bowtie2
bowtie2_index=/home/jhuang/project/chipseq/reference/index/bowtie2/hg19_I/hg19
ls ../clean/*gz | while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
$bin_bowtie2  -p 5  -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done
nohup bash bidui.sh &

比对后得到如下的 *.trimmed.bam 文件
(比对这步好曲折、好曲折呀，记得记录一下)

-rw-rw-r-- 1 jhuang jhuang 1.1G Mar 22 08:05 ChIP-Seq_H3K27Ac_1_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 1.1G Mar 22 08:16 ChIP-Seq_H3K27Ac_2_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 713M Mar 22 08:24 ChIP-Seq_H3K27Ac_3_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 720M Mar 22 08:33 ChIP-Seq_H3K4Me1_1_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 731M Mar 22 08:42 ChIP-Seq_H3K4Me1_2_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 721M Mar 22 08:51 ChIP-Seq_H3K4Me1_3_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 567M Mar 22 08:57 ChIP-Seq_H3K4Me3_1_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 666M Mar 22 09:05 ChIP-Seq_H3K4Me3_2_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 571M Mar 22 09:11 ChIP-Seq_H3K4Me3_3_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 2.5G Mar 22 09:45 ChIP-Seq_None_1_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 970M Mar 22 09:56 ChIP-Seq_None_2_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 1.1G Mar 22 10:10 ChIP-Seq_None_3_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 785M Mar 22 10:20 ChIP-Seq_None_4_trimmed.bam
-rw-rw-r-- 1 jhuang jhuang 873M Mar 22 10:30 ChIP-Seq_Pol2_1_trimmed.bam

对bam文件进行QC

QC之前要先建立index，之后进行统计。统计结果就是我们的指标：比对的好坏、比对成功率.....
建立bam文件的index

cd /home/jhuang/project/chipseq/align
nohup ls  *.bam  |xargs -i samtools index {} &

得到索引文件如下：

-rw-rw-r-- 1 jhuang jhuang 3.1M Mar 22 10:46 ChIP-Seq_H3K27Ac_1_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 3.1M Mar 22 10:46 ChIP-Seq_H3K27Ac_2_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 3.0M Mar 22 10:46 ChIP-Seq_H3K27Ac_3_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.9M Mar 22 10:46 ChIP-Seq_H3K4Me1_1_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.9M Mar 22 10:46 ChIP-Seq_H3K4Me1_2_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.8M Mar 22 10:47 ChIP-Seq_H3K4Me1_3_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 4.2M Mar 22 10:47 ChIP-Seq_H3K4Me3_1_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 3.2M Mar 22 10:47 ChIP-Seq_H3K4Me3_2_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 4.2M Mar 22 10:47 ChIP-Seq_H3K4Me3_3_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.7M Mar 22 10:48 ChIP-Seq_None_1_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.2M Mar 22 10:48 ChIP-Seq_None_2_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.1M Mar 22 10:48 ChIP-Seq_None_3_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.2M Mar 22 10:49 ChIP-Seq_None_4_trimmed.bam.bai
-rw-rw-r-- 1 jhuang jhuang 2.2M Mar 22 10:49 ChIP-Seq_Pol2_1_trimmed.bam.bai

得到索引文件后对bam文件进行QC，得到统计结果

cd /home/jhuang/project/chipseq/align
nohup ls  *.bam  | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done &

得到统计结果文件 .stat

-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me3_3_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_None_4_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_None_2_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_Pol2_1_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me3_1_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_None_3_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me3_2_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me1_3_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K27Ac_3_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me1_2_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:55 ChIP-Seq_H3K4Me1_1_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:56 ChIP-Seq_H3K27Ac_1_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  388 Mar 22 10:56 ChIP-Seq_H3K27Ac_2_trimmed.stat
-rw-rw-r-- 1 jhuang jhuang  390 Mar 22 10:56 ChIP-Seq_None_1_trimmed.stat

查看对比成功率

(chipseq) jhuang 14:18:05 ~/project/chipseq/align 
$ grep % *.stat
ChIP-Seq_H3K27Ac_1_trimmed.stat:31875275 + 0 mapped (98.49% : N/A)
ChIP-Seq_H3K27Ac_2_trimmed.stat:29782443 + 0 mapped (98.44% : N/A)
ChIP-Seq_H3K27Ac_3_trimmed.stat:21414528 + 0 mapped (99.08% : N/A)
ChIP-Seq_H3K4Me1_1_trimmed.stat:21775431 + 0 mapped (98.90% : N/A)
ChIP-Seq_H3K4Me1_2_trimmed.stat:22494412 + 0 mapped (98.63% : N/A)
ChIP-Seq_H3K4Me1_3_trimmed.stat:21940837 + 0 mapped (98.52% : N/A)
ChIP-Seq_H3K4Me3_1_trimmed.stat:22588686 + 0 mapped (97.05% : N/A)
ChIP-Seq_H3K4Me3_2_trimmed.stat:21462983 + 0 mapped (97.54% : N/A)
ChIP-Seq_H3K4Me3_3_trimmed.stat:22116810 + 0 mapped (95.96% : N/A)
ChIP-Seq_None_1_trimmed.stat:123234439 + 0 mapped (97.81% : N/A)
ChIP-Seq_None_2_trimmed.stat:24075754 + 0 mapped (98.10% : N/A)
ChIP-Seq_None_3_trimmed.stat:28739221 + 0 mapped (97.34% : N/A)
ChIP-Seq_None_4_trimmed.stat:22603696 + 0 mapped (96.35% : N/A)
ChIP-Seq_Pol2_1_trimmed.stat:23534897 + 0 mapped (96.14% : N/A)

合并bam文件

因为一个样品分成了多个lane进行测序，所以在进行peaks calling的时候，需要把bam进行合并。
代码如下：

mkdir  /home/jhuang/project/chipseq/mergeBam
cd /home/jhuang/project/chipseq/align
nohup ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ../mergeBam/$id.mergeBam $id*.bam ;done &
 ##这个shell大概的意思就是说：用`sed` 命令把

跑完之后得到合并的bam

-rw-rw-r-- 1 jhuang jhuang 2.7G Mar 22 16:34 ChIP-Seq_H3K27Ac.mergeBam
-rw-rw-r-- 1 jhuang jhuang 2.0G Mar 22 16:39 ChIP-Seq_H3K4Me1.mergeBam
-rw-rw-r-- 1 jhuang jhuang 1.7G Mar 22 16:44 ChIP-Seq_H3K4Me3.mergeBam
-rw-rw-r-- 1 jhuang jhuang 5.1G Mar 22 18:22 ChIP-Seq_None.mergeBam
-rw-rw-r-- 1 jhuang jhuang 873M Mar 22 18:23 ChIP-Seq_Pol2.mergeBam
对合并的bam *.mergeBam 进行QC

QC之前要先建立index
建立 *.mergeBam 文件的index

cd /home/jhuang/project/chipseq/mergeBam
nohup ls  *.mergeBam | xargs -i samtools index {} &

得到索引文件(ChIP-Seq_H3K27Ac.mergeBam.bai)后对*.mergeBam 文件进行QC，得到统计结果 *.mergeBam.stat

cd /home/jhuang/project/chipseq/mergeBam
nohup ls *.mergeBam | while read id ;do (samtools flagstat $id > $id.stat);done &

文章里说Duplicate reads were not considered for peak calling.
（接下来我们面临一个很重要的选择叫做：是否去除PCR重复。如果有PCR重复，那在可视化结果中就可以看出一条reads出现很多次，而且起始终点的位置都是一样的。文章我还没有细看，我们这里假如需要去除PCR）

对合并的bam文件*.mergeBam去除PCR重复
nohup ls *.mergeBam | while read id ;do (samtools markdup -r $id $(basename $id ".mergeBam").rmdupBam);done &

跑完后得到去除PCR的mergbam文件 *.rmdupBam

-rw-rw-r-- 1 jhuang jhuang 2.3G Mar 22 22:01 ChIP-Seq_H3K27Ac.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 1.8G Mar 22 22:04 ChIP-Seq_H3K4Me1.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 1.2G Mar 22 22:07 ChIP-Seq_H3K4Me3.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 3.3G Mar 22 22:15 ChIP-Seq_None.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 824M Mar 22 22:16 ChIP-Seq_Pol2.rmdupBam
对 *.rmdupBam进行QC，同样先建立索引

建立 *.rmdupBam 文件的index

nohup ls *.rmdupBam |xargs -i samtools index {} &

得到索引文件(ChIP-Seq_H3K27Ac.rmdupBam.bai)后对*.rmdupBam 文件进行QC，得到统计结果 *.rmdupBam.stat

nohup ls *.rmdupBam | while read id ;do (nohup samtools flagstat $id > $id.stat );done &

合并的bam文件去除PCR重复前后比较：
去除PCR重复前

-rw-rw-r-- 1 jhuang jhuang 2.7G Mar 22 16:34 ChIP-Seq_H3K27Ac.mergeBam
-rw-rw-r-- 1 jhuang jhuang 2.0G Mar 22 16:39 ChIP-Seq_H3K4Me1.mergeBam
-rw-rw-r-- 1 jhuang jhuang 1.7G Mar 22 16:44 ChIP-Seq_H3K4Me3.mergeBam
-rw-rw-r-- 1 jhuang jhuang 5.1G Mar 22 18:22 ChIP-Seq_None.mergeBam
-rw-rw-r-- 1 jhuang jhuang 873M Mar 22 18:23 ChIP-Seq_Pol2.mergeBam

去除PCR重复后

-rw-rw-r-- 1 jhuang jhuang 2.3G Mar 22 22:01 ChIP-Seq_H3K27Ac.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 1.8G Mar 22 22:04 ChIP-Seq_H3K4Me1.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 1.2G Mar 22 22:07 ChIP-Seq_H3K4Me3.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 3.3G Mar 22 22:15 ChIP-Seq_None.rmdupBam
-rw-rw-r-- 1 jhuang jhuang 824M Mar 22 22:16 ChIP-Seq_Pol2.rmdupBam

使用macs2进行找peaks

macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用实例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

一般而言，我们照葫芦画瓢，按照这个实例替换对应部分就行了，介绍一下各个参数的意义

  • -t: 实验组的输出结果

  • -c: 对照组的输出结果

  • -f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，就是自动检测选择。

  • -g: 基因组大小， 默认提供了hs, mm, ce, dm选项， 不在其中的话，比如说拟南芥，就需要自己提供了。

  • -n: 输出文件的前缀名

  • -B: 会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores

  • -q: q值，也就是最小的PDR阈值， 默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。

  • -p： 这个是p值，指定p值后MACS2就不会用q值了。

  • -m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为你很大概率上不会懂。

    找合并的bam *.mergeBam 的peaks:
    ls *mergeBam |cut -d"." -f 1 |while read id; do (nohup macs2 callpeak -c ChIP-Seq_None.mergeBam -t $id.mergeBam -f BAM -B -g hs -n $id 2> $id.log &); done
    找去除PCR重复 *.rmdupBam 的peaks:
    ls *rmdupBam | cut -d"." -f 1 |while read id; do (nohup macs2 callpeak -c ChIP-Seq_None.rmdupBam -t $id.rmdupBam -f BAM -B -g hs -n $id 2> $id.log &); done

    跑完的结果怎么看，最主要的是生成的 .bed 文件，这就是得到的peaks文件，里面有从哪个坐标到哪里坐标是一个peaks，如下。我们统计下bed文件的行数（wc -l）可以得到一共有多少个peaks。该文件也可以去IGV里查看。

    $ less ChIP-Seq_Pol2_summits.bed
chr1    564595  564596  ChIP-Seq_Pol2_peak_1    21.47882
chr1    564940  564941  ChIP-Seq_Pol2_peak_2    9.23926
chr1    565115  565116  ChIP-Seq_Pol2_peak_3    10.15886
chr1    565348  565349  ChIP-Seq_Pol2_peak_4    11.54637

$ wc -l *.bed
  134684 ChIP-Seq_H3K27Ac_summits.bed
  198510 ChIP-Seq_H3K4Me1_summits.bed
   32796 ChIP-Seq_H3K4Me3_summits.bed
       0 ChIP-Seq_None_summits.bed
    1038 ChIP-Seq_Pol2_summits.bed
  367028 total

    得到peaks后，我们chipseq的上游分析结束。下游分析只是一个开始。

相关文章

网友评论

    本文标题：Chip-Seq数据处理— 上游篇

    本文链接：https://www.haomeiwen.com/subject/aytxuqtx.html

    延伸阅读

    深度阅读

    • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

    栏目导航

    热点阅读

    chip-seq

    生命科学-简书专题

    科研信息学

    Linux-NGS-二代测序分析流程

    关于我们|服务条款|联系我们|Chip-Seq数据处理— 上游篇|投稿指南|网站地图|RSS订阅|排版工具|手机版

    提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

    Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

    备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

    本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

    所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！