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

2019年6月13日 384次阅读来源: 白菜代码小推车

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

注：本文会不断更新……

⛱fasta格式是生信中最为常见也是很容易理解的一种格式。那么使用Perl来对它又可以有哪些操作呢？下面是我在平常会用到的一些操作，现在记录下来，希望对大家有帮助。

fasta文件的介绍
注意
序列信息
3.0 显示出标题头
3.1 统计fasta文件中序列的条数
3.2 统计fasta的总的碱基数目
3.3 统计每一条序列的长度
3.4 统计N50或者N60、N70……
3.4.0 关于N50
3.4.1 程序实现
序列筛选
4.1 按照长度筛选
4.2 按照GC含量筛选
4.3 除掉零长度序列
序列提取
5.1 按照序列名称提取
5.2 按照长度提取
5.3 提取某个区域的序列
5.4 随机获取序列（待写）
序列转换
6.1 转换为反向互补序列
6.2 大小写转换
序列分离
7.1 把包含多条序列的一个fasta文件分为多个
7.2 把一条长序列切成多段
7.3 将单行的极长的序列按照特定长度换行
序列合并
序列名称
9.1 序列重名怎么办？

后记

1. fasta文件的介绍

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

fasta文件含有两类信息

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

详情请见fasta格式

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

3.0 注意

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

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

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

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

3.0 显示出标题头

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

cat 123.fasta | grep "^>"

3.1 统计fasta文件中序列的条数

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

方法1
使用perl

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

输出为

方法2
使用linux命令

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

输出为

3.2 统计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

# 使用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来割割换行符这个尾巴

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

输出为

3.3 统计每一条序列的长度

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

方法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

3.4 统计N50或者N60、N70……

3.4.0 关于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之说。这里说明一下

3.4.1 程序实现

思路：参考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中，最后打印出来，你自己可以试一试哦。

4. 筛选fasta序列

4.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

4.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";
        }
    }
'

输出为

>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTCATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAGGGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAAATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAATCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTATTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCGTTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAATATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATTTCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCATGTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACAA

4.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";
        }
    }
'

5. 提取fasta序列

5.1 按照序列名称提取

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

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

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

第二步：将序列行合一

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单行程序。其实上面的第二三步之间不需要生成中间文件，可以直接管道连接第二步与第三步。

5.2 按照长度提取

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

5.3 提取某个区域的序列

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

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

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

第二步：将序列行合一

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.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 = shift(@$line);
      my @scale = @$line;
      if ($title =~ m/>\b\Q$_->[0]\E\b/){
        for my $s (@scale){
          my @temp = split(/-/,$s);
          my ($start,$end) = ($temp[0],$temp[1]);
          my $len = $end - $start + 1;
          my $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];
    }
  }
'

记得重定向哦

6. 序列转换

6.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都能处理。你也可以尝试一下。

6.2 大小写转换

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

7. 序列分离

7.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

7.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++;
    }
  }

7.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

8. 序列合并

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

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

9. 序列名称

9.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的操作暂时想到了这些，如果有朋友能够想到更多有关的操作，希望能告诉我一下，多谢各位了！

转载希望注明来源：https://www.jianshu.com/p/10da73890ef0

    原文作者：白菜代码小推车
    原文地址: https://www.jianshu.com/p/10da73890ef0
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。

算法网