搜索
查看: 10704|回复: 13

【Panda姐-perl练习题14】根据id列表提取fasta文件中对应的序列

[复制链接]

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
发表于 2016-9-25 12:18:54 | 显示全部楼层 |阅读模式
本帖最后由 Panda姐 于 2016-9-25 12:46 编辑

题目:
编写perl脚本,文件2是id列表,文件1是一个fasta文件,需要从文件1中提取出文件2中对应id的序列并输出(文件3)。

文件1:
>ITG_00002 | Phytophthora infestans T30-4 conserved hypothetical protein (426 nt)
ATGCATCGCTCGGGTTCCGCACGGAAAGCCCAAGGTCTGGGATTACGGGGTGGTGGTCGG
TTACACTTGGAATAACCTCGCAAATTCAGAATCTCTACAGGCTACGTTCGCGGATGGAAC
>ITG_00003 | Phytophthora infestans T30-4 protein kinase (297 nt)
ATGACGGCTGGGGTCGGTACGCCCTACTGGATCGCACCGGAGATTCTTGAAGGCAAACGG
TACACTGAGCAAGCGGATATTTACTCGTTCGGAGTGGTTTTATCCGAGCTGGACACGTGC
AAGATGCCGTTCTCTGACGTCGTTACGGCAGAGGGAAAGAAACCCAAACCAGTTCAGATC
>ITG_00005 | Phytophthora infestans T30-4 protein kinase, putative (1969 nt)
ATGCGCGTGTCTGGTCTCCTTTCAATTCTTGCAGCCACTTTGACCACGGCCCAAGACTAC
文件2:
ITG_00003
ITG_00005
文件3:
>ITG_00003 | Phytophthora infestans T30-4 protein kinase (297 nt)
ATGACGGCTGGGGTCGGTACGCCCTACTGGATCGCACCGGAGATTCTTGAAGGCAAACGG
TACACTGAGCAAGCGGATATTTACTCGTTCGGAGTGGTTTTATCCGAGCTGGACACGTGC
AAGATGCCGTTCTCTGACGTCGTTACGGCAGAGGGAAAGAAACCCAAACCAGTTCAGATC
>ITG_00005 | Phytophthora infestans T30-4 protein kinase, putative (1969 nt)
ATGCGCGTGTCTGGTCTCCTTTCAATTCTTGCAGCCACTTTGACCACGGCCCAAGACTAC

代码记录:
方法1:
[Perl] 纯文本查看 复制代码
 #use/bin/perl
           use Data::Dumper;
           %val=();
           open (IN1,"$ARGV[0]");
           open (IN2,"id.txt");
            while(<IN1>){
       chomp $_;
       @file= split(">",$_);
    $file[0]=~ s/\s*$//g;
    if ($file[0] ne ''){
    $seq = $file[0];
    }
    $file[1]=~ s/\s*$//g;
     if ($file[1] ne ''){
     $id = $file[1];
}
$val{$id} = $seq;
}
    while(<IN2>){
        $one = $_;
        $one =~ s/\s*$//g;
        print "$one\n";
        print "\n";
        #print Dumper $one;
        print  "$val{$one}\n";
            }
# Run this perl file.pl fasta.txt


方法2:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
 
use strict;
use warnings;
 
$ARGV[2] or die "use extractSeq.pl LIST FASTA OUT\n";
 
my $list = shift @ARGV;
my $fasta = shift @ARGV;
my $out = shift @ARGV;
my %select;
open L, "$list" or die;
while (<L>) {
    chomp;
    s/>//g;
    $select{$_} = 1;
}
close L;
 
$/ = "\n>";
open O, ">$out" or die;
open F, "$fasta" or die;
while (<F>) {
    s/>//g;
    my ($id) = split (/\n/, $_);
    print O ">$_" if (defined $select{$id});
}
close F;
close O;
# use extractSeq.pl LIST FASTA OUT








上一篇:生信基础版块有志之士招募录
下一篇:perl可以把字符串按长度分成小片度吗
回复

使用道具 举报

2

主题

10

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-7-24 14:19:51 | 显示全部楼层
本帖最后由 bioamin 于 2017-7-28 19:32 编辑

根据pannda姐给予的方法,实测部分无法实现,自己做了一点修改;
改动1:文件fasta不要进行chomp,防止输出序列时,序列变为一行;
改动2:在对序列输入时,使用点“.”进行连接;改动3:一个完整序列传输完毕后,将$seq清空,准备接收下一个完整序列;
代码如下:
#use/bin/perl%val();
open (IN1,"fasta.txt") ;
#fasta文件句柄
open (IN2,"id.txt") ;#id文件句柄
open(OUT,">out") ;#输出out文件句柄
while(<IN1>)
{
    @file= split(">",$_);#按行输入fasta文件时,按>分割赋予数组,如果是id行,则数组元素0为空,数组元素1为id行(包含id和详细信息);如果是数据行,则数组元素1为空,数组元素0是数据行(只是fasta的一部分数据,不是全部,需要连接)
     if ($file[0] ne '')#数组元素0为数据行
    {
    $seq.= $file[0];#有连接符,此处将一行数据赋予变量,每次叠加,数据行自带换行符
    }
    if ($file[1] ne '')#数组1为id行(包含id和详细信息)
    {
     @idfile =split (" ",$file[1]);#将id行按空格划分赋予数组,则数组元素0为真正id;
     $id=$idfile[0];
    }
   
$val{$id} = $seq;#每个id对应一个序列(序列最初不完整,每次进行覆盖)
if(m/^>/){undef $seq;}#当检测到>号时,表明新的序列开始,将$seq清空

}

while(<IN2>)
{
     chomp($_);#由于构建哈希时,键没有换行符,所以此处要去除换行符
     print OUT "$_\n";#将输出的id输出到文件中;
     print OUT "$val{$_}";#将对应的序列输出;
}
close IN1;close IN2;close out;
回复 支持 3 反对 0

使用道具 举报

0

主题

4

帖子

103

积分

注册会员

Rank: 2

积分
103
发表于 2017-10-23 19:03:21 | 显示全部楼层
#!/usr/bin/perl -w
use strict;

open  IN1, "<", $ARGV[0] or die "Can't open file!";#fasta
open  IN2, "<", $ARGV[1] or die "Can't open file!";#id
open OUT,">",$ARGV[2];#得到所需的序列

my $id;
my $seq;
my %seq;
my @id;
my $i;
my $name;
while (<IN1>) {
        chomp;
        if ($_=~/(>.*)/){
       $name = $1;
        }
        else{
       $seq{$name} .= $_;
      }
}

while (<IN2>) {
        chomp;
        push(@id,$_);
}

for ($i = 0; $i <= $#id; $i++){
        foreach $name (sort keys %seq) {
                if ($name =~ /$id[$i]/){
                        print OUT "$name\n$seq{$name}\n";
                }
        }
}



close IN1;
close IN2;
close OUT;

回复 支持 1 反对 0

使用道具 举报

1

主题

16

帖子

124

积分

注册会员

Rank: 2

积分
124
发表于 2017-7-20 16:01:11 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
open(a,'a.txt');
$/=undef;
$a=<a>;
open(b,'b.txt');
$/="\n";
open(c,">c.txt");
while(<b>){chomp($_);$a=~/(>$_.*?)\n(.*?)>/s or $a=~/(>$_.*?)\n(.*)/s;print c $1."\n".$2};
close a;
close b;
close c;
回复 支持 1 反对 0

使用道具 举报

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
 楼主| 发表于 2016-9-25 12:25:56 | 显示全部楼层
本帖最后由 Panda姐 于 2016-9-25 12:45 编辑

方法3:使用软件samtool。
[Shell] 纯文本查看 复制代码
cut -c 2- EXAMPLE.TXT | xargs -n 1 samtools faidx EXAMPLE.FA


方法4:使用python。
[Python] 纯文本查看 复制代码
f2 = open('accessionids.txt','r')
f1 = open('fasta.txt','r')
f3 = open('fasta_parsed.txt','w')
 
AI_DICT = {}
for line in f2:
    AI_DICT[line[:-1]] = 1
 
skip = 0
for line in f1:
    if line[0] == '>':
        _splitline = line.split('|')
        accessorIDWithArrow = _splitline[0]
        accessorID = accessorIDWithArrow[1:-1]
        # print accessorID
        if accessorID in AI_DICT:
            f3.write(line)
            skip = 0
        else:
            skip = 1
    else:
        if not skip:
            f3.write(line)
 
f1.close()
f2.close()
f3.close()


方法5:使用软件faSomeRecords

[Shell] 纯文本查看 复制代码
---> faSomeRecords
faSomeRecords - Extract multiple fa records
usage:
   faSomeRecords in.fa listFile out.fa
options:
   -exclude - output sequences not in the list file.


方法6:使用Ruby。
[Ruby] 纯文本查看 复制代码
#!/usr/bin/env ruby
 
require 'bio'
require 'trollop'
 
opts = Trollop::options do
  opt :get, "Contigs to get", :type => :string, :required => true
end
 
get = File.readlines(opts[:get]).map { |x| x.chomp }
 
Bio::FlatFile.foreach(ARGF) do |entry|
  id = entry.definition.split.first # no idea why I did not use entry.entry_id
  puts entry if get.include?(id)
end
# For your data you can run it like this:
# ruby ./fasta-get-contigs -g <(cat example.txt | sed 's/^>//') example.fasta > result.fasta


以上方法都是来自biostar,部分代码已测试可用,仅做为参考,请各位同学自行练习。


回复 支持 1 反对 0

使用道具 举报

0

主题

9

帖子

97

积分

注册会员

Rank: 2

积分
97
发表于 2016-10-15 05:20:50 | 显示全部楼层
使用seqtk
回复

使用道具 举报

0

主题

2

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-4-10 10:00:15 | 显示全部楼层
本帖最后由 守夜人 于 2017-4-10 10:02 编辑

[Python] 纯文本查看 复制代码
1 #!/usr/bin/env python
  2[color=#000000] # encoding: utf-8[/color]
  3              
  4 from Bio import SeqIO
  5              
  6 def exact_gene_seq(fasta_file,id_list,out_file):
  7     with open(fasta_file,'r') as ft,open(id_list,'r') as fh,open(out_file,'w') as ofh:
  8         records=SeqIO.parse(ft,'fasta')
  9         id_list=fh.readlines()
 10         print(id_list)
 11         fasta_list=[]
 12         for record in records:
 13             if record.id +"\n" in id_list:                              
 14                 fasta_list.append(record)
 15                 print(fasta_list)
 16         SeqIO.write(fasta_list,ofh,'fasta')
 17              
 18 if __name__ == "__main__":
 19     exact_gene_seq('fasta_file','id_list','out_file')
回复 支持 反对

使用道具 举报

2

主题

10

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-7-23 00:12:45 | 显示全部楼层
请问在fasta文件中序列信息是带有换行符的吧?那在第一个方法中按行读取信息时,一个序列对应的基因信息是无法一次性读入的吧?那既然一次性无法读入,如何构建哈希?
回复 支持 反对

使用道具 举报

0

主题

6

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2017-10-31 15:50:10 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!usr/bin/perl -w
open FASTA, "fasta.txt" ;
open ID, "id.txt" ;
open OUT,">fasta_id.txt";
my (@fasta,$fasta,$j);
my $i = -1;
while (<FASTA>) {
    chomp;
    $fasta = $_ ;
    if ($_ =~ />/) {
        $i++;
        $fasta[$i] = "$fasta\n";
    } else {
        $fasta[$i] .= "$fasta\n";
    }
}
close FASTA;
while (<ID>) {
    chomp;
    foreach $j(@fasta) {
        if ($j =~ /$_/) {
            print OUT $j;
        }
    }
}
close ID;
close OUT;
回复 支持 反对

使用道具 举报

1

主题

55

帖子

828

积分

高级会员

Rank: 4

积分
828
发表于 2018-7-15 11:03:45 | 显示全部楼层
本帖最后由 生信小小菜鸟 于 2018-7-15 11:11 编辑

# !usr/bin/perl -w
# 2018-6-26
# Version 1.01 ###perl test-4.pl fasta.txt id.txt result.fa.
  use strict;

  open IN1,"<",$ARGV[0] or die "Can't open file1";
  open IN2,"<",$ARGV[1] or die "Can't open file2";
  open OUT,">",$ARGV[2];

  my ($chrom,%hash,%hash1);
  while (<IN1>){
                          if (/^>(.*?)\s/){
                                   $chrom = $1;
                                   #print "$1\n";
                          }else{
                                          $hash1{$chrom}.=$_; #must start with id
                                          #print "$hash1{$chrom}\n";
                          }       
  }
  while (<IN2>){
                          chomp;
                          #print "$_\n";
                          $hash{$_} = 1;               
  }
  foreach my $chrom1(sort keys %hash1){
                          if (exists $hash{$chrom1}){
                                                  print OUT ">$chrom1\n$hash1{$chrom1}";       
                          }
  }

  close IN1;
  close IN2;
  close OUT;
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 06:25 , Processed in 0.039760 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.