2013-05-01 72 views
0

我试图从Ensembl FASTA文件中找到蛋白质图案。我已经完成了大部分脚本,比如检索序列ID和序列本身,但是我收到了一些有趣的结果。无法从Emsembl FASTA删除换行符

#!/usr/bin/perl 
use strict; 
use warnings; 
use autodie; 

my $motif1 = qr/(HE(\D)(\D)H(\D{18})E)/x; 
my $motif2 = qr/(AMEN)/x; 
my $input; 
my $output; 
my $count_total  = 0; 
my $count_processed = 0; 
my $total_run  = 0; 
my $id; 
my $seq; 
my $motif1_count = 0; 
my $motif2_count = 0; 
my $motifboth_count = 0; 

############################################################################################################################ 
# FILEHANDLING - INPUT/OUTPUT 
# User input prompting and handling 
print "**********************************************************\n"; 
print "Question 3\n"; 
print "**********************************************************\n"; 

#opens the user input file previously assigned to varible to new variable or kills script. 
open my $fh, '<', "chr2.txt" || die "Error! Cannot open file:$!\n"; 

#Opens and creates output file previously assigned to variable to new variable or kills script 
#open(RESULTS, '>', $output)||die "Error! Cannot create output file:$!\n"; 

# FILE and DATA PROCESSING 
############################################################################################################################ 

while (<$fh>) { 

    if (/^>(\S+)/) { 
     $count_total = ++$count_total; # Plus one to count 
     find_motifs($id, $seq) if $seq; # Passing to subroutine 
     $id = substr($1, 0, 15);   # Taking only the first 16 characters for the id 
     $seq = ''; 
    } 
    else { 
     chomp; 
     $seq .= $_; 
    } 
} 

print "Total proteins: $count_total \n"; 
print "Proteins with both motifs: $motifboth_count \n"; 
print "Proteins with motif 1: $motif1_count \n"; 
print "Proteins with motif 2: $motif2_count \n"; 

exit; 

###################################################################################################################################### 
# SUBROUTINES 
# 
# Takes passed variables from special array 
# Finds the position of motif within seq 
# Checks for motif 1 presence and if found, checks for motif 2. If not found, prints motif 1 results 
# If no motif 1, checks for motif 2 

sub find_motifs { 
    my ($id, $seq) = @_; 
    if ($seq =~ $motif1) { 
     my $motif_position = index $seq, $1; 
     my $motif = $1; 
     if ($seq =~ $motif2) { 
      $motif1_count = ++$motif1_count; 
      $motif2_count = ++$motif2_count; 
      $motifboth_count = ++$motifboth_count; 
      print "$id, $motif_position, \n$motif \n"; 
     } 
     else { 
      $motif1_count = ++$motif1_count; 
      print "$id, $motif_position,\n $motif\n\n"; 
     } 
    } 
    elsif ($seq =~ $motif2) { 
     $motif2_count = ++$motif2_count; 
    } 
} 

正在发生的事情是,如果主题是在一个数据线和下一个的开始结束发现,它会返回母题与数据的换行符。这种篡改数据的方法之前运行良好。

样品结果:

ENSG00000119013, 6, HEHGHHKMELPDYRQWKIEGTPLE (CORRECT!) 

ENSG00000142327, 123, HEVAHSWFGNAVTNATWEEMWLSE (CORRECT!) 

ENSG00000151694, 410, **AECAPNEFGAEHDPDGL** 

这就是问题所在。该主题的比赛,但返回上半年,换行符,然后打印下半年在同一行,以及(这是更大的问题的症状 - 摆脱换行的!)

Total proteins: 13653 
Proteins with both motifs: 1 
Proteins with motif 1: 12 
Proteins with motif 2: 22 

我已经尝试了不同的方法,如@seq =~ s/\r//g或`s \ \ n // g并在脚本中的不同位置。

回答

1

从描述中不清楚,但“在同一行上打印下半部分”听起来像您的输出覆盖在本身上,因为它在最后具有回车符。

如果您在Linux系统上运行,并且您的chomp线路来自Windows,则会发生这种情况。

您应该将chomp替换为s/\s+\z//,这将删除所有尾随的空白。而且由于回车和换行都算作“空白”,它将删除所有可能的终止字符。

顺便提一下,您误导了++运营商的目的。它也修改它应用于的变量的内容,所以你需要的只是++$motif1_count等。你的代码的工作原样,因为运算符也返回增量变量的值,所以$motif1_count = ++$motif1_count首先增加变量,然后分配它自己。

此外,你在你的正则表达式中使用\D。你是否知道这个匹配非数字个字符?这似乎是一个非常模糊的分类有用。