SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
extract FASTA from BAM jwhite Bioinformatics 1 05-15-2014 02:33 PM
Making a Blast Script for FASTA files utagenomics Bioinformatics 8 04-26-2013 11:29 AM
A script for Cuffdiff: extract differential expression gene sequences from genome.fa foxriverlin Bioinformatics 17 04-11-2013 02:46 AM
Any script to format headers in fasta files? Shishir Bioinformatics 2 02-05-2013 06:52 AM
Bioperl Script to convert fasta to fq Lizex Bioinformatics 0 01-26-2012 10:21 PM

Reply
 
Thread Tools
Old 08-03-2014, 04:32 PM   #1
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default Extract fasta script

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.
Shorash is offline   Reply With Quote
Old 08-03-2014, 09:30 PM   #2
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default just use grep

why don't use only grep?

cat filein.fasta | grep '>' --after-context=1 > fileout.fast a
hugorody is offline   Reply With Quote
Old 08-03-2014, 10:04 PM   #3
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by hugorody View Post
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 is offline   Reply With Quote
Old 08-03-2014, 10:09 PM   #4
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by Shorash View Post
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.
Shorash is offline   Reply With Quote
Old 08-03-2014, 10:12 PM   #5
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default

Quote:
Originally Posted by Shorash View Post
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 is offline   Reply With Quote
Old 08-03-2014, 10:15 PM   #6
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default

grep will catch all lines that contains the symbol >, and also one line after context ">".
you can use this for hundreds... thousands
hugorody is offline   Reply With Quote
Old 08-03-2014, 10:17 PM   #7
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by hugorody View Post
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.
Shorash is offline   Reply With Quote
Old 08-03-2014, 10:27 PM   #8
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default

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.

Last edited by hugorody; 08-03-2014 at 10:31 PM.
hugorody is offline   Reply With Quote
Old 08-03-2014, 11:56 PM   #9
rhinoceros
Senior Member
 
Location: sub-surface moon base

Join Date: Apr 2013
Posts: 372
Default

Quote:
Originally Posted by hugorody View Post
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..
__________________
savetherhino.org
rhinoceros is offline   Reply With Quote
Old 08-04-2014, 02:18 AM   #10
SylvainL
Senior Member
 
Location: Geneva

Join Date: Feb 2012
Posts: 179
Default

Hi, just use R and the seqinr package. You will find the read.fasta function...
SylvainL is offline   Reply With Quote
Old 08-04-2014, 11:21 AM   #11
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default

Quote:
Originally Posted by rhinoceros View Post
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
hugorody is offline   Reply With Quote
Old 08-04-2014, 04:49 PM   #12
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by hugorody View Post
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
Shorash is offline   Reply With Quote
Old 08-04-2014, 05:00 PM   #13
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,077
Default

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/

Last edited by GenoMax; 08-04-2014 at 05:10 PM.
GenoMax is offline   Reply With Quote
Old 08-04-2014, 05:04 PM   #14
hugorody
Junior Member
 
Location: Vancouver, BC

Join Date: May 2013
Posts: 9
Default

Quote:
Originally Posted by Shorash View Post
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 ( ' ).
hugorody is offline   Reply With Quote
Old 08-04-2014, 05:14 PM   #15
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by hugorody View Post
which linux distribution do you use?
you should try use double quotes ( " ) instead single quotes ( ' ).
I'm using a portable batch system (PBS).
Shorash is offline   Reply With Quote
Old 08-04-2014, 05:21 PM   #16
Shorash
Member
 
Location: Brisbane

Join Date: Sep 2012
Posts: 18
Default

Quote:
Originally Posted by GenoMax View Post
Thank you very much! Works perfectly and pulls both headers out

Thanks everyone for the help.
Shorash is offline   Reply With Quote
Old 10-08-2014, 09:58 PM   #17
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Sounds like you want someone to do your homework for you.
Brian Bushnell is offline   Reply With Quote
Old 10-09-2014, 05:05 AM   #18
amitm
Member
 
Location: Manchester, UK

Join Date: Feb 2011
Posts: 52
Default

Quote:
Originally Posted by hugorody View Post
why don't use only grep?

cat filein.fasta | grep '>' --after-context=1 > fileout.fast a
Thanks hugorody for mentioning about this parameter in grep. I wasn't aware of it. Saves a lot of scripting!
amitm is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:15 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO