基因組注釋01# | 重復(fù)序列注釋（RepeatModler+RepeatMasker）

寫在前面

進(jìn)行這部分重復(fù)序列提取和注釋，主要是為了屏蔽基因組的重復(fù)序列，以保證后面基因結(jié)構(gòu)注釋的準(zhǔn)確性。
大致的流程如下：
1.首先使用RepeatModler對(duì)組裝好的基因組進(jìn)行重復(fù)序列預(yù)測(cè)，從頭構(gòu)建重復(fù)序列本地庫(kù)(de novo repeat library)。
2.隨后用RepeatMasker對(duì)重復(fù)序列進(jìn)行比對(duì)注釋。注釋的流程又分兩步：

  • 基于dfam和repbase重復(fù)序列數(shù)據(jù)庫(kù)對(duì)基因組進(jìn)行近緣物種（species）的同源注釋，拿到第一輪注釋結(jié)果以及用于下一步的對(duì)重復(fù)序列屏蔽過(guò)的基因組.fasta.masked;
  • 基于RepeatModler拿到的de novo重復(fù)序列本地庫(kù)進(jìn)行注釋，但這次輸入的基因組是上面的經(jīng)過(guò)重復(fù)序列屏蔽過(guò)的.fasta.masked，這樣我們又得到第二輪注釋結(jié)果。

3.后面就是合并兩部分注釋結(jié)果，生成我們需要一些文件用于后續(xù)分析，如.out, .align, .gff, .masked.fasta等。

數(shù)據(jù)準(zhǔn)備

只需要組裝好的基因組文件genome.fasta

##選擇自己需要的基因組下載或上傳服務(wù)器，我這里從NGDC數(shù)據(jù)庫(kù)下載基因組序列文件，并解壓出來(lái)
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主頁(yè)：https://github.com/Dfam-consortium/TETools
這里用一個(gè)官方推薦的容器TETools，里面包含了RepeatModler、RepeatMasker等軟件以及他們的諸多依賴。唯一比較麻煩的是這個(gè)容器里只有Dfam數(shù)據(jù)庫(kù)而不包含RepBase重復(fù)序列數(shù)據(jù)庫(kù)，因?yàn)镽epBase是收費(fèi)的。因此需要重新自定義RepeatMasker的數(shù)據(jù)庫(kù)famDB,把這兩部分?jǐn)?shù)據(jù)庫(kù)都添加進(jìn)去。該容器的gihub網(wǎng)頁(yè)上有相應(yīng)的步驟。
GitHub上用的Docker，我這里用的singularity，根據(jù)其語(yǔ)法改了下，具體如下：

##1.下載容器，可修改容器名字，用singularity pull 去拉取Docker Hub里的容器
singularity pull dfam-tetools-latest .sif docker://dfam/tetools:latest 

##2. 從容器中復(fù)制 RepeatMasker/Libraries/ 目錄到宿主機(jī)文件系統(tǒng)
singularity shell --bind $(pwd):/work dfam-tetools-latest.sif  # 啟動(dòng)交互式Singularity容器并掛載當(dāng)前目錄
ls /work  # 確認(rèn)掛載的當(dāng)前目錄可見(jiàn)
cp -r /opt/RepeatMasker/Libraries/  /work/  # 復(fù)制容器內(nèi)的Libraries/到宿主機(jī)的當(dāng)前目錄
exit  # 退出Singularity容器的交互模式
chown -R $USER ./Libraries/  # 修改宿主機(jī)上的Libraries/ 權(quán)限

##3. 添加Dfam數(shù)據(jù)庫(kù)到宿主機(jī)的Libraries/famdb/里
Dfam數(shù)據(jù)庫(kù)有一個(gè)必需的主分區(qū)partition 0，還有其他可選分區(qū)partition 1-16。里面包含生物界各個(gè)物種分類的重復(fù)序列數(shù)據(jù)庫(kù)(https://www.dfam.org/releases/current/families/FamDB/README.txt)，根據(jù)自己研究物種進(jìn)行選擇，我這里除了partition 0，還選了partition 12 （脊椎動(dòng)物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"  # 運(yùn)行庫(kù)更新腳本，讓RepeatMasker識(shí)別這些文件并重新配置rmlib.config文件

##4. 整合RepBase數(shù)據(jù)庫(kù)到宿主機(jī)的Libraries/famdb/里
RepBase數(shù)據(jù)庫(kù)現(xiàn)在是收費(fèi)的，但我在GitHub上找到一個(gè)別人上傳的版本RepBaseRepeatMaskerEdition-20181026.tar.gz（https://github.com/yjx1217/RMRB)，將其下載上傳到服務(wù)器上。
將RepBase數(shù)據(jù)庫(kù)解壓到Libraries/目錄時(shí)，注意解壓位置，解壓出來(lái)的兩個(gè)文件會(huì)被自動(dòng)解壓到Libraries/目錄，這點(diǎn)一樣要注意，一定得在Libraries/前一級(jí)解壓
./Libraries/               # 主庫(kù)目錄
├─────  famdb/                    # Dfam庫(kù)的h5數(shù)據(jù)
│   ├── dfam39_full.0.h5  
│   ├── dfam39_full.12.h5  
│   └── rmlib.config            
├── RMRBSeqs.embl  # RepBase庫(kù)的EMBL格式數(shù)據(jù)       
├── README.RMRBSeq
└── RepeatMasker.lib          # 主庫(kù)文件

tar -xvzf RepBaseRepeatMaskerEdition-20181026.tar.gz  # 解壓RepBase庫(kù)
singularity exec --bind $(pwd)/Libraries:/opt/RepeatMasker/Libraries dfam-tetools-latest.sif \
    bash -c "cd /opt/RepeatMasker && echo '/opt/RepeatMasker/Libraries' | ./tetoolsDfamUpdate.pl"  #重新更新庫(kù)

##5. 設(shè)置環(huán)境變量
# 將以下內(nèi)容添加到 ~/.bashrc 或 ~/.bash_profile
export LIBDIR="/home/ug2117/wangtao_study/annotation/GenomeAnn_genek/container1/Libraries"
export SINGULARITY_BIND="$LIBDIR:/opt/RepeatMasker/Libraries"
tetools容器包含的軟件

由于以上將repbase庫(kù)整合到TETools容器的內(nèi)容是我后面弄的，在這之前我一直用的genek課程提供的容器TETools202309.sif。所以下面的操作都是用的這個(gè)容器。

RepeatModler預(yù)測(cè)重復(fù)序列，構(gòu)建de novo repeat library

RepeatModler2的工作流程如下圖，先通過(guò)RepeatScout/RECON進(jìn)行1-6輪的denovo檢測(cè)，再利用ltr_retriever進(jìn)行LTR結(jié)構(gòu)預(yù)測(cè)。前面6輪不一定全跑完，視物種的大小和復(fù)雜度而定。

RepeatModler2工作流程 https://doi.org/10.1073/pnas.1921046117
# 1.建庫(kù)BuildDatabase 
singularity exec ../container/TETools202309.sif BuildDatabase \
   -name SlYY \ # 數(shù)據(jù)庫(kù)名稱
   GWHERQF00000000.genome.fasta # 需要預(yù)測(cè)的基因組

# 2.運(yùn)行RepeatModeler
singularity exec ../container/TETools202309.sif RepeatModeler \
   -database SlYY \
   -threads 50 \ # 線程數(shù)
   -LTRStruct \ # 運(yùn)行LTR鑒定
   1>repeatmodeler.log 2>repeatmodeler.err

結(jié)果如下：
.fa文件就包含了預(yù)測(cè)到的所有重復(fù)序列，即de novo repeat library，另外我們這里RepeatModeler只跑了5輪。

RepeatMasker注釋重復(fù)序列

1.基于dfam和repbase數(shù)據(jù)庫(kù)進(jìn)行第一輪同源注釋

singularity exec ../container/TETools202309.sif RepeatMasker \
   -e rmblast \ # 比對(duì)軟件
   -pa 10 \ # 并行任務(wù)數(shù)，實(shí)際占用線程為pa*4
   -gff \ # 輸出gff
   -s \ # 低速模式，0-5%敏感度，比默認(rèn)模式慢2-3倍
   -a \ # 輸出align文件
   -species Teleostei \ # 物種范圍，我這個(gè)物種是魚，所以設(shè)置了硬骨魚類
   -dir ./repeatmasker_out_species \ # 第一輪注釋結(jié)果目錄
   GWHERQF00000000.genome.fasta  \ # 需要注釋的基因組
   1>repeatmasker_spe.log  2>&repeatmasker_spe.err

repeatmasker_out_species目錄下的輸出文件及.tbl文件的統(tǒng)計(jì)結(jié)果如下：

2.基于de novo repeat library進(jìn)行第二輪從頭注釋

singularity exec ../container/TETools202309.sif  RepeatMasker \ 
   -e rmblast \   
   -pa 10 \   
   -gff -s -a  \  
   -lib  SlYY-families.fa \ #  基于自己構(gòu)建的de novo repeat library（SlYY-families.fa）進(jìn)行注釋
   -dir ./repeatmasker_out_denovo  \ # 第二輪注釋結(jié)果目錄
   repeatmasker_out_species/GWHERQF00000000.genome.fasta.masked \ #  輸入為上一輪注釋輸出的masked fasta
  1>repeatmasker_denovo.log 2>repeatmasker_denovo.err

repeatmasker_out_denovo目錄下的輸出文件及.tbl文件的統(tǒng)計(jì)結(jié)果如下：

3.合并兩輪注釋結(jié)果
combineRMFiles.pl腳本只能合并兩輪注釋結(jié)果中的.out和.align文件，其他文件分別用其他幾個(gè)腳本合并。

# 構(gòu)建合并目錄
mkdir  repeatmaskerout_combine
# 合并
singularity exec ../container/TETools202309.sif   \ 
   perl /opt/RepeatMasker/util/combineRMFiles.pl  \ #combineRMFiles.pl腳本合并，注意使用容器內(nèi)路徑
   ./repeatmasker_out_species/GWHERQF00000000.genome.fasta  \ # 第一輪結(jié)果前綴
   ./repeatmasker_out_denovo/GWHERQF00000000.genome.fasta.masked  \ # 第二輪結(jié)果前綴
   repeatmaskerout_combine/SlYYgenome   # 合并結(jié)果附加前綴

合并輸出.out和.align文件：

4.重新生成.gff重復(fù)序列注釋文件
rmOutToGFF3.pl腳本將合并目錄中的.out文件轉(zhuǎn)換成.gff注釋文件。

singularity exec ../container/TETools202309.sif  rmOutToGFF3.pl  \ 
   repeatmaskerout_combine/SlYYgenome.out   \ # 合并結(jié)果中的.out文件作為輸入
   > repeatmaskerout_combine/SlYYgenome.out.gff  # 輸出到合并目錄，命名.out.gff以區(qū)分前面兩輪.gff文件

輸出.out.gff文件：

5.重新生成屏蔽過(guò)重復(fù)序列的基因組注釋文件.masked.fasta
屏蔽基因組重復(fù)序列有兩種方式：

  1. hard mask：就是把基因組中的重復(fù)序列都轉(zhuǎn)換成字母N；
  2. soft mask：就是把基因組中的重復(fù)序列轉(zhuǎn)換成小寫字母（一般基因組文件中的堿基序列都是大寫字母）；

兩種屏蔽方式都是用maskFile.pl腳本來(lái)轉(zhuǎn)換，以原始的基因組文件為基礎(chǔ)，利用合并目錄中的.out文件里的重復(fù)序列信息將其轉(zhuǎn)換成.softmask.fasta或.hardmask.fasta。兩者的區(qū)別就是加不加-softmask參數(shù)，默認(rèn)是hard mask，屏蔽過(guò)的基因組序列用于后續(xù)基因結(jié)構(gòu)注釋。

## soft mask
singularity exec ../container/TETools202309.sif  maskFile.pl  \ 
   -fasta  GWHERQF00000000.genome.fasta  \  # 原始的基因組文件作為基礎(chǔ)
   -annotations repeatmaskerout_combine/SlYYgenome.out  \  # 合并結(jié)果中的.out文件作為輸入
   -softmask  # 該參數(shù)用于指定屏蔽方式
mv ./GWHERQF00000000.genome.fasta.masked  repeatmaskerout_combine/SlYYgenome.softmask.fasta #默認(rèn)會(huì)在當(dāng)前文件夾中生成屏蔽過(guò)的序列文件.fasta.masked，將其挪到合并目錄里并修改序列名

## hard mask
singularity exec ../container/TETools202309.sif  maskFile.pl  \ 
   -fasta  GWHERQF00000000.genome.fasta  \  # 原始的基因組文件作為基礎(chǔ)
   -annotations repeatmaskerout_combine/SlYYgenome.out  \  # 合并結(jié)果中的.out文件作為輸入
mv ./GWHERQF00000000.genome.fasta.masked  repeatmaskerout_combine/SlYYgenome.hardmask.fasta #默認(rèn)會(huì)在當(dāng)前文件夾中生成屏蔽過(guò)的序列文件.fasta.masked，將其挪到合并目錄里并修改序列名

輸出兩種.masked.fasta文件：

注意

至此，我們拿到了用于后續(xù)基因結(jié)構(gòu)注釋所需的.masked.fasta，但是這個(gè)文件是將所有的重復(fù)序列都屏蔽了，包括其中的Low_complexity和Simple_repeat重復(fù)元件，這些元件在基因的UTR，Intro，甚至Exon區(qū)域都有可能出現(xiàn)，盡管概率不高，但是屏蔽這些序列勢(shì)必影響我們的基因結(jié)構(gòu)注釋時(shí)的比對(duì)識(shí)別。因此，我們?cè)谄帘沃貜?fù)序列的時(shí)候，選擇性地不屏蔽Low_complexity和Simple_repeat元件。

通過(guò)grep命令濾掉合并目錄中的.out文件里的Low_complexity和Simple_repeat行，然后再進(jìn)行上述的mask步驟，那樣輸出的.masked.fasta文件就保留了這兩種重復(fù)元件：

## 1. 重新生成.out文件
grep -v "Low_complexity"  ./repeatmaskerout_combine/SlYYgenome.out  |  \  # grep -v過(guò)濾Low_complexity行
  grep -v "Simple_repeat"   \  # grep -v進(jìn)一步過(guò)濾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ù)的基因結(jié)構(gòu)注釋。

6.最后基于木村距離繪制重復(fù)序列景觀圖
首先，使用calcDivergenceFromAlign.pl腳本將合并目錄中的.align文件里每個(gè)重復(fù)元件的木村距離信息Kimura (with divCpGMod)調(diào)出來(lái)，形成.div文件。
隨后，使用createRepeatLandscape.pl腳本將進(jìn)行圖片繪制重復(fù)序列統(tǒng)計(jì)圖。

## 1. 調(diào)取木村距離
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結(jié)果

html 結(jié)果：

寫在后面

最后看一下重復(fù)序列注釋輸出的其他文件格式吧.align、 .out、 .gff ：

.align文件格式

.out文件格式
.gff文件格式
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容