Learn:用Perl处理生物信息数据

# 用perl进行计算
# -e 表示直接执行
perl -e 'print 10*10,"\n";'

# 用perl读取fasta文件的序列id
# 按照顺序一行一行读取文件
# chomp; 把换行符去除
perl -ne 'chomp; my $id = $1 if(/>(\S+)/); print "$id\n" if($id);' Danio_rerio.Zv9.75.dna.toplevel.fa

# 用单行perl计算基因组的长度
# BEGIN{my $sum = 0;} 程序开始读文件前执行，定义了一个全局变量
# END{print "Total: $sum\n";} 程序完成文件读取后执行
perl -ne 'BEGIN{my $sum = 0;} chomp; next if(/^>/); $sum += length($_); END{print "Total: $sum\n";}' Danio_rerio.Zv9.75.dna.toplevel.fa

# 用同样的方法可以完成所有基因的蛋白序列平均长度统计
perl -ne 'BEGIN{my $sum = 0; my $count = 0;} chomp; $count++ if(/^>/); next if(/^>/); $sum += length($_); END{my $avg = $sum/$count; print "PEP average length: $avg\n";}' Danio_rerio.Zv9.75.pep.all.fa

这时候，想一下我们在做生物学研究时，还需要那些统计分析，可以尝试用perl完成。

perl脚本

当然还有一些更加复杂的分析，用perl脚本去实现更加方便。

以下是如何计算zebra fish的染色体长度？

# 新建一个文件名
vim seq_len.pl

vim是linux系统下的文本编程器，当然你可以用gedit。

具体请参考：

Interactive Vim tutorial

以下是seq_len.pl的代码：

#!/usr/bin/perl -w

use strict; #用严格的语法格式，这样不容易出错

my ($fa_f) = @ARGV; # 读入所要处理的fasta文件名

open IN,"$fa_f" || die $!; # 打开fasta文件，如果打开失败的话，输出标准错误信息
$/ = ">"; <IN>; $/ = "\n"; # $/ 变量是perl是存储文件句柄的分隔符
while(<IN>){
	my $id = $1 if(/\S+/); # 获取序列的id
	$/ = ">";              # 分隔符$/ 变量变成">"，可以每次读进属于id的序列
	my $seq = <IN>;
	chomp($seq);           # 去除分隔符">"
	$seq =~ s/\n//g;       # 把换行符"\n"去除
	$/ = "\n";
	print "$id\t".length($seq)."\n";  # 输出结果
}
close IN;

运行perl程序:

perl ./seq_len.pl Danio_rerio.Zv9.75.dna.toplevel.fa

Mix togethter

尝试把perl和shell结合在一起，我们要对每一个基因组都进行DNA和蛋白质序列长度进行统计。

建立并打开basic_stat.sh

vim basic_stat.sh

basic_stat.sh代码如下：

#!/bin/bash

# $1 表示输入的第一个参数，$2代表第二个参数
perl ./seq_len.pl $1.dna.toplevel.fa > $1.ChrLen.lst
perl ./seq_len.pl $1.pep.all.fa > $1.ProteinLen.lst

运行shell脚本：

# sh ./basic_stat.sh <genome_name_prefix>
sh ./basic_stat.sh Danio_rerio.Zv9.75

转载请标明出处：

http://zhoujj2013.github.io/UnCoverIt/bioinformatics/2014/07/17/perl-for-bioinformatics/

bioinformatics 10