提取近缘物种共有重复序列库

近日拜读了一篇文献 【Dynamics and Impacts of Transposable Element Proliferation in the Drosophila nasuta Species Group Radiation】，其中提取了几个近缘果蝇物种的重复序列库。但奈何其方法写的不够详细，本文填充了大量的个人理解，如有不同想法，欢迎讨论！

1. RepeatModeler 预测单一物种重复序列库

    RepeatModeler -LTRStruct -database ZH -pa 30 1>repeatmodeler.log 2>&1

得到结果文件 **-families.fa

2. 合并所有结果文件 并修改序列名字
   注意去掉 () 及 [] 内 内容
   将名字内 "/" 替换为任意字符 例 "@"

    cat *-families.fa > All-families.fa

3. 统计每条序列长度

    python3 01.AllLen.py All-families.fa entry.len

4.运行CD-hit 对序列聚类

    cd-hit-est -i All-families.fa -o All-families.cluster.fa -c 0.95 -n 8,9 -d 0 -M 0 -T 10
    注：-n 8,9    for thresholds 0.90 ~ 0.95
           -c   sequence identity threshold, default 0.9
           -d   length of description in .clstr file, default 20
        if set to 0, it takes the fasta defline and stops at first space
           -M   memory limit (in MB) for the program, default 800; 0 for unlimitted

5.对CD-hit 结果补充序列长度信息

    python3 02.class.clstr.py entry.len All-families.cluster.fa.clstr > All-families.typ.clstr

6.根据序列长度进行筛选 【标准】： 允许最大长度向下取10%， 需要长度大于最大值的90%

    python3 03.filter-entry-len.py All-families.typ.clstr > All-families-filter-len.clstr

7.将cluster分类： 1. 一个cluster只有一条序列 认为该序列可代表该cluster 2. 一个cluster含有多条序列，该cluster继续进行之后过滤

    python3 04.filter-len-split2.py All-families-filter-len.clstr All-families-filter-len-pass.clstr All-families-filter-len-Notpass.clstr

【以下进行self-score过滤】

8.生成符合要求的输入文件

    python3 05.change-geshi.py All-families-filter-len-Notpass.clstr All-families-filter-self.input

9.生成输入的文件夹 即 提取每条cluster序列每个cluster一个文件夹 【提前建好self文件夹】

    python3 06.selfInput.py All-families-filter-self.input All-families.fa

10.生成sh文件【脚本中路径需要更改】并运行

    python3 07.selfRun.py All-families-filter-self.input

11.统计结果并过滤 【标准】：允许最小 identity向上取10%，需要identity小于最小值的110%

    for i in *
    do
    cd $i && cat *.out > $i.final && cd ..
    done
    ll self/ | grep "Cl" | awk '{print$9}' > cluster.name   
    python3 08.selfStat.py cluster.name > All-families-self.stat

12. 将cluster分类： 1. 一个cluster只有一条序列 认为该序列可代表该cluster 2. 一个cluster含有多条序列，该cluster继续进行之后过滤

    python3 09.self-split2.py All-families-self.stat All-families-self-pass.clstr All-families-self-Notpass.clstr

【以下进行hit过滤】

13.生成符合要求的输入文件

    python3 05.change-geshi.py All-families-self-Notpass.clstr All-families-hit.input

14.生成输入的文件夹 即 提取每条cluster序列每个cluster一个文件夹 【提前建好hit文件夹】

    python3 10.hitInput.py All-families-hit.input All-families.fa

15.生成sh文件【脚本中路径需要更改】并运行 【提前对序列建库】

    python3 11.hitRun.py All-families-hit.input

16.统计每条序列的hit情况

    python3 12.hitStat.py All-families-hit.input entry.len > All-families-hit.stat

17.针对hit count 分类结果【只看最大，但会有多个最大值重复情况】

    python3 13.hit-pass.py All-families-hit.stat All-families-hit-pass.clstr All-families-hit-Notpass.clstr

【以下进行手动过滤，综合判断】

    python3 14.hand-filter.py All-families-hit-Notpass.clstr All-families-self.stat All-families-hit.stat entry.len > All-families-hand.stat
    cp All-families-hand.stat All-families-hand-pass.stat 【手动删除All-families-hand-pass.stat，至每个cluster含一条序列】

19. 合并所有pass结果，并检查cluster数是否和All-families.cluster.fa.clstr一致

    cat All-families-filter-len-pass.clstr All-families-hand-pass.stat All-families-hit-pass.clstr All-families-self-pass.clstr > All-families.3.clstr

20. 提取上述每个cluster的序列

    python3 15.extract.seq.py All-families.fa All-families.3.clstr> All-families.round-1.fa