Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Help with script, sort by coverage

    Hello everyone,

    I am not good with with writing optimal scripts, so I would appreciate someones help, either perl or python is fine.

    I have contigs generated by assemblers in the following format:
    NODE_2361_length_509_cov_1.43018_ID_745236

    This is common to velvet and SPAdes. I was wondering how could I extract all contig names that in their name have a coverage above a user specified value. So I could perl selectAboveCoverage.pl 1.0 and get all above 1.0 including 1.0?

    I know I have to first split by >, and then the index 0 will be contig name and index 1 will be the sequence. I need to further split index 0 by _, and index 5 will have the coverage value. I can then somehow sort by that and pick all of those above argv[0], and write them to a new file argv[1].

    I know this from experience in TCL, but scripts on those languages take too long to execute.

    Any help will be appreciated, thank you!

  • #2
    How about

    Code:
    awk '{FS="_"; if($6>=1) print}' file
    savetherhino.org

    Comment


    • #3
      Originally posted by rhinoceros View Post
      How about

      Code:
      awk '{FS="_"; if($6>=1) print}' file

      You sir, are awesome.

      Comment


      • #4
        How would I use a similar script to sort the contigs by lengths, and tell me what the top 10 or 20 lengths are?

        Comment


        • #5
          Could try adding to the above:

          Code:
          sort -r -n file | head
          In short:

          Code:
          sort -r -n file
          sorts the contents of file, -r to use reverse order (highest to lowest), and -n uses numeric sort. Otherwise you get 1, 10, 100, 2, 20, 200, 3, 30, 300.

          head just displays the top 10 of the file, or stdout (lest it dump a pile of numbers into stdout). To adjust the number of lines displayed, use -n. For example,
          Code:
          head -n 20
          displays the first 20 lines of a file.

          Comment


          • #6
            Sorry but the above shows contig names, how will sort know what to sort by? Since it's given the full length of the contig?

            Comment


            • #7
              Oh, whoops.

              Contents of file I wish to process
              Code:
              NODE_2361_length_509_cov_1.43018_ID_745236
              NODE_2361_length_509_cov_5.43018_ID_745236
              NODE_2361_length_509_cov_3.43018_ID_745236
              NODE_2361_length_509_cov_2.43018_ID_745236
              NODE_2361_length_509_cov_21.43018_ID_745236
              The one-liner
              Code:
              sed 's/_/ /g' filename | awk '{FS=" "; print $6}' | sort -r -n | head -n 10
              sed will replace _ with space, then spit out the output to stdout. | to pipe to awk.

              awk will use spaces to separate columns, print $6 instructs awk to display everything from column six (which is the numbers after cov, but before ID).
              You can then pipe that again to sort, with head to display just what you want to display.

              The output
              Code:
              21.43018
              5.43018
              3.43018
              2.43018
              1.43018
              Obviously not the only solution, and there are probably better ones.

              --------------
              Edit: If you're dealing with fasta you will have to separate the headers from the sequences. Standard format is >NODE_...; which can be captured as
              Code:
              grep "^>NODE" contigs.fasta > contignames
              Last edited by winsettz; 10-07-2013, 08:32 AM.

              Comment


              • #8
                I would just pipe the awk output to:

                Code:
                sort -r -g -k 4 -t _
                When you have an idea what tool might work, but don't know exactly how, just google "man sort" or whatever the name of the tool happens to be. Also, with sort, it's sometimes better to use -g (general numerical value) instead of -n. For example, -n (string numerical value) doesn't work when you have exponents (1e-10 would be smaller than 2e-100 so e.g. sorting blast output by e-value would fail). Actually, I don't know if there's any reason why -g shouldn't be used always by default..
                Last edited by rhinoceros; 10-07-2013, 08:46 AM.
                savetherhino.org

                Comment

                Latest Articles

                Collapse

                • seqadmin
                  Latest Developments in Precision Medicine
                  by seqadmin



                  Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

                  Somatic Genomics
                  “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
                  05-24-2024, 01:16 PM
                • seqadmin
                  Recent Advances in Sequencing Analysis Tools
                  by seqadmin


                  The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
                  05-06-2024, 07:48 AM

                ad_right_rmr

                Collapse

                News

                Collapse

                Topics Statistics Last Post
                Started by seqadmin, 05-24-2024, 07:15 AM
                0 responses
                198 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 05-23-2024, 10:28 AM
                0 responses
                220 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 05-23-2024, 07:35 AM
                0 responses
                229 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 05-22-2024, 02:06 PM
                0 responses
                13 views
                0 likes
                Last Post seqadmin  
                Working...
                X