BLAT(BLAST-Like Alignment Tool)是Jim Kent开发的一个类Blast的比对工具,有以下几个优点:1.比对速度快;2.可以把小的比对区域(例如exon)连接成大的比对区域。适合做转录本与基因组的比较以及相近物种的基因同源性分析。
blat database query [-ooc=11.ooc] -minIdentity=85 output.psl
##-minIdentity 指定sequence identity阈值
注意:blat不支持过大的基因组,会报“blat: genoFind.c:1806: clumpHits: Assertion `(hit->tStart >> bucketShift) < bucketCount' failed.”错误,可以把基因组拆开比对再合并比对结果。
比对后的结果文件如下:
各列含义如下:
倒数第4列是比对包含的不连续的block的数目,并不一定都是外显子,也可能是由于其他插入/缺失断开的block。
比对完成后可以利用python的pandas包提取最佳比对结果,代码如下:
import pandas as pd
import sys
infile, outfile = sys.argv[1], sys.argv[2]
inp = pd.read_table(infile, skiprows=5, header= None) # skiprows可用于跳过文件前几个注释行。
inp.columns=["matches", "misMatches", "repMatches", "nCount", "qNumInsert", "qBaseInsert",
"tNumInsert", "tBaseInsert", "strand", "qName", "qSize", "qStart", "qEnd", "tName", "tSize", "tStart", "tEnd",
"blockCount", "blockSize", "qStarts", "tStarts"]
output1 = inp.sort_values(by=['qName']) ##按照query name排序
output2 = output1.sort_values(by=['matches'], ascending= False).drop_duplicates(subset='qName') #对于每个query序列,按照 matches降序排序( ascending= False)),排在第一位的是最佳比对结果; drop_duplicates去除其他比对结果。
output2.to_csv(outfile, sep="\t", index=None)
pandas同样可以用于提取其他比对软件的比对结果,其中minimap2比对结果每行列数不一致,需修改几个参数:
import pandas as pd
import sys
infile, outfile = sys.argv[1], sys.argv[2]
inp = pd.read_csv(infile, sep="\t", usecols=[0,1,2,3,4,5,6,7,8,9,10,11] , skiprows=0, header=None) ##使用usecols提取想要的列,这里我保留了1-12列。
inp.columns=["query", "query_length", "q_start", "q_end", "strand", "target", "target_length", "t_start",
"t_end", "matches", "match_block", "mapping quality"]
output1 = inp.sort_values(by=['query'])
output2 = output1.sort_values(by=['matches'], ascending= False).drop_duplicates(subset='query')
output2.to_csv(outfile, sep="\t", index=None)
参考:
Kent, W. J. (2002). BLAT—the BLAST-like alignment tool. Genome research, 12(4), 656-664.
blat User Guide: https://genome.ucsc.edu/goldenPath/help/blatSpec.html#pslSortUsage
psl文件格式:https://www.ensembl.org/info/website/upload/psl.html
blat使用:https://www.cnblogs.com/adawong/articles/7460300.html
python脚本参考:https://blog.csdn.net/weixin_40099163/article/details/83215747