This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Fasta sequence extraction

I have tab limited file like this

GCF_000707685.1_2840 0 145
GCF_000706885.1_542 0 150
GCF_000593365.1_489 0 156
GCF_000593345.1_3957 256 289
GCF_000593325.1_3041 780 958

I want to extract sequence based on position like 0 to 145 , 256 to 289

>GCF_000707685.1_2840
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
>GCF_000706885.1_542
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

How to do this? awk or perl

blast awk perl genome alignment
• 4,037 views • link
updated by Ram Ram • written by Eva_Maria Eva_Maria

its not working for me (files about 7 gb). is there any programme or awk?

• 0 views • link
updated by Ram Ram • written by Eva_Maria Eva_Maria

Hi, I write a perl script as follows. Usage: perl extr_seq.pl genome.fa list.txt. I am just wondering why bedtools does not work for you. Personally I think it's a good tool.

#!/usr/bin/perl -w
use strict;

my %data;
my $id = '';
my $seq = '';

open my $fasta, '<', $ARGV[0] or die $!;
while (<$fasta>) {
    chomp;
    s/\r//;
    next if /^\s*$/;
    if (/>(\S+)/) {
        $data{$id} = $seq if $seq ne '';
        $id = $1;
        $seq = '';
    } else {
        $seq .= $_;
    }
}
$data{$id} = $seq;
close $fasta;

open my $list, '<', $ARGV[1] or die $!;
while (<$list>) {
    chomp;
    s/\r//;
    next if /^\s*$/;
    my ($gene, $start, $end) = split;
    if (exists $data{$gene}) {
        my $extr_seq = substr($data{$gene}, $start - 1, $end - $start + 1);
        print ">$gene\n$extr_seq\n";
    }
}
close $list;
• 0 views • link
updated by Ram Ram • written by biolab biolab

thanks it works

• 0 views • link
updated by Ram Ram • written by Eva_Maria Eva_Maria

Here is an awk, one line solution to you problem.

awk 'FS="\t" {if(($2==0) && ($3==145) || ($2==256) && ($3 == 289)) print $1}' < foo.file.txt

if this is your file:

GCF_000707685.1_2840 0 145
GCF_000706885.1_542 0 150
GCF_000593365.1_489 0 156
GCF_000593345.1_3957 256 289
GCF_000593325.1_3041 780  958

the output is:

$ awk 'FS="\t" {if(($2==0) && ($3==145) || ($2==256) && ($3 == 289)) print $1}' < foo.file.txt
GCF_000707685.1_2840
• 0 views • link
updated by Ram Ram • written by Phil S. Phil S.

I want to do as a Batch processing because i have 20 mb text file and 256 mb fasta

• 0 views • link
updated by Ram Ram • written by Eva_Maria Eva_Maria

what do you mean by "batch processing"?

• 0 views • link
written by Phil S. Phil S.

0 answers

No answers yet.

Log in to answer this question.

More posts like this