三角褐脂藻的基因组则是经典株CCAP 1055最常用的版本 ASM15095v2,而基因组注释则来自Ensembl Protists,然而目前为止,基因组仍然没有官方的重复序列注释,所以我们自己来进行重头注释。
三角褐脂藻重复序列分析,采用RepeatModeler + RepeatMasker 的方法进行。两个软件通过bioconda,可方便地安装最新版:
conda install -c bioconda RepeatModeler=2.0.1 RepeatMasker=4.1.1
分析方法解析
首先获取基因组序列:
wget ftp://ftp.ensemblgenomes.org/pub/protists/release-48/fasta/phaeodactylum_tricornutum/dna/Phaeodactylum_tricornutum.ASM15095v2.dna_sm.toplevel.fa.gz
gunzip Phaeodactylum_tricornutum.ASM15095v2.dna_sm.toplevel.fa.gz
Ensembl的参考序列,有dna
,dna_sm
,dna_rm
的区分。若为模式物种,dna_sm
命名的序列中,重复区域的碱基会以小写字母表示,而dna_rm
序列中的重复区域,碱基统一被N
取代。然而,三角褐脂藻基因组中并未有被mask的重复序列。
然后,开始分析重复序列,依次输入下面的命令运行RepeatModeler:
Genome=Phaeodactylum_tricornutum.ASM15095v2.dna_sm.nonchromosomal.fa
cpu=48
# RepeatModeler 版本 2.0.1
mkdir ./RM_Ref
BuildDatabase $Genome -name ./RM_Ref/ptr
RepeatModeler -database ./RM_Ref/ptr -pa $cpu $Genome
这段命令需要在48线程的机器上面,运行一小时左右。最终,产生了一个文件夹:./RM_1032.ThuNov120837532020/。而文件夹当中最重要的一个文件,就是consensi.fa.classified,为三角褐脂藻基因组专门产生的一个重复序列模式文件,下面调用RepeatMasker做重复序列分析也要用到它:
## RepeatMasker版本: 4.1.1
RepeatMasker -pa $cpu -e ncbi \
-lib ./RM_1032.ThuNov120837532020/consensi.fa.classified \
-dir ./RepeatMasker -no_is $Genome
其中-e ncbi
设置序列比对方法,我们选用的比对Engine是rmblast
。在RepeatMasker的说明文档中,对该参数进行了详细描述:
-e(ngine) [crossmatch|wublast|abblast|ncbi|rmblast|hmmer]
Use an alternate search engine to the default.
Note: 'ncbi' and 'rmblast' are both aliases for the rmblastn search engine engine. The generic NCBI blastn program is not sensitive enough for use with RepeatMasker at this time.
blast
工具家族中有专门针对重复序列进行比对的rmblast
,而普通的blastn
比对重复序列不够敏感。
RepeatMasker只运行10分钟左右,即得到了结果。
重复序列结果解读
在结果文件夹./RepeatMasker中,产生了重复序列的结果文件。几个比较重要的文件如下:
- $Genome.tbl: 重复元件的Summary。该文件记录了各类主要重复元件的区间大小和比例信息。
- $Genome.mask: 重复序列被hard masked(即重复序列的碱基全部被N代替)的基因组。
- $Genome.out: RepeatMasker产生的重复序列结果的详细记录。在RepeatMasker官网上,主要模式物种的重复序列也是以这个形式记录的。另外,我使用自编脚本,将这个文件转成了gff格式,方便后续分析。
$Genome.tbl
中内容如下:
==================================================
file name: Phaeodactylum_tricornutum.ASM15095v2.dna_sm.toplevel.fa
sequences: 89
total length: 27568093 bp (27135064 bp excl N/X-runs)
GC level: 48.76 %
bases masked: 2248897 bp ( 8.16 %)
==================================================
number of length percentage
elements* occupied of sequence
--------------------------------------------------
Retroelements 1023 1236666 bp 4.49 %
SINEs: 0 0 bp 0.00 %
Penelope 0 0 bp 0.00 %
LINEs: 0 0 bp 0.00 %
CRE/SLACS 0 0 bp 0.00 %
L2/CR1/Rex 0 0 bp 0.00 %
R1/LOA/Jockey 0 0 bp 0.00 %
R2/R4/NeSL 0 0 bp 0.00 %
RTE/Bov-B 0 0 bp 0.00 %
L1/CIN4 0 0 bp 0.00 %
LTR elements: 1023 1236666 bp 4.49 %
BEL/Pao 0 0 bp 0.00 %
Ty1/Copia 1023 1236666 bp 4.49 %
Gypsy/DIRS1 0 0 bp 0.00 %
Retroviral 0 0 bp 0.00 %
DNA transposons 206 156271 bp 0.57 %
hobo-Activator 24 4610 bp 0.02 %
Tc1-IS630-Pogo 0 0 bp 0.00 %
En-Spm 0 0 bp 0.00 %
MuDR-IS905 0 0 bp 0.00 %
PiggyBac 79 84490 bp 0.31 %
Tourist/Harbinger 0 0 bp 0.00 %
Other (Mirage, 0 0 bp 0.00 %
P-element, Transib)
Rolling-circles 0 0 bp 0.00 %
Unclassified: 2350 741230 bp 2.69 %
Total interspersed repeats: 2134167 bp 7.74 %
Small RNA: 0 0 bp 0.00 %
Satellites: 0 0 bp 0.00 %
Simple repeats: 2612 109600 bp 0.40 %
Low complexity: 117 5130 bp 0.02 %
==================================================
* most repeats fragmented by insertions or deletions
have been counted as one element
RepeatMasker version 4.1.1 , default mode
run with rmblastn version 2.9.0+
The query was compared to classified sequences in ".../consensi.fa.classified"
结果中几个关键点:其一,反转录转座子中只有LER/Copia,占全基因组比例为4.49%,其二,DNA转座子占比较低。这部分结论,与2009年的BMC Genomics文章吻合:这篇文章给的LTR占比为5.8%,但是RNA转座子也只有Copia。
接下来工作建议
1、以这部分工作为基础,建立统一的基因组重复序列注释标准,将来分析其他物种;
2、 这部分的结果与先前重复序列的结果进行比较和整合,评估可否升级先前的注释内容。