====== 获取序列内容seq_content.pl ======

===== 序列格式 =====

基因序列文件的格式有多种，目前碰到的有如下几种：

==== FASTA格式 ====

<code>
>I_CSL_PS2507MT02
TTAAACTCAGCGGGTGGCCTCGCCCATCCTGGGGTCGCATTCCAAACGCCATCCATGGATTGACGCATGATTTTGG
ATCTCCCAAAAAGTATGATCTTCTCTAAGGCATACGACACAATAACATAGGCCTTGATTTAATCCACCCCTTATTG
TTCGCAGCCAACGAGATGATGACCTGCTCTTCAGCACACCGCGCCGAGGGGCACGAGGAGCCACTCTGCAC
</code>

这种格式文件的特点是有个“>”开头的注释行。

==== Seq格式 ====

<code>
"Contig 1" (1,717)
  Contig Length:                  717 bases
  Average Length/Sequence:        686 bases
  Total Sequence Length:         1372 bases
  Top Strand:                       1 sequences
  Bottom Strand:                    1 sequences
  Total:                            2 sequences
^^                                                                              
TTGGAAGTAAAAGTCGTAACAAGGTTTCCgTAGGTGAaCcTGCGGAaGGAtCATTGTCGAATCCTGCAATAGCAGAATGA
CCCGCTAACACGTTAAAAATTTAGGCGAGCGTCGGGAGGCCTCGGTCTCCTGTTTGCGAATCCATGGTAGGTGGCCACTC
CCGGGTGGCCACTGGCCTACAAAATCATTCGGGCGCGGTATGCGCCAAGGACCTTAAAATCGAATTGTACGTCCGTATCC
CGTTAGCGGGCATCGGCGTCTTTCCAAAATACAACGACTCTCGACAACGGATATCTCGGCTCTCGCATCGATGAAGAACG
TAGCGAAATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAACCATCGAGTCTTTGAACGCAAGTTGCGCCCGAAGCCA
CTTGGCTGAGGGCACGCCTGCCTGGGTGTCACGCATTGTATTGCCCACAAACCAGTCACACCTGAGAAGTTGTGCCGGTT
TGGGGCGGAAATTGGCCTCCCGTACCTTGTCGTGCGGTTGGCGGAAAAATGAGTCTCCGGCGATGGACGTCGCGACATCG
GTGGTTGTAAAAGACCCTCTTGTCTTGTCGCGTGAATCCTCGTCATCTTAGAGAGCTCCAGGACCCTTAGGCAGCACGTA
CTCTGTGCGCTTCGACTGTGACCCCAGGTCAGGCGGGACtACCCGCTGAGTTTAAGCATATCAATAAGCGGAGGA
</code>

这种格式文件的特点是有一行“^^”。

有时“^^”还不止一个：
<code>
^^:                                  351,652
Contig 44 (1,652)
  Contig Length:                  652 bases
  Average Length/Sequence:        652 bases
  Total Sequence Length:         1304 bases
  Top Strand:                       1 sequences
  Bottom Strand:                    1 sequences
  Total:                            2 sequences                                 
^^
TTGTCGAAACCTGCATAGCAGAACGACCCGTGAACATGTAACAACAATTGGGTGTCCTTGGTATCGGGCTCTTGTTCGAT
TAATTGGATGCCTTGTCAATGTGCGTCTTTGGTCAGCCCTCTGGGTCCTAAGGACGTCACATTGGCGCAACAACAAACCC
CCGGCACGGCATGTGCCAAGGAAAATTAAACTTAAGAAGGGCTTGTATCATGCTTCTCCGTTTACGGGGTTTGAATGAGA
CGTAGCTTCTTTATAATCACAAACGACTCTCGGCAACGGATATCTCGGCTCACGCATCGATGAAGAACGTAGCAAAATGC
GATACTTGGTGTGAATTGCAGAATCCCGTGAACCATCGAGTTTTTGAACGCAAGTTGCGCCCGAAGTCTTTCGGCCGAGG
GCACGTCTGCCTGGGCGTCACACATCGCGTCGCCCCCACCACGCCTCCTCGATGAGGATGCTTGGATGTGGGCGGAGATT
GGTCTCCCGTTCCTATGGTGCGGTTGGCTAAAACAGGAGTCCCCTTTGACGGATGCACGATTAGTGGTGGTTGACAAGAC
CCTCTTATCAAGTTGTGCGTTCTAAGGAGCAAGGAATATCTCTTCAATGACCCCAATGTGTCGTCTTGTGACGATGCTTC
GATCGCGACCCC
</code>

也有的干脆没有注释信息：
<code>
^^
GCTGCTATTGAAGCTCCATCTACAAATGGATAAGATTTCGATCTTTGTGTATACGAGTTTTTGAAAATAACGGAACAATG
CCGATTCTCTTCCAAGAAGTTGGTATTGCTCCGTTATTTATTAGGTTTTTTCTTGAATTTTTTTATTTAGGTCCTTGTCC
TTGTTTTACTTCAACAAAACAAAAAGTATTTTTATGGCTTTTGATTTAGTATCCTATTATTATTATGTGCTAATAATTAA
ATTTCCCTTTAGTTATTTTGGCTTTACAGTCATATTATTGGTTTAGAGTCAATTTAATTAAAAAATTATGGAATTTTCTG
CCGAATTAAGGTAAGATAAAAAAAGAATCATCAAAAAAGATGATCAATGGTAGGAATTGCACTCTTTTTTTTGGTAATTT
TGTAGAGAGTAGGGGGCGGATGTAGCCAA
</code>

==== 类GCG格式 ====

姑且这么称呼，因为GCG格式不完全是这样的：

<code>
        1 ggactctaga gatttacgag ctaaagttct agcgcatgaa agtcgaagta tatactttag
       61 tcgatacaaa gtccgttttt tcgaagatcc actatgataa tgaaaaagat ttctacatat
      121 ccgaccaaat cgatcaagaa tatcccaatc tgataaatcc gtccaaatgg gtttactaat
      181 aggatgcccc gatccagtac aaaattgagc ttttgataag tatcctatga ggagagtagc
      241 ggggactatg gtatcgaatt ttttcattcg agtatctatt agaaatgaat tctccagcat
      301 ttgattcctt actaacaaag tactttttgg tacacttgaa aggtacccca taaaatcgaa
      361 gcaagagttt gctaattggt ttatatggat tcttcgtggc tgagtccaaa aacagaaata
      421 atattgccag aaattgataa ggtagcattt ccatttcttc ttcaaaaaaa aagtgccttt     
      481 tgatgcaaga attccctttc cttgatatcg aacataatgc ataagaggat ccataaagaa
      541 ccatagggtt ttccgagaaa aaccagggta cattatccca aaatgttcca tcttcctaga
      601 aaagtggatt cgttccagaa aagttccaga agatgctaat ggtaagcaag aagattgttt
      661 acgaagaaac aacaaaaaaa attcatattc tgatacataa gagttatata ggaatcgaaa
      721 tagtctttta ttttcttttt gaaaaagaaa aatggatttc attgaagtaa taaaactatt
      781 ccaattcgaa tagtagttaa gaaagaatcg caataaatgc aaagatggaa catcttggat
      841 acggtattga aggagttgaa ccaggatttc caaatggata ggatagggta aatcgtcgac
      901 ctgcaggcat gc
</code>

==== Plain格式 ====

这种格式很明确，就是序列数据本身，没有别的多余信息：

<code>
AGTGACTGCGGACGACATTGTCGAACCTGCCTAGCAGAATGACCCGCGAACTTTGTCTAATACTCTCGGGGAAGCAGAAG
GTTGGTTTTTATGGCCTCCCTTTTGTTCCCTTTGTCGGGTGTGCTCGTGTTGCCCTTTGGGTGACACGCTCATTCCCCGG
CCGAACAACGAACCCCGGCGCGAAACGCGTCAAGGAACTTGAACAAGAACGCAACATCCATGCCCCGTTTCTGGGTGCTC
GTGGTGCTTGCTCTCTCATGAACAAAACGACTCTCGGCAACGGATATCTCGGCTCTCGCATCGATGAAGAACGTAGCGAA
ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAACCATCGAGTTTTTGAA
</code>

===== 程序功能 =====

程序要做的，就是从这多种格式的文件中把其中的序列数据取出来（如同Plain格式那样），以便后续放入结构化的数据库中。

根据上述对文件格式的分析，代码依次如下判断：
  - 如果文件中有一行以“>”打头，则认为是FASTA格式，序列数据从该行后读取；
  - 如果文件中有“^^”行，则认为是Seq格式，序列数据从最后一个“^^”行后读取；
  - 如果文件中有数字，则认为是类GCG格式，序列数据是去掉数字和其他多余字符后剩余的部分；
  - 文件是Plain格式。

===== 代码简述 =====

先把文件内容全部读到变量里（文件都很小）：
<code perl>
my @lines = <>;
</code>

FASTA格式判断：
<code perl>
my $flags = grep { /^>/ } @lines;
if ($flags) {
    $seq_format = "FASTA";
    for (my $i = 0; $i <= $#lines; $i++) {
        if ($lines[$i] =~ /^>/) {
            $begin_line = $i + 1;
            last;
        }
    }
</code>

变量$begin_line记录序列数据开始的行。

Seq格式判断：
<code perl>
} elsif ($flags = grep { /^\^\^/ } @lines) {
    $seq_format = "Seq";
    for (my $i = 0, my $j = 0; $i <= $#lines; $i++) {
        if ($lines[$i] =~ /^\^\^/) {
            $j++;
            if ($j == $flags) {
                $begin_line = $i + 1;
                last;
            }
        }
    }
</code>

这里变量$flags记录了文件中存在“^^”行的个数。

类GCG格式判断：
<code perl>
} elsif ($flags = grep { /[0-9]/ } @lines) {
    $seq_format = "GCG-like";
</code>

其余都认为是Plain格式：
<code perl>
} else {
    $seq_format = "Plain";
}
</code>

读取序列数据，并去掉所有多余字符：
<code perl>
for (my $i = $begin_line; $i <= $#lines; $i++) {
    $seq_content .= $lines[$i];
}
$seq_content =~ s/\r|\n|\s|-|~|[0-9]//g;
</code>

计算序列长度：
<code perl>
$seq_length = length($seq_content);
</code>

输出结果：
<code perl>
print "Format: $seq_format\n";
print "Length: $seq_length\n";
print "Content:\n";
print "\U$seq_content\n";

exit 0;
</code>

完毕。


{{tag>Perl}}