根据类似的行标题拆分文件

是否有人可以帮助我为awk，grep，sed，perl或python中的以下要求编写脚本？根据类似的行标题拆分文件

输入文件“raw.fa”：

>CLocus_1_Sample_61_Locus_1_Allele_0 [JPKM01095229.1, 31450, +] 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>CLocus_1_Sample_67_Locus_1_Allele_0 [JPKM01095229.1, 31450, +] 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>CLocus_1_Sample_107_Locus_1_Allele_0 [JPKM01095229.1, 31450, +] 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>CLocus_1_Sample_107_Locus_1_Allele_1 [JPKM01095229.1, 31450, +] 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGGTTGAAG 
>CLocus_41_Sample_158_Locus_53_Allele_0 [JPKM01105094.1, 1700, +] 
TGCAGGTTATCCAGCTCTATTCTGCACTGGCCATCGTACCAAATAGCAGGAGGGT 
>CLocus_41_Sample_159_Locus_31_Allele_0 [JPKM01105094.1, 1700, +] 
TGCAGGTTATCCAGCTCTATTCTGCACTGGCCATCGTACCAAATAGCAGGAGGGT 
>CLocus_86_Sample_161_Locus_103_Allele_0 [JPKM01106288.1, 770, -] 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACACACTTGGCTCCCATTGGGATGACTCCTTT 
>CLocus_86_Sample_164_Locus_98_Allele_0 [JPKM01106288.1, 770, -] 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACACACTTGGCTCCCATTGGGATGACTCCTTT 
>CLocus_86_Sample_166_Locus_110_Allele_0 [JPKM01106288.1, 770, -] 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACTCACTTGGCTCCCATTGGGATGACTCCTTT 
>CLocus_86_Sample_167_Locus_123_Allele_0 [JPKM01106288.1, 770, -] 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACTCACTTGGCTCCCATTGGGATGACTCCTTT

我想通过轨迹分割上述文件，每个基因座1个文件，保持DNA（第二行）和样品＃从第一行，三个产分别.fa文件：

“locus1.fa”：

>Sample_61 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>Sample_67 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>Sample_107 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGATTGGAG 
>Sample_107 
TGCAGGTGTGTTCTGCAGATCCAAACACAAAGAGGCAGGGGTTGAAG

“locus41.fa”：

>Sample_158 
TGCAGGTTATCCAGCTCTATTCTGCACTGGCCATCGTACCAAATAGCAGGAGGGT 
>Sample_159 
TGCAGGTTATCCAGCTCTATTCTGCACTGGCCATCGTACCAAATAGCAGGAGGGT

“locus86.fa”：

>Sample_161 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACACACTTGGCTCCCATTGGGATGACTCCTTT 
>Sample_164 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACACACTTGGCTCCCATTGGGATGACTCCTTT 
>Sample_166 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACTCACTTGGCTCCCATTGGGATGACTCCTTT 
>Sample_167 
TGCAGGGAACCGTGCTCAGCTCTGGAGTATTCCCACTCACTTGGCTCCCATTGGGATGACTCCTTT

感谢您的帮助！我发现awk代码首次出现分裂，但不是如何拆分相似行的分组（例如，所有具有locus86头部和第二行DNA序列的行）。

克里斯·马丁

来源

2014-10-11 Chris Martin

你有任何编程语言编写的任何代码来尝试自行解决问题呢？ – 2014-10-11 21:04:09

如果你有一小部分的 “轨迹” 的价值观，你可以用grep他们：

grep -A 1 locus86 raw.fa > locus86.tmp

然后重新格式化行：

sed 's/>.*(Sample .*)_Locus_.*/\1/' locus86.tmp > locus86.fa

（显然你可以将它们与管道结合起来，而不是使用中间文件）。

如果你有一个更大或者未知的基因座值集合，像perl这样的脚本是有意义的。需要提醒的是它可以打开/关闭文件贵，这里的一些伪代码：

open(IN, "raw.fa"); 
my $OUT = undef; 
while (<IN>) { 
    if (/>/) { 
     my ($sample, $locus) = ($_ ~= /.*\(Sample_.*\)_\(Locus_.*\)_Allele/); 
     if (defined($OUT)) { 
      close($OUT); 
     } 
     open($OUT, "$locus.fa"); 
     print $OUT ">$sample\n"; 
    } 
    else { 
     print $OUT $_; 
    } 
}

来源

2014-10-11 22:17:12

感谢您的快速响应！我得到这个错误在Ubuntu中运行你的脚本（我也添加了Perl的bash，并改为我的文件名而不是“raw.fa”）：“语法错误在第7行附近”$ _〜“ – 2014-10-11 23:29:06

它应该是'=〜 '而不是'〜='。 – 2014-10-11 23:33:18

再次感谢您的快速帮助!!! – 2014-10-11 23:38:31

#!/usr/bin/perl 
$file = shift; 
my $locusfile; 
open FILE, $file; 
while(<FILE>) 
{ 
    chomp($_); 
    if (/CLocus_([0-9]*)_Sample_([0-9]*)_.*/) 
    { 
    $locusfile="locus" . $1 . ".fa"; 
    $cmd = qq{echo ">Sample_$2" >> $locusfile}; 
    } 
    else 
    { 
    $cmd = qq{echo "$_" >> $locusfile}; 
    } 
    system($cmd); 
} 
close(FILE);

我无法在线 - 年初追加了“>”角色得到了很多的问题当插入到system命令...所以它不是一个完整的答案，但它做你想做的。

[固定]

干杯

来源

2014-10-11 22:32:14

我会承认@Alain，你的回答比较好，但是我使用system（）取而代之获得奖励积分...;） – 2014-10-11 22:34:26

哇 - 工作得很好 - 感谢快速回复！（我作为一个完整的perl新手唯一的挂机是知道用引号中的文件名替换“shift”）。我真的很感激它 - 你已经保存了这个博士后！ – 2014-10-11 23:27:42

如果你必须使用'system'，那么至少得到正确的回应。 ''cmd = qq {echo“> Sample_ $ 2”>> $ locusfile};''和'$ cmd = qq {echo“$ _”>> $ locusfile};'并且将'system'调用放在'if/else'。我不认为使用'system'是个不错的选择，但它确实有效。除此之外，它很昂贵 - 它需要为每一行输入启动一个完整的外壳。输入量为50行时不成问题，但如果有5000万行，则不会出现问题。 – 2014-10-11 23:30:46

根据类似的行标题拆分文件

回答

相关问题