其他
MetaProfile on Transcript
1peaks metaprofile
Immunoprecipitation （IP）免疫共沉淀技术使用特异性抗体来研究蛋白与其相互作用的分子，结合二代测序技术 NGS（Next Generation Sequenceing）可以在全基因组或转录组层面研究蛋白结合的靶标。在 基因组层面 的技术有经典的 ChIP-seq，较新的 CUT & Run 和 CUT & Tag。在 转录组层面 的技术有 iCLIP-seq、PAR-CLIP-seq、MeRIP-seq、eCLIP-seq 等。
这些数据到最后的分析流程经过 callpeak 步骤都会鉴定到所研究的 protein 结合到的很多 peak 位点，对于转录组层面的 peak 位点在转录本上的分布的可视化，在文章里有很多展现：
上图展示了 peaks 在转录本上不同区域的密度分布图。对于 MeRIP-seq 结合的 peak 的分布可视化有两款工具，MetaPlotR 和 Guitar：
下载方式：
$ wget https://github.com/olarerin/metaPlotR/archive/refs/heads/master.zip
# Guitar
BiocManager::install("Guitar")
2cliProfiler 包
今天介绍一款针对于 CLIP-seq 等结合的 peak 的可视化 R 包：cliProfiler ，功能还是挺丰富的。
以下是包的介绍：
In order to make it easier for the researchers to visualize their high-throughput IP experiment data (peaks), we develop cliProfiler package. The cliProfiler includes seven functions which allow users easily make different profile plots.
安装方式
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("cliProfiler")
作者说已经放到 biocondctor 网站上了，但是我找不到，用 BiocManager 也安装不了，比较奇怪。
既然安装不了，就用 devtools 直接从 githup 安装：
library(devtools)
install_github('Codezy99/cliProfiler')
√ checking for file 'C:\Users\admin\AppData\Local\Temp\RtmpUxkLmo\remotesb805d695958\Codezy99-cliProfiler-5317248/DESCRIPTION' ...
- preparing 'cliProfiler':
√ checking DESCRIPTION meta-information ...
- checking for LF line-endings in source and make files and shell scripts
- checking for empty or unneeded directories
- building 'cliProfiler_0.99.0.tar.gz'
...
** testing if installed package keeps a record of temporary installation path
* DONE (cliProfiler)
使用说明
输入文件格式：
1、需要一个 GRanges 格式储存的 peaks 文件，GRanges 格式可以用 GenomicRanges 包来读取 peak 文件。
BiocManager::install("GenomicRanges")
2、需要 gencode 数据库的 GFF3 格式的注释文件，目前只支持 Human 和 Mouse 的。
我们使用包内自带的测试数据。
1、metaGeneProfile
metaGeneProfile 展示 peaks 在基因组区域的定位，在表观学中使用广泛，特别是 m6A-seq（N6-Methyladenosine）。
# 加载R包
library(cliProfiler)
# 载入测试数据
testpath <- system.file("extdata", package = "cliProfiler")
# 加载注释文件
test_gff3 <- file.path(testpath, "annotation_test.gff3")
# 查看数据类型
class(test)
[1] "GRanges"
attr(,"package")
[1] "GenomicRanges"
# 查看数据内容
head(test,4)
GRanges object with 4 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr17 28748198-28748218 +
[2] chr10 118860137-118860157 -
[3] chr2 148684461-148684481 +
[4] chr2 84602546-84602566 -
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
metaGeneProfile ( )函数输出的是一个 list 类型的数据，第一个 list 包含了 GRanges 格式的计算结果。第二个 list 是 ggplot2 类型的绘图结果，我们可以使用 ggplot2 函数添加并修改内容。
# 计算并绘图
meta <- metaGeneProfile(object = test, annotation = test_gff3)
# 查看list1内容
head(meta[[1]],4)
GRanges object with 4 ranges and 4 metadata columns:
seqnames ranges strand | center location Gene_ID Position
<Rle> <IRanges> <Rle> | <integer> <character> <character> <numeric>
[1] chr10 118860137-118860157 - | 118860147 CDS ENSMUSG00000028630.9 0.674444
[2] chr2 84602546-84602566 - | 84602556 UTR3 ENSMUSG00000034101.14 0.122384
[3] chr18 6111874-6111894 - | 6111884 CDS ENSMUSG00000041225.16 0.199836
[4] chr11 33213145-33213165 - | 33213155 UTR3 ENSMUSG00000040594.19 0.159303
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
结果说明：
center ：peak 中间的位置 location ：peak 属于的基于区域 Gene_ID ：peak 属于哪一个基因 id Position ：每个 peak 的相对位置，其值越接近 0，表示离基因区域的 5'端越近，越接近 1，表示离基因区域的 3'端越近，如果等于 5，表示该 peak 没有被注释到。
library(ggplot2)
meta[[2]] + ggtitle("Meta Profile")
meta[[2]] +
geom_density(size=1.2,color='#D44000') +
geom_vline(xintercept = c(1,2),linetype='dashed',size=1) +
theme_prism() +
ylab('Density') +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.ticks.length.x = unit(0.3,'cm'))
# 或者直接提取数据自己画
# meta_data <- meta[[2]]$data
metaGeneProfile ( )函数提供两个方法来计算相对距离。第一种方法计算的相对距离是基因区域中不包含内含子的，第二种方法计算的相对距离是基于完整的转录本的（包括内含子），我们可以使用 include_intron 参数来转换两种方法。如果你的测序数据不是 plolyA 类型的，推荐使用 include_intron = TRUE 参数。
meta <- metaGeneProfile(object = test, annotation = test_gff3,
include_intron = TRUE)
meta[[2]]
如果有不同的分组，可以使用 group 参数：
test$Treat <- c(rep("Treatment 1",50), rep("Treatment 2", 50))
meta <- metaGeneProfile(object = test, annotation = test_gff3,
group = "Treat")
meta[[2]]
此外还可以对 exlevel 和 extranscript_support_level 进行过滤后再绘图，具体可以参考 gencode 注释文件格式说明。
metaGeneProfile(object = test, annotation = test_gff3, exlevel = 3,
extranscript_support_level = c(4,5,6))
作者还提供一个 split 参数可以分布绘制 5'UTR、CDS、3'UTR 的密度分布图，灰色区域为整个基因区域的密度分布图：
library(ggsci)
meta <- metaGeneProfile(object = test, annotation = test_gff3, split = TRUE)
meta[[2]] +
geom_density(size=1.2) +
geom_vline(xintercept = c(1,2),linetype='dashed',size=1) +
theme_prism() +
ylab('Density') +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.ticks.length.x = unit(0.3,'cm')) +
scale_color_lancet()
2、intronProfile
intronProfile（）函数会计算 peaks 在 intron 区域的分布情况，输出结果跟 metaGeneProfile 类似，输出结果保存在 list 中，第一个 list 为计算结果，第二个为绘图结果。
intron <- intronProfile(test, test_gff3)
# 查看结果
head(intron[[1]],4)
GRanges object with 4 ranges and 7 metadata columns:
seqnames ranges strand | Treat center Intron_S Intron_E Intron_length
<Rle> <IRanges> <Rle> | <character> <integer> <numeric> <numeric> <numeric>
[1] chr17 28748198-28748218 + | Treatment 1 28748208 0 0 0
[2] chr2 148684461-148684481 + | Treatment 1 148684471 0 0 0
[3] chr7 5097955-5097975 + | Treatment 1 5097965 0 0 0
[4] chr4 139648373-139648393 + | Treatment 1 139648383 139645102 139653534 8433
Intron_transcript_id Intron_map
<character> <numeric>
[1] NO 3.000000
[2] NO 3.000000
[3] NO 3.000000
[4] ENSMUST00000178644.1 0.389067
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
计算结果说明：
center ：peak 中间的位置 Intron_S and Intron_E ：内含子的 5'和 3'可变剪切位点（splice site） Intron_length ：peak 所属的内含子的长度 Intron_transcript_id ：intron 属于的 transcript_id Intron_map ：peak 的相对位置，其值越接近 0，表示离 5' SS 位点端越近，越接近 1，表示离 3' SS 位点端越近，如果等于 3，表示该 peak 没有被注释到。
绘图：
# 单个绘图
intron[[2]]
# 多组样本绘图
intron <- intronProfile(test, test_gff3, group = "Treat")
intron[[2]]
intronProfile 也能支持分组和对 exlevel 和 extranscript_support_level 进行过滤和筛选，此外还可以对内含子的最大长度和最小长度 maxLength、minLength 进行设置。
intronProfile(test, test_gff3, group = "Treat", maxLength = 10000,
minLength = 50)
3、exonProfile
exonProfile（）函数可以计算 peaks 在外显子上的分布密度并绘图，输出结果类似同上。
exon <- exonProfile(test, test_gff3)
head(exon[[1]],4)
GRanges object with 4 ranges and 7 metadata columns:
seqnames ranges strand | Treat center exon_S exon_E exon_length
<Rle> <IRanges> <Rle> | <character> <integer> <numeric> <numeric> <numeric>
[1] chr17 28748198-28748218 + | Treatment 1 28748208 28746271 28748406 2136
[2] chr2 148684461-148684481 + | Treatment 1 148684471 148683594 148684968 1375
[3] chr7 5097955-5097975 + | Treatment 1 5097965 5097572 5098178 607
[4] chr4 139648373-139648393 + | Treatment 1 139648383 139648158 139649690 1533
exon_transcript_id exon_map
<character> <numeric>
[1] ENSMUST00000004990.13 0.906835
[2] ENSMUST00000028928.7 0.637818
[3] ENSMUST00000098845.9 0.647446
[4] ENSMUST00000039818.9 0.146771
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
计算结果说明：
center ：peak 中间的位置 exon_S and exon_E ：外显子的 5'和 3'可变剪切位点（splice site） exon_length ：peak 所属的外显子的长度 exon_transcript_id ：外显子属于的 transcript_id exon_map ：peak 的相对位置，其值越接近 0，表示离 5' SS 位点端越近，越接近 1，表示离 3' SS 位点端越近，如果等于 3，表示该 peak 没有被注释到。
绘图：
exon[[2]]
该函数用法与 intronProfile 相似，可参考 intronProfile 函数用法。
4、windowProfile
基于以上三个函数只能用于从 gencode 下载的 gff3 注释文件分析 Human 和 Mouse 物种，如果研究者研究的是其他物种，作者开发了 windowProfile（） 这个函数用于分析。需要 GRanges 格式的注释信息，比如你想做一个一个 exon profile 的分析，你只需要从 GRanges 里提取所有的外显子区域，然后运行 windowProfile 函数就行了。
library(rtracklayer)
# 提取所有外显子的注释信息
test_anno <- rtracklayer::import.gff3(test_gff3)
test_anno <- test_anno[test_anno$type == "exon"]
# 计算分析
window_profile <- windowProfile(test, test_anno)
window_profile <- windowProfile(test, test_anno)
# 查看结果
head(window_profile[[1]],4)
GRanges object with 4 ranges and 6 metadata columns:
seqnames ranges strand | Treat center window_S window_E window_length
<Rle> <IRanges> <Rle> | <character> <integer> <numeric> <numeric> <numeric>
[1] chr17 28748198-28748218 + | Treatment 1 28748208 28746271 28748406 2136
[2] chr2 148684461-148684481 + | Treatment 1 148684471 148683594 148684968 1375
[3] chr7 5097955-5097975 + | Treatment 1 5097965 5097572 5098178 607
[4] chr4 139648373-139648393 + | Treatment 1 139648383 139648158 139649690 1533
window_map
<numeric>
[1] 0.906835
[2] 0.637818
[3] 0.647446
[4] 0.146771
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
计算结果说明：
center ：peak 中间的位置 window_S and window_E ：peak 所分布位置的边界 window_length ：peak 所属的基因结构特征的长度 window_map ：peak 的相对位置，其值越接近 0，表示离基因结构特征的起始位点端越近，越接近 1，表示离基因结构特征的终止位点端越近，如果等于 3，表示该 peak 没有被注释到。
绘图：
window_profile[[2]]
5、motifProfile
motifProfile（） 函数可用来做 motif 富集分析，分析特定 motif 在 peaks 中心的富集情况，需要提前安装好相应的基因组文件。fraction 参数能把 y 轴从 frequency 到 frequency。
motif <- motifProfile(test, motif = "DRACH",
genome = "BSgenome.Mmusculus.UCSC.mm10",
fraction = TRUE, title = "Motif Profile",
flanking = 10)
Attaching package: ‘Biostrings’
The following object is masked from ‘package:base’:
strsplit
输出结果为 list 文件，储存了 motif 富集分析的 data.frame 结果和 ggplot2 的绘图结果。
# 查看 motif enrichment结果
head(motif[[1]],4)
Position Fraction
5 -10 0.02
6 -9 0.04
7 -8 0.04
8 -7 0.02
# 绘图
motif[[2]]
6、geneTypeProfile
geneTypeProfile（）函数利用 gff 注释信息对 peaks 进行 genetype 注释，输出结果和类型同上。
geneTP <- geneTypeProfile(test, test_gff3)
# 查看内容
head(geneTP[[1]],4)
GRanges object with 4 ranges and 4 metadata columns:
seqnames ranges strand | Treat center geneType Gene_ID
<Rle> <IRanges> <Rle> | <character> <integer> <character> <character>
[1] chr17 28748198-28748218 + | Treatment 1 28748208 protein_coding ENSMUSG00000053436.15
[2] chr10 118860137-118860157 - | Treatment 1 118860147 protein_coding ENSMUSG00000028630.9
[3] chr2 148684461-148684481 + | Treatment 1 148684471 protein_coding ENSMUSG00000027439.9
[4] chr2 84602546-84602566 - | Treatment 1 84602556 protein_coding ENSMUSG00000034101.14
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
# 绘图
geneTP[[2]]
结果说明：
center ：peak 中间的位置 geneType ：peak 所属的基因类型 Gene_ID ：peak 所属的基因 id
7、spliceSiteProfile
spliceSiteProfile（）以绝对距离提供 peaks 在 5' SS（剪切位点）和 3' SS（剪切位点）的富集程度信息。输出结果和类型同上。
SSprofile <- spliceSiteProfile(test, test_gff3, flanking=200, bin=40)
# 查看结果
head(SSprofile[[1]],4)
GRanges object with 4 ranges and 4 metadata columns:
seqnames ranges strand | Treat center Position5SS Position3SS
<Rle> <IRanges> <Rle> | <character> <integer> <integer> <integer>
[1] chr10 118860137-118860157 - | Treatment 1 118860147 <NA> <NA>
[2] chr2 84602546-84602566 - | Treatment 1 84602556 <NA> <NA>
[3] chr18 6111874-6111894 - | Treatment 1 6111884 -200 <NA>
[4] chr11 33213145-33213165 - | Treatment 1 33213155 <NA> <NA>
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
# 绘图
SSprofile[[2]]
结果说明：
center ：peak 中间的位置 Position5SS ：peak 和 5' SS（剪切位点）的绝对距离，负值表示 peak 位于外显子区域上，正值表示 peak 在内含子区域上 Position3SS ：peak 和 3' SS（剪切位点）的绝对距离，负值表示 peak 位于内含子区域上，正值表示 peak 在外显子区域上
示意图：
以上就是包的大部分函数的介绍了，感兴趣的小伙伴可以去尝试一下。
3尝试自己的数据
从 GEO 数据库里下载了 MeRIP-seq 的两个 peaks 数据，GEO 编号为：GSE132306，物种是人的。
下载 gff 注释文件：
$ wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gff3.gz
# 解压
$ gunzip gencode.v38.annotation.gff3.gz
导入 R 里：
# 加载R包
BiocManager::install("GenomicRanges")
library(GenomicRanges)
library(ggplot2)
# 读取gtf文件
gff3_human <- 'gencode.v38.annotation.gff3'
# 读取peaks文件
m3kd <- read.table('M3KD_Ctl_peak.txt',header = T)
# 使用 GRanges 函数把 peaks 文件变为 GRanges 格式的
m3kd_peak <- GRanges(seqnames = m3kd$chr,
ranges = IRanges(start = m3kd$start,end = m3kd$end,names = m3kd$name),
strand = m3kd$strand)
计算和绘图：
# 计算相对位置
meta <- metaGeneProfile(object = m3kd_peak, annotation = gff3_human)
# 查看计算结果
head(meta[[1]],4)
GRanges object with 4 ranges and 4 metadata columns:
seqnames ranges strand | center location Gene_ID Position
<Rle> <IRanges> <Rle> | <integer> <character> <character> <numeric>
AAMP chr2 218264179-218264228 - | 218264204 UTR3 ENSG00000127837.10 0.873762
AAMP chr2 218269385-218270110 - | 218269748 NO Nan 5.000000
AARS chr16 70252444-70253303 - | 70252874 CDS ENSG00000090861.17 1.094943
AARS chr16 70253736-70253785 - | 70253761 CDS ENSG00000090861.17 0.888889
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
# 绘图
meta[[2]] +
geom_density(size=1.2,color='#D44000') +
geom_vline(xintercept = c(1,2),linetype='dashed',size=1) +
theme_prism() +
ylab('Density') +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.ticks.length.x = unit(0.3,'cm'))
# 分割分别画图，加上 split=TRUE 参数
meta <- metaGeneProfile(object = m3kd_peak, annotation = gff3_human,split = T)
# 绘图
meta[[2]] +
geom_density(size=1.2) +
geom_vline(xintercept = c(1,2),linetype='dashed',size=1) +
theme_prism() +
ylab('Density') +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 18),
axis.ticks.length.x = unit(0.3,'cm')) +
scale_color_lancet()
画个 genetype 的分布：
geneTP <- geneTypeProfile(m3kd_peak, gff3_human)
# 查看内容
head(geneTP[[1]],4)
GRanges object with 4 ranges and 3 metadata columns:
seqnames ranges strand | center geneType Gene_ID
<Rle> <IRanges> <Rle> | <integer> <character> <character>
AACS chr12 125065429-125065626 + | 125065528 protein_coding ENSG00000081760.17
AACS chr12 125065676-125073883 + | 125069780 protein_coding ENSG00000081760.17
AACS chr12 125142186-125142383 + | 125142285 protein_coding ENSG00000081760.17
AACS chr12 125142433-125142482 + | 125142458 protein_coding ENSG00000081760.17
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
# 绘图
geneTP[[2]] +
theme_prism(base_size = 18,border = T,axis_text_angle = 45) +
theme(axis.title.x = element_blank(),
axis.title.y = element_text(size = 18))
多个分组绘图暂时还没想到怎么解决，欢迎小伙伴留言交流！
参考资料
Deep and accurate detection of m6A RNA modifications using miCLIP2 and m6Aboost machine learning: https://www.biorxiv.org/content/10.1101/2020.12.20.423675v1
欢迎小伙伴留言评论！
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，打赏一下吧！