操作记录-2020-11-13：精简代码处理ChIP_seq数据
生活随笔
收集整理的这篇文章主要介绍了
操作记录-2020-11-13：精简代码处理ChIP_seq数据
小编觉得挺不错的,现在分享给大家,帮大家做个参考.
今天准备尝试编写一组精简代码用于处理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_log2.检查数据完整性
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: OK3.根据需要精简文件名称
操作记录如下：
(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.gz4. 在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 bps2. 使用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: 42361873. 使用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数据的全部内容，希望文章能够帮你解决所遇到的问题。
- 上一篇： 格兰杰检验的基本步骤_Toda-Yama
- 下一篇： 自然环境资源数据集分享——资源环境数据云