注册 登录
编程论坛 Perl论坛

[求助]如何用perl计算密码子不同位点的GC含量?

h17785687955 发布于 2019-11-14 18:25, 3068 次点击
如何用perl计算密码子不同位点的GC含量???急用!!!!!
我是想计算每一条CDS密码子三个位点的GC含量,用下面这个脚本但没有数据输出,我的文本是fasta格式的,文件名是M.fasta,fasta结构的文本是以下形式:
>a
ATG GAA GAA GAG AAC GCG CAA CTC CGT ACC GAA TTG
>b
ATG GTT TTG TCA CAA GAA AAT GTT GAT CAC AAT GAT
“  $li=~s/^\w+://;        # 去除每条序列前的序列编号和冒号”怎么修改?
脚本:
#!/usr/bin/perl

sub gc {
  my $x;
  if($_[0] eq "G" or $_[0] eq "C" {
    $x=1;
  }
  $x;
}

my $input=$ARGV[0];
my ($li,$i,$len,$s,$tm,@t,@g);

open TMP, $M.fasta|| die;

while($li=<TMP> {
  my @t;
          my @g;
  chomp $li;
    #print "$li\n";
  $li=~s/^\w+://;        # 去除每条序列前的序列编号和冒号
  $li=~s/\s//g;
  #print "$li\n";
  $len=length $li;
  $i=0;
  while($i<$len) {
    $s=substr($li,$i,1);
        $y=$i%3;        # 判断是三个位点中的哪个
        $t[$y]++;        # 该位点的总碱基数目
        $g[$y]+=gc($s);        # 该位点的GC碱基数目
        $i++;
  }
  $tm=$g[0]/$t[0]*100;
printf "The GC content in site1 is %.3f\%\n",$tm;
$tm=$g[1]/$t[1]*100;
printf "The GC content in site2 is %.3f\%\n",$tm;
$tm=$g[2]/$t[2]*100;
printf "The GC content in site3 is %.3f\%\n",$tm;
}
close TMP;
3 回复
#2
fall_bernana2019-11-29 17:50
回复 楼主 h17785687955
程序代码:
#!/usr/bin/perl

sub gc {
  my $x;
  if($_[0] eq "G" or $_[0] eq "C") {
    $x=1;
  }
  $x;
}

my $input=$ARGV[0];
my ($li,$i,$len,$s,$tm,@t,@g);



while($li=<DATA>) {
  my @t;
    my @g;
  chomp $li;
    print "$li\n";
  $li=~s/^\w+://;        # 去除每条序列前的序列编号和冒号
  $li=~s/\s//g;
  #print "$li\n";
  $len=length $li;
  $i=0;
  while($i<$len) {
    $s=substr($li,$i,1);
        $y=$i%3;        # 判断是三个位点中的哪个
        $t[$y]++;        # 该位点的总碱基数目
        $g[$y]+=gc($s);        # 该位点的GC碱基数目
        $i++;
  }
  $tm=$g[0]/$t[0]*100;
    printf "The GC content in site1 is %.3f\%\n",$tm;
    $tm=$g[1]/$t[1]*100;
    printf "The GC content in site2 is %.3f\%\n",$tm;
    $tm=$g[2]/$t[2]*100;
    printf "The GC content in site3 is %.3f\%\n",$tm;
}


__DATA__
ATG GAA GAA GAG AAC GCG CAA CTC CGT ACC GAA TTG
ATG GTT TTG TCA CAA GAA AAT GTT GAT CAC AAT GAT
#3
TOLSWRF2021-09-23 11:34
回复 2楼 fall_bernana
大神这个编码,只能两行两行的进行运算,有没有多行运算的
#4
fall_bernana2021-11-02 11:22
回复 3楼 TOLSWRF
没懂。什么叫2行运算,什么叫多行运算
1