数据分析-cuttag分析流程分享3-个性化分析内容
原创数据分析-cuttag分析流程分享3-个性化分析内容
原创
在进行了前面两次的流程分析，目前已经得到了bedgarph文件和peak文件，需要在后面对peak文件进行相关的分析，主要有差异peak分析、peak的注释、注释基因的富集分析以及motif分析，我做了几次，发现里面的坑还是很多的。
对bedgraph文件进行可视化以及热图绘制
前面拿到bedgraph文件，就可以在IGV上对峰的结果进行查看了，我还是引的官网上的图片，我感觉我做出来的峰没有官网上这么好看。
通过对全基因组call到的峰来进行特定区域的热图绘制，可以确定自己所研究的蛋白是处于哪个区域，我也没搞明白，为什么我做出来的东西都在gene body上，但是看了igv确实是没差，好奇一枚。
这个画图，只能在linux服务器上，耗的时间还是挺长的，8个线程拟南芥的需要8-9个小时，所以我一般看没人用服务器，直接来20个线程，将工作时间缩短一下。
先对bam文件转换成bw文件。
// ##== linux command ==##
for i in D_rep1 D_rep2 D_rep3 D_rep4 IgG_rep1; do
{
samtools sort -o ${i}.sorted.bam ${i}_bowtie2.mapped.bam
samtools index ${i}.sorted.bam
bamCoverage -b ${i}.sorted.bam -o ${i}_raw.bw
}&
done在进行热图绘制之前，需要对注释文件进行提取来获得基因的注释内容
// ##== linux command ==##
cat TAIR_exons.gtf | cut -f1,4,5,9 | cut -f1 -d";" | awk '{print $1, $2, $3, $5}' | sed -e 's/ /\t/g' | sed -e 's/\"//g' | sed -e 's/transcript\://g' > TAIR_exons.genomic.bed
cores=10
computeMatrix scale-regions -S D_rep1_raw.bw D_rep2_raw.bw D_rep3_raw.bw \
-R /CUTTAG/TAIR_exons.genomic.bed \
--beforeRegionStartLength 3000 --regionBodyLength 5000 \
--afterRegionStartLength 3000 --skipZeros -o LDLDmatrix_gene.mat.gz -p $cores
plotHeatmap -m matrix_gene.mat.gz -out Histone_gene.png --sortUsing sum --missingDataColor '#cc2200'差不多经过一晚，应该就能出来了，其实这个也不是绝对的，需要提前对参考基因组及样本的大小进行估计，像另一个师姐的样本比较多，差不多就在服务器上跑了4d。
加入--missingDataColor '#cc2200'参数主要是因为有的时候bed文件的比对不到注释文件里面，出来的图就会有很多黑线，影响美观，所以可以加入这个参数，给图片美化。
官网还推荐了基于峰的长度的中位来进行热图绘制，我感觉这个比较漂亮，但是大家的文章一般也没有放这个图。
//
for i in D_rep1 D_rep2 D_rep3; do
{
cores = 10
awk '{split($6, summit, ":"); split(summit[2], region, "-"); print summit[1]"\t"region[1]"\t"region[2]}' ${i}_seacr_control.peaks.stringent.bed >${i}_seacr_control.peaks.summitRegion.bed
computeMatrix reference-point -S ${i}_raw.bw \
-R ${i}_seacr_control.peaks.summitRegion.bed \
--skipZeros -o ${i}_SEACR.mat.gz -p $cores -a 3000 -b 3000 --referencePoint center
plotHeatmap -m 4 ${i}_SEACR.mat.gz -out ${i}_SEACR_heatmap.png --sortUsing sum --startLabel "Peak Start" -\
-endLabel "Peak End" --xAxisLabel "" --regionsLabel "Peaks" --samplesLabel "${i}"
}&
done也是写了一个循环，在后台对几个文件进行热图的绘制，一般这个出的一个一个的单图，后面可以做图片的拼接。
差异peak分析
首先是官网推荐的差异peak分析，官网上主要是根据DEseq2来进行差异分析的，但是一直不太理解的就是我改了循环之后，还是表头读不进去我的样本名称。
可以发现在做差异peak的时候主要有两种分析思路，一种先做差异peak的基因注释，来对差异基因进行分析；第二个是先进行差异peak，然后做差异peak的基因注释，拿到基因后在进行下面的分析。因为师姐的实验目的主要是比较不同的时间处理之后是不是基因发生了什么丰度变化，所以我是先做的差异peak的分析。
首先我们看一下官网上推荐的代码，是可以运行的，但是不适合我要用差异peak的位置来进行后面的基因注释，所以放弃了这个分析流程。
// CUTTAG官网推荐差异peak流程
##R语言##
##差异分析，前面已经读了histL及repL，大家可以直接在前面在加上
mPeak = GRanges()
## overlap with bam file to get count
for(hist in histL){
for(rep in repL){
peakRes = read.table(paste0(projPath, "/peakCalling/SEACR/", hist, "_", rep, "_seacr_control.peaks.stringent.bed"), header = FALSE, fill = TRUE)
mPeak = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = "*") %>% append(mPeak, .)
}
}
masterPeak = reduce(mPeak)
### 这个masterPeak存了所有的峰的位置
##获取主峰列表中每个峰的碎片计数
##=== R command ===##
library(DESeq2)
bamDir = paste0(projPath, "/2.cleandata")
countMat = matrix(NA, length(masterPeak), length(histL)*length(repL))
##官网代码有问题，加循环
for(hist in histL){
for(rep in repL){
countMat = matrix(NA,length(masterPeak), length(histL)*length(repL))
}
}
##overlap with bam file to get count
i = 1
for(hist in histL){
for(rep in repL){
bamFile = paste0(bamDir, "/", hist, "_", rep, "_bowtie2.mapped.bam")
fragment_counts <- getCounts(bamFile, masterPeak, paired = TRUE, by_rg = FALSE, format = "bam")
countMat[, i] = counts(fragment_counts)[,1]
i = i + 1
}
}
##测序深度归一化和差异富集峰检测
##=== R command ===##
selectR = which(rowSums(countMat) > 5) ## remove low count genes
dataS = countMat[selectR,]
condition = factor(rep(histL, each = length(repL)))
dds = DESeqDataSetFromMatrix(countData = dataS,
colData = DataFrame(condition),
design = ~ 1)
DDS = DESeq(dds)
normDDS = counts(DDS, normalized = TRUE) ## normalization with respect to the sequencing depth
colnames(normDDS) = paste0(colnames(normDDS), "_norm")
res = results(DDS, independentFiltering = FALSE, altHypothesis = "greaterAbs")
countMatDiff = cbind(dataS, normDDS, res)
head(countMatDiff)
write.table(countMatDiff,"difference.txt",sep="\t",row.names=F)
我的数据运行了这个代码，也有结果，但是在前面说到的表头不太对，改了好几次，还是这样，所以一直觉得这个不靠谱，我就选了下面这个DiffBind进行分析，结果还是可以的。
DiffBind差异peak分析
其实这个包只要安装对了，相对来说后面的就很好走了。如果一直安装失败，就得按照报错的一直补没有的包，要耐心呐，友友们，我是花了一天才装完。
我在分析的时候是直接用的没有排序的bam文件，导致一直在报错，说建不了index，也是折腾了一天，后来在服务器上搞了搞，才明白。
// linux代码
for i in D_rep1 D_rep2 D_rep3;do
{
samtools sort -@ 8 ${i}_bowtie2.mapped.bam -o ${i}_bowtie2.mapped.sorted.bam
samtools index ${i}_bowtie2.mapped.sorted.bam
}&
done也是几个循环后，基本就排好序了。
下面是用DiffBind进行差异peak，需要提前建一个sample.csv文件，把自己要进行差异分析的文件提前读进去，大致可以按照下面这个方式，我主要是参照这个博主的分析流程（https://www.jianshu.com/p/f849bd55ac27），前面给bam文件排好序后，真香，全程无错。
// ##读入文件
library(DiffBind)
dbObj <- dba(sampleSheet="SampleSheet.csv")
##亲和结合矩阵
dbObj <- dba.count(dbObj, bUseSummarizeOverlaps=TRUE)
dba.plotPCA(dbObj, attributes=DBA_FACTOR, label=DBA_ID)
plot(dbObj)
##差异分析
# Establishing a contrast
dbObj <- dba.contrast(dbObj, categories=DBA_FACTOR,minMembers = 2)
dbObj <- dba.analyze(dbObj, method=DBA_ALL_METHODS)
# summary of results
dba.show(dbObj, bContrasts=T)
# overlapping peaks identified by the two different tools (DESeq2 and edgeR)
dba.plotVenn(dbObj,contrast=1,method=DBA_ALL_METHODS)
##提取两个差异分析方法的结果
comp1.edgeR <- dba.report(dbObj, method=DBA_EDGER, contrast = 1, th=1)
comp1.deseq <- dba.report(dbObj, method=DBA_DESEQ2, contrast = 1, th=1)
保存文件
# EdgeR
out <- as.data.frame(comp1.edgeR)
write.table(out, file="D_edgeR.txt", sep="\t", quote=F, col.names = NA)
# DESeq2
out <- as.data.frame(comp1.deseq)
write.table(out, file="D_deseq2.txt", sep="\t", quote=F, col.names = NA)
# Create bed files for each keeping only significant peaks (p < 0.05)
# EdgeR
out <- as.data.frame(comp1.edgeR)
edge.bed <- out[ which(out$p.value < 0.05),
c("seqnames", "start", "end", "strand", "Fold")]
write.table(edge.bed, file="results/Nanog_vs_Pou5f1_edgeR_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)
# DESeq2
out <- as.data.frame(comp1.deseq)
deseq.bed <- out[ which(out$p.value < 0.05),
c("seqnames", "start", "end", "strand", "Fold")]
write.table(deseq.bed, file="results/Nanog_vs_Pou5f1_deseq2_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)我改了博主的FDR参数，主要是因为我师姐的数据FDR不显著，只能换成p.value才能得到一些相应的结果。
差异peak注释
在前面获得了差异peak后，可以对差异peak进行基因注释了，主要是选用的GenomicFeatures、ChIPseeker与ChIPpeakAnno包做的基因注释，需要提前下载基因组的注释文件，才能做接下来的分析。
##R包加载
library("GenomicFeatures")
library("ChIPseeker")
library("ChIPpeakAnno")
##设置工作目录，我为了防止报错，把需要注释的文件都放到下面的这个文件夹下面了
system.file("extdata",package="ChIPpeakAnno")
bed <- system.file("extdata", "D_rep1_deseq2_sig.bed", package="ChIPpeakAnno")
rep1 <- toGRanges(bed, format="BED", header=FALSE)
rep1
##这个可以读gtf或者gff3文件
txdb <-makeTxDbFromGFF(file="E://基因组文件/TAIR_exons.gtf",format="gtf")
annoData <- toGRanges(txdb, feature="gene")
annoData[1:2]
##峰注释
rep1 <- annotatePeak(bed, tssRegion=c(-3000, 3000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=4000)
write.table(as.data.frame(rep1), file="H://ZPP/D_rep1-Anno.xls", sep='\t', quote = F)
#对注释结果进行可视化
plotAnnoPie(rep1)其实如果做的是单个处理的，多重复的，直接就做基因注释就可以了，不用折腾前面的差异peak了。
富集分析
我主要是对GO和KO进行富集分析，因为做的不是模式物种，所以需要提前把ordgb的包hub下来。
// ##R包加载
library(DOSE)
library(topGO)
library(clusterProfiler)
library(pathview)
library(ggplot2)
##没有相关的db包，则需要hub后进行保存进行调用，eg：拟南芥，但是拟南芥是有orgdb的包的
require(AnnotationHub)
hub <- AnnotationHub()
query(hub, "Arabidopsis thaliana")
Ara <- hub[["AH95951"]]
keytypes(Ara)
columns(Ara)
length(keys(Ara))[1]
head(keys(Ara))
saveDb("Ara", file）
##读取已经加载的org.db包
library(AnnotationDbi)
gly=loadDb("D:/R-4.1.3/library/AnnotationDbi/extdata/tair.orgdb")
keytypes(tair)
##基因list文件读取
gene_list <- read.table('H://D-Anno.xls', stringsAsFactors = FALSE,sep='\t',header=T)
head(gene_list)
##转换geneid,有的物种直接用这个就能转换了，有的需要提前转换好
gene.df<-bitr(gene_list$geneId,fromType="TAIR",
toType=c("ENTREZID", "SYMBOL", "REFSEQ"),
OrgDb=Ara)
head(gene.df)
ego = enrichGO(gene_list$GeneID, OrgDb = "gly", keyType="ENTREZID", pvalueCutoff=1, qvalueCutoff=1,ont='all')
head(ego)
write.table(ego,'H://D-go.xls',sep = '\t', quote = FALSE, row.names = FALSE)
##Dotplot visualization
##label_format可以调节字符叠加的问题
plot1=dotplot(ego, split="ONTOLOGY",showCategory = 20,label_format=100) + facet_grid(ONTOLOGY~., scale="free")
ggsave("H://D-ego_dotplot.png",plot=plot1,height= 16,width=14)
plot2= barplot(ego, split="ONTOLOGY",showCategory = 20,label_format=100) + facet_grid(ONTOLOGY~., scale="free")
ggsave("H://D-ego_barplot.png",plot=plot2,height= 16,width=14)
# Multiple samples KEGG analysis
compKEGG = enrichKEGG(gene_list$GeneID, organism="ath",pvalueCutoff=1,keyType='kegg'ath
head(compKEGG)
plot1= dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")
ggsave("H://D-kegg_plot.png",plot=plot1,height= 10,width=14)
plot2= barplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")
ggsave("H://D-kegg_bar.png",plot=plot2,height= 10,width=14)
# Output peakAnnolist file
write.table(compKEGG,'H://D-kegg.xls',sep = '\t', quote = FALSE, row.names = FALSE)可以发现我的富集的阈值设的相对不太严格，主要是害怕有些感兴趣的去掉了，所以在后面可以手动进行筛选。
下面就是motif分析，可以选择自己感兴趣的peak区域、基因区域或者基因启动子区域进行分析，我主要一般是在meme网站进行分析，节约时间。
总结
通过对相关的peak进行分析，可以得到我们杂的蛋白是不是结合到我们想要研究的DNA区域，来对下游的一些功能进行研究。
可以发现整个流程大致可以分为三类 ，第一个就是NGS数据的处理，第二个就是相关结果的可视化，第三个就是根据自己的实验目的去定制自己后续的个性化分析的内容，可以加入其他的公共数据来进行相关的整合，来完善整体的内容，最后加上试验来进行后续的功能验证，完善整个流程。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。