操作记录-2020-11-13：精简代码处理ChIP_seq数据
2024-04-14 18:32:12
今天准备尝试编写一组精简代码用于处理ChIP_seq数据，希望能成功吧。
1.建立相应目录
对新数据建立对应实验人员（lizexing）、测序类型（ChIP_seq）和日期（2020_11_13）的目录。
# 建立后如下：
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13$# 新建对应的目录
mkdir raw_data clean_data bam bam_bw bam_sort sam macs2_bdgdiff macs2_callpeak matrix_reference_point matrix_scale_regions fastqc_report MD5_txt scripts_log
2.检查数据完整性
cat md5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
操纵记录如下
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/clean_data$ cat *.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
Scr-S-L-KQ_FKDL202604880-1a_1.clean.fq.gz: OK
Scr-S-L-KQ_FKDL202604880-1a_2.clean.fq.gz: OK
shTgm1-2-S-L-KQ_FKDL202604881-1a_1.clean.fq.gz: OK
shTgm1-2-S-L-KQ_FKDL202604881-1a_2.clean.fq.gz: OK
shTgm2-1-S-L-KQ_FKDL202604882-1a_1.clean.fq.gz: OK
shTgm2-1-S-L-KQ_FKDL202604882-1a_2.clean.fq.gz: OK
3.根据需要精简文件名称
操作记录如下：
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/clean_data$ ll
total 2.7G
drwxrwxr-x 2 zexing zexing 4.0K 11月 16 11:45 .
drwxrwxr-x 15 zexing zexing 4.0K 11月 13 10:36 ..
-rw-rw-r-- 1 zexing zexing 476 11月 16 11:31 check_md5sum.txt
-rw-rw-r-- 1 zexing zexing 152 11月 16 10:50 MD5_Scr-S-L-KQ_FKDL202604880-1a.txt
-rw-rw-r-- 1 zexing zexing 162 11月 16 11:03 MD5_shTgm1-2-S-L-KQ_FKDL202604881-1a.txt
-rw-rw-r-- 1 zexing zexing 162 11月 16 11:21 MD5_shTgm2-1-S-L-KQ_FKDL202604882-1a.txt
-rw-rw-r-- 1 zexing zexing 401M 11月 16 10:57 Scr_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 410M 11月 16 11:01 Scr_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 495M 11月 16 11:11 shTgm1-2_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 510M 11月 16 11:20 shTgm1-2_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 423M 11月 16 11:26 shTgm2-1_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 437M 11月 16 11:27 shTgm2-1_2.clean.fq.gz
4. 在Linux服务器中对ChIP_seq数据进行处理并常规call peak。
vim新建ChIP_seq_script_1将数据质控、比对、格式转换、排序、生成目录、bamCoverage命令转换文件格式和macs2 callpeak综合在一起。
#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
# This program is used for ChIP-seq data analysis.
# History
# 2020/11/13 zexing First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13 # 用户名称和日期需要更改# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz# 利用for循环进行后续操作
for i in Scr shTgm1-2 shTgm2-1 # 样品名称需要修改
do
# 对数据进行比对
bowtie2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_mm10/bowtie2_index/mm10 -1 ${dir}/clean_data/${i}_1.clean.fq.gz -2 ${dir}/clean_data/${i}_2.clean.fq.gz -S ${dir}/sam/${i}.sam# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam# 对数据生成目录
samtools index -@ 16 ${dir}/bam_sort/${i}.bam.sort # bamCoverage命令转换文件格式
bamCoverage -p 16 -v -b ${dir}/bam_sort/${i}.bam.sort -o ${dir}/bam_bw/${i}.bam.sort.bw# 使用macs2进行常规callpeak
macs2 callpeak -t ${dir}/bam_sort/${i}.bam.sort -f BAM -g mm -B -q 0.05 --outdir ${dir}/macs2_callpeak/ -n ${i}done
在后台运行ChIP_seq_script_1：
nohup bash ChIP_seq_script_1 > ChIP_seq_script_1_log &
5. 使用deeptools软件绘制热图/密度图
1. computeMatrix scale-regions模式计算信号强度并用plotHeatmap/plotProfile作图
scale-regions模式计算的是区域形式，所以指定作图位置的BED或GTF格式文件为macs2 callpeak生成的后缀为peaks.narrowPeak的BED文件。
vim新建ChIP_seq_script_2，脚本如下：
#! /bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
# This program is used for computeMatrix scale-regions.
#History:
# 2020/11/03 zexing First release
# In the scale-regions mode, all regions in the BED file are stretched or shrunken to the length (in bases) indicated by the user.
# 参数-R 指定作图位置的BED或GTF格式文件，可用#标记同一组区域，默认无。
# 参数-S 输入bigwig文件。
# 参数-o 指定输出为文件名用于plotHeatmap, plotProfile
# 参数-b上游（默认0bp），-a下游（默认0bp）设定感兴趣的区域，如果该区域是基因，则为基因TSS上游或TES下游。
# 参数--skipZeros设定0分区域的处理
# 参数-p 设置线程数bed=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/macs2_callpeak
bw=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/bam_bw
results=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/matrix_scale_regionscomputeMatrix scale-regions \
-R ${bed}/Scr_peaks.narrowPeak ${bed}/shTgm1-2_peaks.narrowPeak ${bed}/shTgm2-1_peaks.narrowPeak \
-S ${bw}/Scr.bam.sort.bw ${bw}/shTgm1-2.bam.sort.bw ${bw}/shTgm2-1.bam.sort.bw \
-o ${results}/matrix_scale_threegroups.gz \
-b 1000 -a 1000 \
-p 16# 使用plotHeatmap对结果绘制热图并聚类
plotHeatmap -m ${results}/matrix_scale_threegroups.gz \
-o ${results}/scale_threegroups_heatmap.png \
--dpi 750 \
--whatToShow 'heatmap and colorbar' \
--startLabel "Start" \
--endLabel "End" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1
# 使用plotProfile对结果绘制密度图并聚类
plotProfile -m ${results}/matrix_scale_threegroups.gz \
-o ${results}/scale_threegroups_profile.png \
--dpi 750 \
--legendLocation upper-right \
--startLabel "Start" \
--endLabel "End" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1 \
--perGroup
后台运行ChIP_seq_script_2脚本如下
nohup bash ChIP_seq_script_2 > ChIP_seq_script_2_log &
2. computeMatrix reference-point模式计算信号强度并用plotHeatmap/plotProfile作图
reference-point模式计算的是峰值高点模式，所以指定作图位置的BED或GTF格式文件为macs2 callpeak生成的后缀为summits.bed的BED文件。
vim新建ChIP_seq_script_3，脚本如下：
#! /bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
# This program is used for computeMatrix reference_point.
#History:
# 2020/11/03 zexing First release
# Reference-point refers to a position within a BED region(e.g., the starting point). In this mode, only those genomicpositions before (upstream) and/or after(downstream) of the reference point will be plotted.
# 参数-R 指定作图位置的BED或GTF格式文件，可用#标记同一组区域，默认无。
# 参数-S 输入bigwig文件。
# 参数-o 指定输出为文件名用于plotHeatmap, plotProfile
# 参数-b上游（默认0bp），-a下游（默认0bp）设定感兴趣的区域，如果该区域是基因，则为基因TSS上游或TES下游。
# 参数--skipZeros设定0分区域的处理
# 参数-p 设置线程数dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13
bed=${dir}/macs2_callpeak
bw=${dir}/bam_bw
results=${dir}/matrix_reference_point
computeMatrix reference-point \
-R ${bed}/Scr_summits.bed ${bed}/shTgm1-2_summits.bed ${bed}/shTgm2-1_summits.bed \
-S ${bw}/Scr.bam.sort.bw ${bw}/shTgm1-2.bam.sort.bw ${bw}/shTgm2-1.bam.sort.bw \
-o ${results}/matrix_reference_threegroups.gz \
-b 1000 -a 1000 \
-p 16# 使用plotHeatmap对结果绘制热图并聚类
plotHeatmap -m ${results}/matrix_reference_threegroups.gz \
-o ${results}/reference_threegroups_heatmap.png \
--dpi 750 \
--whatToShow 'heatmap and colorbar' \
--refPointLabel "Peak_center" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1
# 使用plotProfile对结果绘制密度图并聚类
plotProfile -m ${results}/matrix_reference_threegroups.gz \
-o ${results}/reference_threegroups_profile.png \
--dpi 750 \
--legendLocation upper-right \
--refPointLabel "Peak_center" \
--regionsLabel Scr-peak shTgm1-2-peak shTgm2-1-peak \
--samplesLabel Scr shTgm1-2 shTgm2-1 \
--perGroup
后台运行ChIP_seq_script_3脚本如下
nohup bash ChIP_seq_script_3 > ChIP_seq_script_3_log &
6.在Linux服务器对ChIP_seq数据不同样本间的差异peak进行处理
1. 使用macs2 predictd预测插入片段长度
vim新建ChIP_seq_script_4 预测插入片段长度确定均值。
对于插入片段长度，大多数的转录因子chip_seq数据推荐值为200， 大部分组蛋白修饰的chip_seq数据推荐值为147。具体实验可以根据预测来确定均值。
#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
# This program is used for ChIP-seq data analysis.
# History
# 2020/11/13 zexing First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13 # 用户名称和日期需要更改
# 利用for循环进行后续操作
for i in Scr shTgm1-2 shTgm2-1 # 样品名称需要更改
do
macs2 predictd -i ${dir}/bam_sort/${i}.bam.sort
done
后台运行ChIP_seq_script_4脚本如下
nohup bash ChIP_seq_script_4 > ChIP_seq_script_4_log &
在运行结果日志“ChIP_seq_script_4_log”中查看预测片段长度，操作记录如下：
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/scripts_log$ egrep "predicted fragment length is" ChIP_seq_script_4_log
INFO @ Mon, 16 Nov 2020 17:38:09: # predicted fragment length is 275 bps
INFO @ Mon, 16 Nov 2020 17:38:57: # predicted fragment length is 292 bps
INFO @ Mon, 16 Nov 2020 17:39:36: # predicted fragment length is 284 bps
2. 使用macs2 callpeak查看测序深度
vim新建ChIP_seq_script_5 查看测序深度。
#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
# This program is used for ChIP-seq data analysis.
# History
# 2020/11/13 zexing First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13 # 用户名称和日期需要更改
# 利用for循环进行后续操作
for i in Scr shTgm1-2 shTgm2-1 # 样品名称需要更改
do
macs2 callpeak -t ${dir}/bam_sort/${i}.bam.sort --outdir ${dir}/macs2_bdgdiff/ -n ${i} -q 0.05 -g mm -B --nomodel --extsize 260
egrep "tags after filtering in treatment|tags after filtering in control" ${dir}/macs2_bdgdiff/${i}_peaks.xls >> ${dir}/macs2_bdgdiff/${i}_tags.txt
done
后台运行ChIP_seq_script_5脚本如下
nohup bash ChIP_seq_script_5 > ChIP_seq_script_5_log &
查看各样本的测序深度，操作记录如下：
(base) zexing@DNA:~/projects/lizexing/ChIP_seq/2020_11_13/macs2_bdgdiff$ cat Scr_tags.txt
# tags after filtering in treatment: 4236187
3. 使用macs2 bdgdiff提取样品间差异peak
vim新建ChIP_seq_script_6 提取样品间差异peak。
#! /bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
# This program is used for calling differential binding events of ChIP-seq data by macs2.
# History:
# 2020/11/06 zexing First release
# 用法说明：
# usage: macs2 bdgdiff [-h] --t1 T1BDG --t2 T2BDG --c1 C1BDG --c2 C2BDG
# [-C CUTOFF] [-l MINLEN] [-g MAXGAP] [--d1 DEPTH1]
# [--d2 DEPTH2] [--outdir OUTDIR]
# (--o-prefix OPREFIX | -o OFILE OFILE OFILE)
# 可选参数：
# optional arguments:
# 参数 --t1是读取MACS pileup bedGraph for condition 1.
# 参数 --t2是读取MACS pileup bedGraph for condition 2.
# 参数 --c1是读取MACS control lambda bedGraph for condition 1.
# 参数 --c2是读取MACS control lambda bedGraph for condition 2.
# 参数 -g 是Maximu gap to merge nearby differential regions.
# 参数 -l Minimum length of differential region. Try bigger value to remove small regions. DEFAULT: 200
# 参数 --d1 Sequencing depth (# of non-redundant reads in million) for condition 1.
# 参数 --d2 Sequencing depth (# of non-redundant reads in million) for condition 2.
# 参数 --o-prefix diff_c1_vs_c2保存输出文件名。# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/lizexing/ChIP_seq/2020_11_13/macs2_bdgdiffmacs2 bdgdiff --t1 ${dir}/Scr_treat_pileup.bdg --c1 ${dir}/Scr_control_lambda.bdg --d1 4236187 \
--t2 ${dir}/shTgm1-2_treat_pileup.bdg --c2 ${dir}/shTgm1-2_control_lambda.bdg --d2 5103633 \
-g 60 -l 260 --o-prefix ${dir}/diff_Scr_vs_shTgm1-2macs2 bdgdiff --t1 ${dir}/Scr_treat_pileup.bdg --c1 ${dir}/Scr_control_lambda.bdg --d1 4236187 \
--t2 ${dir}/shTgm2-1_treat_pileup.bdg --c2 ${dir}/shTgm2-1_control_lambda.bdg --d2 4491626 \
-g 60 -l 260 --o-prefix ${dir}/diff_Scr_vs_shTgm2-1后台运行ChIP_seq_script_6脚本如下
nohup bash ChIP_seq_script_6 > ChIP_seq_script_6_log &
其中-d1和-d2的值就是第二步运行时输出的reads数目，-o参数指定输出文件的前缀。运行成功后，会产生3个文件
diff_Scr_vs_shTgm1-2_c3.0_cond1.bed # 保存在condition1中上调的peakdiff_Scr_vs_shTgm1-2_c3.0_cond2.bed # 保存了在condition2中上调的peakdiff_Scr_vs_shTgm1-2_c3.0_common.bed # 保存的是没有达到阈值的，非显著差异peak
上述3个文件格式是完全相同的，最后一列的内容为log10 likehood ratio值，用来衡量两个条件之间的差异，默认阈值为3，大于阈值的peak为组间差异显著的peak, 这个阈值可以通过-c参数进行调整。
7.在RStudio中利用ChIPseeker进行基因注释及转换
此部分使用上一步生成的bed文件在本地机Rstudio中，利用ChIPseeker进行操作。
1. 对常规macs2 callpeak的BED文件进行注释
narrowPeak文件是针对一段peak区域注释，summits.bed是针对peak峰值注释。
# 编辑脚本如下：
# This script is used for analysis of daizhongye ChIP-seq data
# History
# Lizexing 2020-11-13 First release
# 原始测序数据经过在服务器上进行bowtie2比对和macs2 callpeak分析后，得到的bed文件，
# 将其下载之本地电脑后进行后续操作
# 安装ChIPseeker包
# BiocManager::install("ChIPseeker")
# 设置工作目录
setwd("G:/daizhongye/ChIP-seq/2020_10_29/macs2_callpeak/")
#加载ChIPseeker包
library(ChIPseeker)
# 加载基因组注释库
# 安装小鼠注释包
# BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
# 安装人的注释包
# BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
# 读取chipseq峰的bed文件
Scr_peak <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_callpeak/Scr_summits.bed")
shTgm1_2_peak <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_callpeak/shTgm1-2_summits.bed")
shTgm2_1_peak <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_callpeak/shTgm1-2_summits.bed")# 注释，TSS的范围可自定义
# 加载小鼠基因组注释包
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
# 对txdb进行指定
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene# 进行注释
Scr_peakAnno <- annotatePeak(Scr_peak, tssRegion = c(-3000, 3000), TxDb = txdb)
shTgm1_2_peakAnno <- annotatePeak(shTgm1_2_peak, tssRegion = c(-3000, 3000), TxDb = txdb)
shTgm2_1_peakAnno <- annotatePeak(shTgm2_1_peak, tssRegion = c(-3000, 3000), TxDb = txdb)# 输出结果
# 设置工作目录
setwd("G:/daizhongye/ChIP-seq/2020_10_29/peakanno/")
write.table(Scr_peakAnno, file = "Scr_peak.txt",sep = '\t', quote = FALSE, row.names = FALSE)
write.table(shTgm1_2_peakAnno, file = "shTgm1_2_peak.txt",sep = '\t', quote = FALSE, row.names = FALSE)
write.table(shTgm2_1_peakAnno, file = "shTgm2_1_peak.txt",sep = '\t', quote = FALSE, row.names = FALSE)# 对Scr数据分布进行绘图
tiff("Scr_peakAnno_1.tiff")
plotAnnoBar(Scr_peakAnno)
dev.off()tiff("Scr_peakAnno_2.tiff")
vennpie(Scr_peakAnno)
dev.off()tiff("Scr_peakAnno_3.tiff")
plotAnnoPie(Scr_peakAnno)
dev.off()tiff("Scr_peakAnno_4.tiff")
plotDistToTSS(Scr_peakAnno)
dev.off()#对一组数据分布进行绘图
tiff("shTgm1_2_peakAnno_1.tiff")
plotAnnoBar(shTgm1_2_peakAnno)
dev.off()tiff("shTgm1_2_peakAnno_2.tiff")
vennpie(shTgm1_2_peakAnno)
dev.off()tiff("shTgm1_2_peakAnno_3.tiff")
plotAnnoPie(shTgm1_2_peakAnno)
dev.off()tiff("shTgm1_2_peakAnno_4.tiff")
plotDistToTSS(shTgm1_2_peakAnno)
dev.off()# 对二组数据分布进行绘图
tiff("shTgm2_1_peakAnno_1.tiff")
plotAnnoBar(shTgm2_1_peakAnno)
dev.off()tiff("shTgm2_1_peakAnno_2.tiff")
vennpie(shTgm2_1_peakAnno)
dev.off()tiff("shTgm2_1_peakAnno_3.tiff")
plotAnnoPie(shTgm2_1_peakAnno)
dev.off()tiff("shTgm2_1_peakAnno_4.tiff")
plotDistToTSS(shTgm2_1_peakAnno)
dev.off()
2. 利用下述脚本对提取出来的cond1.bed和cond2.bed进行注释
参考文章：CHIP-seq流程学习笔记（6）-peak注释软件ChIPseeker
# 编辑脚本如下：
# This script is used for Annotate peaks from macs2_bdgdiff of daizhongye ChIP-seq data。
# History
# Lizexing 2020-11-07 First release
# 利用macs2 bdgdiff命令得到差异化的peak后，得到3种类型的bed文件，将其下载之本地电脑后进行操作。
# 安装ChIPseeker包
# BiocManager::install("ChIPseeker")
# 设置工作目录
setwd("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/")
#加载ChIPseeker包
library(ChIPseeker)
# 加载基因组注释库
# 安装小鼠注释包
# BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
# 安装人的注释包
# BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
# 读取差异化peak的bed文件
Scr_vs_shTgm1_2_cond1 <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/diff_Scr_vs_shTgm1-2_c3.0_cond1.bed")
Scr_vs_shTgm1_2_cond2 <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/diff_Scr_vs_shTgm1-2_c3.0_cond2.bed")
Scr_vs_shTgm2_1_cond1 <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/diff_Scr_vs_shTgm2-1_c3.0_cond1.bed")
Scr_vs_shTgm2_1_cond2 <- readPeakFile("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/diff_Scr_vs_shTgm2-1_c3.0_cond2.bed")# 注释，TSS的范围可自定义
# 加载小鼠基因组注释包
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
# 对txdb进行指定
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene# 进行注释
Scr_vs_shTgm1_2_cond1_peakAnno <- annotatePeak(Scr_vs_shTgm1_2_cond1, tssRegion = c(-3000, 3000), TxDb = txdb)
Scr_vs_shTgm1_2_cond2_peakAnno <- annotatePeak(Scr_vs_shTgm1_2_cond2, tssRegion = c(-3000, 3000), TxDb = txdb)
Scr_vs_shTgm2_1_cond1_peakAnno <- annotatePeak(Scr_vs_shTgm2_1_cond1, tssRegion = c(-3000, 3000), TxDb = txdb)
Scr_vs_shTgm2_1_cond2_peakAnno <- annotatePeak(Scr_vs_shTgm2_1_cond2, tssRegion = c(-3000, 3000), TxDb = txdb)# 保存注释结果
write.table(Scr_vs_shTgm1_2_cond1_peakAnno, file = "Scr_vs_shTgm1_2_cond1_peakAnno.csv",sep = '\t', quote = TRUE, row.names = FALSE)
write.table(Scr_vs_shTgm1_2_cond2_peakAnno, file = "Scr_vs_shTgm1_2_cond2_peakAnno.csv",sep = '\t', quote = TRUE, row.names = FALSE)
write.table(Scr_vs_shTgm2_1_cond1_peakAnno, file = "Scr_vs_shTgm2_1_cond1_peakAnno.csv",sep = '\t', quote = TRUE, row.names = FALSE)
write.table(Scr_vs_shTgm2_1_cond2_peakAnno, file = "Scr_vs_shTgm2_1_cond2_peakAnno.csv",sep = '\t', quote = TRUE, row.names = FALSE)
3. 利用下述脚本将注释文件中的gene_id转化为gene_symbol
参考文章：R中常用函数使用说明
# 构建待处理基因集的向量
setwd("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/")Scr_vs_shTgm1_2_cond1 <- read.csv("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond1_peakAnno.csv", header = TRUE,sep = "\t" )
Scr_vs_shTgm1_2_cond2 <- read.csv("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond2_peakAnno.csv", header = TRUE,sep = "\t" )
Scr_vs_shTgm2_1_cond1 <- read.csv("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond1_peakAnno.csv", header = TRUE,sep = "\t" )
Scr_vs_shTgm2_1_cond2 <- read.csv("G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond2_peakAnno.csv", header = TRUE,sep = "\t" )# 提取待处理基因集中的gene_id一列并转化为向量格式
Scr_vs_shTgm1_2_cond1_V1 <- as.vector(Scr_vs_shTgm1_2_cond1[, 14])
Scr_vs_shTgm1_2_cond2_V1 <- as.vector(Scr_vs_shTgm1_2_cond2[, 14])
Scr_vs_shTgm2_1_cond1_V1 <- as.vector(Scr_vs_shTgm2_1_cond1[, 14])
Scr_vs_shTgm2_1_cond2_V1 <- as.vector(Scr_vs_shTgm2_1_cond2[, 14])# 由鼠的gene_id转化到gene_symbol
library("clusterProfiler")
library("org.Mm.eg.db")
Scr_vs_shTgm1_2_cond1_V2 <- bitr(Scr_vs_shTgm1_2_cond1_V1, # 待转化的文件名fromType = "ENTREZID", # fromType是指你的数据ID类型是属于哪一类的toType = "SYMBOL", # toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种OrgDb = org.Mm.eg.db) # Orgdb是指对应的注释包是哪个Scr_vs_shTgm1_2_cond2_V2 <- bitr(Scr_vs_shTgm1_2_cond2_V1, # 待转化的文件名fromType = "ENTREZID", # fromType是指你的数据ID类型是属于哪一类的toType = "SYMBOL", # toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种OrgDb = org.Mm.eg.db) # Orgdb是指对应的注释包是哪个Scr_vs_shTgm2_1_cond1_V2 <- bitr(Scr_vs_shTgm2_1_cond1_V1, # 待转化的文件名fromType = "ENTREZID", # fromType是指你的数据ID类型是属于哪一类的toType = "SYMBOL", # toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种OrgDb = org.Mm.eg.db) # Orgdb是指对应的注释包是哪个Scr_vs_shTgm2_1_cond2_V2 <- bitr(Scr_vs_shTgm2_1_cond2_V1, # 待转化的文件名fromType = "ENTREZID", # fromType是指你的数据ID类型是属于哪一类的toType = "SYMBOL", # toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种OrgDb = org.Mm.eg.db) # Orgdb是指对应的注释包是哪个# 查看转化后的结果
View(Scr_vs_shTgm1_2_cond1_V2)
View(Scr_vs_shTgm1_2_cond2_V2)
View(Scr_vs_shTgm2_1_cond1_V2)
View(Scr_vs_shTgm2_1_cond2_V2)# 保存差异gene_symbol以便后续处理
write.csv(Scr_vs_shTgm1_2_cond1_V2, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond1_V2.txt", row.names = FALSE)
write.csv(Scr_vs_shTgm1_2_cond2_V2, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond2_V2.txt", row.names = FALSE)
write.csv(Scr_vs_shTgm2_1_cond1_V2, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond1_V2.txt", row.names = FALSE)
write.csv(Scr_vs_shTgm2_1_cond2_V2, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond2_V2.txt", row.names = FALSE)# 更改列名称进行合并文件
colnames(Scr_vs_shTgm1_2_cond1_V2)[1] <- "geneId"
colnames(Scr_vs_shTgm1_2_cond2_V2)[1] <- "geneId"
colnames(Scr_vs_shTgm2_1_cond1_V2)[1] <- "geneId"
colnames(Scr_vs_shTgm2_1_cond2_V2)[1] <- "geneId"# 合并转化后的文件
Scr_vs_shTgm1_2_cond1_V3 <- merge(Scr_vs_shTgm1_2_cond1_V2, Scr_vs_shTgm1_2_cond1, by = "geneId")
Scr_vs_shTgm1_2_cond2_V3 <- merge(Scr_vs_shTgm1_2_cond2_V2, Scr_vs_shTgm1_2_cond2, by = "geneId")
Scr_vs_shTgm2_1_cond1_V3 <- merge(Scr_vs_shTgm2_1_cond1_V2, Scr_vs_shTgm2_1_cond1, by = "geneId")
Scr_vs_shTgm2_1_cond2_V3 <- merge(Scr_vs_shTgm2_1_cond2_V2, Scr_vs_shTgm2_1_cond2, by = "geneId")# 对结果进行输出保存
write.csv(Scr_vs_shTgm1_2_cond1_V3, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond1_V3.csv", row.names = FALSE, quote = TRUE)
write.csv(Scr_vs_shTgm1_2_cond2_V3, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm1_2_cond2_V3.csv", row.names = FALSE, quote = TRUE)
write.csv(Scr_vs_shTgm2_1_cond1_V3, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond1_V3.csv", row.names = FALSE, quote = TRUE)
write.csv(Scr_vs_shTgm2_1_cond2_V3, "G:/daizhongye/ChIP-seq/2020_10_29/macs2_bdgdiff/Scr_vs_shTgm2_1_cond2_V3.csv", row.names = FALSE, quote = TRUE)
4. 利用intersect()函数将ChIP-seq和RNA-seq中相同的基因提取出来
参考文章：R中常用函数使用说明
# 利用intersect函数对RNA_seq和ChIP_seq中的交集gene_symbol进行提取。
# 参考文章：https://blog.csdn.net/woodcorpse/article/details/80494605?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160471856319195264746932%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=160471856319195264746932&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~baidu_landing_v2~default-1-80494605.pc_first_rank_v2_rank_v28&utm_term=R%E4%B8%ADintersect&spm=1018.2118.3001.4449
# 交集intersect：两个向量的交集，集合可以是数字、字符串等# ChIP_seq中的差异peak在上面已经定义，分别是：
# Scr_vs_shTgm1_2_cond1即Scr中上调的peak，即敲除shTgm1_2后降低的peak
ChIP_seq_group_1_down_genes <- as.vector(Scr_vs_shTgm1_2_cond1_V2$SYMBOL)
# Scr_vs_shTgm1_2_cond2即敲除shTgm1_2后升高的peak
ChIP_seq_group_1_up_genes <- as.vector(Scr_vs_shTgm1_2_cond2_V2$SYMBOL)
# Scr_vs_shTgm2_1_cond1即Scr中上调的peak，即敲除shTgm2_1后降低的peak
ChIP_seq_group_2_down_genes <-as.vector(Scr_vs_shTgm2_1_cond1_V2$SYMBOL)
# Scr_vs_shTgm2_1_cond2即敲除shTgm2_1后升高的peak
ChIP_seq_group_2_up_genes <- as.vector(Scr_vs_shTgm2_1_cond2_V2$SYMBOL)# RNA_seq中的差异gene需要再重新读入，分别是：
# 敲除shTgm1_2后降低的genes
RNA_seq_group_1_down_genes <- read.csv("G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/significant_different_genes/significant_pvalue_different_genes_group_1_genecount_down.csv", header = TRUE,sep = ",")
# 敲除shTgm1_2后升高的genes
RNA_seq_group_1_up_genes <- read.csv("G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/significant_different_genes/significant_pvalue_different_genes_group_1_genecount_up.csv", header = TRUE,sep = ",")
# 敲除shTgm2_1后降低的genes
RNA_seq_group_2_down_genes <- read.csv("G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/significant_different_genes/significant_pvalue_different_genes_group_2_genecount_down.csv", header = TRUE,sep = ",")
# 敲除shTgm2_1后升高的genes
RNA_seq_group_2_up_genes <- read.csv("G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/significant_different_genes/significant_pvalue_different_genes_group_2_genecount_up.csv", header = TRUE,sep = ",")# 将RNA_seq中的差异genes对应的一列提取出来，并转化为向量形式
RNA_seq_group_1_down_genes_V1 <- as.vector(RNA_seq_group_1_down_genes$X)
RNA_seq_group_1_up_genes_V1 <- as.vector(RNA_seq_group_1_up_genes$X)
RNA_seq_group_2_down_genes_V1 <- as.vector(RNA_seq_group_2_down_genes$X)
RNA_seq_group_2_up_genes_V1 <- as.vector(RNA_seq_group_2_up_genes$X)# 提取ChIP_seq和RNA_seq中共有的gene交集
group_1_down_genes <- intersect(ChIP_seq_group_1_down_genes,RNA_seq_group_1_down_genes_V1)
group_1_up_genes <- intersect(ChIP_seq_group_1_up_genes,RNA_seq_group_1_up_genes_V1)
group_2_down_genes <- intersect(ChIP_seq_group_2_down_genes,RNA_seq_group_2_down_genes_V1)
group_2_up_genes <- intersect(ChIP_seq_group_2_up_genes,RNA_seq_group_2_up_genes_V1)
操作记录-2020-11-13：精简代码处理ChIP_seq数据相关推荐
- 操作记录-2020-11-08：精简代码处理RNA_seq数据
今天准备尝试编写一组精简代码用于处理RNA_seq数据,希望能成功吧. 1.建立相应目录 对新数据建立对应实验人员(lizexing).测序类型(RNA_seq)和日期(2020_10_20)的目录. ...
- 工作篇-佛山三水恒大-2020.11.13
** 工作篇-佛山三水恒大-2020.11.14 **TAG:此篇文章估计会很长,因为工作的时候变数太多了,预计五千字左右,想看的可以耐心看完,均为个人实战经验.===害,其实是上学期间请假去做的,还 ...
- 人工蜂群算法求解TSP旅行商问题C++（2020.11.13）
ABC算法求解TSP问题的C++实现 1.输入数据文件:bayg29.tsp 2.头文件 3.所需的类 3.1 城市类City 3.2 包含城市的地图类Graph 3.3 蜜蜂类Bee 3.4 人工蜂 ...
- Leetcode每日一题2020.11.13第328题：奇偶链表
328.奇偶链表 题目描述 思路.算法及代码实现 方法:分离节点后合并 如果链表为空,则直接返回链表. 对于原始链表,每个节点都是奇数节点或偶数节点.头节点是奇数节点,头节点的后一个节点是偶数节点,相 ...
- GopherCon 2020技术演讲slide先睹为快 | Gopher Daily (2020.11.13) ʕ◔ϖ◔ʔ
每日一谚:Simplicity is the art of hiding complexity. 1.GopherCon 2020技术演讲slide先睹为快 - 链接: https://pan.bai ...
- 2020/11/13·Java·人脸识别一键登录/注册
Java·人脸识别一键登录/注册 1.eclipse 和 tomcat 服务器的安装与使用 1.1 Eclipse 的安装 1.2 Eclipse 配置 Tomcat 1.3 新建 Web 项目 1. ...
- 【2020.11.13 八上】期中总结~~
目录 浅谈 DAY 1 DAY 2 感想 计划 成绩 浅谈 --愿你以渺小启程,以伟大结束 闲话 这次的期中考准备的非常匆忙,加上CSP J/S,和一些杂七杂八的机构测试,感觉时间上安排的不是很好. ...
- QIIME 2教程. 09数据导入Importing data(2020.11)
文章目录 QIIME 2用户文档. 9数据导入 导入带质量值的FASTQ测序数据 EMP标准混样单端数据 EMP混样双端数据 Casava1.8单端混样数据 Casava 1.8双端拆分后数据 **F ...
- 【不忘初心】Win10_20H2_2009_19042.610_X64_七合一_[纯净精简版][2.83G](2020.11.2)
母版来自MSDN WIN10_20H2.19042.508,集成补到19042.610,20H2相比1909 2004版本要稳定很多,此版修复了上次的一些问题,应微软毒.粉的要求新增一版带Defend ...
- SVN管理代码查看文件（代码）操作记录
两种方法: 一.从开发工具Eclipse查看 1.右键项目--Team--显示资源历史记录 2.此时可在History视图看到代码操作记录 3.选中某一条操作记录右键--比较 4.选择互相比较的两个版 ...
最新文章
- python打完代码怎么运行-Python的代码是如何去进行运行的？
- 耳鼻喉专科服务机构“仁树医疗”完成数千万元A轮融资...
- 特朗普的《AI 倡议》存在一个致命问题
- 试图使用removebg工具的在线网站去除图片背景时遇到的错误
- 一次前后端分离的实践
- 如何使用Instruments诊断App（Swift版）：起步
- python3 urllib3文档_python urllib3
- 光伏巨头“脱轨” 英利确认债务重组
- thinkphp 如何调用百度echarts 数据报表插件
- android开发01 --开发工具
- Zotero使用记录----1 下载与安装
- 谷歌AI平均每天发表2篇论文！Jeff Dean执笔年度汇总：16大方向
- macbook打开网页慢解决办法
- 海致大数据京信_海致网聚提出公安大数据背景下的个人计算新理念
- B. Shifting Sort（rotate函数旋转应用）
- 程序员的520，送给女友的几行漂亮的代码（js版）
- 【FFmpeg杂记】音频解码输出PCM格式数据分析
- 从志愿军“断刀”再论敏捷之道（上篇）
- Android Instant App 介绍
- 一中模拟赛3.15——树上gcd