使用三代测序数据能获得较好的、甚至完整的基因组序列。通过检测基因组序列两端的属于端粒（Telomere）的特定的重复序列，可以知道基因组组装是否得到完整的染色体水平的序列。若在序列中间检测到端粒序列，可以知道基因组组装过程中对Contigs有错误的连接。

在染色体序列首尾存在端粒序列。在人类中，端粒序列由重复单元TTAGGG，串联重复约2500次组成。 不同的物种，端粒的重复单元可能不一样，可以在端粒数据库中查询。 我对一种子囊菌使用PacBio测序数据进行了基因组组装，得到12条序列，发现大部分序列首尾均出现了TTAGGG/CCCTAA的串联重复序列。我对一种昆虫物种PacBio组装基因组序列进行分析，发现序列首尾出现一些长度较长微卫星重复序列，而没有固定的重复单元，这和端粒数据库中的结果一致。

编写Perl程序对全基因组序列的端粒序列进行搜索，查看基因组序列的完整情况：

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 genome.fasta > telomere_info.txt

    大部分物种端粒序列的重复单元是TTAGGG/CCCTAA。本程序能在基因组中寻找端粒重复单元的串联重复序列，并给出位点信息。

    --split-length <int>    default: 100000
    --overlap-length <int>    default: 10000
    程序会将每条序列打断后进行重复单元搜索。这两个参数设置打断的序列长度和相邻两序列之间的重叠长度。

    --repeat-unit <string>    default: TTAGGG
    设置重复单元碱基序列，该重复单元的反向互补也将作为重复单元进行搜索。可以在端粒数据库（http://telomerase.asu.edu/sequences_telomere.html）中寻找目标端粒重复单元。
    vertebrate sp.      TTAGGG
    plants sp.          TTTAGGG
    Pezizomycotina      TTAGGG

    --min-repeat-num <int>    default: 4
    设置重复单元最小重复次数。


USAGE
if (@ARGV==0){die $usage}

my ($splitLength, $overlapLength, $repeatunit, $minRepeatNum);
GetOptions(
    "split-length:i" => \$splitLength,
    "overlap-length:i" => \$overlapLength,
    "repeat-unit:s" => \$repeatunit,
    "min-repeat-num:s" => \$minRepeatNum,
);
$splitLength ||= 100000;
$overlapLength ||= 10000;
$repeatunit ||= "TTAGGG";
$repeatunit = uc($repeatunit);
my $repeatunit_reverse = reverse $repeatunit;
$repeatunit_reverse =~ tr/ATCG/TAGC/;
$minRepeatNum ||= 4;

# 读取基因组序列
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
my (%seq, $seq_id);
while (<IN>) {
    chomp;
    if (m/^>(\S+)/) {
        $seq_id = $1;
    }
    else {
        $_ = uc($_);
        $seq{$seq_id} .= $_;
    }
}
close IN;

# 将基因组序列打断
my (%seq_split, %seq_length);
foreach my $id (keys %seq) {
    my $seq = $seq{$id};
    my $length = length($seq);
    $seq_length{$id} = $length;
    my $pos = 0;
    while ($pos < $length) {
        $seq_split{$id}{$pos} = substr($seq, $pos, $splitLength + $overlapLength);
        $pos += $splitLength;
    }
}

print "SeqID\tSeqLength\tStart\tEnd\tLength\tType\n";
foreach my $id (sort keys %seq_split) {
    foreach my $pos (sort {$a <=> $b} keys %{$seq_split{$id}}) {
        my $seq = $seq_split{$id}{$pos};
        while ($seq =~ m/(($repeatunit){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit\n";
        }
        while ($seq =~ m/(($repeatunit_reverse){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit_reverse\n";
        }
    }
}