WGDI:blockinfo输出结果解读

blockinfo模块输出文件以csv格式进行存放，共23列，可以用EXCEL直接打开。

其中16列非常容易裂解，描述如下

1. id 即共线性的结果的唯一标识

2. chr1 , start1 , end1 即参考基因组（点图的左边）的共线性范围（对应GFF1的位置）

3. chr2 , start2 , end2 即参考基因组（点图的上边）的共线性范围（对应GFF2的位置)

4. pvalue 即共线性结果评估，常常认为小于0.01的更合理些

5. length 即共线性片段中基因对数目

6. ks_median 即共线性片段上所有基因对 ks 的中位数（主要用来评判ks分布的）

7. ks_average 即共线性片段上所有基因对 ks 的平均值

8. block1 , block2 分别为共线性片段上基因 order 的位置。

9. ks 共线性片段上所有基因对的 ks 值

10. density1 , density2 共线性片段的基因分布密集程度。值越小表示稀疏。

最后两列， class1 和 class2 会在 alignment 模块中用到，对应的是两个block分组，默认值是0表示两个block是同一组。这两列后期需要自己根据覆盖率，染色体核型等多个方面进行确定。举个例子，我们可以根据 homo1 的取值范围对class1进行赋值，例如-1~-0.5 是 1，-0.5 ~ 0.5 是2，0.5~1是3，最后在alignment中会就会用三种颜色来展示，例如下图的1,2,3分别对应red,blue,green.

关于class如何应用在自己的实际项目中，可以加入QQ群 966612552 讨论

中间的 homo1 , homo2 , homo3 , homo4 , homo5 并非那么直观，先说结论：

• 这里的homoN(N=1,2,3,4,5) 表示一个基因有N个最佳匹配时的取值

• N由mutiple参数确定，对应点阵图(dotplot)中的红点

• multiple的取值一般取1即可，表示最近一次的WGD可能是一次二倍化事件，因此每个基因只会有一个最佳匹配。如果设置为2，可能是一次3倍化，每个基因由两个最佳匹配。当然实际情况可能会更加复杂，比如说异源四倍体，或者异源六倍体，或者没有多倍化只是小规模的基因复制(small-scale gene duplication) 等情况，也会影响multiple的设置。

• homoN会在后面过滤共线性区块时用到，一般最近的WGD事件所产生的共线性区块会比较接近1，而古老的WGD产生的共线性区块则接近-1.

接着，我们将根据源代码 blast_homo和blast_position 来说明结算过程。

首先需要用到blast_homo函数，用来输出每个基因对在不同最佳匹配情况下的取值（-1,0,1)。

    def blast_homo(self, blast, gff1, gff2, repeat_number):
        index = [group.sort_values(by=11, ascending=False)[:repeat_number].index.tolist()
                 for name, group in blast.groupby([0])]
        blast = blast.loc[np.concatenate(
            np.array([k[:repeat_number] for k in index], dtype=object)), [0, 1]]
        blast = blast.assign(homo1=np.nan, homo2=np.nan,
                             homo3=np.nan, homo4=np.nan, homo5=np.nan)
        for i in range(1, 6):
            bluenum = i+5
            redindex = np.concatenate(
                np.array([k[:i] for k in index], dtype=object))
            blueindex = np.concatenate(
                np.array([k[i:bluenum] for k in index], dtype=object))
            grayindex = np.concatenate(
                np.array([k[bluenum:repeat_number] for k in index], dtype=object))
            blast.loc[redindex, 'homo'+str(i)] = 1
            blast.loc[blueindex, 'homo'+str(i)] = 0
            blast.loc[grayindex, 'homo'+str(i)] = -1
        return blast

for循环前的代码作用是提取每个基因BLAST后的前N个最佳结果。循环的作用基因对进行赋值，主要规则是基因对如果在点图中为 「红色」 ，赋值为1， 「蓝色」 赋值为0， 「灰色」 赋值为-1。