手记

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

慕姐8265434
2018-12-29 23:05:09浏览 3941

1. fasta文件的介绍

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

fasta文件含有两类信息

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

详情请见fasta格式

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

fasta格式是生信中最为常见也是很容易理解的一种格式。那么使用Perl来对它又可以有哪些操作呢？

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
ACCAACGGCACGCTGTGCGCCTAAGTGCGAAAAGACAGTGCAAGCTTTGTTATGCCGAAACCTAAATGTAAAGTCAGCAACACTTCCAAATGCCACTTCTTCCCGTCGCATCCGTCTTCAGGACGATATAGTAACAGGTTTTCACTTCTCAGTGAGTGAAAGATTTTTCCCCGGGTGTACGTTGAAAGCTTCTATCGTAGAACTCATTCGAGAGGGCTTGGTGGTATTAAGAATGGTTCG
GGTGGGGGGTTCTCTTTTTTAA
>atp6
ACGATTACGCCCAACAGCCCACTTGAGCAATTTGCCATTCTCCCATTGATTCCTATGAATATTGGAAAAATTTATTTCTC
ATTCACAAATCCATCTTTGTCTATGCTGCTAACTCTCAGTTTGGTCCTACTTCTGGTTCATTTTGTTACTAAAAACGGAG
GGGGAAACTCAGTACCAAATGCTTGGCAATCCTTGGTAGAGCTTATTCATGATTTCGTGCCGAACCCGGTAAACGAACAA
ATAGGTGGTCTTTCCGGAAATGTTCAACAAAAGTTTTCCCCTCGCATCTCGGTCACTTCTACTTTTTCGTTATTTCGTAA
TCCCCAGGGTATGATACCTTATAGCTTCACAGTCACAAGTCATTTTCTCATTACTTTGGGTCTCTCATTTCCGATTTTTA
TTGGCATTACTATAGTGGGATTTCAAAGAAATGGGCTTCATTTTTTAAGCTTCTCATTACCCGCAGGAGTCCCACTGCCG
TTAGCACCTTTTTTAGTACTCCTTGAGCTAATCCCTCATTGTTTTCGCGCATTAAGCTCAGGAATACGTTTATTTGCTAA
TATGATGGCCGGTCATAGTTCAGTAAAGATTTTAAGTGGGTCCGCTTGGACTATGCTATGTATGAATGATCTTTTTTATT
TCATAGGAGATCCTGGTCCTTTATTTATAGTTCTTGCATTAACCGGTCCGGAATTAGGTGTAGCTATATCACAAGCTCAT
GTTTCTACGATCTCAATCTGTATTTACTTGAATGATGCTACAAATCTCCATCAAAGTGGTTATTTATTTATAATTGAACA
A
>atp8
ATGCCTCAACTGGATAAATTCACTTATTTCACACAATTCTTCTGGTCATGCCTTCTCCTCTTTACTTTTTATATTCCCAT
ATGCAATGATGGAGATGGAGTACTTGGGATCAGCAGAATTCTCAAACTACGGAACCAACTGGTTTCACACCGGGGGAACA
ACATCCGGAGCAACGACCCCAACAGTTTGGAAGATATCTCGAGAAAAGGTTTTAGCACCGGTGTATCCTATATGTACTCA
AGTTTATTCGAAGTATCCCAATGGTGTAACGCCGTCGACTTATTGGGAAAAAGGAGGAAAATCGCTTTGATCTCTTGTTTCGGAGAAATAAGTGGCTCACGAGGAATGGAAAGAAACATTCTATATTTGATCTCGAAGTCCTCATATAGCACTTCTTCCAGTCCTGGATGGGGGATCACTTGTAGGAATGACATAATGCTAATCCATGTTCCACACGGCCAAGGAAGCATCGTTTTTTAA
>atp9
ATGTTAGAAGGTGCAAAATCAATAGGTGCCGGAGCAGCTACAATTGCTTCAGCGGGAGCTGCTGTCGGTATTGGAAACGT
CCTTAGTTCCTCGATCCATTCCGTGGCGCGAAATCCATCATTGGCTAAACAATCATTTGGTTATGCCATTTTGGGCTTTGCTCTAACCGAAGCTATTGCATCGTTTGCCCCAATGATGGCGTCTCTGATCTCATCCGTATTCCGA

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

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

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

  • 方法1
    使用perl

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

输出为

4

  • 方法2
    使用linux命令

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

输出为

4

3.2 统计fasta的总的碱基数目

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

  • 方法1

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

输出为

2088

  • 方法2

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

输出为

2088

  • 方法3

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

输出为

2089

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

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

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

输出为

2088

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



作者：白菜代码小推车
链接：https://www.jianshu.com/p/10da73890ef0


相关标签
Python
0人推荐
随时随地看视频
慕课网APP

相关阅读

python 3.7极速入门教程5循环
Python3序列赋值、序列解包详解（下）
Python实用技法第9篇：从序列中移除重复项且保持元素间顺序不变
leetcode每日一题：376.摆动序列
python有序序列的字典序列推导式运用技巧！