专业背景：计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学
2010-06-24 16:35:49| 分类： perl&bioperl | 标签： |举报 |字号大中小 订阅
在做了大量的测序之后，测序公司给你的可能只是一个又一个的序列文件，而且在大多数情况下这些序列文件的文件名还含有测序公司的一些信息，而这些信息对你来说毫无用处。那么如何才能把有用的信息抽提出来，如何才能把这些测序文件快速地转变成一个fasta文件呢？这正是我攥写本帖所要解决的问题。
我结合一个具体的例子来分享一下我用perl来处理测序文件的方法。测序公司发送给我的文件的如下：
L0408022_taG6_PEGFP-N.txt
L0408022_taG7_PEGFP-N.txt
L0408023_taG8_PEGFP-N.txt
L0408023_taG9_PEGFP-N.txt
L0408025_taY1_PEGFP-N.txt
L0408025_taY2_PEGFP-N.txt
本例字中的所有文件均在附件上，有兴趣的战友可以下载下来自己操作一遍。
我希望得到的结果是：
>taG6
GCTCCTCGCCCTTGCTCACCATGGATCCCGGGCCCGCGGTACCAACTATTGT
>taG7
GCTCCTCGCCCTTGCTCACCATGGATCCCGGGCCCGCGGTACC
……那么如何快速实现？你先下载我的附件按照我的指示一步一步地操作就可以了：
1，安装perl，在windows下，去下载activePerl直接安装就可以。大多数的Linux好像都安装了perl。这一步可以省略。如果没装请自己解决，linux用户解决这一点小问题一定是没有问题的，就不多说什么了。
２，解压附件，到你的文件夹。
３，运行tyhy2fasta.pl脚本文件。Window用户直接双击就可以了，linux用户在终端中运行
perl tyhy2fasta.pl
４，运行完之后你可以看到有一个名字为seqconversed.seq的文件，打开看看。
里面的fasta格式的序列这是我想要的结果。
接下来要实战性地解决你的测序文件问题。用文本编辑器打开tyhy2fasta.pl文件，你可以看到以下代码：
#! /usr/bin/perl;
open(OUTFILE,">seqconversed.seq");
while(defined($nextname=<*.txt> )){
open(FILE,$nextname);
@line=<FILE>;
$seq=join("",@line);
$seq=~s/\W//g;
$fn0=$nextname;
$pst1=index($fn0,"_",5);
#print "$pst1\n";
$pst2=index($fn0,"_",$pst1+1);
#print $pst2,"\n";
$fn=substr($fn0,$pst1+1,$pst2-$pst1-1);
print ">$fn\n$seq\n";
print OUTFILE ">$fn\n$seq\n";
close(FILE);
}
close(OUTFILE);接下来我来解释一下每一行代码的意思，以及如何修改代码使之为你工作。
#! /usr/bin/perl;
#这是文件头，Windows用户可以不管它是什么意思
open(OUTFILE,">seqconversed.seq");
#新建一个命名为seqconversed.seq文件，用来储存fasta格式的序列。
#这个文件在内存里暂时性地把它记忆为OUTFILE。这些名字都是可以随便该的。
while(defined($nextname=<*.txt> )){
#逐个读取当前文件夹中有扩展名为.txt的文件名，把它赋值给$nextname。
#如果测序公司给你的文件是.seq,那么你要把代码中的“.txt“替换成”.seq"。
#注意：如果测序公司给的是.seq文件上一行代码中的seqconversed.seq要改成seqconversed.txt。
#总之不能和测序公司给的文件有相同扩展名。
#注意：当前文件夹最好不要有非序列文件，例如在本例中你不能把“情书.txt“和你的序列放在一起。
open(FILE,$nextname); #打开$nextname文件，把内容存到FILE中去。
@line=<FILE>; #把FILE的每一行作为一个元素，赋值给数组@line
$seq=join("",@line); #把数组@line的所有元素合并，合并后变成一个字符串，并赋值给$seq
$seq=~s/\W//g; #把字符串中杂七杂八的非字母符号如段落符号等删除（实际上是把他们替换为无）
$fn0=$nextname; #序列的文件名赋值给$fn0(意为原始filename)
$pst1=index($fn0,"_",5);
#在原始文件名（本例第一个文件名为L0408022_taG6_PEGFP-N.txt）中查找“_”的位置（position)，
#从第5位开始查起。找到"_"之后把它的位置赋值给$pst1
#注意perl的第一位是第0位，第二位才是第1位，本例中的第5位，"_"的位置是第8位。
#所以$pst1此时的值为8。
#print "$pst1\n";
#把$pst1打印在屏幕上，看看是否计算正确。
#当然要执行这一步要先把前面的这一行最前面的”#“好删除掉，否则这一行是不工作的。
$pst2=index($fn0,"_",$pst1+1); #此行与同前两行类
#print $pst2,"\n";#此行与同前两行类
$fn=substr($fn0,$pst1+1,$pst2-$pst1-1);
#此行是关键！
#它表示从原始文件名（本例第一个文件名为L0408022_taG6_PEGFP-N.txt）的
#第$pst1+1位（在本例中即第９位）开始
#抽出$pst2-$pst1-1（本例中为13-8-1=4)个字符，赋值给$fn
print ">$fn\n$seq\n";#把用fasta格式打印在屏幕上，看看是否正确
print OUTFILE ">$fn\n$seq\n"; #用fasta格式打印到OUTFILE 里面去
close(FILE);#关闭读取的第一个序列文件
}#进入下一个序列文件
close(OUTFILE);#关闭OUTFILE文件。如果你的序列文件是H96-yourname-pcdna3.1-.ab1.seq，那么你就要把上面的部分代码改为：
open(OUTFILE,">seqconversed.txt");
.
.
$pst1=index($fn0,"-",1);
#print "$pst1\n";
$pst2=index($fn0,"-",$pst1+1);
#print $pst2,"\n";
$fn=substr($fn0,$pst1+1,$pst2-$pst1-1);
.
.只要仔细认真地去看看我的注释，你会发现perl原来如此简单和如此亲近……
Tag: perl入门，perl实用技术，dna测序，序列分析，测序文件，转换fasta格式
2fasta.rar (3.39k)
（By tanklao）http://www.dxy.cn/bbs/post/view?bid=73&id=14484479&sty=1&tpg=1&age=0
本文详细出处参考：http://liucheng.name/477/
评论