使用序列 ID 提取 FASTA 序列
extract FASTA sequence using sequence ID
我有两个文件:
文件 1:
84C2_Locus_14_Transcript_1/3_Confidence_0.571_Length_1244
AAACTAGTCAATAGAGAAAATCCAAAGTGGATGAAATTGAAGTGATTGTATGGCACAAGT...so on
>84C2_Locus_14_Transcript_2/3_Confidence_0.857_Length_1961
AAACTAGTCAATAGAGAAAATCCAAAGTGGATGAAATTGAAGTGATTGTATGGCACAAGT...so on
>84C2|Locus_15_Transcript_1/9_Confidence_0.190_Length_757
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA.....so on
>84C2_Locus_15_Transcript_5/9_Confidence_0.333_Length_1841
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA....so on
File2: 只有序列 IDS
84C2_Locus_14_Transcript_1/3_Confidence_0.571_Length_1244
84C2_Locus_14_Transcript_2/3_Confidence_0.857_Length_1961
84C2_Locus_14_Transcript_3/3_Confidence_0.571_Length_1248
84C2_Locus_15_Transcript_1/9_Confidence_0.190_Length_757
............这么多。
我的输出文件应该包含与 header 关联的序列。即匹配序列 id 文件 header 部分与原始 fasta 序列文件,那些序列 header 匹配 fasta 序列 header 存储在另一个包含 header 和 [=25= 的输出文件中] 像这样:
原始输出文件应该是这样的:
>84C2_Locus_15_Transcript_5/9_Confidence_0.333_Length_1841
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA......so on
请给我建议适用于我的问题的 perl 方法。
编辑:新代码,第一次没看懂你的问题。该脚本读取两个文件,将所有 ID 存储在哈希中,然后遍历第一个文件中的所有序列。只有那些 ID 在第二个文件中的序列才会写入输出文件。请注意,输出是用 |
作为分隔符编写的,就像在第一个文件中一样,而不是像在第二个文件中的 ID 那样使用 _
。第一个文件开头缺少的 >
也被正确复制。
#!/usr/bin/perl -w
use strict;
# Check command line arguments
unless ($#ARGV == 1 && -e $ARGV[0] && -e $ARGV[1]) {
print STDERR "Usage: split-fasta.pl DATA IDS\n";
exit 1;
}
my (%ids, $id);
# Read sequence IDs
open(IDS, "<$ARGV[1]") or die "Can't open IDs file: ";
while (<IDS>) {
$ids{$_} = 1;
}
close(IDS);
# Read sequence data and write results to output.fasta
open(DATA, "<$ARGV[0]") or die "Can't open sequences file: ";
open(OUTFILE, ">>output.fasta") or die "Can't open output file: ";
while (<DATA>) {
my $line = $_;
if ($line =~ /^>?(\w+\|.+\n)/) {
$id = ;
$id =~ tr/|/_/;
}
print OUTFILE $line if defined $ids{$id};
}
close(DATA);
close(OUTFILE);
那么,您基本上是想将文件 1 拆分成多个文件,每个文件只包含一个序列?也许这个丑陋的(可能关闭每个文件句柄,而不仅仅是最后一个......)Perl 片段会帮助你。
#!/usr/bin/perl -w
use strict;
while (<>) {
my $line = $_;
if ($line =~ />?([0-9]+|\w+)\//) {
my $file_name = ;
open(OUTFILE, ">>$file_name");
}
print OUTFILE $line;
}
close(OUTFILE);
编辑: 如果您想从第二个文件中输入 ID,只需添加另一个 if
,使其与来自 headers 的 ID 相匹配第一个文件。只要两个文件的顺序相同,这应该就可以工作,因为输入在写入时会立即被丢弃并且不能再被搜索。
我有两个文件: 文件 1:
84C2_Locus_14_Transcript_1/3_Confidence_0.571_Length_1244
AAACTAGTCAATAGAGAAAATCCAAAGTGGATGAAATTGAAGTGATTGTATGGCACAAGT...so on
>84C2_Locus_14_Transcript_2/3_Confidence_0.857_Length_1961
AAACTAGTCAATAGAGAAAATCCAAAGTGGATGAAATTGAAGTGATTGTATGGCACAAGT...so on
>84C2|Locus_15_Transcript_1/9_Confidence_0.190_Length_757
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA.....so on
>84C2_Locus_15_Transcript_5/9_Confidence_0.333_Length_1841
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA....so on
File2: 只有序列 IDS
84C2_Locus_14_Transcript_1/3_Confidence_0.571_Length_1244
84C2_Locus_14_Transcript_2/3_Confidence_0.857_Length_1961
84C2_Locus_14_Transcript_3/3_Confidence_0.571_Length_1248
84C2_Locus_15_Transcript_1/9_Confidence_0.190_Length_757
............这么多。
我的输出文件应该包含与 header 关联的序列。即匹配序列 id 文件 header 部分与原始 fasta 序列文件,那些序列 header 匹配 fasta 序列 header 存储在另一个包含 header 和 [=25= 的输出文件中] 像这样:
原始输出文件应该是这样的:
>84C2_Locus_15_Transcript_5/9_Confidence_0.333_Length_1841
ATTTGCTCGGAAAAACACTTCTCGTGGAACTTGTTAGCGCTGAGCTTGATCCCAAGACGA......so on
请给我建议适用于我的问题的 perl 方法。
编辑:新代码,第一次没看懂你的问题。该脚本读取两个文件,将所有 ID 存储在哈希中,然后遍历第一个文件中的所有序列。只有那些 ID 在第二个文件中的序列才会写入输出文件。请注意,输出是用 |
作为分隔符编写的,就像在第一个文件中一样,而不是像在第二个文件中的 ID 那样使用 _
。第一个文件开头缺少的 >
也被正确复制。
#!/usr/bin/perl -w
use strict;
# Check command line arguments
unless ($#ARGV == 1 && -e $ARGV[0] && -e $ARGV[1]) {
print STDERR "Usage: split-fasta.pl DATA IDS\n";
exit 1;
}
my (%ids, $id);
# Read sequence IDs
open(IDS, "<$ARGV[1]") or die "Can't open IDs file: ";
while (<IDS>) {
$ids{$_} = 1;
}
close(IDS);
# Read sequence data and write results to output.fasta
open(DATA, "<$ARGV[0]") or die "Can't open sequences file: ";
open(OUTFILE, ">>output.fasta") or die "Can't open output file: ";
while (<DATA>) {
my $line = $_;
if ($line =~ /^>?(\w+\|.+\n)/) {
$id = ;
$id =~ tr/|/_/;
}
print OUTFILE $line if defined $ids{$id};
}
close(DATA);
close(OUTFILE);
那么,您基本上是想将文件 1 拆分成多个文件,每个文件只包含一个序列?也许这个丑陋的(可能关闭每个文件句柄,而不仅仅是最后一个......)Perl 片段会帮助你。
#!/usr/bin/perl -w
use strict;
while (<>) {
my $line = $_;
if ($line =~ />?([0-9]+|\w+)\//) {
my $file_name = ;
open(OUTFILE, ">>$file_name");
}
print OUTFILE $line;
}
close(OUTFILE);
编辑: 如果您想从第二个文件中输入 ID,只需添加另一个 if
,使其与来自 headers 的 ID 相匹配第一个文件。只要两个文件的顺序相同,这应该就可以工作,因为输入在写入时会立即被丢弃并且不能再被搜索。