计算FASTA文件中每条序列的长度
输入文件fasta.txt
>con1
GTTCATAACTTACCATTGACTGTTCATGATTTAAACTGTGTCTTCCTGTTCTTTTAGTTGCTTTGGATACTATAAGGTCAGAACTAGAA
>con2
GGGACCAGGTGGGAAGCTGCTGCTTTCTTTTCCCTTTTGGTCTTCTTGGGCCCACCCTTCAGCTTCTGCTTTTCTTCATCTTCTCGGTTTTGAGGCCAGGAGGCAGCCAGTATCCTGGCCGCTTCTGCTTGAGAGCTGGTCCCCTCCTCT
>con3
TCCAGAGCCAGTTCCCTGACGATGCCGAGGTTTACCGGCTCATCGAGGTGACTGGCCTTGCCAGGAGCGAGATCAAGAAGTGGTTCAGTGACCACCGATATCGGTGTCAAAGGGGCATCGTCCACATCACCAG
>con4
TTTTCCTCCAAGTGTACAAGACTGATCTGGAGAAGGACATTATTTCGGACACATCTGGTGACTTCCGCA
>con5
GTTTGCATCGTCATCCAATTTTTTTTCATATGCCCCAAACCCATTCAATTTCTGATTGTGTTGTGTTCCCTGGTGTAAGATATCTCCCAGGAGAGGGCCACACATTCCCCAGAGGTGGACCTTCTTGGTACATACACC
自己写脚本
#!/usr/bin/perl
my ($f,$out)=@ARGV;
open(DATA,$f) or die "infile 文件无法打开,$f";
open(OUT,">$out")||die $!;
$/=">";<DATA>; #设置输入记录分隔符为">",并去除第一个">"
while(<DATA>){
$line=$_;
my $id=$1 if($line=~/^(\S+)/); #获取序列ID
chomp $line; #去掉末尾的换行
$line=~s/^.+?\n//; #删除第一行
$line=~s/\s//g; #删除序列中的空白字符
my $gc = () = /G/g;
my $cc = () = /C/g;
my $len=length($line); #计算序列长度
my $gcc=($gc+$cc)/$len;
print OUT "$id\t$gcc\n"; #输出结果到输出文件中
}
close IN;
close OUT;
网上其他答案
#!/usr/bin/perl -w
use strict;
unless (@ARGV==2){ #@ARGV 传给脚本的命令行参数列表
die"Usage: perl $0 <input.fa><out.gc>\n"; #当命令行参数不是2的时候输出使用说明
}
my ($infile,$outfile)=@ARGV; #把命名行参数赋值给输入、输出文件
open IN,$infile||die"error:can't open infile:$infile"; #打开输入文件句炳IN
open OUT,">$outfile"||die $!; #打开输出文件句柄OUT
$/=">";<IN>; #设置输入记录分隔符为 ">",并去除第一个”>"
while (<IN>){ #$_=<IN>,把序列ID行和序列赋值给$_,$_=可以省略不写
my $id=$1 if(^/(\S+)/); #获取序列ID
chomp; #去掉末尾换行符
s/^.+?\n//; #删除第一行
s/\s//g; # 删除序列中的空白字符
my $GC=(tr/GC/GC/);#计算G或C碱基个数
my $AT=(tr/AT/AT/); #计算A或T碱基个数
my $len=$GC+$AT; # 计算序列非N长度
my $gc_cont = $len ? $GC/$len :0 ; #计算GC含量,如果长度为0,GC含量为0
print OUT "$id\t$gc_cont\n"; #输出结果到输出文件
}
close IN;
close OUT;
自己写的还是会有欠缺,没有考虑GC为0的情况,还分别计算了G 和C的个数,然后相加,完全可以用(tr/GC/GC/)来计算。