Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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.

  • #2
    just use grep

    why don't use only grep?

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

    Comment


    • #3
      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?

      Comment


      • #4
        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.

        Comment


        • #5
          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?

          Comment


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

            Comment


            • #7
              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.

              Comment


              • #8
                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, 10:31 PM.

                Comment


                • #9
                  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

                  Comment


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

                    Comment


                    • #11
                      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

                      Comment


                      • #12
                        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

                        Comment


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

                          Discussion of next-gen sequencing related bioinformatics: resources, algorithms, open source efforts, etc




                          Last edited by GenoMax; 08-04-2014, 05:10 PM.

                          Comment


                          • #14
                            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 ( ' ).

                            Comment


                            • #15
                              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).

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                The Impact of AI in Genomic Medicine
                                by seqadmin



                                Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
                                02-26-2024, 02:07 PM
                              • seqadmin
                                Multiomics Techniques Advancing Disease Research
                                by seqadmin


                                New and advanced multiomics tools and technologies have opened new avenues of research and markedly enhanced various disciplines such as disease research and precision medicine1. The practice of merging diverse data from various ‘omes increasingly provides a more holistic understanding of biological systems. As Maddison Masaeli, Co-Founder and CEO at Deepcell, aptly noted, “You can't explain biology in its complex form with one modality.”

                                A major leap in the field has
                                ...
                                02-08-2024, 06:33 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, 02-28-2024, 06:12 AM
                              0 responses
                              28 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 02-23-2024, 04:11 PM
                              0 responses
                              74 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 02-21-2024, 08:52 AM
                              0 responses
                              85 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 02-20-2024, 08:57 AM
                              0 responses
                              69 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X