重复序列分析文档

1 重复序列分析
重复序列广泛存在于真核生物基因组中，这些重复序列或集中成簇，或分散在基因之间，根据分布把重复序列分为分散重复序列和串联重复序列。

分散重复序列分为四种：LTR、LINE、SINE、和DNA转座子、
LTR，长末端重复转座子（long terminal repeat），是由RNA反转录而成的元件，它在两端有长大数百碱基对的LTR。

Length： 1.5-10kbp Encode reverse transcriptase Flanked by 300-1000bps terminal repeats
LINE，长散在重复序列（long interspersed nuclear elements），意为散在分布的长细胞核因子，是散在分布在哺乳动物基因组中的一类重复，这种重复序列比较长，平均长度大于1000bp，平均间隔3500-5000bp，如：rRNA，tRNA基因，形成基因家族。

SINE 为短散在重复序列（short interspersed nuclear elements）。

SINE是非自主转座的反转录转座子，来源于RNA聚合酶III的转录物，它的平均长度约为300bp，平均间隔1000bp，如：Alu家族，Hinf家族序列。

DNA 转座子: single intron-less open reading frame Encode transposase Two short inverted repeat sequences flanking the reading frame。

串联重复序列根据重复序列的重复单位的长度可分为卫星DNA、小卫星DNA 和微卫星DNA。

微卫星DNA又称为串联重复序列（short Tandem Repeat. STR）
●Simple Sequence Repeats (SSR)+Satellites
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG (G)
ATATATATATATATATATATATATATATATATATATAT (AT)n
●Lower complexity region（低复杂性区段）
TTTTTTATTTTTTGTTTTTTTTTT
（1）研究表明一些简单的重复序列与许多疾病有关。

（2）STR是存在人类基因DNA中的一类具有长度多态性的DNA序列，其多态性成为法医物证检验个人识别和亲子鉴定的丰富来源。

1.1 TRF
简介
Trf（Tandem Repeat Finder）是用来搜寻DNA序列中的串联重复序列（相临的重复两次或者多次特定核酸序列模式的重复序列）。

重复单元可以从1bp到500bp，DNA查询序列大小可以超过５Ｍ。

GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGCA GAGTA GAGTA GAGTA GAGTA
■使用
Path/trf File Match Mismatch Delta PM PI Minscore MaxPeriod [options]
例如 trf sequence.fa 2 7 7 80 10 50 800 –f –d –m
■重要参数
File = 输入序列，fasta格式。

Match = 匹配的权重，默认为2。

Mismatch = 没有匹配的权重，默认为7。

Delta = 插入删除的匹配的权重，默认为7。

PM = 匹配上的概率，默认为80。

PI = 插入的概率，默认为10。

Minscore = 最小的匹配分值，默认为50。

MaxPeriod = 最大的重复单元bp数，默认为500bp。

[options]：
-m 产生一个将串联重复序列屏蔽为Ｎ的序列文件。

-f 输出每一个串联重复序列两侧200bp 的侧翼序列，输出到比对文件中。

-d 产生一个屏蔽文件记录了与列表文件一样的信息，及比对信息，可用于后续程序的处理。

Match ,Mismatch, Delta为匹配上、没匹配上、插入的权重值。

低的权重值将允许更多的“没匹配上”、“插入”情况。

匹配的权重值“2”已被证明对“没匹配上”、“插入”的罚分权重值在3-7 范围内是有效的。

“没匹配上”、“插入”的罚分权重将被自动解释为负值。

“3”就比较宽松，“7”就比较严格。

对Match，Mismatch，Delta的推荐缺省值分别为2,7,7。

■实例
Path/trf SAC.fa 2 7 7 80 10 50 500 -f -d –m
■结果说明
程序运行完，会在目录下生成以下几个文件
*.dat -d 参数产生的屏蔽的的串联重复序列信息文件。

Tandem Repeats Finder Program writen by:
Gary Benson
Program in Bioinformatics 软件的信息
Boston University
Version 4.04
Sequence: Scaffold1
Parameters: 2 7 7 80 10 50 500 输出的参数
8891 8923 17 1.9 17 100 0 66 36 0 9 54 1.32 TTTTTAATGTTAATAAG TTTTTAATGTTAATAAG TTTTTAATGTTAATAA 串联重复
其中：8891 8923 是串联重复在输入序列里的起始、终止位置。

17 串联重复的大小bp。

1.9 重复的次数。

17 库里的重复序列得大小。

100 与邻近的重复匹配的百分比。

0 与邻近的重复匹配的插入删除百分比。

66 匹配的分值。

36 核苷酸“A”在串联重复单元里的百分比。

0 核苷酸“C”在串联重复单元里的百分比。

9 核苷酸“G”在串联重复单元里的百分比。

54 核苷酸“T”在串联重复单元里的百分比。

1.32 重复序列包含的信息量，即熵值。

*.mask -m 参数产生的串联重复序列被屏蔽为Ｎ的序列文件TTACATTCAATATTTTACATTTTACTCTACTTTTTCAC NNNNNNNNNNNNNNNNNNNNNN 重复序列处NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN N屏蔽NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN CCATG
*.html 记录串联重复序列信息的文件
*.txt.html 记录相关串联重复序列比对信息的文件
■流程说明
perl /nas/MG01/FUNGUS/PGAP/repeat/bin/repeat_pipeline.pl –trf <seq.file> 由于程序默认参数为 2 7 7 80 10 50 2000 –d –h
我们TRF重复序列预测流程只输出*.dat 文件。

然后将*.dat 文件转换成*.gff文件。

1.2 Repbase（RepeatMasker）
RepeatMasker 是一种屏蔽DNA序列中转座子重复序列和低复杂度序列的程，有Arian Smit 和Robert Hubley 开发，它将输序列中已知的重复序列都屏蔽为Ｎ或Ｘ，并给出相应的重复序列统计列表。

RepeatMasker可以选择cross_match 或wu-blast做为比对的搜索引擎。

RepeatMasker 以Repbase为默认的转座子库。

RepeatMasker 是在碱基水平上预测转座子。

■使用RepeatMasker [-options] <seqfiles(s) in fasta format>
当不带任何参数时，缺省设置是屏蔽灵长类动物所有类型的重复许乐
■重要参数
-w(ublast) **deprecated** 使用WU-blast做为搜索引擎
-de(cypher) **deprecated** 使用DeCypher做为搜索引擎
-e(ngine) [crossmatch|wublast|abblast|ncbi|decypher] 选择一个默认的搜索引擎
-pa(rallel) [number]
The number of processors to use in parallel (only works for batch
files or sequences over 50 kb)
-s 慢速搜索，敏感度提高0-5%，速度慢2-3倍
-q 快速搜索，敏感度降低5-10%，速度快2-5倍
-qq 急速搜索，敏感度降低10%，速度快4-10倍
-nolow /-low 不屏蔽低复杂度DNA 或简单重复序列
-noint /-int 只屏蔽低复杂度DNA 或简单重复序列
-norna 不屏蔽小RNA基因（假基因）
-alu 只屏蔽Alus 序列（7SLRNA，SVA and LTR5）
-lib [filename] 自定义的重复序列
Contamination options（污染参数）
-no_is 跳过细菌插入因子的检查
■结果说明
程序执行之后会生成如下几个文件：
1 *. masked 重复序列被屏蔽之后的文件
>scaffold1 0
TATTACAGGCAGCGGTATAATATTTATATTACTTAGTCTAAGTTTTAGCG TTTTGCTAAAAAAAAAAAAAACTATAATAATNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN CTTTAAAAAACTTGTAGGTAGCAAAAGCTAAAAAAAAATAACTTTATTA 重复序列被屏蔽为N
2 *.out 被比对上重复序列的说明文件
以第一行为例，其代表的意义是：
1306 = 比对上的Smith-Waterman 分值
15.6 = % 比上区间与共有序列相比的替代率
6.2 =% 在查询序列中的碱基缺失的百分率（删除碱基）
0.0 =% 在repeat 库序列中碱基缺失的百分率（插入碱基）
HSU08988 = 查询序列的名称
6563 = 比上区间在查询序列中的起始位置
7714 =比上区间在查询序列中的终止位置
（22462）= 在查询序列中超出比上区域的碱基数
C = 比上了库中重复序列得互补链
MER7A =比对上重复序列名称
DNA/MER2_type = 比对上重复序列得类型，这里是属于一个ＤＮＡ转座子ＭＥＲ２类（0）= 重复序列在比对区间之前（是互补链的情况）的碱基数（所以“0”表示比对上的区间已经是重复序列得末端）
336 =比上区间在重复序列中的起始位置
103=比上区间在重复序列中的终止位置
3 *.tbl 文件。

此文件在使用自定义的重复序列库时不产生
4 *.cat 文件此文件内容同*.out
对于真菌基因组，我们现有流程及具体参数为
perl /nas/MG01/FUNGUS/PGAP/repeat/bin/repeat_pipeline.pl –repbase –prefix <seq.file>
以上程序会自动调用RepeatMasker 程序参数默认如下：
-nolow -no_is -norna -engine wublast -parallel 1 –lib path/lib
现有流程的结果只需要*.out 文件。

然后将所有*.out 文件cat到一起，在转换为GFF 文件，最后统计相关结果。

1.3 RepeatProteinMasker
RepeatProteinMasker软件将输入的序列与RepeatProteinMasker自带的转座子蛋白库比较，找到转座子序列。

■使用
RepeatProteinMask [-pvalue #] [-minscore #] [-wordsize #] [-maxAADist] [-noLowSimple] [-noTRF] [-queryStatLen #] <fasta file>
■重要参数
-pvalue 可以接受匹配的临界值。

匹配必须有一个概率值，但概率值要比设置的pvalue值小，默认没有pvalue临界值。

但以前设定的值是0.0001
-minscore 最小分值的临界值，注意没有默认值。

因此所有比对的结果都会返回，除非设定了最小分值。

-wordsize 使用wublastx 搜索的字长，默认为3
-querystatlen query序列的有效长度用来统计计算
-maxaadist 两个blastx hit 相同时，它们之间的最大距离，默认333
-noLowSimple 关闭屏蔽/注释最终输出结果中的低复杂度和简单重复。

低复杂度和简单重复序列仍旧会出现，当搜寻匹配的蛋白转座子数据库
-noTRF 关闭输入序列中的串联重复的屏蔽或注释
■输出结果
■ 1 *. masked重复序列被屏蔽之后的文件
>scaffold164
ACTAACTAAACTAAAATAAAATAGTTAATAGTATTAAGTTATATTTAACT ATACTAAAGTAAAAAATATAAAAATACTTAAAATAGCTAGCTTTTTATAT TTAACTAATAAATATTAAGTAAAAATAAAAAATAAAGTAAAAATATAGTA GAAATAAAGTAAAATAGTAGAGTTTGTGGCAGNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTGTCTCACATTT ACCACGTATTAGATTTTTGAAAACAAGCTACTTAACTTAAAGTTGGTG
2 *.annot RepeatProteinMask软件预测的最初结果文件。

文件格式如下：
Pvalue Score Method SeqID Begin End Strand Repeat Type Begin End
1.70e-06 35 WUBlastX scaffold1 7083 7379 + PYRET LTR/Gypsy 940 1038
2.30e-29 156 WUBlastX scaffold1 8661 9689 + CFT1_pol LTR/Gypsy 806 1135
3.80e-16 219 WUBlastX scaffold1 811109 811312 - BEL1_SM LTR/Pao 183 249
对于真菌基因组RepeatProteinMask 转座子预测，现有流程的主要流程及其参数如下：
perl /nas/MG01/FUNGUS/PGAP/repeat/bin/repeat_pipeline.pl –proteinmasker –prefix <seq.file> 以上程序会自动调用RepeatProteinMask 程序。

RepeatProteinMask的参数设置如下：-noLowSimple -pvalue 1e-04
现有流程的结果只需要*.annot 文件。

然后将所有*.annot文件cat到一起，在转换为GFF 文件，最后统计相关结果。

1. 4 Denovo 预测转座子
■原理
通过buildXDFDatabase软件以自身序列建立数据库，根据数据库用RepeatModeler 建立转座子模型，利用所建模型用Repeatmasker软件找转座子。

■使用
（1）为RepeatModeler建立一个XDF（eXtended Database Format）数据库
<RepeatModelerPath>/BuildXDFDatabase -name mydb <seq.file>
-name： XDF数据库的名字
（2）运行RepeatModeler
<RepeatModelerPath>/RepeatModeler –database mydb
RepeatModeler会生成大量的输出文件。

这些输出结果会在一个以RM_<PID>.<DATE>命名的目录下，例如"RM_5098.MonMar141305172005" 。

在这个目录下，找到一个以"consensi.fa.classified"命名的文件，这个文件就是最终的转座子模型文件。

（3）利用第二步的结果进行转座子预测
<RepeatMaskerPath>/RepeatMasker -lib consensi.fa.classified <seq.file>
这里RepeatMasker预测转座子的参数与1.2 Repbase相同。

2、总结
2.1程序路径及脚本
对于以上四种重复序列的预测方法，我们已经有整合好的流程程序。

分步流程程序路径：
Path ：/nas/MG01/FUNGUS/PGAP/script/find_repeat.pl
perl /nas/MG01/FUNGUS/PGAP/script/find_repeat.pl [options] genome.fa --all run all program
--denovo run Denovo Repeat Finding
--repeatmodeler run RepeatModeler
--proteinmasker run RepeatProteinMasker
--repbase run RepeatMask by Repbase database
--trf run TRF
--outdir <str> set the output directory, default .
--prefix <str> set a prefix name for the gene ID in gff3
--cutf <int> set the number of cutted files
--run <str> set the parallel type, qsub or multi, default=qsub
--cpu <int> set the cpu number to use in parallel, default=3
--help output help information to screen
perl /nas/MG01/FUNGUS/PGAP/script/find_repeat.pl –proteinmasker –repbase –trf --denovo -repeatmodeler genome.fa
perl /nas/MG01/FUNGUS/PGAP/script/find_repeat.pl --all genome.fa
2.2重复序列结果说明
流程程序运行完之后会自动生成下表中的六个文件。

这些文件的结果都是冗余的。

在撰写结题报告的时候，我们会对GFF文件的结果进行统计。

表1 | 重复序列结果文件列表
表2 | 重复序列统计表
表 2 第一列是预测重复序列的方法；第二列是重复序列的总长；第三列是重复序列占基因组的百分比。

Total是三种方法找到的重复序列去冗余后的总的结果。

表3 |转座子分类信息统计表
表3 第一列指的是转座子的类型；第二列是用Repbase数据库预测出的转座子的结果统计；第三列是用RepeatProteinMasker的预测结果；第四列是综合两种方法的去冗余后结果。

Total是几类转座子去冗余后的综合结果。