注释得到的基因集中，可能某些基因存在被转座子插入的情况，该基因会在后续功能注释的时候被注释上，但实际在基因组中该基因可能已经被插入失活。因此在基因组的功能注释前，需要用检测转座子软件（如TransposonPSI、TEsorter等）将含有转座子的基因找出并去除。
这里记录下TEsorter软件筛选预测基因的方法。
1. TEsorter安装
TEsorter原本是用于调用LTR_retriever鉴别长末端重复序列反转座子（LTR-RTs），也可以用于其他类型TE的鉴别，其鉴定原理为将待测序列与数据库REXdb（整合viridiplantae_v3.0 + metazoa_v3）的TE序列进行比对。
也可以使用GyDB数据库进行比对，官网上有具体的参数用法。
zhangrengang/TEsorter: TEsorter: an accurate and fast method to classify LTR-retrotransposons in plant genomes (github.com)
TEsorter提供conda安装，但是我没有安装成功，这里还是新建conda环境后手动安装各种依赖：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| conda create -n "TEs"
conda activate TEs
conda install python==3.11 # 官方要求python版本高于3，否则运行会报错
conda install biopython
conda install xopen
conda install hmmer
conda install blast
git clone https://ghproxy.com/https://github.com/zhangrengang/TEsorter.git # 从github镜像网站下载
cd TEsorter
python setup.py install
TEsorter TEsorter/test/rice6.9.5.liban # 测试
|
2. TEsorter运行
这里有一个问题，Braker预测基因有gtf和gff3两种格式的输出结果，但是两者的行数不一样。
行数不一致是因为BRAKER v3.0.2存在bug，GeneMark.hmm3预测结果在转gtf到gff的过程中丢失，已于BRAKER v3.0.3版本修复（来自半年后的注释）。
方便起见我这里处理了gtf结果文件。
TEsorter软件可以输入基因序列或者蛋白序列，以基因序列为例，简单编写脚本如下：
1
2
3
4
5
| #!/bin/bash
#SBATCH -n 8
#SBATCH -t 7200
TEsorter /public/home/wlxie/baima_pre_mydb/braker.codingseq -eval 1e-6 -p 8
|
大约十几分钟运行完毕。
3. 结果文件处理
结果文件如下：
.tsv后缀的文件中以列表形式列出了所有预测的TE类型，一个基因可能有多种类型的TE插入，因此需要处理结果文件，统计含有TE的基因，并在gtf结果文件中将该基因去除。
1
2
3
4
5
6
7
8
| # 筛选含有TE的基因
grep -v "^#" braker.codingseq.rexdb.cls.tsv | cut -f1 | sort | uniq | cut -f1 -d "_" | sort | uniq > TE-genes.txt
# 去除含有TE的基因序列
grep -Fvf /public/home/wlxie/biosoft/TEsorter/baima_mydb/TE-genes.txt braker.gtf | awk '$3 ~ /gene/' > baima_gene_only.gtf
# 去除含有TE的转录本序列
grep -Fvf /public/home/wlxie/biosoft/TEsorter/baima_mydb/TE-genes.txt braker.gtf | awk '$3 ~ /transcript/' > baima_transcript.gtf
|
可以看看去除TE序列后的转录本和基因数量：
1
2
3
4
5
6
| (base) wlxie 17:03:52 ~/baima_pre_mydb
$ cat baima_transcript.gtf | wc -l
23716
(base) wlxie 17:04:05 ~/baima_pre_mydb
$ cat baima_gene_only.gtf | wc -l
20742
|
本想通过gffread软件根据处理后的gtf文件重新提取基因组的蛋白序列，但是运行过程中总是报错no genomic sequence available，原因暂时未知（可能是因为braker预测结果braker.gtf不是标准的gtf文件格式，同样用gffread做gtf2gff转换的时候会有部分信息丢失）。
可以直接写一个脚本处理braker.aa文件，根据前面筛选的TE-genes.txt文件，去除含有TE的蛋白序列，这里后续用到再做更新。
2023.10.13 更新
gffread报错的原因归根到底还是braker的gtf不是标准gtf文件：
gff文件格式没有问题，但是内容存在缺失，需要重新转换一次gff文件：
1
2
3
| cat braker.gtf | perl -ne 'if(m/\tAUGUSTUS\t/ or m/\tGeneMark.hmm\t/ or m/\tGeneMark.hmm3\t/ or m/\tgmst\t/) {print $_;}' | perl gtf2gff.pl --gff3 --out=final.gff3
# 其实第二列一共就只有AUGUSTUS、GeneMark.hmm3和gmst三种
# perl文件来自于singularity容器的/usr/share/augustus/scripts/gtf2gff.pl
|
作者说该bug已经在v3.0.3版本修复，我用的v3.0.2，只能通过这种方式重新转换。
之后筛选TE和提取蛋白：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| # 筛选含有TE的基因
grep -v "^#" braker.codingseq.rexdb.cls.tsv | cut -f1 | sort | uniq | cut -f1 -d "_" | sort | uniq > TE-genes.txt
# 去除含有TE的序列
grep -Fvf /public/home/wlxie/biosoft/TEsorter/Ap_mydb/TE-genes.txt final.gff3 > Ap_rmTE.gff3
# 从基因组重新提取去除了TE的CDS序列和蛋白序列
gffread Ap_rmTE.gff3 -g ~/Genome/Ap.fasta -x Ap_rmTE.codingseq
gffread Ap_rmTE.gff3 -g ~/Genome/Ap.fasta -y Ap_rmTE.aa
# 计算基因总长度，基因数量和平均基因长度
cat Ap_rmTE.gff3 | awk '{if ($3=="gene") print $0}' | awk '{sum+=($5 -$4)} END {print sum, NR, sum/NR}'
67389881 20094 3353.73
|
