☰

在 fasta 文件上查找具有位置、行和列信息的引物

tag-icon • tag-icon grep bash-scripting
在 fasta 文件上查找具有位置、行和列信息的引物

我有一些fasta 文件每行恰好包含 80 列。忘记开头的定义行>，考虑其余大约有 1250000 行。我需要找到一些入门知识，例如CAATCGCCGT它们的位置、行和列。我将所有行转换为一行，然后使用来grep查找它们在 Bash 脚本中的位置。

grep -o -P -b -n CAATCGCCGT input.fasta

3:3206721:CAATCGCCGT以此格式在原始文件上输出。

我现在需要的是显示原始 fasta 文件中的行号和列号，例如

 3:3206721[40084:1]:CAATCGCCGT      #that is just `div 80` and `mod 80` see below

 3:3206721[3206721 / 80:3206721 % 80]:CAATCGCCGT

样本快拍；

ATCCATTTGTCTTTCCCTCAATCGCCGTGCTCCTTATAAACAATCGCCTTCGGTGTACCCCTGTTAGGCG
CGTACGAATGTATCTGCGGTGAGTCGAAAGTGAGAGCTCCGCCTGCCAGAAGCCTTCCTGCTGCAACCAA
TCGCCGTGAGTATGACAAATTCATTTTCAAACAACCTGCCCCCTACGATCCGACTCCGACTCTCATTCCC
AATTCCGGCAAATCCCTGTGTGCATGTCAACCCCCGGAGGGAAGCACACATTTCCTTGCGGAAGAAAGTT
AGCTTTGGCCGGCCGTCGTCTTCTTGCCTTCCGGGATTTGCCCGTCCCCGGTGGGGATGCCCTAACGGTA
TCAAATTGGTGACGTTTCAAATCGTTTGCAATCAATCGCCGTCCAGCTGCATAACTTGGGGCGCGTTTCG

我不能简单地得到输出grep并将其除以 80 并取模。我需要一个解决方案，而直接在 fasta 文件上工作的解决方案更好，因为不需要将输入转换为一行。两者都很好。

样本输入中，有一个引物占据两条线，其他的引物没有分开。

答案1

我看过了，最简单的解决办法是映泰这是用 Perl 实现的。

use strict;
use warnings;
use Bio::SeqIO;

my $usage = "perl dnamotif.pl <fasta file> <primers file>";
my $fasta_filename = shift(@ARGV) or die("Usage: $usage $!");
my $pfile = shift(@ARGV) or die("Usage: $usage $!");
my $start = 0;
my $motifCounter=1; 

my $fasta_parser = Bio::SeqIO->new(-file => $fasta_filename, -format => 'Fasta');



while(my $seq_obj = $fasta_parser->next_seq())
{
    printf("Searching sequence '%s'...\n", $seq_obj->id);

    open my $info, $pfile or die "Could not open primers $pfile: $!";

    while( my $motif = <$info>)  {   
        chomp $motif;
        printf("\n[%2s Looking motif [%s]]\n", $motifCounter, $motif);    
        $start = 0;

        while((my $pos = index($seq_obj->seq(), $motif, $start)) != -1) {

            printf("\nmotif found at position %8d ", $pos + 1);
            printf("[%8d,%8d]", ((($pos)/70))+1, (($pos) % 70)+1);
            $start = $pos + 1;
        
        }
        
        $motifCounter++;
        printf("\n");
    
    }
    close $info;

}

我找到了基本代码所以并根据我的需要对其进行了修改（我需要搜索 80 个引物）。我将其放在这里以满足任何人的需求。实际文件的部分输出是；

Searching sequence 'ref|NW_024108998.1|:1-6077679'...

[ 1 Looking motif [CAATCGCCGT]]


[ 2 Looking motif [TTCCGAACCC]]

motif found at position   660830 [    9441,      30]
motif found at position   788016 [   11258,      26]
motif found at position  1349822 [   19284,      12]

相关内容