一、简单了解原理

主要原理：构建各种重复序列的库，然后把所有的库合并在一起，最后RepeatMasker运行

有以下库：

1. LTR和EDTA鉴定的除LTR以外的非冗余的TE库。
2. repeatmodeler已知分类重复数据库
3. repeatmodeler未知分类重复数据库
4. homology-based重复库

小培老师的文档：基因组注释.note

鹏哥的文档：https://www.kdocs.cn/l/cjJbtLYZ0TWI

鹏哥打包好的镜像：repeat_mark3b ；93上的容器名：gx_repeat

二、重复序列同源基因注释

一些用到的路径
/root/miniconda3/share/RepeatMasker/famdb.py
/root/miniconda3/share/RepeatMasker/Libraries/RepeatMaskerLib.h5
/root/miniconda3/share/RepeatMasker/util/buildRMLibFromEMBL.pl
/root/miniconda3/bin/gt

python /root/miniconda3/share/RepeatMasker/famdb.py -i /root/miniconda3/share/RepeatMasker/Libraries/RepeatMaskerLib.h5 lineage -ad 物种名 #查看你的物种的所属类群
python /root/miniconda3/share/RepeatMasker/famdb.py -i /root/miniconda3/share/RepeatMasker/Libraries/RepeatMaskerLib.h5 info #查看版本信息

/root/miniconda3/share/RepeatMasker/famdb.py -i /root/miniconda3/share/RepeatMasker/Libraries/RepeatMaskerLib.h5 families -f embl -a -d Rosaceae > Rosaceae_ad.embl #由于李是蔷薇科 所以选择Rosaceae，香榧换成 Pinopsida
/root/miniconda3/share/RepeatMasker/util/buildRMLibFromEMBL.pl Rosaceae_ad.embl> Rosaceae_ad.fasta #转化格式，以便后续库之间合并在一起

同源基因重复序列的注释部分参考鹏哥的文档：https://www.kdocs.cn/l/cn2JSCKl0rg7

三、重复序列从头注释

1. 构建基因组索引

/root/miniconda3/bin/gt suffixerator -db genome.fa -indexname Ps.genome -tis -suf -lcp -des -ssp -sds -dna 

1. gt ltrharvest

nohup /root/miniconda3/bin/gt ltrharvest -index Ps.genome -similar 85 -vic 10 -seed 20 -seqids yes -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1 > ctg.harvest.scn &

1. LTR_FINDER

nohup ltr_finder -C -L 7000 -M 0.85 genome.fa > ctg.finder.scn &

1. 合并预测结果

nohup /home/soft/LTR_retriever/LTR_retriever -genome genome.fa -inharvest ctg.harvest.scn -infinder ctg.finder.scn -threads 30 &

1. EDTA

这一步要换一个容器跑，93上的镜像名：quay.io/biocontainers/edta:1.9.6--hdfd78af_2

进入你的基因组所在的路径下，运行以下命令

docker run \
-v $PWD:/in \
-w /in \
quay.io/biocontainers/edta:1.9.6--hdfd78af_2 \
EDTA.pl --genome genome.fa \
--species others \
--sensitive 0 \
--anno 0 \
--threads 32 \
> ctg.edta.log 2>&1 &

查看后台运行进程

ps -ef | grep EDTA

四、repeatmodeler鉴定其余重复序列

grep '>' genome.fa.mod.EDTA.TElib.fa | grep -v 'LTR' | wc -l #统计LTR 元件数量，槜李2042个
grep '>' genome.fa.mod.EDTA.TElib.fa | grep -v 'LTR' | sed 's/^>//' | seqkit grep -f - genome.fa.mod.EDTA.TElib.fa > nonLTR_sequences.fa #使用 sekit根据上一步获得的非LTR的序列的id提取对应的序列，得到没有 LTR 的 DNA 转座子库。数量也正好能够对上。
cat genome.fa.mod.EDTA.TElib.fa nonLTR_sequences.fa > all.LTR_noLTR.lib.fa #合并 LTR 库与 EDTA 的 DNA 转座子库

nohup RepeatMasker -norna -pa 10 -lib all.LTR_noLTR.lib.fa -dir .  \
genome.fa > ctg.round1.RepeatMasker.log 2>&1 & #使用 RepeatMasker 进行重复序列屏蔽

使用 RepeatModeler 鉴定其余的重复序列

tr -d 'N' < genome.fa.masked | seqkit seq -w 0 > rmmasked.genome.fasta
BuildDatabase -name Ps rmmasked.genome.fasta 
nohup RepeatModeler -threads 32 -database Ps -LTRStruct 1>repeatmodeler.log 2>&1 & # RepeatModeler 进行转座子重复序列鉴定

perl /home/soft/CRL_Scripts1.0/repeatmodeler_parse.pl \
  --fastafile Ps-families.fa \
  --unknowns repeatmodeler_unknowns.fasta \
  --identities repeatmodeler_identities.fasta

makeblastdb \
  -in /home/data/Tpases020812 \
  -dbtype prot \
  -out Tpases020812

/root/miniconda3/bin/blastx \
  -query repeatmodeler_unknowns.fasta \
  -db Tpases020812 \
  -evalue 1e-10 \
  -num_descriptions 10 \
  -out modelerunknown_blast_results.txt #blastx搜索未知序列

perl /home/soft/CRL_Scripts1.0/transposon_blast_parse.pl \
  --blastx modelerunknown_blast_results.txt \
  --modelerunknown repeatmodeler_unknowns.fasta #解析blastx的结果

mv unknown_elements.txt ModelerUnknown.lib
cat identified_elements.txt repeatmodeler_identities.fasta > ModelerID.lib #修改文件名

五、合并重复序列屏蔽

1. 合并所有的库

cat \
all.LTR_noLTR.lib.fa \
ModelerID.lib \
ModelerUnknown.lib \
Rosaceae_ad.fasta \
> Ps_allrepeats.lib

1. 硬屏蔽

mkdir Ps.result
cd Ps.result
nohup RepeatMasker \
  -norna \
  -pa 10 \
  -lib ../Ps_allrepeats.lib \
  -gff \
  -dir . \
  ../genome.fa &

1. 软屏蔽

跑甲基化与TE联合需要软屏蔽

nohup RepeatMasker \
  -xsmall \
  -norna \
  -pa 10 \
  -lib ../Ps_allrepeats.lib \
  -gff \
  -dir . \
  /home/gengxin/genome.fa &