请问用Linux或Perl如何从蛋白质序列文件中取出存储在另一个文件中的10个蛋白质ID的Fasta格式序列

#!/usr/bin/perl -w

use strict

#contact perlcoder weixin

my %hash

while(<DATA>)

{

    if($_=~/\S+/)

    {

        

    }

    else

    {

       next    

    }

    chomp

    my ($key)=$_=~/>(\S+)/

    my $value=<DATA>

    chomp($value)

    $hash{$key}=$value    

}

print $hash{'Gh_A04G0739'},"\n"

__DATA__

>Gh_A12G2559              

ATGGCTACTTTCTTTGGCTATTTTAC

>Gh_A12G2554              

ATGGCGGGAACTATCCAATCCCTAAT

>Gh_A04G0755              

ATGGCCTCCGATCAGACCCTTTTTCA

>Gh_A04G0739              

ATGGAGAAAGCAAAACCTGAAGCACC

>Gh_A04G0738              

ATGGCGCAGGTTTTAGACGACGCTGA

>Gh_A07G1346              

ATGGTTCATTGTCTGCCAAAGGTTTC

>Gh_A07G1331              

ATGGCAGACAAGGATTCTTCAAGGCC

>Gh_A07G1334              

ATGACGACGCCAACTCGAGATGCAAT

>Gh_A07G1335              

ATGGCCGCAACTAGATTCCTCTCTCA

>Gh_A08G1510              

ATGGCTACTGCACCGATAAAGTCTCA

>Gh_A08G1433              

ATGGGTAAAACACCTACTGGCAAGGA

在SPDE中，查找最长转录本是用来为转录组服务的。在转录组测序完成后，测序结果中往往会有本次测序的基因序列。正基于此，将fasta格式的序列文件放入①，在保存位置②设定保存位置，然后单击按钮即可实现对fasta文件中各基因最长转录本的查找

结果如下图所示：

在得到ID后，就可以根据ID从序列中提取最长转录本的序列了。

这里需要稍微补充一句，该功能是针对转录组序列的，或者是不同序列是通过 ".1"、".2"等以及"_1"、"_2"等进行区分的（序列ID的结尾是以它们结尾的），这个时候该功能可以发挥作用。如下例子所示：

这应该是转录组序列文件命名的基本规则。除去这两种方式，如果同学们确实需要对你文件的不同转录本进行判断，这个时候似乎并没有什么太好的办法，可以考虑使用本地blast的相应功能，SPDE的界面如下：

其基本使用在该专题的其他文章中已有提到。在blast结束后，可以通过比较相似序列的比对情况来判断是否是同一个基因的不同转录本。同一个基因的不同转录本间的序列是极为相似的，甚至短的那个可能只是长的那个的一部分，从而再找出最长转录本。

1. 提取fasta文件abc.fas中序列>LG02的第164~202碱基之间序列，另存为abc_LG02_164_202.fas：

  $ more abc.fas|awk 'BEGIN{ORS=""}{if(/>/)print "\n"$0"@"else print $0}'|grep LG02|sed 's/@/\n/g'|grep -v '>'|awk '{print substr ($0,164+1,202-164-1)}' >abc_LG02_164_202.fas

在文件abc_LG02_164_202.fas第一行插入“>abc_LG02_164_202”

  $ sed  -i '1i >abc_LG02_164_202' abc_LG02_164_202.fas

2. 截取fasta文件中序列XXXX及其前后n个碱基的序列：

  $ egrep  -o ‘.{n}XXXX.{n}’  abc.fas 

3. 从fasta文件中提取列表list中的序列（list格式：>a\n>b\n）：

  $ awk 'NR==FNR{a[$0]next} /^>/{b=($1 in a)} b'  list abc.fas >list.fas

4. 在fasta格式的氨基酸序列文件中，将序列中字符“.”去掉：

  $ more abc.fas|awk '{if(!/>/) gsub (/\./,"")print}' > abc.fasta

5. 去除空行,去除wei行

  $  sed '/^\s*$/d'

  $  sed -i '$d'  a.txt 

6. 批量重命名

  $  ls |grep fastq.gz|awk -F "." '{print$1}'|xargs -I {} mv {}.fastq.gz {}.fq.gz

7. 拆分文件，可将sh脚本拆分多个脚本同时后台运行，避免运行任务过多内存不足

   $  split -l 10 abc.sh -d -a 4 abc_

8. 提取rsem分析结果的FPKM

ls|grep genes.results|awk '{print "awk \047{print $NF }\047 "$0 " >"$0".list"}'|sh

ls |grep results.list|awk 'BEGIN{ORS="\t"print "paste\tID"}{print$0}END{print" >Expression.FPKM\n"}'|sh

ls |grep results.list|sed 's/.genes.results.list//g'|awk 'BEGIN{ORS="\t"print "sed -i \0471i ID"}{print $0}END{print"\047 Expression.FPKM\n"}'|sh

more abc.list|awk '{print"grep \047\\b"$0"\\b\047 abc_Expression.FPKM >>abc.FPKM"}'|sh

---