LTR-RTs的鉴定,LAI值的计算

2020-12-24 本文已影响0人 wo_monic

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.