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 位点在转录本上的分布的可视化，在文章里有很多展现：

Deep and accurate detection of m6A RNA modifications using miCLIP2 and m6Aboost machine learning^[1]

上图展示了 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:UsersadminAppDataLocalTempRtmpUxkLmoremotesb805d695958Codezy99-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))

多个分组绘图暂时还没想到怎么解决，欢迎小伙伴留言交流！

参考资料