-
Notifications
You must be signed in to change notification settings - Fork 0
/
reverse_ATGC_TACG
52 lines (46 loc) · 2.1 KB
/
reverse_ATGC_TACG
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# usr/bin/perl -w
use strict;
unless (@ARGV==2) { #@ARGV传给脚本的命令行参数列表
die"Usage: perl $0 <input:fasta> <output:ID,len,%GC,seq>\n"; #当命令行参数不是2的时候输出使用说明
}
my ($infile,$outfile) = @ARGV; #把命令行参数赋值给输入文件和输出文件
open(IN,"<$infile") || die "error: can't open infile: $infile"; #打开输入文件句柄IN
open(OUT,">$outfile") || die "error: can't make outfile: $outfile"; #打开输出文件句柄OUT
my $line=<IN>; #读取第一行,>cog1
if($line=~m/^>/) #~m//与^>(\S+)进行匹配,^>表示匹配以>开头的
{
chomp($line); #去除末尾换行符
print OUT "$line\t"; #输出第一行,>cog1
}
my $seq; #申明全局变量:$seq
while($line=<IN>) #while为一行一行遍历输入信息
{
if($line !~ /^>/) #不以>开头
{
$line=~s/[\n\s+]//g; #对于序列,去除换行符,[]表示匹配符合[]内的内容,然后~s/a/b/用b替换a
$seq.=$line; #$seq=$seq.$line 起连接的作用
}elsif($line=~/^>/) #判断是否匹配以>开头
{
chomp($line); #去除末尾换行符
my $GC = $seq =~ tr/GC/GC/; #计算G或C碱基个数
my $len=length($seq); #获取$seq的长度
my $gc_cont = $len ? $GC / $len : 0; #计算GC含量,如果长度为0,GC含量算0
my $seq_reverse=reverse $seq; #将序列反向
$seq_reverse=~ tr/TACG/ATGC/; #将序列进行互补
#输出结果到输出文件
print OUT $len;
print OUT "\t$gc_cont";
print OUT "\n$seq_reverse";
$seq=""; #每次取完序列长度后清空变量
print OUT "\n$line\t";
}
}
#最后一个序列后没有>号了要单独输出
my $GC = $seq =~ tr/GC/GC/;
my $len=length($seq);
my $gc_cont = $len ? $GC / $len : 0;
my $seq_reverse=reverse $seq;
$seq_reverse=~ tr/TACG/ATGC/;
print OUT "$len\t$gc_cont\n$seq_reverse";
close IN;
close OUT;