Day59-用Perl和Python脚本提取FASTA中最长转录

    这段时间简书因为一些不可描述的文章被整顿了，所以没有做学习笔记。加上最近开学，乱七八糟的事情很多，目测要下个礼拜开始上课才能够回归正轨。
    最近做了一个事情就是，注释出来的初始基因组pep文件会存在许多个转录本，很多冗余的氨基酸序列，我们需要去除这些冗余，取最长的那个转录本。
    对于序列只有一行的很容易可以用grep办到，但是还是那个问题，fasta格式存在自动换行问题，所以我们最好还是写脚本提取最长转录本比较稳妥，这里分享一下所用的脚本。
    原数据格式基本如下：

>Gene1.m1 gene=Gene1
MSVAADSPIHSSSSDDFIAYLDDALAASSPDASSDKEVENQDELESGRIKRCKFESAEETEESTSEGIVK
QNLEEYVCTHPGSFGDMCIRCGQKLDGESGVTFGYIHKGLRLHDEEISRLRNTDVKNLLIRKKLYLILDL
...

    我写的Perl脚本修改了原文件格式和顺序（哈希表的sort真是个谜），以后有时间我会考虑一下如何保留所有原格式输出，暂时将就着用吧。

#use strict;
#use warnings;

my %pep=();

my $trans = "";
my $gene = "";

open IN,"$ARGV[0]" or die "fail to open!";
open OUT,">$ARGV[1]" or die "fail to create!";
while (<IN>){
    chomp;
    my $line = $_;
    if ($line =~ /^>/){
        my @items = split /\s+/,$line;
            $gene = $items[1]; 
            $gene =~ s/gene=//g;
            $trans = $items[0];
            $trans =~ s/>//g;
    }else{
        $line =~ s/\.//g;
        $pep{$gene}{$trans}=$pep{$gene}{$trans}.$line;
    }
}
for my $key1(sort{$pep{$a} <=> $pep{$b}} keys %pep){
    my $hash2 = $pep{$key1};
    my $maxlen = 0;
    my $pepfin;
    my $transfin;
    for my $key2(sort{$hash2->{$a} <=> $hash2->{$b}} keys %$hash2){
#       print OUT $key1."\t".$key2."\t".$hash2->{$key2}."\n";
        my $len = length($hash2->{$key2});
        if ($len > $maxlen){
            $maxlen = $len;
            $pepfin = $hash2->{$key2};
            $transfin = $key2;
        }
    }
    print OUT ">".$transfin." gene=".$key1."\n";
    print OUT $pepfin."\n";
}
close IN;
close OUT;

    另外同学依据我的需求也写了一个对应功能的Python脚本，不得不承认Python v3.0的字典在保留原顺序输出方面的能力就很强，这个脚本几乎对原文件没有改动，很实用。有时间的话会对它进行全方位的注释学习，希望自己能在Python脚本书写能力上有所进步【狗头苦笑】！

import sys,getopt

def usage():
    print('usage:python3 removeRedundantProteins.py -i <in_fasta> -o <out_fasta> <-h>')
    return

def removeRedundant(in_file,out_file):
    gene_dic = {}
    flag = ''
    with open (in_file) as in_fasta:
        for line in in_fasta:
            if '>' in line:
                line1 = line.strip('>\n')
                line2 = line1.split('.')
                li = line2[0]
                flag = li
                try:
                    gene_dic[li]
                except KeyError:
                    gene_dic[li] = [line]
                else:
                    gene_dic[li].append(line)
            else:
                gene_dic[flag][-1] += line
    with open (out_file,'w') as out_fasta:
        for k,v in gene_dic.items():
            if len(v) == 1:
                out_fasta.write(gene_dic[k][0])
            else:
                trans_max = ''
                for trans in gene_dic[k]:
                    a = len(list(trans))
                    b = len(list(trans_max))
                    if a > b:
                        trans_max = trans
                out_fasta.write(trans_max)

def main(argv):
    
    try:
        opts, args = getopt.getopt(argv,'hi:o:')
    except getopt.GetoptError:
        usage()
        sys.exit()

    for opt, arg in opts:
        if opt == '-h':
            usage()
            sys.exit()
        elif opt == '-i':
            in_fasta_name = arg
        elif opt == '-o':
            outfile_name = arg
    try:
        removeRedundant(in_fasta_name,outfile_name)
    except UnboundLocalError:
        usage()

    return

if __name__ == '__main__':
    main(sys.argv[1:])