perl实战练习-计算fastq文件中每条序列的GC含量

计算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/）来计算。