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.zip
export 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 subprocess
cgitb.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运行后返回如下界面，点击对应的行，就可以跳转到对应联配中。
BLAT返回结果
网友评论