Question: fasta files: keep only lines containing certain characters + the respective sequences
0
gravatar for Tmatamatamas
3.0 years ago by
Tmatamatamas0 wrote:

Hi. I have a fasta file containing 18S and 16S sequences. Can I use sed or grep to find all headers containing 18S and print a new fasta file that includes only these headers and the respective sequences?

Edit: What I have:

>GQ118297.1 Tricorythus sp. BYU IGCEP206 16S ribosomal RNA gene, partial sequence; mitochondrial
TTGTCTCTTGGAAATATAGGAGATAGGGCCTGCCCAATGAAATTTCAATGGCCGCAGTAATTTGACTGTG
CAAAGGTAGCATAATCATTAGTTTTTTAATTGAGGACTGGTATGAAAGGCATAATGAGGTACTTGTTTTC
TTAAATAAAAGAATAAAATTTTACTTTTTAGTTAAAAGGCTAAAGTAATATAAAGGGACGAGAAGACCCT
ATAGAGTTTTATAAATTAAATTAATTTATTTTAGTAAAATTAAAGAATTGATGGAATTTATTTAGTTGGG
GAGATTTTGTAATAAAACTTATAATTATATAAACATTTATATATGATTTTAAGATCCATAAATTGATTAA
AAAATTAAATTACCTTAGGGATAACAGCGTTATTTTCTTGGAGAGTTCTTATCAATAGGAAAGTTTGCGA
CCTCGATGTTGGATTAAGAAAATAGTTAAATGAAGCCGTTTAATTAATAGGTCTGTTCGACCTTTAAAAT
CTTACA
>GQ118277.1 Teloganella sp. BYU IGCEP161 18S ribosomal RNA gene, partial sequence
GCTGTCTCAGTGCAAGCCTAATTAAAGTGAAACCGCAAATGGCTCATTAAATCAGTTATGGTTCATTAGA
TGATACATTATTTTACTTGGATAACTGTGGTAATTCTAGAGCTAATACATGCAAATTGAGTTCCCATCGG
TGACGGTAGGAACGCTTTTATTAGATCAAAACCAATACGTTTGCTTCGGCATACGATTTCAATGGTGATT
CTGAATAACTTTTTGATGATCGTACGGTCCTTGTATCGACGACAAATCTTTCAAATGTCTGCCTTATCAA
CTGTCGATGGTAGGCTCTGCGCCTACCATGGTTGTAACGGGTAACGGGGAATCAGGGTTCGATTCCGGAG
AGGGAGCCTGAGAAACGGCTACCACATCCAAGGAAGGCAGCAGGCGCGCAAATTACCCACTCCCGGCACG
GGGAGGTAGTGACGAAAAATAACGATACGGGACTCATCCGAGGCCCCGTAATCGGAATGAGMACACTTTA
AATMYTTTAACRAGTAYCYAWTGGAGGGCAAGTCTGGTGCCAGCAGCCGCGGTAATTCCAGCTCCATTGG
CGTATATTAAAGTTGTTGCGGTTAAAAAGCTCGGTAGTTGGTATCATGTGTMTCGGACGGTCGGTTCGCC
NCNCNCGGTGTTCAACTGACCGGTCCGGACGTCCTGCCGGTGGGACCCGGTTCGCGCCGGGCCCCGT

What I would like to have after processing by identifying all sequences wit "18S" in the header and printing them into a new file:

>GQ118277.1 Teloganella sp. BYU IGCEP161 18S ribosomal RNA gene, partial sequence
GCTGTCTCAGTGCAAGCCTAATTAAAGTGAAACCGCAAATGGCTCATTAAATCAGTTATGGTTCATTAGA
TGATACATTATTTTACTTGGATAACTGTGGTAATTCTAGAGCTAATACATGCAAATTGAGTTCCCATCGG
TGACGGTAGGAACGCTTTTATTAGATCAAAACCAATACGTTTGCTTCGGCATACGATTTCAATGGTGATT
CTGAATAACTTTTTGATGATCGTACGGTCCTTGTATCGACGACAAATCTTTCAAATGTCTGCCTTATCAA
CTGTCGATGGTAGGCTCTGCGCCTACCATGGTTGTAACGGGTAACGGGGAATCAGGGTTCGATTCCGGAG
AGGGAGCCTGAGAAACGGCTACCACATCCAAGGAAGGCAGCAGGCGCGCAAATTACCCACTCCCGGCACG
GGGAGGTAGTGACGAAAAATAACGATACGGGACTCATCCGAGGCCCCGTAATCGGAATGAGMACACTTTA
AATMYTTTAACRAGTAYCYAWTGGAGGGCAAGTCTGGTGCCAGCAGCCGCGGTAATTCCAGCTCCATTGG
CGTATATTAAAGTTGTTGCGGTTAAAAAGCTCGGTAGTTGGTATCATGTGTMTCGGACGGTCGGTTCGCC
NCNCNCGGTGTTCAACTGACCGGTCCGGACGTCCTGCCGGTGGGACCCGGTTCGCGCCGGGCCCCGT
sed grep fasta • 1.3k views
ADD COMMENTlink modified 3.0 years ago by Paul1.4k • written 3.0 years ago by Tmatamatamas0
1

Can I use sed or grep to find all headers containing 18S and print a new fasta file that includes only these headers and the respective sequences?

no, you need something like awk.

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by Pierre Lindenbaum129k
1

Whether or not you can do this depends on the contents of the file, and how precise you need to be. In the general case, no, it is not possible. I suggest you post the headers of some of your sequences, along with as much other data as possible. And clarify whether you wish to separate sequences based on their header or contents.

ADD REPLYlink written 3.0 years ago by Brian Bushnell17k
1

if you know a little bit of python, Biopython SeqIO would be an easy and intuitive solution.

ADD REPLYlink written 3.0 years ago by grant.hovhannisyan2.0k

Unfortunately, I don't (yet), so Python does not help me at the moment

ADD REPLYlink written 3.0 years ago by Tmatamatamas0

Agreed, BioPython is totally the way to go for this

ADD REPLYlink written 3.0 years ago by Joe17k

Tmatamatamas : You have multiple options of correct answers below. Please "accept" (green check mark) as many as you wish to provide closure to this thread.

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by genomax85k

Hi, if you want to use grep, just try:

grep -w "18S" -A 1 input.fasta > 18S.fasta

parameter A means, that grep print one lines of trailing context. Parameter w match only whole word -in our case 18S.

EDIT: This works only if you have in header 18S and no-where else. Also you need correct format of FASTA file - first line started with > is header and second line nucleotides.

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by Paul1.4k

OP wants to keep the entire sequence where the header contains 18S word.

Edit: Moved to a comment since it is not an appropriate answer for the original question.

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by genomax85k

This is a Downvote:)

ADD REPLYlink written 3.0 years ago by grant.hovhannisyan2.0k

This would work if the fasta was a single-line fasta and not multi-line like the OP has.

ADD REPLYlink written 3.0 years ago by st.ph.n2.5k

st.ph.n, thank you for the clarity on the above answer. Is there any way to do this with grep?

ADD REPLYlink written 3.0 years ago by vinayjrao170
1

Unless the fasta's are linearised, as in the above comment, then no not really. I might be possible, but would be ugly and clunky I'm sure.

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by Joe17k
3
gravatar for Pierre Lindenbaum
3.0 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum129k wrote:
awk '/^>/ { ok=index($0,"18S")!=0;} {if(ok) print;}'  input.fa
ADD COMMENTlink written 3.0 years ago by Pierre Lindenbaum129k

That's it, thanks a lot!

ADD REPLYlink written 3.0 years ago by Tmatamatamas0
3
gravatar for Frédéric Mahé
3.0 years ago by
France, Montpellier, CIRAD
Frédéric Mahé3.0k wrote:

Please Pierre Lindenbaum, don't underestimate sed :-)

sed -n '/^>.*18S/,/^>/ {/^>.*18S/p ; /^>/! p}' input.fa
ADD COMMENTlink modified 3.0 years ago • written 3.0 years ago by Frédéric Mahé3.0k

excellent ! :-)

ADD REPLYlink written 3.0 years ago by Pierre Lindenbaum129k
0
gravatar for cpad0112
3.0 years ago by
cpad011213k
India
cpad011213k wrote:

my fasta:

 cat test2.fa 
>GQ118297.1 Tricorythus sp. BYU IGCEP206 16S ribosomal RNA gene, partial sequence; mitochondrial
TTGTCTCTTGGAAATATAGGAGATAGGGCCTGCCCAATGAAATTTCAATGGCCGCAGTAATTTGACTGTG
CAAAGGTAGCATAATCATTAGTTTTTTAATTGAGGACTGGTATGAAAGGCATAATGAGGTACTTGTTTTC
TTAAATAAAAGAATAAAATTTTACTTTTTAGTTAAAAGGCTAAAGTAATATAAAGGGACGAGAAGACCCT
ATAGAGTTTTATAAATTAAATTAATTTATTTTAGTAAAATTAAAGAATTGATGGAATTTATTTAGTTGGG
GAGATTTTGTAATAAAACTTATAATTATATAAACATTTATATATGATTTTAAGATCCATAAATTGATTAA
AAAATTAAATTACCTTAGGGATAACAGCGTTATTTTCTTGGAGAGTTCTTATCAATAGGAAAGTTTGCGA
CCTCGATGTTGGATTAAGAAAATAGTTAAATGAAGCCGTTTAATTAATAGGTCTGTTCGACCTTTAAAAT
CTTACA
>GQ118277.1 Teloganella sp. BYU IGCEP161 18S ribosomal RNA gene, partial sequence
GCTGTCTCAGTGCAAGCCTAATTAAAGTGAAACCGCAAATGGCTCATTAAATCAGTTATGGTTCATTAGA
TGATACATTATTTTACTTGGATAACTGTGGTAATTCTAGAGCTAATACATGCAAATTGAGTTCCCATCGG
TGACGGTAGGAACGCTTTTATTAGATCAAAACCAATACGTTTGCTTCGGCATACGATTTCAATGGTGATT
CTGAATAACTTTTTGATGATCGTACGGTCCTTGTATCGACGACAAATCTTTCAAATGTCTGCCTTATCAA
CTGTCGATGGTAGGCTCTGCGCCTACCATGGTTGTAACGGGTAACGGGGAATCAGGGTTCGATTCCGGAG
AGGGAGCCTGAGAAACGGCTACCACATCCAAGGAAGGCAGCAGGCGCGCAAATTACCCACTCCCGGCACG
GGGAGGTAGTGACGAAAAATAACGATACGGGACTCATCCGAGGCCCCGTAATCGGAATGAGMACACTTTA
AATMYTTTAACRAGTAYCYAWTGGAGGGCAAGTCTGGTGCCAGCAGCCGCGGTAATTCCAGCTCCATTGG
CGTATATTAAAGTTGTTGCGGTTAAAAAGCTCGGTAGTTGGTATCATGTGTMTCGGACGGTCGGTTCGCC
NCNCNCGGTGTTCAACTGACCGGTCCGGACGTCCTGCCGGTGGGACCCGGTTCGCGCCGGGCCCCGT

Command:

seqkit grep -r -n -i -p  18s test2.fa

output:

 seqkit grep -r -n -i -p  18s test2.fa 
>GQ118277.1 Teloganella sp. BYU IGCEP161 18S ribosomal RNA gene, partial sequence
GCTGTCTCAGTGCAAGCCTAATTAAAGTGAAACCGCAAATGGCTCATTAAATCAGTTATG
GTTCATTAGATGATACATTATTTTACTTGGATAACTGTGGTAATTCTAGAGCTAATACAT
GCAAATTGAGTTCCCATCGGTGACGGTAGGAACGCTTTTATTAGATCAAAACCAATACGT
TTGCTTCGGCATACGATTTCAATGGTGATTCTGAATAACTTTTTGATGATCGTACGGTCC
TTGTATCGACGACAAATCTTTCAAATGTCTGCCTTATCAACTGTCGATGGTAGGCTCTGC
GCCTACCATGGTTGTAACGGGTAACGGGGAATCAGGGTTCGATTCCGGAGAGGGAGCCTG
AGAAACGGCTACCACATCCAAGGAAGGCAGCAGGCGCGCAAATTACCCACTCCCGGCACG
GGGAGGTAGTGACGAAAAATAACGATACGGGACTCATCCGAGGCCCCGTAATCGGAATGA
GMACACTTTAAATMYTTTAACRAGTAYCYAWTGGAGGGCAAGTCTGGTGCCAGCAGCCGC
GGTAATTCCAGCTCCATTGGCGTATATTAAAGTTGTTGCGGTTAAAAAGCTCGGTAGTTG
GTATCATGTGTMTCGGACGGTCGGTTCGCCNCNCNCGGTGTTCAACTGACCGGTCCGGAC
GTCCTGCCGGTGGGACCCGGTTCGCGCCGGGCCCCGT
ADD COMMENTlink written 3.0 years ago by cpad011213k

seqkit can be found here.

ADD REPLYlink written 3.0 years ago by genomax85k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1525 users visited in the last hour