寫在前面
進行這部分重復序列提取和注釋，主要是為了屏蔽基因組的重復序列，以保證后面基因結構注釋的準確性。
大致的流程如下：
1.首先使用RepeatModler對組裝好的基因組進行重復序列預測，從頭構建重復序列本地庫(de novo repeat library)。
2.隨后用RepeatMasker對重復序列進行比對注釋。注釋的流程又分兩步：
- 基于dfam和repbase重復序列數(shù)據(jù)庫對基因組進行近緣物種（species）的同源注釋，拿到第一輪注釋結果以及用于下一步的對重復序列屏蔽過的基因組
.fasta.masked; - 基于RepeatModler拿到的de novo重復序列本地庫進行注釋，但這次輸入的基因組是上面的經(jīng)過重復序列屏蔽過的
.fasta.masked，這樣我們又得到第二輪注釋結果。
3.后面就是合并兩部分注釋結果，生成我們需要一些文件用于后續(xù)分析，如.out, .align, .gff, .masked.fasta等。
數(shù)據(jù)準備
只需要組裝好的基因組文件genome.fasta
##選擇自己需要的基因組下載或上傳服務器，我這里從NGDC數(shù)據(jù)庫下載基因組序列文件，并解壓出來
mkdir RepeatAnnoWorkdir
cd RepeatAnnoWorkdir
wget https://download.cncb.ac.cn/gwh/Animals/Silurus_lanzhouensis_Sl_YY_GWHERQF00000000/GWHERQF00000000.genome.fasta.gz
guzip GWHERQF00000000.genome.fasta.gz
軟件
TETools主頁：https://github.com/Dfam-consortium/TETools
這里用一個官方推薦的容器TETools，里面包含了RepeatModler、RepeatMasker等軟件以及他們的諸多依賴。唯一比較麻煩的是這個容器里只有Dfam數(shù)據(jù)庫而不包含RepBase重復序列數(shù)據(jù)庫，因為RepBase是收費的。因此需要重新自定義RepeatMasker的數(shù)據(jù)庫famDB,把這兩部分數(shù)據(jù)庫都添加進去。該容器的gihub網(wǎng)頁上有相應的步驟。
GitHub上用的Docker，我這里用的singularity，根據(jù)其語法改了下，具體如下：
##1.下載容器，可修改容器名字，用singularity pull 去拉取Docker Hub里的容器
singularity pull dfam-tetools-latest .sif docker://dfam/tetools:latest
##2. 從容器中復制 RepeatMasker/Libraries/ 目錄到宿主機文件系統(tǒng)
singularity shell --bind $(pwd):/work dfam-tetools-latest.sif # 啟動交互式Singularity容器并掛載當前目錄
ls /work # 確認掛載的當前目錄可見
cp -r /opt/RepeatMasker/Libraries/ /work/ # 復制容器內的Libraries/到宿主機的當前目錄
exit # 退出Singularity容器的交互模式
chown -R $USER ./Libraries/ # 修改宿主機上的Libraries/ 權限
##3. 添加Dfam數(shù)據(jù)庫到宿主機的Libraries/famdb/里
Dfam數(shù)據(jù)庫有一個必需的主分區(qū)partition 0，還有其他可選分區(qū)partition 1-16。里面包含生物界各個物種分類的重復序列數(shù)據(jù)庫(https://www.dfam.org/releases/current/families/FamDB/README.txt)，根據(jù)自己研究物種進行選擇，我這里除了partition 0，還選了partition 12 （脊椎動物Vertebrata）。
wget https://www.dfam.org/releases/current/families/FamDB/dfam39_full.0.h5.gz # 下載主分區(qū)
wget https://www.dfam.org/releases/current/families/FamDB/dfam39_full.12.h5.gz # 下載可選分區(qū)
gunzip -c dfam39_full.0.h5.gz > Libraries/famdb/dfam39_full.0.h5 # 解壓到Libraries/famdb/目錄里
gunzip -c dfam39_full.12.h5.gz > Libraries/famdb/dfam39_full.12.h5
singularity exec --bind ./Libraries:/opt/RepeatMasker/Libraries dfam-tetools-latest.sif \
bash -c "cd /opt/RepeatMasker && ./tetoolsDfamUpdate.pl" # 運行庫更新腳本，讓RepeatMasker識別這些文件并重新配置rmlib.config文件
##4. 整合RepBase數(shù)據(jù)庫到宿主機的Libraries/famdb/里
RepBase數(shù)據(jù)庫現(xiàn)在是收費的，但我在GitHub上找到一個別人上傳的版本RepBaseRepeatMaskerEdition-20181026.tar.gz（https://github.com/yjx1217/RMRB)，將其下載上傳到服務器上。
將RepBase數(shù)據(jù)庫解壓到Libraries/目錄時，注意解壓位置，解壓出來的兩個文件會被自動解壓到Libraries/目錄，這點一樣要注意，一定得在Libraries/前一級解壓
./Libraries/ # 主庫目錄
├───── famdb/ # Dfam庫的h5數(shù)據(jù)
│ ├── dfam39_full.0.h5
│ ├── dfam39_full.12.h5
│ └── rmlib.config
├── RMRBSeqs.embl # RepBase庫的EMBL格式數(shù)據(jù)
├── README.RMRBSeq
└── RepeatMasker.lib # 主庫文件
tar -xvzf RepBaseRepeatMaskerEdition-20181026.tar.gz # 解壓RepBase庫
singularity exec --bind $(pwd)/Libraries:/opt/RepeatMasker/Libraries dfam-tetools-latest.sif \
bash -c "cd /opt/RepeatMasker && echo '/opt/RepeatMasker/Libraries' | ./tetoolsDfamUpdate.pl" #重新更新庫
##5. 設置環(huán)境變量
# 將以下內容添加到 ~/.bashrc 或 ~/.bash_profile
export LIBDIR="/home/ug2117/wangtao_study/annotation/GenomeAnn_genek/container1/Libraries"
export SINGULARITY_BIND="$LIBDIR:/opt/RepeatMasker/Libraries"
由于以上將repbase庫整合到TETools容器的內容是我后面弄的，在這之前我一直用的genek課程提供的容器TETools202309.sif。所以下面的操作都是用的這個容器。
RepeatModler預測重復序列，構建de novo repeat library
RepeatModler2的工作流程如下圖，先通過RepeatScout/RECON進行1-6輪的denovo檢測，再利用ltr_retriever進行LTR結構預測。前面6輪不一定全跑完，視物種的大小和復雜度而定。
# 1.建庫BuildDatabase
singularity exec ../container/TETools202309.sif BuildDatabase \
-name SlYY \ # 數(shù)據(jù)庫名稱
GWHERQF00000000.genome.fasta # 需要預測的基因組
# 2.運行RepeatModeler
singularity exec ../container/TETools202309.sif RepeatModeler \
-database SlYY \
-threads 50 \ # 線程數(shù)
-LTRStruct \ # 運行LTR鑒定
1>repeatmodeler.log 2>repeatmodeler.err
結果如下：
.fa文件就包含了預測到的所有重復序列，即de novo repeat library，另外我們這里RepeatModeler只跑了5輪。
RepeatMasker注釋重復序列
1.基于dfam和repbase數(shù)據(jù)庫進行第一輪同源注釋
singularity exec ../container/TETools202309.sif RepeatMasker \
-e rmblast \ # 比對軟件
-pa 10 \ # 并行任務數(shù)，實際占用線程為pa*4
-gff \ # 輸出gff
-s \ # 低速模式，0-5%敏感度，比默認模式慢2-3倍
-a \ # 輸出align文件
-species Teleostei \ # 物種范圍，我這個物種是魚，所以設置了硬骨魚類
-dir ./repeatmasker_out_species \ # 第一輪注釋結果目錄
GWHERQF00000000.genome.fasta \ # 需要注釋的基因組
1>repeatmasker_spe.log 2>&repeatmasker_spe.err
repeatmasker_out_species目錄下的輸出文件及.tbl文件的統(tǒng)計結果如下：
2.基于de novo repeat library進行第二輪從頭注釋
singularity exec ../container/TETools202309.sif RepeatMasker \
-e rmblast \
-pa 10 \
-gff -s -a \
-lib SlYY-families.fa \ # 基于自己構建的de novo repeat library（SlYY-families.fa）進行注釋
-dir ./repeatmasker_out_denovo \ # 第二輪注釋結果目錄
repeatmasker_out_species/GWHERQF00000000.genome.fasta.masked \ # 輸入為上一輪注釋輸出的masked fasta
1>repeatmasker_denovo.log 2>repeatmasker_denovo.err
repeatmasker_out_denovo目錄下的輸出文件及.tbl文件的統(tǒng)計結果如下：
3.合并兩輪注釋結果
combineRMFiles.pl腳本只能合并兩輪注釋結果中的.out和.align文件，其他文件分別用其他幾個腳本合并。
# 構建合并目錄
mkdir repeatmaskerout_combine
# 合并
singularity exec ../container/TETools202309.sif \
perl /opt/RepeatMasker/util/combineRMFiles.pl \ #combineRMFiles.pl腳本合并，注意使用容器內路徑
./repeatmasker_out_species/GWHERQF00000000.genome.fasta \ # 第一輪結果前綴
./repeatmasker_out_denovo/GWHERQF00000000.genome.fasta.masked \ # 第二輪結果前綴
repeatmaskerout_combine/SlYYgenome # 合并結果附加前綴
合并輸出.out和.align文件：
4.重新生成.gff重復序列注釋文件
rmOutToGFF3.pl腳本將合并目錄中的.out文件轉換成.gff注釋文件。
singularity exec ../container/TETools202309.sif rmOutToGFF3.pl \
repeatmaskerout_combine/SlYYgenome.out \ # 合并結果中的.out文件作為輸入
> repeatmaskerout_combine/SlYYgenome.out.gff # 輸出到合并目錄，命名.out.gff以區(qū)分前面兩輪.gff文件
輸出.out.gff文件：
5.重新生成屏蔽過重復序列的基因組注釋文件.masked.fasta
屏蔽基因組重復序列有兩種方式：
- hard mask：就是把基因組中的重復序列都轉換成字母N；
- soft mask：就是把基因組中的重復序列轉換成小寫字母（一般基因組文件中的堿基序列都是大寫字母）；
兩種屏蔽方式都是用maskFile.pl腳本來轉換，以原始的基因組文件為基礎，利用合并目錄中的.out文件里的重復序列信息將其轉換成.softmask.fasta或.hardmask.fasta。兩者的區(qū)別就是加不加-softmask參數(shù)，默認是hard mask，屏蔽過的基因組序列用于后續(xù)基因結構注釋。
## soft mask
singularity exec ../container/TETools202309.sif maskFile.pl \
-fasta GWHERQF00000000.genome.fasta \ # 原始的基因組文件作為基礎
-annotations repeatmaskerout_combine/SlYYgenome.out \ # 合并結果中的.out文件作為輸入
-softmask # 該參數(shù)用于指定屏蔽方式
mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.softmask.fasta #默認會在當前文件夾中生成屏蔽過的序列文件.fasta.masked，將其挪到合并目錄里并修改序列名
## hard mask
singularity exec ../container/TETools202309.sif maskFile.pl \
-fasta GWHERQF00000000.genome.fasta \ # 原始的基因組文件作為基礎
-annotations repeatmaskerout_combine/SlYYgenome.out \ # 合并結果中的.out文件作為輸入
mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.hardmask.fasta #默認會在當前文件夾中生成屏蔽過的序列文件.fasta.masked，將其挪到合并目錄里并修改序列名
輸出兩種.masked.fasta文件：
注意
至此，我們拿到了用于后續(xù)基因結構注釋所需的
.masked.fasta，但是這個文件是將所有的重復序列都屏蔽了，包括其中的Low_complexity和Simple_repeat重復元件，這些元件在基因的UTR，Intro，甚至Exon區(qū)域都有可能出現(xiàn)，盡管概率不高，但是屏蔽這些序列勢必影響我們的基因結構注釋時的比對識別。因此，我們在屏蔽重復序列的時候，選擇性地不屏蔽Low_complexity和Simple_repeat元件。通過
grep命令濾掉合并目錄中的.out文件里的Low_complexity和Simple_repeat行，然后再進行上述的mask步驟，那樣輸出的.masked.fasta文件就保留了這兩種重復元件：## 1. 重新生成.out文件 grep -v "Low_complexity" ./repeatmaskerout_combine/SlYYgenome.out | \ # grep -v過濾Low_complexity行 grep -v "Simple_repeat" \ # grep -v進一步過濾Simple_repeat行 >./repeatmaskerout_combine/SlYYgenome.filter.out # 重新輸出不帶Low_complexity和Simple_repeat行的.filter.out文件 ## 2. 重新mask基因組 singularity exec ../container/TETools202309.sif maskFile.pl \ -fasta GWHERQF00000000.genome.fasta \ -annotations repeatmaskerout_combine/SlYYgenome.filter.out \ # 重新生成的.filter.out文件作為輸入 mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.filter_hardmask.fasta最后生成的
.filter_hardmask.fasta基因組序列文件即可用于后續(xù)的基因結構注釋。
6.最后基于木村距離繪制重復序列景觀圖
首先，使用calcDivergenceFromAlign.pl腳本將合并目錄中的.align文件里每個重復元件的木村距離信息Kimura (with divCpGMod)調出來，形成.div文件。
隨后，使用createRepeatLandscape.pl腳本將進行圖片繪制重復序列統(tǒng)計圖。
## 1. 調取木村距離
singularity exec ../container/TETools202309.sif calcDivergenceFromAlign.pl \
-s SlYY-repeat.div \ # 輸出div文件
repeatmaskerout_combine/SlYYgenome.align # 輸入合并目錄里的.align文件
## 2. 作圖
singularity exec ../container/TETools202309.sif createRepeatLandscape.pl \
-div SlYY-repeat.div \ # 輸入div文件
-g 776800000 \ # 基因組大小
> SlYY-repeat.div.html # 輸出html結果
html 結果：
寫在后面
最后看一下重復序列注釋輸出的其他文件格式吧.align、 .out、 .gff ：
