Extract fasta script

[Shorash]

Hi all, I'm trying to extract nucleotide sequences from a fasta file using the following script which allows multiple extractions:

perl -ne 'if(/^>(\S+)/){$c=grep{/^$1$/}qw()}print if $c' file.fasta

My fasta file "file.fasta" contains contigs with two different headers from two different assemblers which I've merged. It doesn't seem to extract sequences from the second assembly headers but works fine for the first. Here's an example of each header:

Header 1 (working) - comp89447_c0_seq4
Header 2 (won't extract) - BN2_l1_1_(paired)_merged_contig_14067

I can manually search and find the "BN2_l1_1_(paired)_merged_contig_14067" but I cannot extract it using the script. Any help would be greatly appreciated.

[hugorody]

just use grep

why don't use only grep?

cat filein.fasta | grep '>' --after-context=1 > fileout.fast a

[Shorash]

why don't use only grep?

cat filein.fasta | grep '>' --after-context=1 > fileout.fast a

Sorry what exactly can I use this for?

[Shorash]

Sorry what exactly can I use this for?

I don't think I could use this for multiple extractions (up to 200-300)? It also doesn't seem to extract the full sequence following the contig name.

[hugorody]

Sorry what exactly can I use this for?

could you explore a bit more what exactly you need to do?

based on your post I think you just want to capture fasta sequences from a file. is that right?

[hugorody]

grep will catch all lines that contains the symbol >, and also one line after context ">".
you can use this for hundreds... thousands

[Shorash]

could you explore a bit more what exactly you need to do?

based on your post I think you just want to capture fasta sequences from a file. is that right?

Ok, I have a fasta file containing nucleotide sequences with two different headings, eg:

Heading 1:
>comp35107_c0_seq1
GTGGCAGGGACCAGAGCAAGCAGTTCTTCACAGACTTGTGGGAGTTCAGCCTGAAGGACC
TTGAGTGGAAGGACAAGAGTCAACTCATCATTAGTGATGTGGCAGGCATGGTGCCCAGTG
GCCGAAGTGGCGCCTCCACATGGGTCGGAAAGGATCAAGCACTCTACATGTTTGGTGGAA
ACACTGTGGTCCGCACAGACTCAGGTCTCCGCAGCGGCATTGGATATGGAGCTGATCTCT
GGAGGATGTCCACAAACAACCACAGCTGGCAGCTTTTGTCAGGCACTACAAAACCTGGGA
CTCCAGCCAAGTTTGGTCGCCTTGGGGAGTACACTATAATGAGTCAGCCTGGCAGTCGGT
GTGGGGCCATCACCTGGGTGGACACAGCCGGCAACCTGTGGATGTTTGGCGGTGATGGCA
CAGACACAAGTCTTCCTTCTCCCTACCACGCATCACTGCTGCTCTCTGACCT

Heading 2:
>BN2_l1_1_(paired)_merged_contig_20016
TCTCTCTCTCTCTCTCTGTGCCTATTCACATATCTCTTTTTTGTGCGTCTTCTCCTCTAA
ACCACTGCAATAAAACTGTCCGAGTGCAGTCTCTGTCGGGACGCTGATGAAGGGAGGCTG
GGGAGATGGAGAGAGGAGATGACACCCCCAGGTCCTGATTAAGCTGAGAGCTATTGCCGT
AATGGACTAAAAGCACACGGGCGCCGTATTTCCGCTCNNCCGCTCAGACTCCATCCGCTT
TATTCGGGACTTCGATGAGATGAAAGTCCTCGTTGCATTACGCCAATTTGATTACGGCAC
TGATTTGACCCTGCAAACGAACCCCTGCAACTTCAGGAGTGCTCGCCCAATTGGGGTTGG
CACGCTGTGGAACGCTCGAGGCACCGTGGGCAGCCGCCAGACCTTCGGTCTCCAATCTGC
AACGCCGTGGCAGGTGGAATTACAAGGAAATGGACACTCGAACCTCTTTGTGTCAGGAGC
AGATTGCTTGCGGCTGTGGGATTTATTGTAGG

I require a script to extract multiple of these sequences using the headers I have in bold above. So basically I will copy a large number of the headers above from Excel and using a script I want to extract all the corresponding sequences. My current script which I've outlined in the opening post, for some reason, only extracts the sequences from header 1.

Hope that clarifies a little.

[hugorody]

OK. are u using Linux right?

So, create a list.txt file with all headers you need separated by enter:

header1
header2
...

then type on shell:

$cat name_your_fasta_file.fasta | grep --file=list.txt --after-context=1 > my_sequences.fasta

and that's it.
you don't need a script.

[rhinoceros]

OK. are u using Linux right?

So, create a list.txt file with all headers you need separated by enter:

header1
header2
...

then type on shell:

$cat name_your_fasta_file.fasta | grep --file=list.txt --after-context=1 > my_sequences.fasta

and that's it.
you don't need a script.

It doesn't work if there are linebreaks in the seqs like OP posted. You first have to deal with them..

[SylvainL]

Hi, just use R and the seqinr package. You will find the read.fasta function...

[hugorody]

It doesn't work if there are linebreaks in the seqs like OP posted. You first have to deal with them..

To remove the linebreaks:

$cat your_file_with_linebreaks.fasta | sed 's/ //g' | sed 's/\(>.*\)/\1 /g' | sed ':a;N;s/\n//g;ta' | sed 's/>/\n>/g' | sed 's/ /\n/g' > file_without_linebreaks.fasta

[Shorash]

To remove the linebreaks:

$cat your_file_with_linebreaks.fasta | sed 's/ //g' | sed 's/\(>.*\)/\1 /g' | sed ':a;N;s/\n//g;ta' | sed 's/>/\n>/g' | sed 's/ /\n/g' > file_without_linebreaks.fasta

Hi there,

I receive the following error when attempting to remove the linebreaks:

:~/Extract_fasta/Orthologs> $cat BN_clc.fasta | sed 's/ //g' | sed 's/\(>.*\)/\1 /g' | sed ':a;N;s/\n//g;ta' | sed 's/>/\n>/g' | sed 's/ /\n/g' > BN_Fixed.fasta
./BN_clc.fasta: line 1: syntax error near unexpected token `('
'/BN_clc.fasta: line 1: `>BN2_l1_1_(paired)_merged_contig_4

[GenoMax]

faSomeRecords from Kent utilities would be the simplest/fast solution (http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/)

http://seqanswers.com/forums/showthread.php?t=9498

https://www.biostars.org/p/2822/

https://www.biostars.org/p/1195/

[hugorody]

Hi there,

I receive the following error when attempting to remove the linebreaks:

:~/Extract_fasta/Orthologs> $cat BN_clc.fasta | sed 's/ //g' | sed 's/\(>.*\)/\1 /g' | sed ':a;N;s/\n//g;ta' | sed 's/>/\n>/g' | sed 's/ /\n/g' > BN_Fixed.fasta
./BN_clc.fasta: line 1: syntax error near unexpected token `('
'/BN_clc.fasta: line 1: `>BN2_l1_1_(paired)_merged_contig_4

which linux distribution do you use?
you should try use double quotes ( " ) instead single quotes ( ' ).

[Shorash]

which linux distribution do you use?
you should try use double quotes ( " ) instead single quotes ( ' ).

I'm using a portable batch system (PBS).