perl常用数据分析

这次讲讲perl里跟模式匹配或叫正则表达式有关的东西。

比如说，给出一个序列文件，里面都是Fasta格式的序列。 然后序列里面有一些NNNNNN的连续字符。

问题就是要得出这些NNNN的一段字符在该序列的具体位置。（就是匹配某字符串）
例子文件：seq.fasta
>CMM_00532
CGCGCGCTGTGCTACGCAGGCCTCTTCCAGGCCCATCTCCCGGCGGCGTGCACCACTACC
AGGATGGTGTGCGTGGGCGGGGGCGCCGCCGAGCTGGTCGCCTTTGCCAGCTTCTTGGGC
GACGACGACGACGACGACGGGGCGCACAGCAGCAGGCGCGGGGAGCTGACGCTNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGGTCGACGCGGCGCCGTGGGAGGGCGT
CGCGGACACGGTGCTGCGCGCGCTCACGACGCCGCTGCCGCTGTCCCCAGTCCGGAGCAG
>CMM_00589
ACGGGCGTGTTCCTGGCGTACGGCGGCAGCGACGATGCGCTGCCGGAGGCGGGCCTCGCG
GTGCGCATGAACGACGGGCCTTCGGGCCCTGCGTTTTGGCCGCAGCCGCGCCTGCGGCTC
ATGGAGATGCTGCTGCCGTACCTCGACCAGCACCGCTTCGCGGCCGGCGATATNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGCCTGGCATCTGCGGCAGCAGTGGGAC
GTGCCGCGGACGCACGCGTACTACGTGCCGCCCGGCGCCGTGCGGACGGCCGCGCCGCTG
CTGCTCATGGCGGCGACGCGCGACCCCGTGACGCCGTATGCGGCGGCGCGCGCGGCGCTC
>CMM_00662
GCCGTACTCTCCCAGAACGACTTGGCCTCTGCCCGTACCCTCTTTAAAGACAACCTCAAC
CTGACGCCCTATATTGCCTCGACCGAGTGCAGCGGCGTGTGGGCGCGCCGNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGGGCCGCCGGACGAGGAGCGCGGCATGGTC
GAGGTCGGGTACGGGATCGACCCGGCGTGCCGGCGGCGGGGACACGCGCGGGCGGCGCTG
>CMM_00942
CTCAACCTGCGCGACGCCGGCGCCGTGGCGGGCAGCGCGATCCCCGCCGGGCGCGTGTAC
CGCTGCGGCACGCTCGAGTACGCGGCCGCCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNCTCGCCGACTTTGCCGAGGCCGGGGCGTCGCCGCCTGGCGCACGCAGTAC
CTCCACGTCGCGCTGGCGTATGCGCCCACGTTCCGCGCCGTCCTGGAGCATGTGCGCGAC只取了一小部分数据。

分析：
◦1，CMM_开头的是该序列的ID
2，序列有几行。所以要去掉换行符，先变成一行。这样才能得出正确的位置。
代码例子：temp.pl
$/ = ">";

print "id/tstart/tend/tlen/n"; #先输出Title
while(<>){
if($_ =~ /(CMM.*?)/s(.*)>/ms){ #第一个括号匹配ID。第二个匹配序列
$id = $1;
$seq = $2;
$seq =~ s//s//g; #把序列里的换行去掉。变成一列
while ($seq =~ m/(N+)/g) { #匹配一个N或以上的字符
$len = length($1); #返回这段匹配的长度
$end = pos($seq); #用pos函数返回该匹配的终止位置
$start = $end - $len + 1; #计算出起始位置
print "$id/t$start/t$end/t$len/n"; #输出结果
}
}
}运行 perl temp.pl seq.fasta >output输出output文件。

~完



#######################################

如果不懂fasta文件，可以再看一下解释：
或是查看: Fasta格式的详细说明

>cel-miR-1 MIMAT0000003 Caenorhabditis elegans miR-1
UGGAAUGUAAAGAAGUAUGUA
>cel-miR-2 MIMAT0000004 Caenorhabditis elegans miR-2
UAUCACAGCCAGCUUUGAUGUGCUAUCACAGCCAGCUUUG
UAUCACAGCCAGCUUUGAUGUGC
……其中标识符就是大于号'>'。按两个为一组分成若干个文件。大意上是这样。

分割fasta文件
#!/usr/bin/perl

#fasta2.pl
#Usage:perl fasta2.pl in_fasta out_file

open(IN,"<$ARGV[0]");
$i = 0;
$j = 1;
while(<IN>){
if(/^>/){
$i++;
}
if($i == 3){
$j++;
$i = 1;
}
open(OUT,">>$ARGV[1]_$j");
print OUT $_;
}注： in_fasta是指要处理的fasta文件。

out_file是指输出的文件。（如命名为out, 则生成的文件名为out_1, out_2, out_3等）

主要是利用循环嘛，第一步是按大于号'>'来统计个数。再用$j来循环输出文件名

###########################################

这里简单介绍一下用Perl来实现抓好取网页的源代码，以及用POST的方法来提交表格，并返回结果。难的讲不来，讲讲简单的。

这里讲到的Perl模块有：

use LWP::Simple;use LWP::UserAgent;用perldoc查看详细的用法。
1，用perl抓取网页
如果只是要拿到某个网页，那使用 LWP::Simple 里的函数是最简单的。通过调用 get($url) 函数，就可以得到相关网址的内容。

my $url = 'http://freshair.npr.org/dayFA.cfm?todayDate=current'

use LWP::Simple;
my $content = get $url;
die "Couldn't get $url" unless defined $content;

# $content 里是网页内容，下面是对此内容作些分析：

if($content =~ m/jazz/i) {
print "They're talking about jazz today on Fresh Air!/n";
} else {
print "Fresh Air is apparently jazzless today./n";
}非常简单易懂。拿网页内容是容易的，难的是用正则过滤需要的内容。

2，通过 POST提交表格
部分HTML表格使用HTML POST 向服务器提交数据，在这里你可以这样：

$response = $browser->post( $url,
[
formkey1 => value1,
formkey2 => value2,
...
],
);实例分析：例如在(http://www.enzim.hu/hmmtop/html/submit.html)提交一段序列并返回结果，用perl来实现。代码如下：

#!/usr/bin/perl

use LWP::UserAgent;
my $browser = LWP::UserAgent->new;
$protein = "MSSSTPFDPYALSEHDEERPQNVQSKSRTAELQAEIDDTVGIMRDNINKVAERGERLTSI";
my $SUSUI_URL = "http://www.enzim.hu/hmmtop/server/hmmtop.cgi";
my $response = $browser->post( $SUSUI_URL, [ 'if' => $protein, ] );

if ($response->is_success) {
print $response->content;
} else {
print "Bad luck this time/n";
}通过分析http://www.enzim.hu/hmmtop/html/submit.html的页面可知，这个要提交的input只有一个，就是name="if"。$protein就是要提交的序列。$response->content就是返回结果



###########################

两种方法查看文件的行数

对于我们所操作的文件或是数据，行数是一个最常用的值。最后的统计结果当中，这个行数也是差不多作为一个必需项出现的，因为行数在大部分情况下，就是代表着总数。

我每天工作都是在接触两种系统：XP和Linux。所以介绍两种我常用的计算行数的方法，Excel的方法及linux命令的方法。
1，用Excel查看行数
打开Excel，在右下角的地方点右键，有“平均值”、“计数”、“计数值”、“最大值”、“求和”等。选择“计数”，因为“计算”是最常用的。这里的“计数”可以计算行数，也可以计算列数。计数的范围就是鼠标选择的范围。

“计数”与“计数值”的区别：“计数值”是指计算仅是数字的数值，不包括字符串。“计数”就没有这限制。

缺点：只有一个值时无法计算。这在数据大时会引起一些问题。例如误删等情况。需要特别留意。具体看图：

只有一个值时无法计数

2，用Linux的wc命令
在Linux下用wc进行计数。返回文件的行数、字数、字节数等。

看个例子：

wc wc1.txt
3 5 16 wc1.txt
输出信息依次是：行数 字数 字节数 文件名称。再具体点，单个统计。

wc -m filename：显示一个文件的字符数
wc -l filename：显示一个文件的行数
wc -L filename：显示一个文件中的最长行的长度
wc -w filename：显示一个文件的字数需要留意的：貌似wc统计的行算是用换行符来确定的。就是说最后一行要有换行符，最后wc的行数才是正确的，否则将会少一行。

为了说明这个问题，看一个perl的测试：

perl -e 'print "a"'|wc
0 1 1
perl -e 'print "a/n"'|wc
1 1 2够清楚了吧



############################

Linux下大文件的排序和去重复

Linux下我们用 sort 与 uniq 的命令来实现去重复行。

去重复行
简单的用法如下，如一个文件名：happybirthday.txt

cat happybirthday.txt (显示文件内容)

Happy Birthday to You!
Happy Birthday to You!
Happy Birthday Dear Tux!
Happy Birthday to You!

cat happybirthday.txt|sort （排序）

Happy Birthday Dear Tux!
Happy Birthday to You!
Happy Birthday to You!
Happy Birthday to You!

cat happybirthday.txt|sort|uniq (去重复行)

Happy Birthday Dear Tux!
Happy Birthday to You!
去大文件重复行
但有时碰到一个大文件时（例如G级的文件），用上面的命令时报错，提示空间不足。我尝试了一下，最后是用 split 命令把大文件分割为几个小文件，单独排完序后再合并 uniq 。

split -b 200m happybirthday.big Prefix_

用-b参数切割happybirthday.big，小文件为200M。切割后的文件名前缀是Prefix_切割后的文件名如

Prefix_aa

Prefix_ab再分别sort

sort Prefix_aa >Prefix_aa.sort

sort Prefix_ab >Prefix_ab.sort再用 sort -m合并，再 uniq

cat Prefix_aa.sort Prefix_ab.sort |sort -m |uniq这是好早前碰到的一个问题了。没记错的话应该是这么回事。~

sort 与 uniq 命令还有许多有用的参数，如sort -m、uniq -u、uniq -d等。sort 与 uniq的组合是很强大的。

~完

################################

这里介绍两种办法，去掉重复的数据。说之前来复习一下我喜欢的一句话：柳城博客(Lc.)， 努力在数据的海洋里畅游。

1，用Excle，适合不算太大量的数据
如果是用Excle，太大的数据打开会有问题的。打开十几M的大小的Excle都够吃力的。如果电脑内存差些，那更加惨。不过，这种情况是适合大部分人的。





2，用Linux，sort与uniq命令
假设数据放在一个文件，取名file.txt。

cat file.txt | sort | uniq >newfile.txt这样就是去掉重复数据，并输出到一个新的文件newfile.txt

简单吧。