LTR注释

1、安装LTR_harvest/LTR_Finder+LTR_retriever+LTR_digest
环境需要同时配置Perl、RepeatMasker、rmblast、genometools（LTR_harvest+LTR_digest）、cd-hit
conda create -n ltr_env
conda activate ltr_env
2、安装genometools（不能用conda，版本太低）
cd /mnt/e/Data/RT/software
wget https://github.com/genometools/genometools/releases/download/v1.6.6/gt-1.6.6-Linux_x86_64-64bit-barebone.tar.gz
tar -xzf gt-1.6.6-Linux_x86_64-64bit-barebone.tar.gz
cd gt-1.6.6-Linux_x86_64-64bit-barebone
设置环境变量
export PATH=$PWD/bin:$PATH
确认GenomeTools安装成功
gt -help
3、安装其他软件
conda install -c bioconda ltr_finder_parallel hmmer cd-hit blast perl trf tesorter
conda install -c bioconda -c conda-forge ltr_retriever
检查是否安装完全

点击查看代码

echo "==== CORE CHECK ===="

echo "LTR tools:"
which LTR_retriever
which LTR_FINDER_parallel
which ltr_finder

echo "Genome tools:"
which gt

echo "BLAST:"
which blastn
which makeblastdb

echo "HMMER:"
which hmmsearch

echo "TE tools:"
which RepeatMasker
which cd-hit-est
which TEsorter

4、下载玉米基因组，进入玉米基因组文件夹
cd /mnt/e/Data/RT/Zea/Ensembl-Zm-B73-dna.fa
5、识别LTR
（1）LTR_Finder:
LTR_FINDER_parallel -seq Ensembl-Zm-B73-dna.fa -threads 4 -harvest_out -size 1000000 -time 300
开新终端检查运行
ps -ef | grep LTR_FINDER
（2）LTR_harvest:
建立索引
mkdir -p index
gt suffixerator -db Ensembl-Zm-B73-dna.fa -indexname index/sample -tis -suf -lcp -des -ssp -sds -dna
需要17G内存，内存不够，换轻量级
gt suffixerator -db Ensembl-Zm-B73-dna.fa -indexname index/sample -tis -suf -lcp -des -ssp -sds -dna -memlimit 2GB
LTR_harvest预测
gt ltrharvest -index index/sample -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 2 -similar 80 -vic 10 -seed 20 -seqids yes > Ensembl-Zm-B73-dna.harvest.scn
6、LTR_retriever对前面预测结果进行精确识别、去冗余、分类注释
LTR_retriever -genome Ensembl-Zm-B73-dna.fa -inharvest Ensembl-Zm-B73-dna.harvest.scn -infinder Ensembl-Zm-B73-dna.fa.finder.combine.scn -threads 4 -u 1.3e-8 -out retriever
7、发现没有生成LTR组装指数(LTR Assembly Index, LAI，是用完整LTR-RTs在所有LTR-RTs的占比来评估基因组组装连贯性的一个指数），检查，发现在whole-genome annotation（RepeatMasker步骤）这一步Terminated，可能是内存不够，用服务器
8、上传文件到服务器
scp "E:\Data\RT\software\genometools-1.6.6.tar.gz" rhizo@100.70.93.1:"D:\zhongyunan\software"
scp "E:\Data\RT\Zea\Ensembl-Zm-B73-dna.fa\Ensembl-Zm-B73-dna.fa" rhizo@100.70.93.1:"D:\zhongyunan\RT\Zea\Ensembl-Zm-B73-dna.fa"
scp "E:\Data\RT\Zea\Ensembl-Zm-B73-dna.fa\Ensembl-Zm-B73-dna.harvest.scn" rhizo@100.70.93.1:"D:\zhongyunan\RT\Zea\Ensembl-Zm-B73-dna.fa"
scp "E:\Data\RT\Zea\Ensembl-Zm-B73-dna.fa\Ensembl-Zm-B73-dna.fa.finder.combine.scn" rhizo@100.70.93.1:"D:\zhongyunan\RT\Zea\Ensembl-Zm-B73-dna.fa"
服务器ltr_env环境中只用LTR_retriever
cd /mnt/d/zhongyunan/RT/Zea/Ensembl-Zm-B73-dna.fa/
mkdir retriever2
cd retriever2
nohup LTR_retriever -genome ../Ensembl-Zm-B73-dna.fa -inharvest ../Ensembl-Zm-B73-dna.harvest.scn -infinder ../Ensembl-Zm-B73-dna.fa.finder.combine.scn -threads 14 -u 1.3e-8 > ../retriever2.log 2>&1 &
安装TEsorter
conda activate ltr_env
conda install -c bioconda tesorter
mkdir te1
cd te1
TEsorter ../Ensembl-Zm-B73-dna.fa.mod.LTRlib.fa -db rexdb -p 10

9、提取序列（此内容在自己电脑上做的）
LTR_retriever的输出包括：有坐标和结构信息的完整基因信息[汇总表 (.pass.list)、GFF3 格式输出 (.pass.list.gff3)]、所有非冗余LTR-RT的fasta序列 (.LTRlib.fa)
提取完整LTR-RT序列
awk 'BEGIN{OFS="\t"} $3 ~ /LTR_retrotransposon/ { id=$9;sub(/^ID=/,"",id);sub(/;.*/,"",id);class = $9; sub(/.*classification=/, "", class); sub(/;.*/, "", class);name = id "|" class;print $1,$4-1,$5,name,".",$7 }' Ensembl-Zm-B73-dna.fa.mod.pass.list.gff3 > test/1LTR_RT.bed
conda activate quast_env
bedtools getfasta \ -fi Ensembl-Zm-B73-dna.fa \ -bed test/1LTR_RT.bed \ -s \ -name \ -fo test/2LTR_RT.fa
【没啥大用的步骤，可以省略跳过
注释蛋白保守结构域
conda activate ltr_env
TEsorter test/2LTR_RT.fa -p 8 -pre test/3LTR_RT
对其中含有RT结构域的DNA序列提取
grep "gene=RT" test/3LTR_RT.dom.gff3 |awk -F'\t' '{print $9}' |sed 's/.*ID=//' |sed 's/|Class.*//' |sed 's/_LTR\//|LTR\//' > test/4RT_ID.txt
conda activate tree_env
seqkit grep -n -f test/4RT_ID.txt test/2LTR_RT.fa > test/5LTR_RT_positive.fa】
直接对LTR_retriever输出的完整LTR的DNA序列进行ORF预测
conda activate emboss_env
getorf -sequence test/2LTR_RT.fa -outseq test/8LTR_RT_ORF.fa -find 1 -minsize 300 > orf_summary.txt
但是getorf运行太慢了，换成ORFipy
conda create -n orfipy_env
conda activate orfipy_env
conda install -c bioconda orfipy
orfipy --help
orfipy test/2LTR_RT.fa --pep 9LTRRT.ORF.pep.fa --dna 9LTRRT.ORF.dna.fa --bed 9LTRRT.ORF.bed --min 300 --procs 8 --outdir test --ignore-case
用HMMMER定位RT位置
ltr_env里面有hmmsearch工具
conda activate ltr_env
下载Pfam数据库
cd /mnt/e/Data/RT/software/HMM_Pfam
wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
gunzip Pfam-A.hmm.gz
hmmpress Pfam-A.hmm
cd /mnt/e/Data/RT/Zea/Ensembl-Zm-B73-dna.fa/retriever2/test
hmmsearch --cpu 4 --domtblout 10RT.domtblout /mnt/e/Data/RT/software/HMM_Pfam/Pfam-A.hmm 9LTRRT.ORF.pep.fa > 10LTRRT.hmmersearch.out
提取匹配到RVT1和RVT2的蛋白序列
seqkit grep -n -f <(grep -E "RVT_1|RVT_2" 10RT.domtblout | awk '{print $1}') 9LTRRT.ORF.pep.fa > 11RVT_ORF_candidates.pep.fa