基因组注释（1）——串联重复序列注释

本系列笔记开始记录如何对组装的植物基因组进行注释。前面通过一系列组装过程，我们拿到了组装好的基因组草图，而这个基因组草图只是研究的开始，我们关注的是基因组中有哪些我们感兴趣的功能基因或者结构基因，以及怎么用这些基因阐述生物学问题等等，这个时候一个高准确度的基因组注释结果就非常重要了。

基因组注释可以分为两部分：基因组的结构注释（重复序列识别、非编码基因预测、编码基因预测）和基因功能注释，结构注释是功能注释的基础。

这里先从结构注释中的重复序列注释开始。我们知道植物基因组多倍化频繁，且基因组中存在大量的重复序列（有的植物基因组中重复序列甚至能达到80%），这些重复序列控制植物表型调控中有非常重要的作用。基因组中的重复序列可以分为以下几种：

1. 串联重复序列注释

串联重复（Tandem Repeat, TR）指DNA中的一个或多个核苷酸前后相连接的重复。串联重复又分为卫星DNA(Satellite DNA)、小卫星（Minisatellite）、微卫星（Microsatellite）。微卫星在植物中一般称为SSR（Simple Sequence Repeats）SSR在植物基因组常被用做遗传标记使用。下面我用两款软件跑下串联重复序列注释。

1.1 GMATA

这个软件主要用来搜索重复单元较短的简单重复序列，也就是微卫星SSR序列。这软件运行速度比较快，而且可以同时设计SSR引物，还可以预测elect-PCR结果，或者将预测结果显示在基因组浏览器上，可以在github上找到项目地址：XuewenWangUGA/GMATA: software GMATA (github.com)

需要注意如果是在linux系统上跑一键流程的话，需要单独安装primer3和e-pcr（可以直接用conda安装），分别是设计SSR引物和模拟PCR的时候需要调用。如果没有这方面需要，可以在设置文件default_cfg.txt修改为不启用后面的模块。我这里统一用linux系统演示，只演示命令行操作，这个软件在windows上运行有UI界面，还是比较直观的。

#!/bin/bash
#SBATCH -n 50
#SBATCH -t 7200

perl gmata.pl -c default_cfg.txt -i /public/home/wlxie/NextPolish/luobuma_rundir/genome.nextpolish.fasta

我这里直接用了一键流程，修改默认的设置文件中三个模块[set]:doprimer_smt、[set]:elctPCR和[set]:mk2gff3的ModulRun = N，虽然可以批量设计引物，但是我这里用不着…..

预测的SSR结果在原fasta文件路径下，以.ssr和.ssr.sat2为后缀：

在sat结果文件中，最终结果以4个表格的方式呈现，分别统计motif k-mer、motif和成对的motif信息以及最后每个contig的SSR统计信息。以上是其中两个表格。

1.2 TRF

这个软件和上面软件类似，可以统计整个基因组上的串联重复序列，在上面一个软件输出结果上稍微有些不同。

#!/bin/bash
#SBATCH -n 50
#SBATCH -t 7200

trf /public/home/wlxie/NextPolish/luobuma_rundir/genome.nextpolish.fasta 2 7 7 80 10 50 500 -f -d -m -r -h

说明一下这个软件使用过程中传参的一堆数字代表什么：

Match  = matching weight		# 匹配上的权重，缺省值2
Mismatch  = mismatching penalty	# 未匹配上的权重，缺省值7
Delta = indel penalty		# 插入罚分，缺省值7
PM = match probability (whole number)		# 比对上的概率，可选值为80和75
PI = indel probability (whole number)		# 插入的概率，可选值为10和20
Minscore = minimum alignment score to report		# 串联重复序列的比对必须达到或超过要报告的比对分数
MaxPeriod = maximum period size to report		# 最大重复单元的bp数，不指定的话从1到2000

其他主要可选参数（列一部分）:
-m 输出屏蔽重复序列后的基因组
-f 记录每个重复序列侧翼的500个核苷酸，主要用于PCR引物设计
-d 生成屏蔽数据文件，与汇总表有相同的信息，不包含标签，主要方便做其他处理
-h 不生成html结果文件（contig数量多的话建议使用，否则有大量的文件生成）

运行结束后可以生成.mask后缀的屏蔽后的序列文件，还有一个.dat后缀的结果文件，包含了重复序列的详细信息。

主要讲解下这个dat文件，先less -S打开看看：

上面记录的参数和其他信息就不说了，主要是底下的数据，每一行是一个重复序列的信息，每行分为15列：

第1列和第2列是预测到的重复序列的起始和结束位点；

第3列是重复单元的长度；

第4列是重复单元的拷贝数，不一定是整数，因为可能存在插入缺失；

第5列是一致性序列的长度；

第6列是匹配的百分率；

第7列是插入缺失的百分率；

第8列是TRF软件给的分值，越高越可靠；

第9-12列分别为ACGT碱基的个数；

第13列表示比对的熵值；

第14列是一致性序列的具体碱基排列；

第15列是整个重复序列的具体碱基排列顺序。

对于结果文件的处理有两种方法，一种是将.dat后缀的结果文件转换成标准的.gff3文件格式，然后用bedtools提取trf特征。转化gff的的方法github上有不少开源的项目，这里推荐一个Adamtaranto/TRF2GFF: Convert Tandem Repeat Finder dat file output into gff3 format (github.com)

还有一种方法就是自己写个脚本，可以看到同一位点处可能有多条预测的串联重复序列，也就是说这些串联重复序列之间可能存在交叠，思路是将同一位点预测的重复序列只保留最短的一条（起始位点相同保留前一条，结束位点相同保留后一条），然后统计第3列重复序列k-mer数量和类型，根据第1列和第2列计算长度，统计总长度和占比即可。

'''
自编TFR dat格式结果文件统计脚本
2023.3.12
'''
from collections import Counter

# 数据过滤
loci_start = []
loci_finish = []
total_line = []
with open('./genome.baima.fasta.2.7.7.80.10.50.500.dat', 'r') as input:
    for i in input.readlines()[15:]:
        if i.find('Sequence') != -1 or i.find('Parameters') != -1:
            continue
        else:
            lst = i.strip().split(' ')
            if len(lst) < 15:
                continue
            if len(loci_start) < 1 and len(loci_finish) < 1:    # 处理列表为空的情况
                loci_start.append(lst[0])
                loci_finish.append(lst[1])
                total_line.append(lst)
            else:
                if lst[0] != loci_start[-1] and lst[1] != loci_finish[-1]:  # 开始结束位点都不同，则记录数据
                    loci_start.append(lst[0])
                    loci_finish.append(lst[1])
                    total_line.append(lst)
                elif lst[0] == loci_start[-1]:        # 开始位点相同，跳过
                    continue
                elif lst[1] == loci_finish[-1]:       # 结束位点相同，删除列表最后一个元素并加入新元素
                    del loci_start[-1]
                    del loci_finish[-1]
                    del total_line[-1]
                    loci_start.append(lst[0])
                    loci_finish.append(lst[1])
                    total_line.append(lst)

# 提取motif长度和重复序列长度
motif_lst = []
leng_lst = []
for i in total_line:
    motif_lst.append(i[2])
    leng = int(i[1]) - int(i[0]) + 1
    leng_lst.append(leng)

# 统计相同motif的总长度和所有重复序列总长度
total_leng = {}
motif_sum = 0
for i in range(len(motif_lst)):
    item = motif_lst[i]
    motif_sum += leng_lst[i]
    if item in total_leng:
        total_leng[item] += leng_lst[i]
    else:
        total_leng[item] = leng_lst[i]

# 统计motif-mer数量，总数，占比
count_motif = Counter(motif_lst)
count_lst = list(count_motif.items())
count_lst.sort(key = lambda x : x[1], reverse = True)
lst_ = []
hit_num = 0
for i in count_lst:
    hit_num += i[1] 
for i in count_lst:
    ls = i[0]
    lst1 = list(i)
    if ls in total_leng:
        lst1.append(total_leng[ls])
    precentage = '%.2f%%'%(100 * i[1] / hit_num)
    lst1.append(precentage)
    lst_.append(lst1)

# 统计结果过滤（取前二十）
lst_filted = []
hit_ = 0
motif_ = 0
for i in range(1, 20):
    lst_filted.append(lst_[i])
    hit_ += lst_[i][1]
    motif_ += lst_[i][2]
lst_filted.append(['Other', int(hit_num - hit_), int(motif_sum - motif_), '%.2f%%'%(100 - 100 * hit_ / hit_num)])
lst_filted.append(['Total', int(hit_num), int(motif_sum), '100%'])

# 输出结果
with open('./stastics.xls', 'w') as output:
    output.write('Motif(-mer)\tNumber\tLength(bp)\tPrecentage\n')
    for i in lst_filted:
        motif = i[0]
        number = i[1]
        length = i[2]
        pre = i[3]
        output.write(motif + '\t' + str(number) + '\t' + str(length) + '\t' + str(pre) + '\n')