python脚本根据抗病基因ID和所有基因的bed文件鉴定抗病基因簇

植物里的抗病基因更倾向于成簇存在，分析抗病基因家族通常也会分析成簇存在的或者单个存在的抗病基因的比例，之前想自己写脚本统计这个数据，但是怎么写代码一直没有想明白，最近看论文

Variation in abundance of predicted resistance genes in the Brassica oleracea pangenome

这个论文里提供了一个python脚本

脚本的链接 https://github.com/AppliedBioinformatics/B_oleracea_R_genes_supplementary/blob/master/makeRGeneClusterAnalysis.py

首先是使用 RGAugury 这个流程鉴定抗病基因类似物，获得抗病基因的id列表，然后根据基因组的gff格式注释文件可以获得所有基因的bed文件。有了这两个文件就可以获得最终的结果。

这篇论文里基因簇的定义：Resistance gene candidates were merged into RGA-gene-rich clusters if there was at least one other resistance gene within 10 upstream or 10 downstream genes using a Python 3 script

某个抗病基因的上游或者下游10个基因如果存在其他抗病基因，那么就是一个抗病基因簇，这个定义也不是固定的，不同论文里定义基因簇的方法也不太一样

这个python脚本里面获取某个基因上下游的基因用到的是通过python的os模块调用grep命令，windows下好像没有这个命令，这个脚本应该是只能在linux系统下用，不确定mac是否能用

所有基因的bed文件要根据位置的从大到小的顺序排好

这个脚本里定义的第一个函数还是没有看懂是什么意思

def merge(lsts):
  # from https://stackoverflow.com/a/9112588
  sets = [set(lst) for lst in lsts if lst]
  merged = 1
  while merged:
    merged = 0
    results = []
    while sets:
      common, rest = sets[0], sets[1:]
      sets = []
      for x in rest:
        if x.isdisjoint(common):
          sets.append(x)
        else:
          merged = 1
          common |= x
      results.append(common)
    sets = results
  return sets

这段代码里有一个符号 |= 查了一下暂时也没看懂是什么意思

目前的状态是能够简单修改脚本，换成自己的数据也能跑

一个简单的小例子

python makeRGeneClusterAnalysis.py RGA.lst gene.bed

RGA.lst 是抗病基因的id列表，一行一个 gene.bed文件是所有基因的bed文件 （这两个数据都是我自己随便构造的） 运行输出