如何在IGV上使用BLAT搜索非模式物种

IGV提供了BLAT，用于进行序列搜索，但可惜我一直用不上，因为它默认是调用了UCSC的CGI工具，将我们的输入序列发送到https://genome.ucsc.edu/cgi-bin/hgBlat处理，处理后返回JSON文件用于展示。因此，除非我们自己搭建一个UCSC类似的网站，否则，无法用到IGV的这个功能。

我觉得肯定不只是我一个人有这个问题，所以我就去谷歌上用关键词 “custom genome BLAT IGV” 进行检索, 果然发现有很多人都有类似的需求，我找到最早一条是2016年，但是距今6年了，仍旧没有一个直观有条理的文章介绍应该怎么做。于是就有了这篇文章。

之所以没有一篇合适的教程，是因为很多人（包括我在内）都冲着自建UCSC的网页服务，然而这条路非常的艰难。但是经过我的检索和测试后，我发现，我们实际上只需要做到向服务器发送一条请求，然后得到包含所需信息的JSON即可。

本篇教程要求

Windows10的WSL或者Linux虚拟机或者有管理员权限的Linux服务器(ubuntu系统)
IGV ≥ 2.10.3

配置BLAT

我们可以用其他工具代替BLAT，只要保证返回的输出结果是PSL格式或者能转成PSL格式即可。

首先，安装BLAT

wget https://users.soe.ucsc.edu/~kent/src/blatSrc35.zip
unzip blatSrc35.zipexport MACHTYPE=x86_64
mkdir -p ~/bin/$MACHTYPE# 软件安装在$HOME/bin/x86_64下
make -j

接着，将fa转成2bit格式，方便blat调用

# 家目录
mkdir genomes
cd genoems
mkdir twoBit
# 设置参考基因组组, 例如Athaliana
ref=Athaliana.fa
$HOME/bin/x86_64/faToTwoBit $ref twoBit/Athaliana.2bit

配置Apache CGI

在安装了blat之后，下一个问题就是如何让服务器调用。这里我使用的是Apache的CGI服务。CGI全称是Common Gateway Interface，公共网关接口，根据维基百科的定义，CGI就是用于让服务器执行外部程序，特别是处理用户的请求

In computing, Common Gateway Interface (CGI) is an interface specification that enables web servers to execute an external program, typically to process user requests.[1]

首先，我们需要安装Apache。Ubuntu就是运行如下的命令

sudo apt install apache2

默认apache的配置文件在 /etc/apache2下，没有启用CGI模块，因此我们需要先启动这个模块

#该目录记录的是哪些模块被启用了
cd /etc/apache2/mods-enabled
# 添加cgi.load
sudo ln -s ../mods-available/cgi.load .

重新加载即可

sudo service apache2 reload

为了测试CGI服务器是否打开，我们创建一个 hw.sh 进行测试，将其复制到 /usr/lib/cgi-bin下

#!/bin/bash
printf "Content-type: text/html\n\n"
printf "Hello World!\n"

通过curl向我们的网页服务器发送一个请求，会返回 Hello World!.

# 启用cgi-bin
$ curl http://127.0.0.1/cgi-bin/hw.sh
Hello World!

在确认CGI服务器启动之后，那么就可以将我用Python编写用于处理BLAT请求的CGI脚本，我将其命名为 myBlat ，复制到/usr/lib/cgi-bin下。

该代码比较简陋，没有考虑多个参考基因组的情况，也没有各种报错检测，但至少能运行了

#!/usr/bin/env python3
import os
import json
import cgitb
import cgi
import tempfile
import subprocesscgitb.enable(display=0, logdir="/tmp/logdir")form = cgi.FieldStorage()
seq = form.getvalue("userSeq")
db = form.getvalue("db")# 需要改成你实际的blat路径
blat = "/home/xzg/bin/x86_64/blat"
# 需要改成你实际的参考位置
ref = "/home/xzg/genomes/twoBit/Athaliana/Athaliana.2bit"fd1, fa_file = tempfile.mkstemp(suffix=".fa")
fd2, psl_file = tempfile.mkstemp(suffix=".psl")with open(fa_file,"w") as f:f.write(f">query\n{seq}\n")process = subprocess.Popen([blat, ref, fa_file, psl_file],stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr = process.communicate()psl_ret = [line for line in open(psl_file, "r") ]
psl_rec = [line.strip().split('\t') for line in psl_ret[5:] ]print("Context-type: application/json\n")data =  {"track": "blat","genome": "ath","fields": ["matches", "misMatches", "repMatches", "nCount", "qNumInsert", "qBaseInsert", \"tNumInsert", "tBaseInsert", "strand", "qName", "qSize", "qStart", "qEnd", "tName", \"tSize", "tEnd", "blockCount", "blockSizes", "qStarts", "tStarts"],"blat" : psl_rec
}print(json.dumps(data))

配置IGV

最后，我们回到IGV，选择View里的Preferences

然后选择Adanved, 修改其中的Blat url

http://127.0.0.1/cgi-bin/myBlat?userSeq=$SEQUENCE&type=DNA&db=$DB&output=json

需要注意的是 127.0.0.1 , 因为我用的是WSL，所以网页服务架设在localhost。如果使用的是虚拟机，或者用的是内网服务器，那么该地址，应该是虚拟机或者内网服务器的IP地址。

此外，我们可以注意到这里面有两个$开头的变量，SEQUENCE和DB，分别对应我们输入的序列，以及我们当前参考基因组的名字，Athaliana.fa.

BLAT运行后返回如下界面，点击对应的行，就可以跳转到对应联配中。