Perl脚本用于对多个DNA序列进行分组

问题描述 投票：1回答：2

我有一个约500个DNA序列的FASTA文件，每个DNA序列都具有我所知的单核苷酸多态性（SNP）的目标位置。

对于文件中的每个条目，我在每行上都有一个单独的制表符分隔文本文件

  1. FASTA序列名称
  2. 起始位置
  3. 结束位置
  4. SNP的位置

文本文件中的序列和位置的顺序相同。

虚拟FASTA文件是：

>AOS-94_S25_L002_R1_001_trimmed_contig_767
GACACACACTGATTGTTAGTGGTGTACAGACATTGCTTCAAACTGCA
>AOS-94_S25_L002_R1_001_trimmed_contig_2199
TAGGTTTTCTTTCCCATGTCCCCTGAATAACATGGGATTCCCTGTGACTGTGGGGACCCCTGAGAGCCTGGT
>AOS-94_S25_L002_R1_001_trimmed_contig_2585
GATAAGGAGCTCACAGCAACCCACATGAGTTGTCC

并且虚拟位置文件是

AOS-94_S25_L002_R1_001_trimmed_contig_767   5   15  10
AOS-94_S25_L002_R1_001_trimmed_contig_2199  8   19  11
AOS-94_S25_L002_R1_001_trimmed_contig_2585  4   20  18

这是我编写和试过的脚本

use warnings;
use strict; 

# Read in the complete FASTA file:
print "What is the name of the fasta contig file?\n";
my $fasta = <STDIN>;
chomp $fasta;

# Read in file of contig name, start pos, stop pos, SNP pos in tab delimited 
text:
print "Name of text file with contig name and SNP position info? \n";
my $text = <STDIN>;
chomp $text;

# Output file
print "What are we calling the output? \n";
my $out= <STDIN>;
chomp $out;

local $/ = "\n>"; #Read by fasta record

my $seq1 = (); 

open(FASTA,$fasta) || die "\n Unable to open the file!\n";
open(POS,$text) || die "\n Unable to open the file! \n";
my @fields = <POS>;
    while (my $seq = <FASTA>){
        chomp $seq;
        my @seq = split(/\n/,$seq);
            if($seq[0] =~ /^>/){
                $seq1 = $seq[0];
            }elsif($seq[0] =~ /[^>]/){ #matches any character except the >
                $seq1 = ">".$seq[0];
            }
    for my $pos (@fields){
        chomp $pos;
        my @field = split(/\t/,$pos);
    open(OUTFILE,">>$out");
    print OUTFILE "$seq1";
    my $subseq = substr $seq[1], $field[1] -1, $field[2] - $field[1]; 
    print OUTFILE "$subseq\n";
    }   
}
close FASTA;
close POS;
close OUTFILE;

这就是我的结果，这就是我想要的：

>AOS-94_S25_L002_R1_001_trimmed_contig_767
CACACTGATT
>AOS-94_S25_L002_R1_001_trimmed_contig_2199
TTTTCTTTCC
>AOS-94_S25_L002_R1_001_trimmed_contig_2585
AGGAGCTCAC

但是，我还需要在序列名称之后打印SNP位置（第4列），例如，

>AOS-94_S25_L002_R1_001_trimmed_contig_767
pos=10
CACACTGATT
>AOS-94_S25_L002_R1_001_trimmed_contig_2199
pos=11
TTTTCTTTCC
>AOS-94_S25_L002_R1_001_trimmed_contig_2585
pos=18
AGGAGCTCAC

我尝试插入print OUTFILE "pos= $field[3]\n";after print OUTFILE "$seq1";，我得到以下内容：

>AOS-94_S25_L002_R1_001_trimmed_contig_767
10
AOS-94_S25_L002_R1_001_trimmed_contig_2199
CACACTGATT
>AOS-94_S25_L002_R1_001_trimmed_contig_2199
10
AOS-94_S25_L002_R1_001_trimmed_contig_2199
TTTTCTTTCC
>AOS-94_S25_L002_R1_001_trimmed_contig_2585
10
AOS-94_S25_L002_R1_001_trimmed_contig_2199
AGGAGCTCAC

显然我搞乱了我的循环，可能还有一些chomp命令。

例如，当我print "$seq1"到一个文件时，为什么它不需要包含在打印字符串中的"\n"？字符串中必须有一个硬回车？

我知道我错过了一些关于这种结构的基础知识，但到目前为止我还无法弄清楚如何解决我的错误。有人可以提供任何建议吗？

Update

Perl代码重新格式化以便易读

use warnings;
use strict;

# Read in the complete FASTA file:
print "What is the name of the fasta contig file?\n";
my $fasta = <STDIN>;
chomp $fasta;

# Read in file of contig name, start pos, stop pos, SNP pos in tab delimited
text:
print "Name of text file with contig name and SNP position info? \n";
my $text = <STDIN>;
chomp $text;

#Output file
print "What are we calling the output? \n";
my $out = <STDIN>;
chomp $out;

local $/ = "\n>";    # Read by FASTA record

my $seq1 = ();

open( FASTA, $fasta ) || die "\n Unable to open the file!\n";
open( POS,   $text )  || die "\n Unable to open the file! \n";

my @fields = <POS>;

while ( my $seq = <FASTA> ) {

    chomp $seq;
    my @seq = split( /\n/, $seq );

    if ( $seq[0] =~ /^>/ ) {
        $seq1 = $seq[0];
    }
    elsif ( $seq[0] =~ /[^>]/ ) {    # matches any character except the >
        $seq1 = ">" . $seq[0];
    }

    for my $pos ( @fields ) {
        chomp $pos;
        my @field = split( /\t/, $pos );

        open( OUTFILE, ">>$out" );
        print OUTFILE "$seq1";

        my $subseq = substr $seq[1], $field[1] - 1, $field[2] - $field[1];
        print OUTFILE "$subseq\n";
    }
}

close FASTA;
close POS;
close OUTFILE;
perl
2个回答
2
投票

您的代码存在许多问题

  • 您的评论与代码不符。例如，当代码从STDIN接受文件名并修剪它时，你有Read in the complete FASTA file。通常最好用精心选择的标识符编写干净的代码;这样程序解释了自己
  • 您正在使用open和全局文件句柄的双参数形式。您也没有die字符串失败的原因，并且最后有一个换行符，这将阻止Perl为您提供发生错误的源文件名和行号 就像是 open( FASTA, $fasta ) || die "\n Unable to open the file!\n" 应该 open my $fasta_fh, '<', $fasta_file or die qq{Unable to open "$fasta_file" for input: $!} 和 open( OUTFILE, ">>$out" ); 应该 open my $out_fh, '>>', $output_file or die qq{Unable to open "$output_file" for appending: $!}
  • 你应该避免在变量名称周围加上引号。 print OUTFILE "$seq1" 应该 print OUTFILE $seq1
  • 您将输入记录分隔符设置为"\n>"。这意味着每次调用<FASTA> Perl都将读取该字符串的下一次出现。这也意味着chomp将从行的末尾删除该字符串，如果它在那里

最大的问题是你在从$/读取之前从未重置POS。请记住，它的设置会影响每个readline（或<>）和每个chomp。因为你的$text文件可能在一行开头没有包含>字符，你将一次读取整个文件

这就是为什么你在输出中看到换行而不需要它们的原因。你已经阅读了整个文件以及所有嵌入的换行符，并且chomp在这里没用，因为你修改了它删除的字符串

local因为某种原因被命名。它会临时和本地更改当前范围的值。但是您的“当前范围”是文件其余部分的全部内容，您正在使用修改后的终结符读取这两个文件

使用一些大括号{ ... }来限制local修改的范围。或者，因为更新版本的Perl中的文件句柄表现为IO::Handle对象，所以您可以编写

$fasta_fh->input_record_separator("\n>")

并且更改将仅适用于该文件句柄，并且根本不需要本地化$/

这是您的程序的修订版本，它还解决了一些不良的标识符选择以及其他一些问题。请注意，此代码未经测试。我目前正在火车上工作，只能查看我在写作的内容

请注意，当没有指定循环变量时，像while ( <$fasta_fh> )和for ( @pos_records )这样的东西使用默认变量$_。同样，当缺少相应的参数时，chomp和split等运算符将适用于$_。这样就不需要明确地提及任何变量，并且它会导致更简洁和可读的代码。 $_相当于英语

我鼓励你理解你写的东西实际上是做什么的。通常的做法是从互联网的一个部分复制代码并将其提供给其他地方的某些人，以使其适合您。这不是“学习编程”，除非你学习语言并把它放在心上，否则你什么都不会理解

请详细说明您的代码。我希望你能看到我对你的问题所做的编辑，以及我的解决方案中的代码，比你发布的程序更容易阅读？虽然你可以随心所欲地让自己的工作变得尴尬，但是向那些你要求免费编程帮助的陌生人世界提供这样的混乱是不公平和不礼貌的。一个很好的中间行是在按下Tab键时改变你的编辑器使用四个空格的缩进。切勿在源代码中使用制表符！

use strict;
use warnings 'all';

print "Name of the FASTA contig file: ";
chomp( my $fasta_file = <STDIN> );

print "Name file with SNP position info: ";
chomp( my $pos_file = <STDIN> );

print "Name of the output file: ";
chomp( my $out_file = <STDIN> );

open my $out_fh, '>', $out_file die qq{Unable to open "$out_file" for output: $!};

my @pos_records = do {
    open $pos_, '<', $pos_file or die qq{Unable to open "$pos_file" for input: $!};
    <$pos_fh>;
};
chomp @pos_records; # Remove all newlines

{
    open my $fasta_fh, '<', $fasta_file or die qq{Unable to open "$fasta_file" for input: $!};

    local $/ = "\n>"; # Reading FASTA format now

    while ( <$fasta_fh> ) {

        chomp;    # Remove "">\n" from the end

        my ( $header, $seq ) = split /\n/; # Separate the two lines

        $header =~ s/^>?/>/; # Replace any chomped >

        for ( @pos_records ) {

            my ( $name, $beg, $end, $pos ) = split /\t/;
            my $subseq = substr $seq, $beg-1, $end-$beg;

            print $out_fh "$header\n";
            print $out_fh "pos=$pos\n";
            print $out_fh "$subseq\n";
        }
    }
} # local $/ expires here

close $out_fh or die $!;
0
投票

好的，通过几个非常小的编辑，您的代码完美无缺。这是对我有用的解决方案：

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

print "Name of the FASTA contig file: ";
chomp( my $fasta_file = <STDIN> );

print "Name file with SNP position info: ";
chomp( my $pos_file = <STDIN> );

print "Name of the output file: ";
chomp( my $out_file = <STDIN> );

open my $out_fh, '>', $out_file or die qq{Unable to open "out_file" for output: $!};


my @pos_records = do {
    open my $pos_, '<' , $pos_file or die qq{Unable to open "$pos_file" for input: $!};
    <$pos_>;
};
chomp @pos_records; #remove all newlines  

{
     open my $fasta_fh, '<', $fasta_file or die qq{Unable to open "$fasta_file" for input: $!};

     local $/ = "\n>"; #Reading FASTA format now

     for ( <$fasta_fh> ) {

         chomp; #Remove ">\n" from the end

         my ( $header, $seq) = split /\n/; #separate the two lines

         $header = ">$header" unless $header =~ /^>/; # Replace any chomped >


     for ( @pos_records ) {

             my ($name,$beg,$end,$pos) = split /\t/;
             my $subseq = substr $seq, $beg-1, $end-$beg;
             my $final_SNP = $end - $pos; 

             if($header =~ /$name/){

               print $out_fh "$header\n";
               print $out_fh "pos=$final_SNP\n";
               print $out_fh "$subseq\n";
     }
    } 
  }
} #local expires here

close $out_fh or die $!;

我改变的唯一实质性的事情是增加了一个if语句。没有它，每个fasta序列被写三次，每一个都有一个具有三个SNP位置之一。我也略微改变了我在做什么来表示SNP位置，在切除序列之后，实际上是$ end - $ pos而不仅仅是$ pos。

再一次，我不能够感谢你，因为很明显你花了很多时间来帮助我。为了它的价值，我真诚地感激它。你的解决方案将作为我未来努力的模板（这可能是对fasta文件的类似操作），你的解释帮助我更好地理解像豌豆大脑可以理解的本地所做的事情。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.