LTR-RTs的鉴定,LAI值的计算
LAI是一种新的评估基因组组装质量的评价标准。
转载自https://www.jianshu.com/p/f962d5c40fdf
LTR_retreiver
是一个整合工具,整合其他的LTR鉴定工具的结果,并且给出LAI值。
LTR-RTs 长末端重复反转录转座子的鉴定 目前推荐的是LTR_Finder
和LTR_harvest
组合鉴定,之后使用LTR_retreiver
整合两者的结果。
安装软件
Install LTR_Finder
git clone https://github.com/xzhub/LTR_Finder.git
cd LTR_Finder/source/
make
Install LTR_harvest
axel -n 16 http://genometools.org/pub/genometools-1.6.1.tar.gz
tar -zxvf genometools-1.6.1.tar.gz
cd genometools-1.6.1
make -j8 install prefix=/share/home/software/genometools/ cairo=no
pathadd /share/home/software/genometools/
source ~/.bashrc
Install LTR_retriever
git clone https://github.com/oushujun/LTR_retriever.git
或者使用conda安装LTR_retriever
conda create -n LTR_retriever
conda activate LTR_retriever
conda install -c conda-forge perl perl-text-soundex
conda install -c bioconda cd-hit repeatmasker
git clone https://github.com/oushujun/LTR_retriever.git
./LTR_retriever/LTR_retriever -h
开始分析
输入文件
基因组文件 genome.fa
输出文件
species.finder.scn 和species.harvest.scn
运行脚本
LTR.find.sh
内容如下
#species="Pco"
#genome="/share/home/database/Pco/Pco.genome.fa"
if [ $# -eq 0 ] || [ $# -eq 1 ];then
echo "Usage:
bash LTR.find.sh Pco /share/home/database/Pco/Pco.genome.fa "
exit 1
fi
species=$1
genome=$2
#LTRharvest
gt suffixerator \
-db ${genome} \
-indexname ${species} \
-tis -suf -lcp -des -ssp -sds -dna
gt ltrharvest \
-index ${species} \
-similar 90 -vic 10 -seed 20 -seqids yes \
-minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 \
-motif TGCA -motifmis 1 > ${species}.harvest.scn
# LTR_FINDER
ltr_finder -D 15000 -L 7000 -C -M 0.9 ${genome} > ${species}.finder.scn
# LTR_retriever 合并前两次的分析
LTR_retriever -genome $genome -inharvest ${species}.harvest.scn -infinder ${species}.finder.scn -threads 16
运行方法:bash LTR.find.sh 物种名 物种的参考基因组文件
例如:bash LTR.find.sh TAIR /home/genome/TAIR10.fa
物种名可以是缩写,主要用于输出文件的前缀。
程序运行比较慢,主要是ltr_finder
非常慢,没有找到多线程的方法。
LTR_retriever
最后一行的进程数,可以根据服务器性能自行修改。
注意:
LTR_retriever
要求输入的基因组文件的chr
的字符串长度不超过15,如果包含scaffold,请修改为15字以内,或者是可以舍去scffold,只留下chr。
输出结果讲解
最终会输出一个genome.fa.out.LAI
依据你输入的genome决定,例如:TaIR10.fa 输出结果就是TaIR10.fa.out.LAI .
cat TaIR10.fa.out.LAI |head
Chr From To Intact Total raw_LAI LAI
whole_genome 1 523245110 0.0967 0.4732 20.43 18.32
Chr01 1 3000000 0.1917 0.7771 24.67 21.56
Chr01 300001 3300000 0.1973 0.7971 24.76 21.65
Chr01 600001 3600000 0.1976 0.8249 23.96 20.85
Chr01 900001 3900000 0.2032 0.8621 23.57 20.46
Chr01 1200001 4200000 0.1981 0.8769 22.59 19.48
第7列是每个位置的LAI的值,第2行最后一个就是整个基因组的LAI值。这个基因组是18.32可以看出质量不错。个人觉得LAI>15都挺好。0.0967代表完整LTR-RT在基因组中占比9.67% , 0.4732代表LTR在基因组中比例为47.32%.
LTR_retriever
最后过滤后的LTR_RTs文件是genome.fa.pass.list
。这里面是过滤重复后通过的LTR_RTs.可以看出里面分为Copia和Gypsy,还有unknown。
head genome.fa.pass.list
#LTR_loc Category Motif TSD 5_TSD 3_TSD Internal Identity Strand SuperFamily TE_type Insertion_Time
Chr01:109718..119470 pass motif:TGCA TSD:AAAAC 109713..109717 119471..119475 IN:111581..117607 0.9973 - Copia LTR 103354
Chr01:151156..160815 pass motif:TGCA TSD:ACCTT 151151..151155 160816..160820 IN:152959..159011 0.9945 ? unknown NA 214113
Chr01:259702..269481 pass motif:TGCA TSD:CGTTG 259697..259701 269482..269486 IN:261573..267609 0.9963 + Copia LTR 144256
Chr01:282588..292147 pass motif:TGCA TSD:CAATG 282583..282587 292148..292152 IN:284331..290404 0.9966 ? unknown NA 132702
Chr01:406690..416424 pass motif:TGCA TSD:AAGGG 406685..406689 416425..416429 IN:408590..414524 0.9979 + Copia LTR 81086
Chr01:434887..439473 pass motif:TGCA TSD:GGCAA 434882..434886 439474..439478 IN:435372..438989 0.9959 - Gypsy LTR 159372
species.LTR.Insertion_Time.csv是genome.fa.pass.list最后一列
可视化物种的LTR的年代的密度图
library("ggplot2")
#读入文件
dat<-read.csv("species.LTR.Insertion_Time.csv",header=FALSE)
#除以100万
VAF<-dat$V1 / 1000000
#画图(出图结果中x轴是以mya(百万年)为单位的。)
ggplot(dat,aes(x=VAF))+geom_density(color = "red")+xlab('Mya')+ylab('Density')+
scale_x_continuous(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0))+ #设置横纵坐标从0开始
theme_classic()
ggsave('Speies.LTR.density.pdf',dpi = 300)
#求峰值
y_peak <- which.max(density(VAF)$y)#找y值最大的
x_peak <- density(VAF)$x[y_peak]#找出最大的y所对应的x
#LTR.typetime.csv是genome.fa.pass.list最后一列和倒数第三列。
#读入文件
dat1<-read.csv("LTR.typetime.csv",header=FALSE)
#预处理
names(dat1) <- c("VAF1","Type")
dat1$VAF1<-dat1$VAF1 / 1000000
#画图
ggplot(dat1,aes(x=VAF1))+geom_density(aes(color = Type))+xlab('Mya')+ylab('Density')+
scale_x_continuous(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0))+ #设置横纵坐标从0开始
theme_classic()
ggsave('SpO_LTR_type.density.pdf',dpi = 300)
By default the mutation rate is 1.3e-8 per
bp per year (rice), so the unit is calendar year. You may specify a
different rate by providing -u/-miu.默认的变异速率是1.3E-8.