Fasta sequence extraction

Question

[Eva_Maria]

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

[biolab]

Hi,

The getfasta in bedtools can accomplish this job. Please follow the link http://bedtools.readthedocs.org/en/latest/content/tools/getfasta.html for details.

[Eva_Maria]

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

[biolab]

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;

[Eva_Maria]

thanks it works

[Matt Shirley]

If bedtools is not working for you then try https://github.com/mdshw5/pyfaidx#cli-script-faidx

I think bedtools might read the entire FASTA into memory before extracting, but not sure.

[Phil S.]

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

[Eva_Maria]

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

[Phil S.]

what do you mean by "batch processing"?