LaTex的安装并不需要管理员权限

You do not need to be root (administrator on Windows) to install, use, or manage TeX Live. Indeed, we recommend installing it as a normal user, except perhaps on MacOSX, where it’s conventional to install as administrator.

下面是Ubuntu的安装方法

第一步: 下载ISO文件(清华镜像源)

下载地址: https://mirrors.tuna.tsinghua.edu.cn/CTAN/systems/texlive/Images/texlive2020-20200406.iso

其他镜像站点，可以在 https://ctan.org/mirrors 进行选择

第二步: 在Windows上解压缩ISO文件，然后上传到服务器, 我们假设上传到家目录下的texlive

第三步: 运行安装命令

1
2

cd texlive
perl install-t 

输出信息如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

======================> TeX Live installation procedure <=====================

======>   Letters/digits in <angle brackets> indicate   <=======
======>   menu items for actions or customizations      <=======

 Detected platform: GNU/Linux on x86_64
 
 <B> set binary platforms: 1 out of 6

 <S> set installation scheme: scheme-full

 <C> set installation collections:
     40 collections out of 41, disk space required: 6516 MB

 <D> set directories:
   TEXDIR (the main TeX directory):
     !! default location: /usr/local/texlive/2020
     !! is not writable or not allowed, please select a different one!
   TEXMFLOCAL (directory for site-wide local files):
     /usr/local/texlive/texmf-local
   TEXMFSYSVAR (directory for variable and automatically generated data):
     /usr/local/texlive/2020/texmf-var
   TEXMFSYSCONFIG (directory for local config):
     /usr/local/texlive/2020/texmf-config
   TEXMFVAR (personal directory for variable and automatically generated data):
     ~/.texlive2020/texmf-var
   TEXMFCONFIG (personal directory for local config):
     ~/.texlive2020/texmf-config
   TEXMFHOME (directory for user-specific files):
     ~/texmf

 <O> options:
   [ ] use letter size instead of A4 by default
   [X] allow execution of restricted list of programs via \write18
   [X] create all format files
   [X] install macro/font doc tree
   [X] install macro/font source tree
   [ ] create symlinks to standard directories
   [X] after install, set CTAN as source for package updates

 <V> set up for portable installation

Actions:
 <I> start installation to hard disk
 <P> save installation profile to 'texlive.profile' and exit
 <H> help
 <Q> quit

Enter command: 

我们需要输入D进入更改界面，然后输入1，输入我们目标安装路径，输入R回到主目录，最后输入I进行安装。

最后安装结果之后会有如下提示

1
2
3
4
5
6

Add /home/xzg/opt/texlive/texmf-dist/doc/man to MANPATH.
Add /home/xzg/opt/texlive/texmf-dist/doc/info to INFOPATH.
Most importantly, add /home/xzg/opt/texlive/bin/x86_64-linux
to your PATH for current and future sessions.
Logfile: /home/xzg/opt/texlive/install-tl.log

即在 ~/.bashrc 或 ~/.zshrc 这些配置文件下添加如下内容

1
2
3

export MANPATH=$MANPATH:/home/xzg/opt/texlive/texmf-dist/doc/man
export INFOPATH=$INFOPATH:/home/xzg/opt/texlive/texmf-dist/doc/info
export PATH=$PATH:/home/xzg/opt/texlive/bin/x86_64-linux

准备阶段

为了实现该爬虫我们需要用到如下工具

• Chrome浏览器
• Python 3 语法知识
• Python的Requests库

此外，这个爬取程序利用的是微信公众号后台编辑素材界面。原理是，当我们在插入超链接时，微信会调用专门的API（见下图），以获取指定公众号的文章列表。因此，我们还需要有一个公众号。

正式开始

我们需要登录微信公众号，点击素材管理，点击新建图文消息，然后点击上方的超链接。

接着，按F12，打开Chrome的开发者工具，选择Network

此时在之前的超链接界面中，点击「选择其他公众号」，输入你需要爬取的公众号（例如中国移动）

此时之前的Network就会刷新出一些链接，其中以”appmsg”开头的便是我们需要分析的内容

我们解析请求的URL

1

https://mp.weixin.qq.com/cgi-bin/appmsg?action=list_ex&begin=0&count=5&fakeid=MzI1MjU5MjMzNA==&type=9&query=&token=143406284&lang=zh_CN&f=json&ajax=1

它分为三个部分

• https://mp.weixin.qq.com/cgi-bin/appmsg: 请求的基础部分
• ?action=list_ex: 常用于动态网站，实现不同的参数值而生成不同的页面或者返回不同的结果
• &begin=0&count=5&fakeid: 用于设置?里的参数，即begin=0, count=5

通过不断的浏览下一页，我们发现每次只有begin会发生变动，每次增加5，也就是count的值。

接着，我们通过Python来获取同样的资源，但直接运行如下代码是无法获取资源的

1
2
3
4

import requests
url = "https://mp.weixin.qq.com/cgi-bin/appmsg?action=list_ex&begin=0&count=5&fakeid=MzI1MjU5MjMzNA==&type=9&query=&token=1957521839&lang=zh_CN&f=json&ajax=1"
requests.get(url).json() 
# {'base_resp': {'ret': 200003, 'err_msg': 'invalid session'}}

我们之所以能在浏览器上获取资源，是因为我们登录了微信公众号后端。而Python并没有我们的登录信息，所以请求是无效的。我们需要在requests中设置headers参数，在其中传入Cookie和User-Agent，来模拟登陆

由于每次头信息内容都会变动，因此我将这些内容放入在单独的文件中，即”wechat.yaml”，信息如下

1
2

cookie:  ua_id=wuzWM9FKE14...
user_agent: Mozilla/5.0...

之后只需要读取即可

1
2
3
4
5
6
7
8
9
10
11
12

# 读取cookie和user_agent
import yaml
with open("wechat.yaml", "r") as file:
    file_data = file.read()
config = yaml.safe_load(file_data) 

headers = {
    "Cookie": config['cookie'],
    "User-Agent": config['user_agent'] 
}

requests.get(url, headers=headers, verify=False).json()

在返回的JSON中，我们就看到了每个文章的标题(title), 摘要(digest), 链接(link), 推送时间(update_time)和封面地址(cover)等信息。

appmsgid是每一次推送的唯一标识符，aid则是每篇推文的唯一标识符。

实际上，除了Cookie外，URL中的token参数也会用来限制爬虫，因此上述代码很有可能输出会是{'base_resp': {'ret': 200040, 'err_msg': 'invalid csrf token'}}

接着我们写一个循环，获取所有文章的JSON，并进行保存。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

import json
import requests
import time
import random

import yaml
with open("wechat.yaml", "r") as file:
    file_data = file.read()
config = yaml.safe_load(file_data) 

headers = {
    "Cookie": config['cookie'],
    "User-Agent": config['user_agent'] 
}

# 请求参数
url = "https://mp.weixin.qq.com/cgi-bin/appmsg"
begin = "0"
params = {
    "action": "list_ex",
    "begin": begin,
    "count": "5",
    "fakeid": config['fakeid'],
    "type": "9",
    "token": config['token'],
    "lang": "zh_CN",
    "f": "json",
    "ajax": "1"
}

# 存放结果
app_msg_list = []
# 在不知道公众号有多少文章的情况下，使用while语句
# 也方便重新运行时设置页数
i = 0
while True:
    begin = i * 5
    params["begin"] = str(begin)
    # 随机暂停几秒，避免过快的请求导致过快的被查到
    time.sleep(random.randint(1,10))
    resp = requests.get(url, headers=headers, params = params, verify=False)
    # 微信流量控制, 退出
    if resp.json()['base_resp']['ret'] == 200013:
        print("frequencey control, stop at {}".format(str(begin)))
        break
    
    # 如果返回的内容中为空则结束
    if len(resp.json()['app_msg_list']) == 0:
        print("all ariticle parsed")
        break
        
    app_msg_list.append(resp.json())
    # 翻页
    i += 1

在上面代码中，我将fakeid和token也存放在了”wechat.yaml”文件中，这是因为fakeid是每个公众号都特有的标识符，而token则会经常性变动，该信息既可以通过解析URL获取，也可以从开发者工具中查看

在爬取一段时间后，就会遇到如下的问题

1

{'base_resp': {'err_msg': 'freq control', 'ret': 200013}}

此时你在公众号后台尝试插入超链接时就能遇到如下这个提示

这是公众号的流量限制，通常需要等上30-60分钟才能继续。为了完美处理这个问题，你可能需要申请多个公众号，可能需要和微信公众号的登录系统斗智斗勇，或许还需要设置代理池。

但是我并不需要一个工业级别的爬虫，只想爬取自己公众号的信息，因此等个一小时，重新登录公众号，获取cookie和token，然后运行即可。我可不想用自己的兴趣挑战别人的饭碗。

最后将结果以JSON格式保存。

1
2
3
4

# 保存结果为JSON
json_name = "mp_data_{}.json".format(str(begin))
with open(json_name, "w") as file:
    file.write(json.dumps(app_msg_list, indent=2, ensure_ascii=False))

或者提取文章标识符，标题，URL，发布时间这四列信息，保存成CSV。

1
2
3
4
5
6
7
8
9

info_list = []
for msg in app_msg_list:
    if "app_msg_list" in msg:
        for item in msg["app_msg_list"]:
            info = '"{}","{}","{}","{}"'.format(str(item["aid"]), item['title'], item['link'], str(item['create_time']))
            info_list.append(info)
# save as csv
with open("app_msg_list.csv", "w") as file:
    file.writelines("\n".join(info_list))         

下一篇，将介绍如何根据每个文章的连接地址，来获取每篇文章的阅读量信息。

参考资料

• https://blog.csdn.net/kindred_joe/article/details/99289890
• https://blog.csdn.net/qq_28804275/article/details/82150874

最终代码如下，使用方法为python wechat_parser.py wechat.yaml

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

import json
import requests
import time
import random
import os
import yaml
import sys

if len(sys.argv) < 2:
    print("too few arguments")
    sys.exit(1)

yaml_file = sys.argv[1]
if not os.path.exists(yaml_file):
    print("yaml_file is not exists")
    sys.exit(1)
    

with open(yaml_file, "r") as file:
    file_data = file.read()
config = yaml.safe_load(file_data)

headers = {
    "Cookie": config['cookie'],
    "User-Agent": config['user_agent'] 
}

# 请求参数
url = "https://mp.weixin.qq.com/cgi-bin/appmsg"
begin = "0"
params = {
    "action": "list_ex",
    "begin": begin,
    "count": "5",
    "fakeid": config['fakeid'],
    "type": "9",
    "token": config['token'],
    "lang": "zh_CN",
    "f": "json",
    "ajax": "1"
}

# 存放结果
if os.path.exists("mp_data.json"):
    with open("mp_data.json", "r") as file:
        app_msg_list = json.load(file)
else:
    app_msg_list = []
# 在不知道公众号有多少文章的情况下，使用while语句
# 也方便重新运行时设置页数
i = len(app_msg_list) // 5
while True:
    begin = i * 5
    params["begin"] = str(begin)
    # 随机暂停几秒，避免过快的请求导致过快的被查到
    time.sleep(random.randint(1,10))
    resp = requests.get(url, headers=headers, params = params, verify=False)
    # 微信流量控制, 退出
    if resp.json()['base_resp']['ret'] == 200013:
        print("frequencey control, stop at {}".format(str(begin)))
        break
    
    # 如果返回的内容中为空则结束
    if len(resp.json()['app_msg_list']) == 0:
        print("all ariticle parsed")
        break
        
    app_msg_list.append(resp.json())
    # 翻页
    i += 1

# 保存结果为JSON
json_name = "mp_data.json"
with open(json_name, "w") as file:
    file.write(json.dumps(app_msg_list, indent=2, ensure_ascii=False))

尽管目前已测序的物种已经很多了，但是对于一些动辄几个G的复杂基因组，还没有某个课题组有那么大的经费去测序，所以仍旧缺少高质量的完整基因组，那么这个时候一个高质量的转录组还是能够凑合用的。

二代测序的组装结果只能是差强人意，最好的结果就是不要组装，直把原模原样的转录组给你是最好的。PacBio Iso-Seq 做的就是这件事情。只不过Iso-Seq测得是转录本，由于有些基因存在可变剪切现象，所以所有将一个基因的所有转录本放在一起看，搞出一个尽量完整的编码区。

Cogent能使用高质量的全长转录组测序数据对基因组编码区进行重建，示意图如下：

软件安装

虽然说软件可以直接通过conda进行安装，但是根据官方文档的流程，感觉还是很麻烦

下面操作默认你安装好了miniconda3, 且miniconda3的位置在家目录下

安装方法如下（更新与2020/7/30）

1
2
3
4
5
6
7
8
9
10
11
12
13

conda create -y -n cogent python=3.7 
conda activate cogent
conda install -y  -c bcbio bx-python
conda install -y  -c conda-forge pulp
pip install networkx==2.0
conda install -y -c bioconda parasail-python

# 安装Github的Cogent
cd ~/miniconda3/envs/cogent
git clone https://github.com/Magdoll/Cogent.git
cd Cogent
python setup.py build
python setup.py install

经过上面一波操作，请用下面这行命令验证是否安装成功

1
2

run_mash.py --version
run_mash.py Cogent 6.0.0

此外还需要安装另外两个软件，分别是Minimap2和Mash

1
2
3
4
5

conda install minimap2
# 目前mash升级了
wget  https://github.com/marbl/Mash/releases/download/v2.2/mash-Linux64-v2.2.tar
tar xf mash-Linux64-v2.2.tar
mv mash-Linux64-v2.2/mash  $HOME/miniconda3/envs/cogent/bin

对待大数据集，你还需要安装Cupcake

1
2
3
4
5
6
7

conda activate cogent 
conda install cython -y
git clone https://github.com/Magdoll/cDNA_Cupcake.git
cd cDNA_Cupcake
git checkout
python setup.py build
python setup.py install

关于这个Cupcake, 如果python3环境出现报错，就新建一个Python2环境进行安装。

创建伪基因组

让我们先下载测试所用的数据集，

1
2
3

mkdir test
cd test
wget https://raw.githubusercontent.com/Magdoll/Cogent/master/test_data/test_human.fa

第一步： 从数据集中搜索同一个基因簇(gene family)的序列

这一步分为超过20,000 条序列的大数据集和低于20,000条序列这两种情况, 虽然无论那种情况，在这里我们都只用刚下载的测试数据集

小数据集

第一步：从输入中计算k-mer谱和配对距离

1
2
3

run_mash.py -k 30 --cpus=12 test_human.fa
# -k, k-mer大小
# --cpus， 进程数

你一定要保证你的输入是fasta格式，因为该工具目前无法自动判断输入是否是fasta格式，所以当你提供诡异的输入时，它会报错，然后继续在后台折腾。

上面这行命令做的事情是：

• 将输入的fasta文件进行拆分，你可以用--chunk_size指定每个分块的大小
• 对于每个分块计算k-mer谱
• 对于每个配对的分块，计算k-mer相似度（距离）
• 将分块合并成单个输出文件

第二步：处理距离文件，创建基因簇

1

process_kmer_to_graph.py  test_human.fa test_human.fa.s1000k30.dist test_human human

这里的test_human是输出文件夹，human表示输出文件名前缀。此外如果你有输入文件中每个转录亚型的丰度，那么你可以用-c参数指定该文件。

这一步会得到输出日志test_human.partition.txt，以及test_human下有每个基因family的次文件夹，文件里包含着每个基因簇的相似序列。对于不属于任何基因family的序列，会在日志文件种专门说明，这里是human.partition.txt

大数据集

如果是超过20,000点大数据集，分析就需要用到Minimap2和Cpucake了。分为如下几个步骤：

• 使用minimap2对数据进行粗分组
• 对于每个分组，使用上面提到的精细的寻找family工具
• 最后将每个分组的结果进行合并

第一步：使用minimap进行分组

run_preCluster.py要求输入文件名为isoseq_flnc.fasta, 所以需要先进行软连接

1
2

ln -s test_human.fa isoseq_flnc.fasta 
run_preCluster.py --cpus=20

为每个分组构建基因簇寻找命令

1

generate_batch_cmd_for_Cogent_family_finding.py --cpus=12 --cmd_filename=cmd preCluster.cluster_info.csv preCluster_out test_human

得到的是 cmd 文件，这个cmd可以直接用bash cmd运行，也可以投递到任务调度系统。

最后将结果进行合并

1
2

printf "Partition\tSize\tMembers\n" > final.partition.txt
ls preCluster_out/*/*partition.txt | xargs -n1 -i sed '1d; $d' {} | cat >> final.partition.txt

第二步：重建编码基因组

上一步得到每个基因簇, 可以分别重构编码基因组，所用的命令是reconstruct_contig.py

1
2
3
4

使用方法: reconstruct_contig.py [-h] [-k KMER_SIZE]
                             [-p OUTPUT_PREFIX] [-G GENOME_FASTA_MMI]
                             [-S SPECIES_NAME]
                             dirname

如果你手头有一个质量不是很高的基因组，可以使用-G GENOME_FASTA_MMI和-S SPECIES_NAME提供参考基因组的信息。毕竟有一些内含子是所有转录本都缺少的，提供基因组信息，可以补充这部分缺失。

由于有多个基因簇，所以还需要批量运行命令 reconstruct_contig.py

1
2

generate_batch_cmd_for_Cogent_reconstruction.py test_human > batch_recont.sh
bash batch_recont.sh

第三步: 创建伪基因组

首先获取未分配的序列， 这里用到get_seqs_from_list.py脚本来自于Cupcake, 你需要将cDNA_Cupcake/sequence添加到的环境变量PATH中。

1
2

tail -n 1 human.partition.txt | tr ',' '\n'  > unassigned.list
get_seqs_from_list.py hq_isoforms.fasta unassigned.list > unassigned.fasta

这里的测试数据集里并没有未分配的序列，所以上面这一步可省去。

然后将未分配的序列和Cogent contigs进行合并

1
2
3

mkdir collected
cd collected
cat ../test_human/*/cogent2.renamed.fasta ../unassigned.fasta > cogent.fake_genome.fasta

最后，将冗余的转录亚型进行合并。 这一步我们需要用Minimap2或者GMAP将Iso-seq分许得到的hq_isoforms.fasta回帖到我们的伪基因组上。 关于参数选择，请阅读Best practice for aligning Iso Seq to reference genome: minimap2, GMAP, STAR, BLAT

1
2
3
4

ln -s ../test_human.fa hq_isoforms.fasta
minimap2 -ax splice -t 30 -uf --secondary=no \
  cogent.fake_genome.fasta hq_isoforms.fasta > \
   hq_isoforms.fasta.sam

然后可以根据collapse tutorial from Cupcake将冗余的转录亚型合并

1
2
3
4

sort -k 3,3 -k 4,4n hq_isoforms.fasta.sam > hq_isoforms.fasta.sorted.sam
collapse_isoforms_by_sam.py --input hq_isoforms.fasta -s hq_isoforms.fasta.sorted.sam \
           -c 0.95 -i 0.85 --dun-merge-5-shorter \
           -o hq_isoforms.fasta.no5merge

由于这里没有每个转录亚型的丰度文件cluster_report.csv，所以下面的命令不用运行, 最终结果就是hq_isoforms.fasta.no5merge.collapsed.rep.fa

1
2

get_abundance_post_collapse.py hq_isoforms.fasta.no5merge.collapsed cluster_report.csv
filter_away_subset.py hq_isoforms.fasta.no5merge.collapsed

如果运行了上面这行代码，那么最终文件就应是hq_isoforms.fasta.no5merge.collapsed.filtered.*

第13章: ArchR的ChromVAR偏离富集分析

正如之前章节所提到的，TF motif富集可以帮助我们预测在我们感兴趣的细胞类型中哪些调控元件最为活跃。只不过这些富集既不是根据每个计算进行计算，也没有考虑到Tn5转座酶的序列偏好性。GreenLeaf Lab开发的 chromVARR包就是为了解决这些问题。chromVAR的设计目标就是要根据每个细胞的稀疏染色质开放数据来预测TF活跃富集情况。chromVAR的两个主要输出为

1. “deviations”: deviation(偏离)是偏好校正值，它根据所有细胞或样本的预期开放度评估给定特征(例如motif)在每个细胞中的偏离程度。
2. “z-score”: z-score也称之为”deviation score”, 是所有细胞中每个偏好纠正偏离值的z-score。偏离得分的绝对值和每个细胞的read深度相关。这是因为read越多，对于给定特征(例如motif)在每个细胞中相对于预期值存在的差异，你会更有把握，觉得它不可能是随机事情。

chromVAR的一个主要局限在于它是为早期scATAC-seq数据开发的，那个时候的实验只有上百个细胞。对于目前的数据规模，chromVAR很难将所有的cell-by-peak矩阵读取到内存中快速计算TF偏离。并且，当前的实验会输出成千上万个细胞，产生的cell-by-peak矩阵很难加载到内存中。于是，即便是中等大小数据(50,000个细胞)，它的运行速度和内存占用都会急剧增加。

为了规避这些局限，ArchR通过分析独立的分析样本的每个子表达矩阵来实现相同的chromVAR分析流程。

首先，ArchR读取每个子样本中所有细胞的全局开放性。然后，对于每个peak，ArchR根据GC含量和开放性确定一组背景peak。其次，对于每个样本，ArchR单独使用chromVAR根据背景peak集和全局开放性计算偏好校正后的离值。这种计算方式每次只会加载5,000-10,000细胞到内存中，因此降低了内存消耗。最终，我们能够在大规模数据中应用chromVAR并提升了运行性能。

13.1 Motif偏离

首先，我们要确保我们已经在ArchRProject中添加了motif注释

1
2
3

if("Motif" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")
}

此外，我们还需要添加一组背景peak用于计算偏离。背景peak通过chromVAR::getBackgroundPeaks()函数进行选择，该函数根据根据GC含量相似性和样本中的fragment数计算马氏距离然后对peak进行抽样。

1

projHeme5 <- addBgdPeaks(projHeme5)

接下来，就可以使用addDeviatonMatrix()函数根据所有的motif注释计算每个细胞的偏离值。该函数有一个可选参数matrixName，用于定义该偏离值矩阵在Arrow文件里的名字。在下面的例子，函数会以”peakAnnotation”里设置的参数为基础，额外在后面添加字符串”Matrix”，因此下面函数运行结束后会为每个Arrow文件都创建了一个”MotifMatrix”的偏离值矩阵

1
2
3
4
5

projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

我们用getVarDeviations()函数提取这些偏离值矩阵。假如我们需要它返回一个ggplot对象，那么只需要设置plot=TRUE即可，函数会返回一个DataFrame对象。函数运行后，会默认展示该DataFrame对象的前几行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

plotVarDev <- getVarDeviations(projHeme5, name = "MotifMatrix", plot = TRUE)
# DataFrame with 6 rows and 6 columns
# seqnames idx name combinedVars combinedMeans
#
# f388 z 388 GATA2_388 11.9292478607949 -0.034894792575792
# f155 z 155 CEBPA_155 11.8070700579364 -0.174087405321135
# f383 z 383 GATA1_383 11.8045825337775 -0.0378306234562619
# f336 z 336 SPIB_336 11.3432739583017 -0.0819836042460723
# f385 z 385 GATA5_385 10.8828679211543 -0.036867577013264
# f651 z 651 SMARCC1_651 10.2885493109675 -0.131812047523969
# rank
#
# f388 1
# f155 2
# f383 3
# f336 4
# f385 5
# f651 6

从上面DataFrame的输出信息中，你会发现MotifMatrix的seqnames并不是chromosome(染色体名)。通常而言，TileMatrix, PeakMatrix, GeneScoreMatrix，我们都是在seqnames中记录染色体信息。MotifMatrix并没有任何对应的位置信息，而是会在相同的矩阵里记录chromVAR输出的”devations”和”z-scores”信息，即deviations和z。如果你后续想在getMarkerFeatures()这类函数中使用MotifMatrix(Sparse.Assays.Matrix类)的话，那么这些信息就非常重要了。在这类操作中，ArchR会希望你从MotifMatrix提取其中一个seqnames(例如，选择z-scores或deviations执行计算)

我们先绘制这些偏离值

1

plotVarDev

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(plotVarDev, name = "Variable-Motif-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

假如我们想提取部分motif用于下游分析，那该怎么做呢？这就需要用到getFeatures()函数。下面的paste(motifs, collapse="|")语句会以”逻辑或”连接所有motifs里的值，用于选择给定的motif。

1
2
3
4
5
6
7
8

motifs <- c("GATA1", "CEBPA", "EBF1", "IRF4", "TBX21", "PAX5")
markerMotifs <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "MotifMatrix")
markerMotifs
# [1] “z:TBX21_780” “z:PAX5_709” “z:IRF4_632”
# [4] “z:GATA1_383” “z:CEBPA_155” “z:EBF1_67”
# [7] “z:SREBF1_22” “deviations:TBX21_780” “deviations:PAX5_709”
# [10] “deviations:IRF4_632” “deviations:GATA1_383” “deviations:CEBPA_155”
# [13] “deviations:EBF1_67” “deviations:SREBF1_22”

正如之前所提到的，MotifMatrix的seqnames包含z-scores(z:)和deviations(deviations:)。为了只提取对应特征的z-scores， 我们需要用到grep。此外，在之前的选择中由于”EBF1”会匹配到”SREBF1”，而后者并不是我们所需要的，因此我们还需要一步过滤。ArchR提供了%ni表达式，它是R提供的%ni%的反义词，表示反向选择。

1
2
3

markerMotifs <- grep("z:", markerMotifs, value = TRUE)
markerMotifs <- markerMotifs[markerMotifs %ni% "z:SREBF1_22"]
markerMotifs

既然，我们已经有了我们感兴趣的特征，我们可以为每个cluster绘制chromVAR偏离得分。注，我们提供的是之前基因得分分析里计算的推断权重。考虑到scATAC-seq数据的稀疏性，推断权重利用邻近细胞对信号进行平滑处理。

1
2
3
4
5
6

p <- plotGroups(ArchRProj = projHeme5, 
  groupBy = "Clusters2", 
  colorBy = "MotifMatrix", 
  name = markerMotifs,
  imputeWeights = getImputeWeights(projHeme5)
)

我们使用cowplot将不同的moitfs的分布组合在一张图中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

p2 <- lapply(seq_along(p), function(x){
  if(x != 1){
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }else{
    p[[x]] + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6) +
    theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")) +
    theme(
        axis.ticks.y=element_blank(),
        axis.title.y=element_blank()
    ) + ylab("")
  }
})
do.call(cowplot::plot_grid, c(list(nrow = 1, rel_widths = c(2, rep(1, length(p2) - 1))),p2))

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(p, name = "Plot-Groups-Deviations-w-Imputation", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

除了检查z-scores的分布，我们也可以和之前展示基因得分一样将z-scores在UMAP嵌入图中进行展示

1
2
3
4
5
6
7

p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "MotifMatrix", 
    name = sort(markerMotifs), 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)

我们可以使用cowplot将motif UMAP放在一张图上展示

1
2
3
4
5
6
7
8
9
10
11
12

p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

为了比较TF deviation z-scores和根据对应TF基因的基因得分推断的基因表达量，我们可以把这两者画在同一个UMAP中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

markerRNA <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "GeneScoreMatrix")
markerRNA <- markerRNA[markerRNA %ni% c("SREBF1","CEBPA-DT")]
# 获取GeneScoreMatrix
p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "GeneScoreMatrix", 
    name = sort(markerRNA), 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)
# 同时处理的多个图
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
#拼图
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

同样的，我们之前将对应的scRNA-seq数据和scATAC-seq数据进行了关联，我们也可以在UMAP图上绘制每个TF对应的基因表达量.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

markerRNA <- getFeatures(projHeme5, select = paste(motifs, collapse="|"), useMatrix = "GeneIntegrationMatrix")
markerRNA <- markerRNA[markerRNA %ni% c("SREBF1","CEBPA-DT")]
p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "GeneIntegrationMatrix", 
    name = sort(markerRNA), 
    embedding = "UMAP",
    continuousSet = "blueYellow",
    imputeWeights = getImputeWeights(projHeme5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

13.2 ArchR和自定义偏离

在Peak Annotation Enrichment一章中，我们介绍了如何根据任意的基因组区域创建peak注释。这包括 1)ArchR支持的区域集，例如来自于ENCODE人工审核过的TF结合位点和混合ATAC-seq;2) 用户自定义的区域集。如果你没有阅读该章节，我们建议先去阅读一遍，这样能很好的理解peak注释是如何工作的。

这些peak注释能和motif一样用于计算偏离。这里，我们提供了一些例子关于这类分析，由于代码和之前motif分析类似，因此我们不会对每一步的代码做详细的解释。一旦你在Arrow文件中创建了偏离矩阵，其他都是相同的。

13.2.1 Encode TFBS

同样，我们要确保我们已经在ArchRProject中添加了”Encode TFBS”注释矩阵

1
2
3

if("EncodeTFBS" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "EncodeTFBS")
}

接着，我们创建偏离矩阵，以”Encode TFBS”作为peakAnnotation参数的输入

1
2
3
4
5

projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "EncodeTFBS",
  force = TRUE
)

我们就可以绘制排序deviations的点图

1
2

plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "EncodeTFBSMatrix")
plotVarDev

保存为PDF

1

plotPDF(plotVarDev, name = "Variable-EncodeTFBS-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

或者，我们可以提取和特定motif相关的TF结合位点，然后在UMAP嵌入中绘制每个细胞的deviation z-scores。代码和上一节类似。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

tfs <- c("GATA_1", "CEBPB", "EBF1", "IRF4", "TBX21", "PAX5")
markerTFs <- getFeatures(projHeme5, select = paste(tfs, collapse="|"), useMatrix = "EncodeTFBSMatrix")
markerTFs <- sort(grep("z:", markerTFs, value = TRUE))
TFnames <- stringr::str_split(stringr::str_split(markerTFs, pattern = "\\.", simplify=TRUE)[,2], pattern = "-", simplify = TRUE)[,1]
markerTFs <- markerTFs[!duplicated(TFnames)]
p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "EncodeTFBSMatrix", 
    name = markerTFs, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

13.2.2 混池ATAC-seq

类似的，我们可以使用ArchR审核过的混池ATAC-seq peak集来计算motif偏离。下面的代码用来增加ATAC的注释

1
2
3

if("ATAC" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "ATAC")
}

接着创建偏离矩阵

1
2
3
4
5

projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "ATAC",
  force = TRUE
)

然后画图

1
2

plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "ATACMatrix")
plotVarDev

保存为PDF

1

plotPDF(plotVarDev, name = "Variable-ATAC-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

另外，我们还可以在UMAP图上绘制deviation z-scores和每个细胞的peak集

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

ATACPeaks <- c("Heme_HSC", "Heme_LMPP", "Heme_Ery", "Heme_Mono", "Heme_CD4", "Heme_CD8", "Heme_B", "Heme_NK", "IAtlas_DC_Plasmacytoid")
markerATAC <- getFeatures(projHeme5, select = paste(ATACPeaks, collapse="|"), useMatrix = "ATACMatrix")
markerATAC <- sort(grep("z:", markerATAC, value = TRUE))
p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "ATACMatrix", 
    name = markerATAC, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

13.2.3 自定义偏离

除了使用ArchR审核过的区域集，我们也能够提供自己区域集作为peak注释。这些注释和ArchR审核过的注释使用方法相同。

首先，如果你没有在之前的章节中创建”EncodePeaks”, 我们需要先下载一些ENCODE peak，然后调用addPeakAnnotations()函数

1
2
3
4
5
6
7
8
9

EncodePeaks <- c(
  Encode_K562_GATA1 = "https://www.encodeproject.org/files/ENCFF632NQI/@@download/ENCFF632NQI.bed.gz",
  Encode_GM12878_CEBPB = "https://www.encodeproject.org/files/ENCFF761MGJ/@@download/ENCFF761MGJ.bed.gz",
  Encode_K562_Ebf1 = "https://www.encodeproject.org/files/ENCFF868VSY/@@download/ENCFF868VSY.bed.gz",
  Encode_K562_Pax5 = "https://www.encodeproject.org/files/ENCFF339KUO/@@download/ENCFF339KUO.bed.gz"
)
if("ChIP" %ni% names(projHeme5@peakAnnotation)){
    projHeme5 <- addPeakAnnotations(ArchRProj = projHeme5, regions = EncodePeaks, name = "ChIP")
}

接着，我们从这些peak注释中创建偏离值矩阵

1
2
3
4
5

projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "ChIP",
  force = TRUE
)

后续的分析就和之前的一模一样

1
2

plotVarDev <- getVarDeviations(projHeme5, plot = TRUE, name = "ChIPMatrix")
plotVarDev

保存PDF

1

plotPDF(plotVarDev, name = "Variable-ChIP-Deviation-Scores", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

在UMAP上同时展示结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

markerChIP <- getFeatures(projHeme5, useMatrix = "ChIPMatrix")
markerChIP <- sort(grep("z:", markerChIP, value = TRUE))
p <- plotEmbedding(
    ArchRProj = projHeme5, 
    colorBy = "ChIPMatrix", 
    name = markerChIP, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme5)
)
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 2),p2))

第12章 使用ArchR进行motif和特征富集分析

在鉴定到可靠的peak集之后，我们也会想预测有哪些转录因子参与了结合事件(binding events)，从而产生了这些染色质开放位点。在分析标记peak或差异peak时，这能帮助我们更好的理解为什么某组的peak会富集某一类转录因子的结合位点。举个例子，我们想在细胞特异的染色质开放区域中找到定义谱系的关键转录因子。同样，我们也想根据其他已知特征对不同的peak进行富集分析。比如说，我们像知道是否细胞类型A的细胞特异性ATAC-seq peak对于另一组基因组区域（如ChIP-seq peak）也富集。这一章会详细介绍ArchR中的富集分析原理。

12.1 差异peak中的motif富集

继续上一章的差异peak分析，我们可以寻找在不同类型细胞富集的peak中的motif。我们需要先将motif的注释信息加入到我们的ArchRProject中。我们调用addMotifAnnotations()函数分析ArchRProject的peak中是否存在motif。运行结束后会在ArchRProject对象中加入一个新的二值矩阵，用于判断peak是否包括motif。

1

projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif")

接着我们使用上一章差异检验得到的markerTest分析motif的富集情况，这是一个SummarizedExperiment对象。我们用peakAnnoEnrichments()函数分析这些差异开放peak是否富集某一类moitf。可以设置cutOff来过滤peak，例如PDR <= 0.1 & Log2FC >=0.5记录的是”Erythroid”比”Progenitor”更开放的peak。

注: peakAnnoEnrichment()能用于多种差异富集检验，在后续章节还会介绍。

1
2
3
4
5
6

motifsUp <- peakAnnoEnrichment(
    seMarker = markerTest,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

输出的peakAnnoEnrichment()是一个SummarizedExperiment对象，里面存放着多个assays, 记录着超几何检验的富集结果。

1
2
3
4
5
6
7
8
9

motifsUp
# class: SummarizedExperiment
# dim: 870 1
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(870): TFAP2B_1 TFAP2D_2 … TBX18_869 TBX22_870
# rowData names(0):
# colnames(1): Erythroid
# colData names(0):

然后，我们创建一个data.frame对象用于ggplot作图，包括motif名，矫正的p值和显著性排序。

1
2
3

df <- data.frame(TF = rownames(motifsUp), mlog10Padj = assay(motifsUp)[,1])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))

正如我们所预期的那样，”Erythroid”里开放的peak富集的motif主要是GATA转录因子，符合以往研究中”GATA1”在erythroid分化中发挥的作用。

1

head(df)

使用ggplot展示结果，以ggrepel来标识每个TF motif名。

1
2
3
4
5
6
7
8
9
10
11
12
13

ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
  geom_point(size = 1) +
  ggrepel::geom_label_repel(
        data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
        size = 1.5,
        nudge_x = 2,
        color = "black"
  ) + theme_ArchR() + 
  ylab("-log10(P-adj) Motif Enrichment") + 
  xlab("Rank Sorted TFs Enriched") +
  scale_color_gradientn(colors = paletteContinuous(set = "comet"))

ggUp

通过设置Log2FC <= 0.5我们可以挑选出在”Progenitor”里更加开放的peak，然后分析其中富集的motif。

1
2
3
4
5
6
7

motifsDo <- peakAnnoEnrichment(
    seMarker = markerTest,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC <= -0.5"
  )
motifsDo

准备绘图所需数据框

1
2
3

df <- data.frame(TF = rownames(motifsDo), mlog10Padj = assay(motifsDo)[,1])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))

此时，我们会发现在”Progenitor”细胞更加开放的peak中，更多富集RUNX, ELF和CBFB。

1
2
3
4
5
6
7
8

head(df)
# TF mlog10Padj rank
# 326 ELF2_326 88.68056 1
# 733 RUNX1_733 64.00586 2
# 801 CBFB_801 53.55426 3
# 732 RUNX2_732 53.14766 4
# 734 ENSG00000250096_734 53.14766 5
# 336 SPIB_336 52.79666 6

使用ggplot展示结果。

1
2
3
4
5
6
7
8
9
10
11
12
13

ggDo <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
  geom_point(size = 1) +
  ggrepel::geom_label_repel(
        data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
        size = 1.5,
        nudge_x = 2,
        color = "black"
  ) + theme_ArchR() + 
  ylab("-log10(FDR) Motif Enrichment") +
  xlab("Rank Sorted TFs Enriched") +
  scale_color_gradientn(colors = paletteContinuous(set = "comet"))

ggDo

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(ggUp, ggDo, name = "Erythroid-vs-Progenitor-Markers-Motifs-Enriched", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

12.2 标记Peak的motif富集分析

和之前利用差异peak的motif富集分析类似，我们同样能用getMarkerFeatures()分析标记peak里富集的motif。

我们向函数peakAnnotationEnrichment()传入存放标记peak的SummarizedExperiment对象，即markersPeaks

1
2
3
4
5
6

enrichMotifs <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

输出的peakAnnoEnrichment()是一个SummarizedExperiment对象，里面存放着多个assays, 记录着超几何检验的富集结果。

1
2
3
4
5
6
7
8
9

enrichMotifs
# class: SummarizedExperiment
# dim: 870 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(870): TFAP2B_1 TFAP2D_2 … TBX18_869 TBX22_870
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

直接用plotEnrichHeatmap()函数绘制不同细胞组的富集的motif。通过设置参数n限制每个细胞分组中展示的motif。

1

heatmapEM <- plotEnrichHeatmap(enrichMotifs, n = 7, transpose = TRUE)

使用ComplexHeatmap::draw()函数展示结果

1

ComplexHeatmap::draw(heatmapEM, heatmap_legend_side = "bot", annotation_legend_side = "bot")

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(heatmapEM, name = "Motifs-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3 ArchR富集分析

除了分析peak中富集motif, ArchR还能进行个性化的富集分析。为了方便这类数据探索，我们已人工确定了一些特征数据集，它们能比较容易地在我们感兴趣的peak区间进行检验。我们接下来将会逐个介绍这些特征数据集。该分析最初受LOLA启发。

12.3.1 Encode TF 结合位点

ENCODE协会已经将TF结合位点(TFBS)匹配到多种细胞类型和因子中。我们可以利用这些TFBS去更好地理解聚类结果。例如，我们可以根据富集结果去判断未知细胞类型的可能类型。为了能够使用ENCODE TFBS特征集进行分析，我们需要调用addArchRAnnotations()函数，设置collection = "EncodeTFBS". 和使用addPeakAnnotations()类似，这会创建一个二值矩阵，记录我们的标记peak是否和ENCODE TFBS有重叠。

1

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "EncodeTFBS")

我们接着使用peakAnnoEnrichment()函数分析这些 ENCODE TFBS是否在我们的peak中富集。

1
2
3
4
5
6

enrichEncode <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "EncodeTFBS",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

和之前一样，该函数返回一个SummarizedExperiment对象。

1
2
3
4
5
6
7
8
9
10

enrichEncode
# class: SummarizedExperiment
# dim: 689 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(689): 1.CTCF-Dnd41… 2.EZH2_39-Dnd41… …
# 688.CTCF-WERI_Rb_1… 689.CTCF-WI_38…
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

我们可以使用plotEnrichHeatmap函数从富集结果中创建热图。

1

heatmapEncode <- plotEnrichHeatmap(enrichEncode, n = 7, transpose = TRUE)

然后用ComplexHeatmap::draw()绘制热图

1

ComplexHeatmap::draw(heatmapEncode, heatmap_legend_side = "bot", annotation_legend_side = "bot")

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(heatmapEncode, name = "EncodeTFBS-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3.2 混池ATAC-seq

和ENCODE TFBS类似，我们还可以使用混池ATAC-seq实验鉴定的peak，分析两者的重叠情况。通过设置collection="ATAC"来调用混池ATAC-seqpeak数据集。

1

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "ATAC")

接着通过设置peakAnnotation = "ATAC"检验我们的标记peak是否富集了混池ATAC-seq的peak。

1
2
3
4
5
6

enrichATAC <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "ATAC",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

和之前一样，该函数会输出SummarizedExperiment对象，记录着富集结果

1
2
3
4
5
6
7
8
9
10

enrichATAC
# class: SummarizedExperiment
# dim: 96 11
# metadata(0):
# assays(10): mlog10Padj mlog10p … CompareFrequency feature
# rownames(96): Brain_Astrocytes Brain_Excitatory_neurons … Heme_MPP
# Heme_NK
# rowData names(0):
# colnames(11): B CD4.M … PreB Progenitor
# colData names(0):

我们用plotEnrichHeatmap()函数基于SummarizedExperiment绘制富集热图

1

heatmapATAC <- plotEnrichHeatmap(enrichATAC, n = 7, transpose = TRUE)

使用ComplexHeatmap::draw()绘制结果

1

ComplexHeatmap::draw(heatmapATAC, heatmap_legend_side = "bot", annotation_legend_side = "bot")

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(heatmapATAC, name = "ATAC-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.3.3 Codex TFBS

相同类型的分析还能用于 CODEX TFBS，只要设置collection = "Codex"即可

1
2
3
4
5
6
7
8
9

projHeme5 <- addArchRAnnotations(ArchRProj = projHeme5, collection = "Codex")
enrichCodex <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "Codex",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )
heatmapCodex <- plotEnrichHeatmap(enrichCodex, n = 7, transpose = TRUE)
ComplexHeatmap::draw(heatmapCodex, heatmap_legend_side = "bot", annotation_legend_side = "bot")

于是我们就保存图片了

1

plotPDF(heatmapCodex, name = "Codex-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

12.4 自定义富集

除了之前这些经过人工审核的注释数据集，ArchR还能处理用户自定义注释信息来执行富集分析。接下来，我们会介绍如何根据ENCODE ChIP-seq实验来创建自定义的注释信息。

首先，我们先提供后续将被使用并下载的数据集，也可以提供本地文件。

1
2
3
4
5
6

EncodePeaks <- c(
  Encode_K562_GATA1 = "https://www.encodeproject.org/files/ENCFF632NQI/@@download/ENCFF632NQI.bed.gz",
  Encode_GM12878_CEBPB = "https://www.encodeproject.org/files/ENCFF761MGJ/@@download/ENCFF761MGJ.bed.gz",
  Encode_K562_Ebf1 = "https://www.encodeproject.org/files/ENCFF868VSY/@@download/ENCFF868VSY.bed.gz",
  Encode_K562_Pax5 = "https://www.encodeproject.org/files/ENCFF339KUO/@@download/ENCFF339KUO.bed.gz"
)

然后，我们用addPeakAnnotation()函数在ArchRProject函数中增加自定义注释。我们这里将其命名为”ChIP”

1

projHeme5 <- addPeakAnnotations(ArchRProj = projHeme5, regions = EncodePeaks, name = "ChIP")

和之前一样，我们使用peakAnnoEnrichment()函数根据自定义的注释信息执行peak注释富集分析

1
2
3
4
5
6

enrichRegions <- peakAnnoEnrichment(
    seMarker = markersPeaks,
    ArchRProj = projHeme5,
    peakAnnotation = "ChIP",
    cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
  )

并以相同的步骤生成注释热图。

1
2

heatmapRegions <- plotEnrichHeatmap(enrichRegions, n = 7, transpose = TRUE)
ComplexHeatmap::draw(heatmapRegions, heatmap_legend_side = "bot", annotation_legend_side = "bot")

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(heatmapRegions, name = "Regions-Enriched-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

第11章: 使用ArchR鉴定标记Peak

在讨论基因得分(gene score)这一章中，我们已经介绍了标记特征的鉴定。相同的函数getMakerFeatures()也能够用于从ArchRProject任意矩阵中鉴定标记特征。所谓的标记特征指的是相对于其他细胞分组唯一的特征。这些特征能帮助我们理解类群或者细胞类型特异的生物学现象。在这一章中，我们会讨论如何使用该功能鉴定标记peak。

11.1: 使用ArchR鉴定标记Peak

通常而言，我们想知道哪些peak是某个聚类或者某一些聚类所特有的。在ArchR中，这可以通过设置addMarkFeatures()函数的useMatrix="PeakMatrix"来实现（无需监督）。

首先，我们需要再看一眼projHeme5中有哪些细胞类型，以及它们的相对比例

1

table(projHeme5$Cluster2)

现在，让我们调用getMarkerFeatures参数，并设置useMatrix="PeakMatrix". 此外，为了降低不同细胞组之间的数据质量对结果的影响，我们可以设置bias参数，其中bias = c("TSSEnrichment", "log10(nFrags)")就是用来避免TSS富集和每个细胞的fragment数对结果的影响。

1
2
3
4
5
6
7

markersPeaks <- getMarkerFeatures(
    ArchRProj = projHeme5, 
    useMatrix = "PeakMatrix", 
    groupBy = "Clusters2",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)

getMarkerFeatures()函数返回一个SummarizedExperiment对象，该对象包含一些不同的assays

1

markerPeaks

接着，我们可以用getMarkers函数从输出的SummarizedExperiment对象中提取我们感兴趣的部分。默认情况下，它会返回一个包含多个DataFrame的列表，不同的DataFrame表示来自不同的细胞分组。

1
2

markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")
markerList

如果我们对特定的一个细胞分组感兴趣，我们可以用$进行提取。

1

markerList$Erythroid

除了返回一个包含多个DataFrame的列表外，我们还可以用getMarkers()返回一个GRangesList，只要设置returnGR=TRUE即可。

1
2

markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1", returnGR = TRUE)
markerList

这个GRangesList同样可以用$提取特定细胞组的结果，返回的是一个GRanges对象

1

markerList$Erythroid

11.2 在ArchR中绘制Marker Peaks

ArchR提供了许多绘图函数用于getMarkerfeatuers()返回的SummarizedExperiment对象的可视化。

11.2.1 Marker Peak Heatmap

markerHeatmap能以热图的形式展示标记Marker Peak（或者其他getMarkerFeatures()输出的特征）

1
2
3
4
5

heatmapPeaks <- markerHeatmap(
  seMarker = markersPeaks, 
  cutOff = "FDR <= 0.1 & Log2FC >= 0.5",
  transpose = TRUE
)

使用draw函数绘制结果

1

plotPDF(heatmapPeaks, name = "Peak-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

plotFDF()函数能够以可编辑的矢量版本保存图片

1

plotPDF(heatmapPeaks, name = "Peak-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme5, addDOC = FALSE)

11.2.2 Marker Peak MA和火山图

除了绘制热图，我们也可以为每个细胞分组绘制MA或者火山图(volcano)。这些图可以用markerPlot()函数绘制。对于MA图，需要设置参数plotAs="MA". 我们以”Erythroid”细胞分组为例，设置参数name = "Erythroid"

1
2

pma <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "MA")
pma

同样的，只要设置plotAs="Volcano"就可以绘制火山图

1
2

pv <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "Volcano")
pv

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(pma, pv, name = "Erythroid-Markers-MA-Volcano", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

11.2.3 Browser Tracks的Marker Peak

此外，我们在基因组浏览器上检查这些peak区域，只需要为plotBrowserTrack()函数的features参数传入 相应的peak区间。这会额外在我们的ArchR track图的下方以BED形式展示marker peak区域。

1
2
3
4
5
6
7
8

p <- plotBrowserTrack(
    ArchRProj = projHeme5, 
    groupBy = "Clusters2", 
    geneSymbol = c("GATA1"),
    features =  getMarkers(markersPeaks, cutOff = "FDR <= 0.1 & Log2FC >= 1", returnGR = TRUE)["Erythroid"],
    upstream = 50000,
    downstream = 50000
)

我们使用grid::grid.draw()绘制结果

1
2

grid::grid.newpage()
grid::grid.draw(p$GATA1)

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(p, name = "Plot-Tracks-With-Features", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

11.3 组间配对检验

标记特征鉴定是一种特别的差异表达检验。此外，使用相同的getMarkerFeatures()函数也能实现标准化的差异分析。我们只需要设置useGroup为一组细胞，然后设置bgdGroup为另一组细胞即可。这就可以对给定两组进行差异分析。在这些差异分析中，在useGroups比较高的peak的倍数变化值是正数，在bgdGroups比较高的peak则是由负的倍数变化值。

这里，我们对”Erythroid”与”Progenitor”细胞组进行配对检验。

1
2
3
4
5
6
7
8
9

markerTest <- getMarkerFeatures(
  ArchRProj = projHeme5, 
  useMatrix = "PeakMatrix",
  groupBy = "Clusters2",
  testMethod = "wilcoxon",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  useGroups = "Erythroid",
  bgdGroups = "Progenitor"
)

使用markerPlot()函数可以绘制MA或者火山图。MA图需要设置plotAs="MA"

1
2

pma <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "MA")
pma

火山图需要设置plotAs="Volcano"

1
2

pv <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "Volcano")
pv

plotFDF()函数能够以可编辑的矢量版本保存图片。

1

plotPDF(pma, pv, name = "Erythroid-vs-Progenitor-Markers-MA-Volcano", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)

在后续章节中，我们还会进行差异分析，因为会在我们的差异开放的peak中搜索富集的motif。

第10章: 使用ArchR识别Peak

识别peak(Peak calling)是ATAC-seq数据分析的常规操作之一。因为每个细胞的scATAC-seq数据都只有两种状态（开放或不开放），所以我们不能从单个细胞中识别peak。因此，我们在之前的章节中先对细胞进行了分组（细胞聚类）。此外，我们还构建了拟混池重复(pseduo-bulk replicates)，用于评估被识别的peak的可重复性。

10.1 迭代重叠Peak合并流程

我们先介绍2018年的迭代重叠Peak合并流程，然后逐一介绍其他策略（这些策略都存在着一些问题）。

10.1.1 等宽和不定宽peak

我们之所以选择501-bp的等宽peak，一方面是下游处理不再需要根据peak宽度进行标准化从而简化了计算，另一方面是ATAC-seq的大部分peak宽度都低于501-bp。如果使用不定宽peak，后续不同样本的peak合并就会变得特别复杂。而且，我们也没有觉得使用不定宽peak带来的潜在好处值得我们花费那么多精力。更进一步的说，无论你选择哪种，大部分的分析结果也都不会受到影响。

下面，我们会以几个含有不同peak的细胞类型为例，用来讲解不同peak合并方法的差异。

10.1.2: Raw Peak Overlap

raw peak overlap就是分析不同细胞类型的peak是否有重叠，如果有将其合并成单个较大的peak。在如下的图解中，我们会发现尽管细胞类型A和C的前两个peak没有直接相连，但是由于细胞类型B中第一个peak和细胞类型A和C的前两个peak相互重叠，最终使得三个peak合并成一个更大peak。除了这个问题外，如果你想跟踪peak的顶点，你要么记录每一个合并后的新peak的顶点，要么记录每个合并后peak对应的原peak的顶点。

最后，这种类型的peak合并方法通常用bedtools merge命令实现。

10.1.3: Clustered Overlap

clustered overlap方法将peak进行聚类，然后根据规则挑选其中一个胜者。我们可以使用bedtools cluster命令对peak进行聚类，然后挑选每个peak聚类中最显著的peak。根据我们的经验，它会漏掉附近较小的peak，导致总体peak数减少。

10.1.4 Iterative Overlap

iterative overlap尽量避免了以上提到的问题。首先根据显著性对peak进行排序。接着我们保留其中最显著的peak，并移除掉所有于最显著peak发生重叠的peak。然后，在所有余下的peak中，我们会重复上述步骤，直到没有peak为止。这个能够避免10.1.2提到的问题，同时还能使用固定宽度peak。

10.1.5 不同方法对比

我们很容易从最终的peak集中看到以上这些方法的区别。我们认为，iterative overlap策略能够较少的代价得到最好的peak集。

10.1.6 ArchR的工作方式

iterative overlap这种以分级形式对数据进行合并，尽可能保留所有细胞类型特异的peak。

想象一下，你有3种细胞类型，A，B和C，每一种细胞类型多有3个拟混池重复。ArchR使用函数addReproduciblePeakSet()实现iterative overlap的peak合并流程。首先，ArchR单独为每个拟混池重复鉴定peak。然后ArchR同时分析每个细胞类型的所有拟混池重复，进行第一次的iterative overlap过滤。需要注意的是，ArchR使用标准化的peak显著性分值矩阵去比较不同样本的peak的显著性。因为有报道称MACS2的显著性和测序深度成比例关系，所以无法直接比较不同样本的peak显著性。在第一轮的iterative overlap过滤后，ArchR检查所有拟混池重复中每个peak的可重复性，只保留符合参数reproducibility的peak。最后，3种细胞类型(A, B, C)每一个都有单独的合并后的peak集。

接着，我们重复以上步骤，用于合并细胞类型A, B 和C的peak集。在这一步，我们还会根据不同的细胞类型对peak的显著性进行标准化，然后执行iterative overlap过滤。

最终我们会得到单个固定宽度的peak数据集。

10.1.7 如何更改策略

加入你不想用iterative overlap合并策略，那么你可以跳过addReproduciblePeakSet()，或者使用ArchRProj <- addPeakSet()添加自己的peak数据集

10.2 使用w/MACS2鉴定peak

如上所述，我们使用ArchR的addReproduciblePeakSet()得到可重复的peak。默认情况下，ArchR会使用MACS2鉴定peak，此外ArchR也实现它原生的peak鉴定工具，用在MACS2无法安装的情况（例如，我们无法在Windows上安装MACS2），该可选peak鉴定方法会在下一节介绍。

为了使用MACS2鉴定peak，ArchR必须能够找到MACS2的执行路径。ArchR会先在你的PATH环境变量进行查找。如果找不到，ArchR会尝试确认你是否用pip或pip3安装过MACS2. 如果以上都没有陈宫，ArchR就会放弃并输出报错信息。如果你安装了MACS2，但是ArchR找不到它，那么你需要在addReproduciblePeakSet()函数中设置pathToMacs2参数。

1

pathToMacs2 <- findMacs2()

这是成功的信息

1
2

## Searching For MACS2..
## Found with $path!

如下是失败信息

1
2
3
4
5
6

Searching For MACS2..
Not Found in $PATH
Not Found with pip
Not Found with pip3
Error in findMacs2() : 
  Could Not Find Macs2! Please install w/ pip, add to your $PATH, or just supply the macs2 path directly and avoid this function!

当你有MACS2的位置信息后，我们接着就可以构建可重复的合并peak数据集(约5到10分钟)。为了避免细胞过少的拟混池重复的影响，我们可以通过设置peaksPerCell参数来限制每个细胞的peak数上限。这能避免细胞数过少的peak给最终合并的peak集贡献低质量的peak。此外，addReproduciblePeakSet()还有许多其他的参数可以调整，用?addReproduciblePeakSet可以看到更多信息。

每个ArchRProject项目只能包括一个peak集。我们最后将addReproduciblePeakSet()的输出赋值给我们想要的ArchRProject()。如果你想要测试不同的peak集，你必须先保存一份ArchRProject()并拷贝Arrow文件。尽管这会消耗更多的硬盘存储，但这是Arrow文件的会保存peak矩阵信息，因此无法避免。

1
2
3
4
5

projHeme4 <- addReproduciblePeakSet(
    ArchRProj = projHeme4, 
    groupBy = "Clusters2", 
    pathToMacs2 = pathToMacs2
)

我们使用getPeakSet()函数以GRanges对象提取peak集。该peak集对象包括每个peak最初来源信息。然而这并不是意味着该peak只是在那个组中出现，而是意味着那组中的该peak拥有最高的标准化后显著性。

1

getPeakSet(projHeme4)

10.3鉴定Peaks w/TileMatrix

如上所述，ArchR提供了原生的peak鉴定工具。尽管经过测试，该工具效果和MACS2表现差不多，但除非实在是没有办法，否则都不建议使用这个原生方法。

ArchR的原生peak鉴定工具基于500-bp的TileMatrix鉴定peak，可以通过设置addReproduciblePeakSet()函数的peakMethod参数来进行调用。注意，我们没有将输出保存在projHeme4对象中，因为我们这里只是测试而已，所以不打算保留这个peak集。将其保存在ArchRProject对象会覆盖之前已经存放在projHeme4的peak集。

1
2
3
4
5
6

projHemeTmp <- addReproduciblePeakSet(
    ArchRProj = projHeme4, 
    groupBy = "Clusters2",
    peakMethod = "Tiles",
    method = "p"
)

我们同样可以检查这个合并后的peak集

1

getPeakSet(projHemeTmp)

10.3.1 比较两种方法

我们可以根据重叠peak的百分比等指标来比较MACS2和ArchR自带的TileMatrix方法的差异。

1，我们检查MACS2鉴定的peak和TileMatrix鉴定的peak的重叠情况

1
2

length(subsetByOverlaps(getPeakSet(projHeme4), getPeakSet(projHemeTmp))) / length(getPeakSet(projHeme4))
### 0.9627246

2，我们检查TileMatrix鉴定的peak和MACS2鉴定的peak的重叠情况。我们发现这个重叠比例并不高。

1
2

length(subsetByOverlaps(getPeakSet(projHemeTmp), getPeakSet(projHeme4))) / length(getPeakSet(projHemeTmp))
### 0.7533365

如果我们增加peak的边界(从500-bp提高到1000)，MACS2鉴定的peak和TileMatrix鉴定的peak的重叠比例几乎没变化

1
2

length(subsetByOverlaps(resize(getPeakSet(projHeme4), 1000, "center"), getPeakSet(projHemeTmp))) / length(getPeakSet(projHeme4))
### 0.9676687

但是增加了TileMatrix的peak和MACS2找到的peak的重叠比例。

1
2

length(subsetByOverlaps(getPeakSet(projHemeTmp), resize(getPeakSet(projHeme4), 1000, "center"))) / length(getPeakSet(projHemeTmp))
### 0.8287639

来自译者的话: 我们从中可以得出，TileMatrix能鉴定的peak，MACS2也差不多都能鉴定到，但是MACS2鉴定出来的peak，未必TileMatrix也能鉴定。此外，TileMatrix找到的一些peak和MACS2找到的peak也比较近，因此通过延长peak的宽度，就能提高比例。

在10.1这个这个章节中，作者用到了菊花链拓扑(daisy chain)这个名词用来描述不同peak相互连接的情况。考虑到这个菊花链并非是一个常用名词，因此我通过意译来翻译这些名词。

10.4 增加Peak矩阵

我们可以用saveArchRProject()函数保存原来的projHeme4对象。该ArchRProject包含来源于MACS2的合并后peak集。

1

saveArchRProject(ArchRProj = projHeme4, outputDirectory = "Save-ProjHeme4", load = FALSE)

为了方便下游分析，我们可以创建新的ArchRProject命名为projHeme5，并增加新的矩阵，该矩阵记录这新的合并后peak集的insertion计数。

1

projHeme5 <- addPeakMatrix(projHeme4)

于是，projHeme5现在又多了一个新的矩阵，名为”Peakmatrix”，这是和”GeneScoreMatrix”和”TileMatrix”类似的专用矩阵变量名。正如之前所说，每个ArchRProject对象只能有一个peak集和一个PeakMatrix。当然，你可以不同命名构建无限个自定义的特征矩阵，但是PeakMatrix只能有一个，它就是用来保存来自于peak集的insertion计数矩阵。

1
2
3

getAvailableMatrices(projHeme5)
## [1] "GeneIntegrationMatrix" "GeneScoreMatrix" "PeakMatrix"
## [4] "TileMatrix"

MAKER深入篇-如何避免重复运算

通常而言，我们会运行不只一轮的MAKER。如果参考组序列没有变化，那么有一些计算只需要做一次就行了，例如将EST, Repeat和Protein序列比对到参考基因组，得到它们对应的位置。

我们有三种方法可以避免不必要的运算，第一种方法是直接修改配置文件，让MAKER重复利用之前的运行结果；第二种方式是利用之前输出的GFF文件，通过配置”Re-annotation Using MAKER Derived GFF3”里的选项来跳过对应的计算；第三种方法于是利用之前输出的GFF文件，从中提取EST/Repeat/Protein的位置信息保存为GFF文件，通过配置”est_gff”, “protein_gff”, “rm_gff”来避免重新计算位置信息。

后续分析建立MAKER高级篇-SNAP模型训练基础上，也就是通过protein和est序列直接输出基因模型，然后训练出初步的HMM模型

方法1

方法1最为简单，我们只需要修改之前的maker_opts.ctl里的参数，然后重新运行即可。运行时会输出如下的警告信息。

注意: MAKER是通过对比maker_opt.ctl里的配置信息和自己运行时记录的maker_opts.log来判断哪些参数发生了改变。因此，如果SNAP第二次训练生成的文件，要是和上一次命名相同，那么它会认为你这次输入的模型文件和上次相同，就会跳过SNAP预测这一步。

实际运行时，MAKER会跳过BLAST步骤，但是依旧会调用”exonerate”来处理BLAST结果。

方法2

如果你不小心把maker的输出文件删掉了，但是你保留着之前gff3_merge默认参数输出的文件，那么你可以使用该文件来跳过BLAST和Exonerate运算。

Step1: 配置”Re-annotation Using MAKER Derived GFF3”里的参数

1
2
3
4
5
6
7
8
9

#-----Re-annotation Using MAKER Derived GFF3
maker_gff=round1.gff #MAKER derived GFF3 file
est_pass=1 #use ESTs in maker_gff: 1 = yes, 0 = no
altest_pass=1 #use alternate organism ESTs in maker_gff: 1 = yes, 0 = no
protein_pass=1 #use protein alignments in maker_gff: 1 = yes, 0 = no
rm_pass=1 #use repeats in maker_gff: 1 = yes, 0 = no
model_pass=1 #use gene models in maker_gff: 1 = yes, 0 = no
pred_pass=1 #use ab-initio predictions in maker_gff: 1 = yes, 0 = no
other_pass=1 #passthrough anyything else in maker_gff: 1 = yes, 0 = no

此处的round1.gff通过gff3_merge从上一论的maker输出中提取，代码如下

1

gff3_merge -d genome.maker.output/genome_master_datastore_index.log -o round1.gff

Step2: 将”EST Evidence”和”Protein Homology Evidence”里的配置清空，如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

#-----EST Evidence (for best results provide a file for at least one)
est= #set of ESTs or assembled mRNA-seq in fasta format
altest= #EST/cDNA sequence file in fasta format from an alternate organismest_gff= #aligned ESTs or mRNA-seq from an external GFF3 file
altest_gff= #aligned ESTs from a closly relate species in GFF3 format

#-----Protein Homology Evidence (for best results provide a file for at least one)
protein=  #protein sequence file in fasta format (i.e. from mutiple organisms)
protein_gff=  #aligned protein homology evidence from an external GFF3 file

#-----Repeat Masking (leave values blank to skip repeat masking)
model_org= #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff= #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)

Step3: 配置”Gene Prediction”，例如SNAP, 同时将”est2genome”和”protein2genome”设置为0

1
2
3
4
5
6
7
8

#-----Gene Prediction
snaphmm=snap.hmm #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species= #Augustus gene prediction species model
# 略过其他参数
est2genome=0 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=0 #infer predictions from protein homology, 1 = yes, 0 = no
# 略过其他参数

会跳过exonerate步骤，直接从snap预测开始。

方法3

我们还可以通过设置est2gff, protein_gff和rm_gff，来避免重复序列屏蔽和BLAST+Exonerate运算

Step1: 从之前的MAKER输出的GFF文件种提取EST/Protein/Repeat的位置信息

1
2
3
4
5
6

# transcript alignment
awk '{ if ($2 ~ "est") print $0 }' round1.gff > est.gff
# protein alignments
awk '{ if ($2 == "protein2genome") print $0 }' round1.gff > protein2genome.gff
# repeat alignments
awk '{ if ($2 ~ "repeat") print $0 }' round1.gff > repeats.gff

Step2: 修改EST Evidence / rotein Homology Evidence /Repeat Masking里的配置参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

#-----EST Evidence (for best results provide a file for at least one)
est= #set of ESTs or assembled mRNA-seq in fasta format
altest= #EST/cDNA sequence file in fasta format from an alternate organismest_gff=est.gff #aligned ESTs or mRNA-seq from an external GFF3 file
est_gff=./est.gff
altest_gff= #aligned ESTs from a closly relate species in GFF3 format

#-----Protein Homology Evidence (for best results provide a file for at least one)
protein=  #protein sequence file in fasta format (i.e. from mutiple organisms)
protein_gff=protein2genome.gff  #aligned protein homology evidence from an external GFF3 file

#-----Repeat Masking (leave values blank to skip repeat masking)
model_org= #select a model organism for RepBase masking in RepeatMasker
rmlib= #provide an organism specific repeat library in fasta format for RepeatMasker
repeat_protein= #provide a fasta file of transposable element proteins for RepeatRunner
rm_gff=repeats.gff #pre-identified repeat elements from an external GFF3 file
prok_rm=0 #forces MAKER to repeatmask prokaryotes (no reason to change this), 1 = yes, 0 = no
softmask=1 #use soft-masking rather than hard-masking in BLAST (i.e. seg and dust filtering)

Step3: 配置”Gene Prediction”，例如SNAP, 同时将”est2genome”和”protein2genome”设置为0

1
2
3
4
5
6
7
8

#-----Gene Prediction
snaphmm=snap.hmm #SNAP HMM file
gmhmm= #GeneMark HMM file
augustus_species= #Augustus gene prediction species model
# 略过其他参数
est2genome=0 #infer gene predictions directly from ESTs, 1 = yes, 0 = no
protein2genome=0 #infer predictions from protein homology, 1 = yes, 0 = no
# 略过其他参数

同样也会跳过exonerate步骤，直接从snap开始。

结果比较

对于这三种方法，从运行日志中看，三者都会跳过重复序列屏蔽，将EST和蛋白序列回帖到参考基因组的步骤，然而最终预测的基因数却不一致。

分析方法2和方法1的输出GFF文件时，发现方法2输出包括exonerate_protein2genome-gene和exonerate_est2genome-gene。推测其原因在第二种方法的model_pass, pred_passs参数在设置为1时会使用之前est2genome和protein2genome输出的基因模型，而由于模型本身就来自于EST和Protein，就变成自我验证，于是输出结果就变多了。当设置model_pass, pred_passs参数为0时，最终保证方法2和方法1输出结果一致。

之后设置model_pass, pred_passs参数为0，然后比较方法1，方法2和方法3输出的GFF。我发现方法1和方法2的第二列信息完全相同，是blastn, blastx, est2genome, maker, protein2genome, repeatmasker, snap_masked, 而方法3的第二列为est_gff:est2genome, maker, protein_gff:protein2genome, repeat_gff:repeatmasker, snap_masked. 目前只能推测是MAKER对这些证据使用方式不同引起了最终输出结果的差异，但具体的原理我没有分析清楚，不过不妨碍使用。

最后，三种方法使用优先级分别是方法1 > 方法2 > 方法3，其中方法2要注意设置model_pass, pred_passs的设置。

准备阶段

训练SNAP模型，需要准备三个文件，分别是参考基因组序列，组装的转录本序列和同源蛋白序列。

对于参考基因组序列，我们要保证N50需要超过预期基因长度的中位数，否则注释效果不好。不过目前的基因组在三代测序的加持下，基本上都不是问题。

对于同源蛋白， 建议只用UniProt/Sprot的人工检查过的高质量蛋白序列，而不是盲目参考文献，使用近源物种的所有蛋白。除非你的近源物种都是模式物种，蛋白序列可靠性高，否则用错误的输入进行训练，数据越多反而错的越多。

我们可以在ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/选择合适的uniprot_sprot数据, 然后将其输出为fasta格式。以植物为例

1
2
3
4

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz
zcat uniprot_sprot_plants.dat.gz |\
    awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' |\
    sed 's/;$//'> protein.fa

对于转录本，我们通常会测一些转录组数据，有三种策略可以得到转录本。（这里暂时不考虑三代全长转录本）

1. Trinity重头组装转录本
2. 使用STAR + Trinity 获取转录本
3. 使用STAR + StringTie + gffread 获取转录本

对于这三种策略，不推荐策略一，因为在有参考基因组的情况下，策略二不但计算效率高，而且能避免组装错误（多倍体等位基因之间相似度高）。对于策略二和策略三，我会推荐策略三。因为对于靠的比较近的基因，Trinity很可能会把这两个基因装成一个。

并且利用该转录本作为输入训练SNAP模型，之后以SNAP模型作为输入，将转录组和同源蛋白作为证据而不是直接用作模型，我们再检查maker的结果, 也会发现使用StringTie进行组装的结果才是对的。

因此使用策略二不但计算量大，而且有些情况下还会导致过近的转录本错误融合，反而影响了最终效果，因此我最终推荐策略三。当然，这是我定性通过IGV浏览结果得出的结论，样本小，结论未必可靠，仅供参考。

训练阶段

假设我们准备的三个文件分别命名为, genome.fa, Trinity-GG.fasta 和 protein.fa

接着使用maker -CTL新建配置文件, 设置如下选项

1
2
3
4
5

genome=genome.fa
est=Trinity-GG.fasta
protein=protein.fa
est2genome=1
protein2genome=1

然后使用mpiexec -n 线程数 maker &> run.log运行程序。

处理结果后，我们新建一个snap目录训练模型

1
2

mkdir snap && cd snap
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log

使用makerzff构建输入文件

1

maker2zff -c 0.8 -e 0.8 -o 0.8 -x 0.2 genome.all.gff

maker2zff会根据AED(-x)和QI值进行过滤，其中QI值一共有9项，每一项的含义如下

1. Length of the 5’ UTR
2. Fraction of splice sites confirmed by an EST alignment (-c)
3. Fraction of exons that overlap an EST alignmetn(-e)
4. Fraction of exons that overlap EST or Protein alignments(-o)
5. Fraction of splice site confrimed by a ab-initio prediction(-a)
6. Fraction of exons that overlap a ab-initio prediction(-t)
7. Number of exons in the mRNA
8. length of the 3’ UTR
9. Length of the protein sequence produced by the mRNA (-l)

如果QI值第二项为-1，表示没有支持该剪切位点的EST证据.

接着构建模型

1
2
3
4

fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap.hmm

修改配置，然后重新运行。MAKER会自动处理冲突的部分，避免重复序列屏蔽等的一些重复计算。

1
2
3
4
5
6

genome=genome.fa
est=Trinity-GG.fasta
protein=protein.fa
snap=snap.hmm
est2genome=0
protein2genome=0

根据输出结果再一次训练模型

1
2
3
4
5
6

mkdir snap2 && cd snap2
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log
fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap2.hmm

通常迭代2-3次就够了，毕竟我们可能还会训练AUGUSTUS和GeneMark模型，通过比较多个模型来得到最终结果。

参考资料: http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018#Training_ab_initio_Gene_Predictors

本文翻译自http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/The_MAKER_control_files_explained

MAKER会生成三个配置文件，如下

• maker_opts.ctl: 控制MAKER分析行为的主要配置文件
• maker_bopts.ctl: BLAST和Exnerate的过滤阈值
• maker_exe.ctl: MAKER运行过程中所依赖软件的路径

maker_opts.ctl

基因组

用于设置被注释的基因组序列的位置和物种类型，包括genome和organisam_type两项

• genome: FASTA序列路径
• organism_type: eukaryotic,prokaryotic二选一

需要注意的是，基因组序列的N50需要超过预期基因长度的中位数，否则注释效果不好。另外最好保证基因组序列只包括A,T,C,G,N, 对于其他类型兼并碱基可以都改成N.

使用MAKER得到GFF3进行重注释

这一项基本上我们用不上，它是在当你把MAKER的中间输出文件都删除了，仅保留了输出的GFF3文件时，你可以用之前相同的输入设置重新运行流程得到相同的输出。

1
2
3
4
5
6
7
8
9

#-----Re-annotation Using MAKER Derived GFF3
maker_gff= #MAKER derived GFF3 file
est_pass=0 #use ESTs in maker_gff: 1 = yes, 0 = no
altest_pass=0 #use alternate organism ESTs in maker_gff: 1 = yes, 0 = no
protein_pass=0 #use protein alignments in maker_gff: 1 = yes, 0 = no
rm_pass=0 #use repeats in maker_gff: 1 = yes, 0 = no
model_pass=0 #use gene models in maker_gff: 1 = yes, 0 = no
pred_pass=0 #use ab-initio predictions in maker_gff: 1 = yes, 0 = no
other_pass=0 #passthrough anyything else in maker_gff: 1 = yes, 0 = no