perl 小练习

最新推荐文章于 2022-11-27 20:29:36 发布

爱笑的小牙

最新推荐文章于 2022-11-27 20:29:36 发布

阅读量1k

点赞数 1

分类专栏： perl

本文链接：https://blog.csdn.net/Cassiel60/article/details/104812804

版权

本文介绍了一个使用Perl进行生物信息学分析的练习，包括计算scaffold.fasta文件的GC含量、N50和N90，以及从基因组scaffold和gff文件中提取CDs区域并以特定格式输出，最后将CDs翻译成蛋白质序列并按行输出。

摘要由CSDN通过智能技术生成

1.以scaffold.fasta作为输入文件，计算GC含量以及N50和N90

#基因组装中N50和N90分别表示把片段按照从大到小的顺序排列，然后累加，分别超过总长度的50%和90%那条序列的长度

2.根据给定的基因组scaffold.fasta文件和相对用的基因注释gff文件提取基因的cds区域，并以每行60个碱基的格式输出到cds.fasta文件中

3.以cds.fasta作为输出文件，将其翻译成蛋白质序列并以每行60个氨基酸的格式输出到pep.fasta文件

#!usr/bin/perl -w
use strict;

if (@ARGV<1){
	die "Usage:perl $0 <genome.fa>\n";
}

my $input = shift @ARGV; #表示取输入数字的第一个 也可以使用$input =$ARGV[0] 

my ($sum,$G_num,$C_num,$N_num) = (0,0,0,0);
my %seq = ();
my $id;
open IN,"<$input" or die $!;
while (my $line= <IN>){
	chomp $line;
	if ($line =~ />([^\s]+)/){       #查找>开头的作为哈希的键，一个()表示$1
		$id = $1;
	}else {
		$sum +=length($line);
		$seq{$id} += length($line);
		$G_num += ($line =~ tr/Gg/Gg/);  #用str返回替换次数，表示Gg个数
		$C_num += ($line =~ tr/Cc/Cc/);
		$N_num +