Perl 命令行实战1 - fasta文件的相关操作

之前的5篇简书小文已经说明了常用的参数的使用方法，学习了当然要致用啊！这里我们就来操作一下fasta文件。

⛱fasta格式是生信中最为常见也是很容易理解的一种格式。那么使用Perl来对它又可以有哪些操作呢？下面是我在平常会用到的一些操作，现在记录下来，希望对大家有帮助。这一篇你应该就能开始慢慢体会到perl单行的威力了！

fasta文件的介绍

>gene1
ATGAGCTGGCGATGCTGACTGTGATCTGATGCT
GTGACTGACTGACGTATGCGAGCTCAGCTGACG
TGTTAA
>gene2
ATGGCAGGCTGCAGCGATGTAGAGTCGACTTAC
GACTGTGATCTGATGCTTAGAGTCGACTTAAAA
AGTGTGGGTTGA
>gene3
ATGGCAGGCTGTGATGCTTATGTAGAGTCGAAT
GACTTTAGAGTCGACTGATGCTTAGAGTCGACT
AGTGTGGGTTGGTGTTGA

fasta文件含有两类信息

第一类是以>符号开头的，是标题头信息，记录了基因名称，有时候下载的fasta文件含有更多的信息（序列说明、编号、版本号、物种来源等等）

详情请见fasta格式

第二类是序列信息，就是跟着>打头的标题头信息的这一行的第二行开始直到遇到下一个>或者到达文件末尾就为这个标题头对应的序列了。

注意

由于我使用的windows系统进行演示的，在文件的一行结束的位置除了一个\n换行符之外，还有\r回车符这样的字符存在，而使用perl中的chomp方法不能除去\r回车符，所以下面代码中，在Mac或者Linux中可以直接写为chomp我换成了s/\r?\n//。

因为一般fasta文件的序列都是以80个字符一行，就是说序列被分成了多行。
不能直接说第1行是序列名称，那么第2行就是序列，第3行是另外一个序列名称。这是错误的！
判断的依据就是碰到下一个>或者到达文件末尾

获取fasta文件的信息

首先新建一个测试文件123.fasta
它的内容为

>atp4
ATGAGATTTAGTTCACGGGATATGCAGGATAGAAAGATGCTATTTGCTGCTATTCCATCTATTTGTGCATCAAGTCCGAA
GAAGATCTCAATCTATAATGAAGAAATGATAGTAGCTCGTCGTTTTATAGGCTTTATCATATTCAGTCGGAAGAGTTTAG
GTAAGACTTTCAAAGTGACTCTCGACGGGAGAATCGAGGCTATTCAGGAAGAATCGCAGCAATTCCCCAATCCTAACGAA
GTAGTTCCTCCGGAATCCAATGAACAACAACGATTACTTAGGATCAGCTTGCGAATTTGTGGCACCGTAGTAGAATCATT
ACCAACGGCACGCTGTGCGCCTAAGTGCGAAAAGACAGTGCAAGCTTTGTTATGCCGAAACCTAAATGTAAAGTCAGCAA
CACTTCCAAATGCCACTTCTTCCCGTCGCATCCGTCTTCAGGACGATATAGTAACAGGTTTTCACTTCTCAGTGAGTGAA
AGATTTTTCCCCGGGTGTACGTTGAAAGCTTCTATCGTAGAACTCATTCGAGAGGGCTTGGTGGTATTAAGAATGGTTCG
GGTGGGGGGTTCTCTTTTTTAA
>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTC
ATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAG
GGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAA
ATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAA
TCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTA
TTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCG
TTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAA
TATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATT
TCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCAT
GTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACA
A
>atp8
ATGCCTCAACTGGATAAATTCACTTATTTCACACAATTCTTCTGGTCATGCCTTCTCCTCTTTACTTTTTATATTCCCAT
ATGCAATGATGGAGATGGAGTACTTGGGATCAGCAGAATTCTCAAACTACGGAACCAACTGGTTTCACACCGGGGGAACA
ACATCCGGAGCAACGACCCCAACAGTTTGGAAGATATCTCGAGAAAAGGTTTTAGCACCGGTGTATCCTATATGTACTCA
AGTTTATTCGAAGTATCCCAATGGTGTAACGCCGTCGACTTATTGGGAAAAAGGAGGAAAATCGCTTTGATCTCTTGTTT
CGGAGAAATAAGTGGCTCACGAGGAATGGAAAGAAACATTCTATATTTGATCTCGAAGTCCTCATATAGCACTTCTTCCA
GTCCTGGATGGGGGATCACTTGTAGGAATGACATAATGCTAATCCATGTTCCACACGGCCAAGGAAGCATCGTTTTTTAA
>atp9
ATGTTAGAAGGTGCAAAATCAATAGGTGCCGGAGCAGCTACAATTGCTTCAGCGGGAGCTGCTGTCGGTATTGGAAACGT
CCTTAGTTCCTCGATCCATTCCGTGGCGCGAAATCCATCATTGGCTAAACAATCATTTGGTTATGCCATTTTGGGCTTTG
CTCTAACCGAAGCTATTGCATCGTTTGCCCCAATGATGGCGTCTCTGATCTCATCCGTATTCCGA

由于是实战，所以上面的序列信息尽量使用真实的信息，这些线粒体基因序列是从网上下载的。

1. 显示出标题头

思路：>打头的就是标题头

1	cat 123.fasta \| grep "^>"

2. 统计fasta文件中序列的条数

思路：有多少个>打头的行就有多少条序列

方法1
使用perl

cat 123.fasta | perl -n -e '
    # 根据fasta文件的特点，每一个以>开头的就为一个序列
    if(m/^>/){$n++;}
    # 由于只有一行指令，所以开闭括号写在一行上
    END{print "$n"}
'

输出为：

方法2
使用linux命令

1	cat 123.fasta \| grep "^>" \| wc -l

输出为：

3. 统计fasta的总的碱基数目

思路：排除>打头的行，将其他的所有字符进行数量统计

方法1

# 使用grep来排除序列名称那一行，只剩下序列
cat 123.fasta | grep -v "^>" | perl -n -e '
    # 除去末尾的回车符、换行符
    s/\r?\n//;
    $n += length($_);
    END{print $n}
'

输出为

方法2

cat 123.fasta | perl -n -e '
    s/\r?\n//;
    # 排除序列名称
    if(m/^>/){next;}
    $n += length($_);
    END{print $n}
'

输出为：

方法3

1
2
3

# 使用sed命令把行末的回车符和换行符除去
# 使用wc 的 -c 参数获得字符的计数
cat 123.fasta | grep -v "^>" | sed "s/\r//" | sed ":a;N;s/\n//g;ba" | wc -c

输出为

不知道为什么数值是2089，与其他方法比多出来了一个，不知道怎么回事，希望有朋友能够告知。
其实这种方法我是拒绝的，因为涉及到使用sed去除行末的换行符，这个问题不太好处理。

方法3.0
不过也可以用perl来割割换行符这个尾巴

1	cat 123.fasta \| grep -v "^>" \| perl -p -e ' s/\r?\n//' \| wc -c

输出为：

4. 统计每一条序列的长度

思路：遇到>打头的行就到了一条新的序列

方法1

cat 123.fasta | perl -n -e '
    s/\r?\n//;
    # 得到序列的名称
    if(m/^>(.+?)\s*$/){
        $title = $1;
    }elsif(defined $title){
    # 将这条序列的长度进行累加，直到遇到>或者文件尾
        $title_len{$title} += length($_);
    }
    # 最后打印出信息来
    # 你也可以个性化的输出
    END{
        # 
        # for my $title (sort {$title_len{$b} <=> $title_len{$a}} keys %title_len){
        for my $title (sort keys %title_len){
            print "$title","\t","$title_len{$title}","\n";
        }
    }
'

输出为

atp4    582
atp6    801
atp8    480
atp9    225

方法2
这个方法就是之前讲$/和$\这两个变量的时候说到过。

cat 123.fasta | perl -n -e '
    # 首先不要将换行符去掉，我们用来作为一个标识
    BEGIN{
        # 首先设置输入分隔符 $/
        $/ = ">";
    }
    # 正则表达式分解
    # (.+?) ： 非贪婪匹配，为了匹配序列名称
    # \s*   ： 如果序列名称后面有空格，用来匹配空格
    # \r?\n ： 匹配第一个换行符
    # 在 > 符号 与 第一个换行符 之间那么肯定是序列名称
    if(m/(.+?)\s*\r?\n/){
        $title = $1;
        # $&为匹配到的序列长度，那么之后的就是序列了。
        $sequence = (substr($_,length($&)) =~ s/\r?\n//rg);
        # 由于是以 > 作为“行”的标示符，所以末尾一般都有>。去除
        $sequence =~ s/>$//;
        $title_len{$title} = length($sequence);
    }
    END{
        for my $title (sort {$title_len{$b} <=> $title_len{$a}} keys %title_len){
            print "$title","\t","$title_len{$title}","\n";
        }
    }
'

输出为：

atp4    582
atp6    801
atp8    480
atp9    225

5. 统计N50或者N60、N70…

关于N50

N      :   10  20  30  40  50  60  70  80  90  100 
Mark   :   v   v   v   v   v   v   v   v   v   v
all    :========================================
contig1:===========        |   |   |
contig2:           ========|   |   |
contig3:                   ======= |
contig4:                          ======
contig5:                                =====
contig6:                                     ===

按照上图的，将所有的contig按照长度从大到小排列起来，首尾相连得到总的长度。
当在序列50%的位置处取一点，这一点对应的组成这个位置的contig，它的长度即为N50。例如上面的是contig3所对应的长度7。
同样的，N70就是区总序列的70%的位置，在这个位置上对应的contig的长度就是N70。例如上面的是contig4所对应的长度6。
上面为了演示，我故意写了一个100，但是实际上没有N100之说。这里说明一下

程序实现

思路：参考2.3节

# 这次借用一下List::Util模块中的求和sum方法
cat 123.fasta | perl -M'List::Util' -n -e '
    #==================# 
    #=第一部分：收集信息=#
    #==================#
    s/\r?\n//;
    # 得到序列的名称
    if(m/^>(.+?)\s*$/){
        $title = $1;
    }elsif(defined $title){
    # 将这条序列的长度进行累加，直到遇到>或者文件尾
        $title_len{$title} += length($_);
    }
    #==================# 
    #=第二部分：综合分析=#
    #==================#
    END{
        # 定义需要求得N值
        my $n = 0.5;
        # 可以将这里的代码与上图进行对比看。
        
        # 将数值进行按照从大到小的排序
        my @lengths = sort {$b <=> $a} values %title_len;
        # 求所有数值的和
        my $all     = List::Util::sum(@lengths);
        # 用来累积数值，以与总的长度进行比较
        my $accumulation = 0;
        # 遍历列表
        for my $len (@lengths){
            $accumulation += $len;
            # 如果累积值达到$all的值的一半以上
            if($accumulation > $all * $n){
                print "N50:$len";
                # 结束循环
                last;
            }
        }
    }
'

输出：

N50:582

验证一下结果

atp4    582
atp6    801
atp8    480
atp9    225

# 按照长度排序
801 582 480 225
# 总长度（2.2节）
2088
# 总长度一半
1044
# 递增
801 < 1044
801 + 582 > 1044
# 结果是
582

除了计算N50，也可以计算其他值，只需要改一下那个$n值就可以了。
然后如果想同时计算多个N值。可以将这些N值也按照从小大的顺序排列然后进行判断，将对应的N值存到Hash中，最后打印出来，你自己可以试一试哦。

筛选fasta序列

1. 按照长度筛选

思路：可以利用上述2.3节的方法。

cat 123.fasta | perl -n -e '
    s/\r?\n//;
    # 得到序列的名称
    if(m/^>(.+?)\s*$/){
        # 之前说过了，一个序列结束的标志是遇到一个>符号打头的行
        # 这个时候先不对$title进行赋值
        # 因为title和sequence是对应的
        if(defined $sequence){ # 如果是第一次进行title的赋值，那么自然就没有sequence，那么这个语句就不会执行
            # 比如我这里限制长度为500
            if(length($sequence) >= 500){
                print ">$title\n$sequence\n";
            }
        }
        # 这个时候再进行赋值
        $title = $1;
        # 由于遇到新的序列名称了，同时需要清空$sequence
        $sequence = "";
    }elsif(defined $title){
    # 将这条序列进行累加，直到遇到>或者文件尾
        $sequence .= $_;
    }
    # 最后打印出信息来
    # 你也可以个性化的输出
    END{
        # 之前说到过，除了遇到>符号打头的行之外，还有就是遇到文件尾也是序列结束的标志。这两个标志是互斥的
        # 所以最后我们还需要判断一下最后一条序列是否符合规范
        if(length($sequence) >= 500){
            print ">$title\n$sequence\n";
        }
    }
'

输出：

>atp4
ATGAGATTTAGTTCACGGGATATGCAGGATAGAAAGATGCTATTTGCTGCTATTCCATCTATTTGTGCATCAAGTCCGAAGAAGATCTCAATCTATAATGAAGAAATGATAGTAGCTCGTCGTTTTATAGGCTTTATCATATTCAGTCGGAAGAGTTTAGGTAAGACTTTCAAAGTGACTCTCGACGGGAGAATCGAGGCTATTCAGGAAGAATCGCAGCAATTCCCCAATCCTAACGAAGTAGTTCCTCCGGAATCCAATGAACAACAACGATTACTTAGGATCAGCTTGCGAATTTGTGGCACCGTAGTAGAATCATTACCAACGGCACGCTGTGCGCCTAAGTGCGAAAAGACAGTGCAAGCTTTGTTATGCCGAAACCTAAATGTAAAGTCAGCAACACTTCCAAATGCCACTTCTTCCCGTCGCATCCGTCTTCAGGACGATATAGTAACAGGTTTTCACTTCTCAGTGAGTGAAAGATTTTTCCCCGGGTGTACGTTGAAAGCTTCTATCGTAGAACTCATTCGAGAGGGCTTGGTGGTATTAAGAATGGTTCGGGTGGGGGGTTCTCTTTTTTAA
>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTCATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAGGGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAAATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAATCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTATTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATTTCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACAA

可以看到序列最后是输出到一行，这里能不能与之前一样80个碱基一行呢？

来改一下程序：

cat 123.fasta | perl -n -e '
    s/\r?\n//;
    # 得到序列的名称
    if(m/^>(.+?)\s*$/){
        # 之前说过了，一个序列结束的标志是遇到一个>符号打头的行
        # 这个时候先不对$title进行赋值
        # 因为title和sequence是对应的
        if(defined $sequence){ # 如果是第一次进行title的赋值，那么自然就没有sequence，那么这个语句就不会执行
            # 比如我这里限制长度为500
            if(length($sequence) >= 500){
                print ">$title\n$sequence\n";
            }
        }
        # 这个时候再进行赋值
        $title = $1;
        # 同时需要清空$sequence
        $sequence = "";
    }elsif(defined $title){
    # 将这条序列进行累加，直到遇到>或者文件尾
        $sequence .= $_;
    }
    # 最后打印出信息来
    # 你也可以个性化的输出
    END{
        # 之前说到过，除了遇到>符号打头的行之外，还有就是遇到文件尾也是序列结束的标志。
        # 所以最后我们还需要判断一下最后一条序列是否符合规范
        if(length($sequence) >= 500){
            print ">$title\n$sequence\n";
        }
    }
' | perl -n -e '
    chomp;
    if(m/^>/){
        print $_,"\n";
    }else{
        s/(\w{80})/$1\n/g;
        # 因为保不齐恰巧有的序列是80的倍数，如果是那样最后一行序列会存在换行符
        if(substr($_,-1,1) =~ m/\n/){
            print $_;
        }else{
            print $_,"\n";
        }
    }
'

可以看到在最后再一次通过管道，再运行一个perl命令，使用正则表达式来进行分隔，得到想要的结果。

输出为：

>atp4
ATGAGATTTAGTTCACGGGATATGCAGGATAGAAAGATGCTATTTGCTGCTATTCCATCTATTTGTGCATCAAGTCCGAA
GAAGATCTCAATCTATAATGAAGAAATGATAGTAGCTCGTCGTTTTATAGGCTTTATCATATTCAGTCGGAAGAGTTTAG
GTAAGACTTTCAAAGTGACTCTCGACGGGAGAATCGAGGCTATTCAGGAAGAATCGCAGCAATTCCCCAATCCTAACGAA
GTAGTTCCTCCGGAATCCAATGAACAACAACGATTACTTAGGATCAGCTTGCGAATTTGTGGCACCGTAGTAGAATCATT
ACCAACGGCACGCTGTGCGCCTAAGTGCGAAAAGACAGTGCAAGCTTTGTTATGCCGAAACCTAAATGTAAAGTCAGCAA
CACTTCCAAATGCCACTTCTTCCCGTCGCATCCGTCTTCAGGACGATATAGTAACAGGTTTTCACTTCTCAGTGAGTGAA
AGATTTTTCCCCGGGTGTACGTTGAAAGCTTCTATCGTAGAACTCATTCGAGAGGGCTTGGTGGTATTAAGAATGGTTCG
GGTGGGGGGTTCTCTTTTTTAA
>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTC
ATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAG
GGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAA
ATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAA
TCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTA
TTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCG
TTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAA
TATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATT
TCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCAT
GTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACA
A

2. 按照GC含量筛选

思路：前半部分与上面按照长度筛选相似

为了看起来方便，我将之前的注释都给去掉

cat 123.fasta | perl -n -e '
    BEGIN{
        # 首先可以定义一个求序列GC含量的
        sub statistic_GC_base{ # 统计序列的GC含量和数量
            my $sequence = shift;
            my $len      = length($sequence);
            my $num      = ($sequence =~ tr/GCgc/GCgc/);
            my $GC       = sprintf("%.2f",$num * 100 / $len); # 返回的是百分数（不带百分号）
            return $GC;
        }
        # 先定义一个GC含量的限制
        $GC_threshold = 40;
    }
    s/\r?\n//;
    if(m/^>(.+?)\s*$/){
        if(defined $sequence){
            # 比如我这里限制GC含量为50%，小于这个值才能通过
            if(statistic_GC_base($sequence) <= $GC_threshold){
                print ">$title\n$sequence\n";
            }
        }
        # 这个时候再进行赋值
        $title = $1;
        # 同时需要清空$sequence
        $sequence = "";
    }elsif(defined $title){
    # 将这条序列进行累加，直到遇到>或者文件尾
        $sequence .= $_;
    }
    # 最后打印出信息来
    # 你也可以个性化的输出
    END{
        # 之前说到过，除了遇到>符号打头的行之外，还有就是遇到文件尾也是序列结束的标志。
        # 所以最后我们还需要判断一下最后一条序列是否符合规范
        if(statistic_GC_base($sequence) <= $GC_threshold){
            print ">$title\n$sequence\n";
        }
    }
'

输出为：

1
2

>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTCATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAGGGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAAATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAATCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTATTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATTTCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACAA

3. 除掉零长度序列

有时候我们用程序批量得到的一些fasta序列可能只包含序列名称却不包含序列，在blast或者其他工具的使用的时候会报错而导致后续无法进行，那么就除掉它！

例如

>seq1
>seq2
>seq3
ATCGATGCTGATCGT

这种情况下seq1、seq2是不包含序列信息的，seq3是有序列的，所以这里需要将seq1、seq2除掉。这里其实就是把4.1里面的500换成了1，异曲同工。

cat 123.fasta | perl -n -e '
    s/\r?\n//;
    if(m/^>(.+?)\s*$/){
        if(defined $sequence){
            if(length($sequence) >= 1){
                print ">$title\n$sequence\n";
            }
        }
        $title = $1;
        $sequence = "";
    }elsif(defined $title){
        $sequence .= $_;
    }
    END{
        if(length($sequence) >= 1){
            print ">$title\n$sequence\n";
        }
    }
'

提取fasta序列

1. 按照序列名称提取

事实上这里适合写perl脚本然后保存脚本文件。
如果需要提取的列表多，放在文件里面的话，那么事先可以读取文件。这里比如我要提取123.fasta文件中的atp4和atp8这两个的序列

第一步：新建一个文件123.list存放这两个序列名

1 2	echo atp4 >> 123.list.txt echo atp8 >> 123.list.txt

第二步：将序列行合一

1	cat 123.fasta \| perl -p -e 's/\r?\n//;s/^>(.+)$/>$1\n/;s/^>/\n>/' > 123.merge.fasta

第三步：提取序列

cat 123.merge.fasta | perl -n -e '
  BEGIN{
    # 由于需要使用到两个文件，这里需要使用到文件句柄
    open my $file_fh,"<","./123.list.txt" or die $!;
    while(<$file_fh>){
      s/\r\n//;
      if(m/^>?(.+?)[\s\t]*$/){
        push @lookup,$1;
      }
    }
    close $file_fh;
  }
  if(m/^>/){
    my $title = $_;
    my $sequence = <>;
    if(grep {$title =~ m/\b\Q$_\E\b/} @lookup){
      printf "%s%s",$title,$sequence;
    }
  }
' > 123.extract.fasta

这里使用临时文件，为了让第三步的代码更加简洁，其实前面为了一步就达到目标写一个perl单行程序有时候不是很好，这里先将序列合并到一行之上一方面加快程序运行速度，另外一方面也可以方便编写perl单行程序。其实上面的第二三步之间不需要生成中间文件，可以直接管道连接第二步与第三步。

2. 按照长度提取

这里与上面的3.1 按照长度筛选是一样的原理，详细可以参考这一小节。

3. 提取某个区域的序列

与上面的5.1步骤有点相似，这里我们分三步走

第一步：新建一个文件123.scale.list存放这两个序列名以及需要提取的区域

1 2	echo "atp4 200-400 500-600" > 123.scale.txt echo "atp8 20-45 68-90" >> 123.scale.txt

第二步：将序列行合一

1	cat 123.fasta \| perl -p -e 's/\r?\n//;s/^>(.+)$/>$1\n/;s/^>/\n>/' > 123.merge.fasta

第三步：提取序列

cat 123.merge.fasta | perl -n -e '
      BEGIN{
        sub seq_comp_rev{
          my $r_c_seq = &seq_com(&seq_rev(shift));
          return $r_c_seq;
        }
        sub seq_com{
          return shift =~ tr/AGTCagtc/TCAGtcag/r;
        }
        sub seq_rev{
          my $temp = reverse shift;
          return $temp;
        }
        open my $file_fh,"<","./123.scale.txt" or die $!;
        while(<$file_fh>){
          s/\r?\n//;
          my @temp = split(/\s+/,$_);
          push @scales,[@temp];
        }
        close $file_fh;
      }
      s/\r?\n//;
      if(m/^>/){
        my $title = $_;
        my $sequence = <>;
        for my $line (@scales) {
          my $name = $line->[0];
          my @scale = @{$line}[1..scalar(@$line) - 1];
          if ($title =~ m/>\b\Q$name\E\b/){
            for my $s (@scale){
              my @temp = split(/-/,$s);
              my ($start,$end) = ($temp[0],$temp[1]);
              my $len = abs($end - $start) + 1;
              my $seg;
              if ($start > $end){
                  $seg = substr($sequence,$end - 1,$len);
                  $seg = seq_comp_rev($seg);
              }else{
                  $seg = substr($sequence,$start - 1,$len);
              }
              push @total , ["$title:$start-$end",$seg];
            }
          }
        }
      }
      END{
        for my $line (@total){
          printf "%s\n%s\n",$line->[0],$line->[1];
        }
      }
    '

记得重定向哦

序列转换

1. 转换为反向互补序列

将序列转换为反向互补的序列。例如ATGC转变为GCAT。实际上这个就是一个序列反转，字符替换的问题。

cat 123.fasta | perl -n -e '
  BEGIN{
    sub seq_comp_rev{
      my $r_c_seq = &seq_com(&seq_rev(shift));
      return $r_c_seq;
    }
    sub seq_com{
      return shift =~ tr/AGTCagtc/TCAGtcag/r;
    }
    sub seq_rev{
      my $temp = reverse shift;
      return $temp;
    }
  }
  s/\r?\n//;
  if(m/^>(.+?)\s*$/){
    $n++;
    if(defined $sequence){
      printf ">%s\n%s\n",$title,&seq_comp_rev($sequence);
    }
    $title = $1;
    $sequence = "";
  }elsif(defined $title){
    $sequence .= $_;
  }
  END{
    if($n >= 1){
      printf ">%s\n%s\n",$title,&seq_comp_rev($sequence);
    }
    }
' > 123.com_rev.fasta

其实这里不一定偏要什么用cat读取文件然后用perl来处理。就是说这里不管数据从cat读取文件来还是从别的东西输出通过管道来perl都能处理。你也可以尝试一下。

2 大小写转换

cat 123.fasta | perl -n -e '
    if(m/^>/){
        print;
    }else{
        # 转大写
        print uc($_);
    }
    '

序列分离

1. 把包含多条序列的一个fasta文件分为多个

mkdir temp
cat 123.fasta \
| sed '/^$/d' \
| perl -nl -e '
  BEGIN{
    $/ = ">"
  }
  s/>$//;
  if($_ ne ""){
    $title = $1 if m/^(.+?)\r?\n/;
    open my $write, ">", "./temp/$title.fasta" or die;
    print $write ">$_";
  }
'

最后得到：

./temp/atp4.fasta
./temp/atp6.fasta
./temp/atp8.fasta
./temp/atp9.fasta

2. 把一条长序列切成多段

拿一条序列来做实验，例如>atp4，把它放到一个新的文件456.fasta中

cat 456.fasta \
| perl -p -e 's/\r?\n//;s/^>(.+)$/>$1\n/;s/^>/\n>/' \
| sed '/^$/d' \
| perl -n -e '
  if(m/^>(.+)\r?\n/){
    $title = $1;
    $seq   = <>;
    # 设定长度
    my $len = 100;
    my $total = 0;
    my $count = 1;
    while(my $temp = substr($seq,0,$len,"")){
       $before = $total + 1;
       $total  = $before + length($temp) - 1;
       open my $write ,">", "${title}_$count.fasta" or die;
       print $write ">${title}:$before-$total\n$temp\n";
       close $write;
       $count++;
    }
  }

3. 将单行的极长的序列按照特定长度换行

有的时候我们拿到的序列都是排布在一行之上，特别的长，这种序列在某些软件分析的时候会出错。而且也不方便复制。
比如123.fasta：

>atp4
ATGAGATTTAGTTCACGGGATATGCAGGATAGAAAGATGCTATTTGCTGCTATTCCATCTATTTGTGCATCAAGTCCGAAGAAGATCTCAATCTATAATGAAGAAATGATAGTAGCTCGTCGTTTTATAGGCTTTATCATATTCAGTCGGAAGAGTTTAGGTAAGACTTTCAAAGTGACTCTCGACGGGAGAATCGAGGCTATTCAGGAAGAATCGCAGCAATTCCCCAATCCTAACGAAGTAGTTCCTCCGGAATCCAATGAACAACAACGATTACTTAGGATCAGCTTGCGAATTTGTGGCACCGTAGTAGAATCATTACCAACGGCACGCTGTGCGCCTAAGTGCGAAAAGACAGTGCAAGCTTTGTTATGCCGAAACCTAAATGTAAAGTCAGCAACACTTCCAAATGCCACTTCTTCCCGTCGCATCCGTCTTCAGGACGATATAGTAACAGGTTTTCACTTCTCAGTGAGTGAAAGATTTTTCCCCGGGTGTACGTTGAAAGCTTCTATCGTAGAACTCATTCGAGAGGGCTTGGTGGTATTAAGAATGGTTCGGGTGGGGGGTTCTCTTTTTTAA
>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTCATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAGGGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAAATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAATCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTATTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATTTCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACAA
>atp8
ATGCCTCAACTGGATAAATTCACTTATTTCACACAATTCTTCTGGTCATGCCTTCTCCTCTTTACTTTTTATATTCCCATATGCAATGATGGAGATGGAGTACTTGGGATCAGCAGAATTCTCAAACTACGGAACCAACTGGTTTCACACCGGGGGAACAACATCCGGAGCAACGACCCCAACAGTTTGGAAGATATCTCGAGAAAAGGTTTTAGCACCGGTGTATCCTATATGTACTCAAGTTTATTCGAAGTATCCCAATGGTGTAACGCCGTCGACTTATTGGGAAAAAGGAGGAAAATCGCTTTGATCTCTTGTTTCGGAGAAATAAGTGGCTCACGAGGAATGGAAAGAAACATTCTATATTTGATCTCGAAGTCCTCATATAGCACTTCTTCCAGTCCTGGATGGGGGATCACTTGTAGGAATGACATAATGCTAATCCATGTTCCACACGGCCAAGGAAGCATCGTTTTTTAA
>atp9
ATGTTAGAAGGTGCAAAATCAATAGGTGCCGGAGCAGCTACAATTGCTTCAGCGGGAGCTGCTGTCGGTATTGGAAACGTCCTTAGTTCCTCGATCCATTCCGTGGCGCGAAATCCATCATTGGCTAAACAATCATTTGGTTATGCCATTTTGGGCTTTGCTCTAACCGAAGCTATTGCATCGTTTGCCCCAATGATGGCGTCTCTGATCTCATCCGTATTCCGA

将这种一行的序列分割到多行上，我这里设置的是80，你也可以自己更改一下。

cat 123.fasta | perl -n -e '
    BEGIN{
        $word_num = 80;
    }
    chomp;
    if(m/^>/){
        print $_,"\n";
    }else{
        s/(\w{$word_num})/$1\n/g;
        if(substr($_,-1,1) =~ m/\n/){
            print $_;
        }else{
            print $_,"\n";
        }
    }
' > 123.multi_line.fasta

序列合并

额，这里不使用perl脚本来执行，说实话不是很合适而且也没有那个必要。

# 如果是.fa，那改为.fa即可
for i in $(ls *.fasta);
do
  cat ${i}
  echo
done > total.fasta

序列名称

1. 序列重名怎么办？

有的时候fasta文件中有重名的，比如有一文件test.fasta：

>mouse
CGTAGCTGGATG
>dog
CGATGCTGACGT
>mouse
GCTACGTACGTG
>mouse
GTACGTGAGCGT

这种fasta文件用于后续分析某些软件可能会报错，怎么做呢？有多种方法处理，这里最为推荐的是在序列名后面加上编号：

cat test.fasta | perl -n -e '
  if(m/^>/){
    chomp;
    my $temp = $_;
    my $n = 1;
    while ( exists $hash{$temp} ){
      $temp = $_ . "_" . $n;
      $n++;
    }
    $hash{$temp}++;
    print "$temp\n";
  }else{
    print;
  }
'

输出为：

>mouse
CGTAGCTGGATG
>dog
CGATGCTGACGT
>mouse_1
GCTACGTACGTG
>mouse_2
GTACGTGAGCGT

后记

这一次相较于之前的perl命令行参数的介绍，代码量明显增加了一些，但是所有的知识都是与之前相互关联的，如果需要可以查看之前的短文。

perl One-Liners | perl命令行学习5 $/和$ \变量 - 介绍$/变量
perl One-Liners | perl命令行学习4 -M与-I参数 - 介绍模块的导入
perl One-Liners | perl命令行学习3 -a和-F参数 - 介绍BEGIN{}块和 END{}块

目前有关fasta的操作暂时想到了这些，如果有朋友能够想到更多有关的操作，希望能告诉我一下，多谢各位了！