CHIPseek笔记记录

2022-08-07 09:30:58 阅读：308 来源： 互联网

标签：记录 CHIPseek NC 笔记 sed pkg packages CHIPseeker TxDb

有人在进行reads比对的时候，参考基因组用的是RefSeq的基因组(如斑马鱼的基因组)，可以发现比对出来，染色体的编号是如图所示：但据我了解(知道的朋友可以交流)，CHIPseeker仅识别染色体编号为chrX（X代表数字，且字母要严格小写）,因此要顺利分析进行，必须进行编号改变

由下图可以发现上染色体编号的对应关系(所以直接替换就行)

# linux 批量修改文件的对应内容
for i in `ls *.fas`; do echo "sed 's/NC_007112.7/chr1/g' $i |sed 's/NC_007113.7/chr2/g' |sed 's/NC_007114.7/chr3/g' |sed 's/NC_007115.7/chr4/g' |sed 's/NC_007116.7/chr5/g' |sed 's/NC_007118.7/chr7/g' |sed 's/NC_007119.7/chr8/g' |sed 's/NC_007120.7/chr9/g' |sed 's/NC_007121.7/chr10/g' |sed 's/NC_007122.7/chr11/g' |sed 's/NC_007123.7/chr12/g' |sed 's/NC_007124.7/chr13/g' |sed 's/NC_007125.7/chr14/g' |sed 's/NC_007126.7/chr15/g' |sed 's/NC_007127.7/chr16/g' |sed 's/NC_007128.7/chr17/g' |sed 's/NC_007129.7/chr18/g' |sed 's/NC_007130.7/chr19/g' |sed 's/NC_007131.7/chr20/g' |sed 's/NC_007132.7/chr21/g' |sed 's/NC_007133.7/chr22/g' |sed 's/NC_007134.7/chr23/g' |sed 's/NC_007135.7/chr24/g' |sed 's/NC_007136.7/chr25/g' |sed 's/NC_002333.2/chrMT/g' > $i.re"; done >> run.sh

结果：

对不同物种进行分析，需要用到不同的注释库，如下方式可以得到对应的库(OrgDb和TxDb库)

分析代码:

# CHIPseeker_zebrafish
# ly-2021-07-29

# org.Xx.eg.db系列包概述: http://www.bio-info-trainee.com/788.html

cran_packages <- c("BiocManager",
                   "reshape",
                   "Vennerable",
                   "RMariaDB",
                   "UpSetR",
                   "optparse",
                   "stringr",
                   "ggplot2")

Biocdouctor_packages <- c("RBGL",
                          "graph",                          
                          "TxDb.Hsapiens.UCSC.hg19.knownGene", #CHIPseeker的依赖
                          "TxDb.Hsapiens.UCSC.hg38.knownGene",
                           "CHIPseeker",
　　　　　　　　　　　　　　　　"TxDb.Drerio.UCSC.danRer11.refGene",
　　　　　　　　　　　　　　　　"clusterProfiler",
　　　　　　　　　　　　　　　　"org.Dr.eg.db")

# 加载或安装必须包
for(pkg in cran_packages){
  if(!require(pkg, character.only = T)){
    install.packages(pkg, ask=F, update=F)
    require(pkg, character.only = T)
  }
}
for(pkg in Biocdouctor_packages){
  if(!require(pkg, character.only = T)){
    BiocManager::install(pkg, ask = F, update = F)
    require(pkg, character.only = T)
  }
}

# 加载org.Dr.eg.db包可能会报错:$ operator is invalid for atomic vectors
options(connectionObserver = NULL)
library("org.Dr.eg.db")  #不同物种对应的库有差别

# 安装Vennerable可能报错
install.packages("Vennerable", repos="http:/ /R-Forge.R-project.org")
library("Vennerable")

# 确定是否能加载所有包
for (pkg in c(cran_packages, Biocdouctor_packages)) {
  require(pkg, character.only = T)
}
library(CHIPseeker) #前面加在可能会出错
search() #[43]

# 开始分析
files <- getSampleFiles() # 读取例子文件
print(files) 
peak <- readPeakFile(files[[4]]) # peak calling
peak

## 自己的bed文件（注意染色体名字要求）
#peak <- readPeakFile("ctrl_1_summits.bed")

# 首先可视化自己的bed文件(支持多个bed文件同时可视化：扩展阅读2)
# ?covplot
covplot(peak, weightCol = "V5",
        chrs = c("chr17", "chr18"), # 默认整个
        xlim = c(4.5e7, 5e7)) #weightCol展现表达量

# peak的注释  
TxDb <- TxDb.Hsapiens.UCSC.hg38.knownGene #这里究竟该用什么版本，一定要跟上有分析对应
promter <- getPromoters(TxDb = TxDb, upstream = 3000, downstream = 3000) # 默认的是TxDb.Hsapiens.UCSC.hg19.knownGene
# filelist <- c("test")
# temp = get(filelist)
PeakAnno <- annotatePeak(peak,tssRegion = c(-3000, 3000),annoDb = "org.Hs.eg.db",TxDb=TxDb) #一定要在bed文件的染色体编号前加chr,否则会出现如下报错： Error in (function (classes, fdef, mtable)  : unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’

#自己定义promoter区域，上下游3000bp
promoter <- getPromoters(upstream=3000, downstream=3000)
#不理解这个函数也没关系，是为下一步做热图提供matrix
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

#一键绘图，效果同上
# peakHeatmap(peak, TxDb=txdb, upstream=3000, downstream=3000, color="red")

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
#加一个置信区间
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)


peak
#共计1331个peak
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000))
peakAnno


class(peakAnno)
peakAnno.df <- as.data.frame(peakAnno)
peakAnno.gr <- as.GRanges(peakAnno)
head(peakAnno.gr, 3)

plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
vennpie(peakAnno)

install.packages("magick")
install.packages("Rcpp")
library(Rcpp)
library(magick)
ChIPseeker::upsetplot(peakAnno, vennpie=TRUE)
upsetplot(peakAnno, vennpie=TRUE) #也行

其他都没问题，就在最后一步用upsetplot(peakAnno, vennpie=TRUE)画图时有点问题(运行就一直卡着不动，但又不报错)：使用该函数时，需要依赖包magick和Rcpp，第一次安装之后，发现也不能加载，但不会报错(一直卡着不动)，最后的解决办法是重新安装这两个包

在linux上安装CHIPseeker软件报错：

安装CHIPseeker时，会依赖TxDb.Hsapiens.UCSC.hg19.knownGene包，但安装时出现如图报错。起初以为是网络的问题，其实不是。解决办法为：Win+R，输入inetcpl.cpl 直接打开Internet选项。打开后，在高级中勾选使用TSL 3.0、使用TLS 1.0、使用TLS 1.1、使用TLS 1.2。

扩展阅读:

1、bioconductor之ChIPseeker学习

2、学习一遍ChIPseeker的使用

标签：记录,CHIPseek,NC,笔记,sed,pkg,packages,CHIPseeker,TxDb
来源： https://www.cnblogs.com/ly-zy/p/15077007.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

CHIPseek笔记记录

分析代码: