Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • #16
    The perl stuff is faster, at least if you have a few 100,000 of contigs

    Sven

    Comment


    • #17
      There is a program called getN50 in the amos package.



      Comment


      • #18
        A single awk, no pipes (except between sort and awk), and thereby somewhat shorter. Note sort -n, not sort -rn

        sort -n contig_lengths.txt | awk '{len[i++]=$1;sum+=$1} END {for (j=0;j<i+1;j++) {csum+=len[j]; if (csum>sum/2) {print len[j];break}}}'

        Comment


        • #19
          Perl clearly 10 times faster...

          macbook:~$ time perl -ne 'chomp(); push(@contigs,$_);$total+=$_;END{foreach(sort{$b<=>$a}@contigs){$sum+=$_;$L=$_;if($sum>=$total*0.5){print "TOTAL: $total\nN50 : $L\n";exit;} ;}}' contigs_100.lines
          TOTAL: 59789620
          N50 : 212

          real 0m0.444s
          user 0m0.348s
          sys 0m0.038s
          macbook:~ $ time sort -n contigs_100.lines | awk '{len[i++]=$1;sum+=$1} END {for (j=0;j<i+1;j++) {csum+=len[j]; if (csum>sum/2) {print len[j];break}}}'
          212

          real 0m5.502s
          user 0m4.055s
          sys 0m0.560s
          --------------------------------------
          Elia Stupka
          Co-Director and Head of Unit
          Center for Translational Genomics and Bioinformatics
          San Raffaele Scientific Institute
          Via Olgettina 58
          20132 Milano
          Italy
          ---------------------------------------

          Comment


          • #20
            Cool. Well, consider it an exercise in awk rather than an attempt to beat perl...

            Comment


            • #21
              maubp's code on Python requires large amount of memory and CPU if numbers in the list are huge (it creates X copies of every number X).

              I suggests to use this faster function in Python for calculating N50 based on this definition http://seqanswers.com/forums/showpos...6&postcount=4:

              PHP Code:
              def N50(numlist):
                
              """
                Abstract: Returns the N50 value of the passed list of numbers.
                Usage: N50(numlist)

                Based on the definition from this SEQanswers post
                http://seqanswers.com/forums/showpost.php?p=7496&postcount=4
                (modified Broad Institute's definition
                https://www.broad.harvard.edu/crd/wiki/index.php/N50)
                
                See SEQanswers threads for details:
                http://seqanswers.com/forums/showthread.php?t=2857
                http://seqanswers.com/forums/showthread.php?t=2332
                """
                
              numlist.sort(reverse True)
                
              sum(numlist)
                
              limit 0.5
                
              for l in numlist:
                  
              -= l
                  
              if <= limit:
                    return 

              Originally posted by maubp View Post
              OK - so the stdin is one integer per line. How about a python script like this,
              see also http://seqanswers.com/forums/showthread.php?t=2332

              Code:
              #!/usr/bin/python
              import sys
              
              def N50(numlist):
                  """
                  Abstract: Returns the N50 value of the passed list of numbers. 
                  Usage:    N50(numlist)
              
                  Based on the Broad Institute definition:
                  https://www.broad.harvard.edu/crd/wiki/index.php/N50
                  """
                  numlist.sort()
                  newlist = []
                  for x in numlist :
                      newlist += [x]*x
                  # take the mean of the two middle elements if there are an even number
                  # of elements.  otherwise, take the middle element
                  if len(newlist) % 2 == 0:
                      medianpos = len(newlist)/2  
                      return float(newlist[medianpos] + newlist[medianpos-1]) /2
                  else:
                      medianpos = len(newlist)/2
                      return newlist[medianpos]
              
              assert N50([2, 2, 2, 3, 3, 4, 8, 8]) == 6
              
              lengths = []
              for line in sys.stdin :
                  lengths.append(int(line))
              print N50(lengths)
              Then at the Unix command line, you could use it like this:
              Code:
              $ grep "^>" 454AllContigs.fna | cut -d"=" -f2 | cut -d" " -f1 | ./stdin_N50.py 
              386

              Comment


              • #22
                Originally posted by flxlex View Post
                A single awk, no pipes (except between sort and awk), and thereby somewhat shorter. Note sort -n, not sort -rn

                sort -n contig_lengths.txt | awk '{len[i++]=$1;sum+=$1} END {for (j=0;j<i+1;j++) {csum+=len[j]; if (csum>sum/2) {print len[j];break}}}'

                So most scripts here actually take the scaffold length at which the sum is bigger than half the genome size. But should it not be equal AND/OR bigger ?


                The N50 of an assembly is a weighted median of the lengths of the sequences it contains, equal to the length of the longest sequence s, such that the sum of the lengths of sequences greater than or equal in length to s is greater than or equal to half the length of the genome being assembled.
                from the Assemblathon paper

                Comment


                • #23
                  Originally posted by ebioman View Post
                  But should it not be equal AND/OR bigger?
                  Hehe, you're right, that is a mistake. The correct version is

                  Code:
                  sort -n contig_lengths.txt | awk '{len[i++]=$1;sum+=$1} END {for (j=0;j<i+1;j++) {csum+=len[j]; if (csum>=sum/2) {print len[j];break}}}'

                  Comment

                  Latest Articles

                  Collapse

                  • seqadmin
                    Strategies for Sequencing Challenging Samples
                    by seqadmin


                    Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                    03-22-2024, 06:39 AM
                  • seqadmin
                    Techniques and Challenges in Conservation Genomics
                    by seqadmin



                    The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                    Avian Conservation
                    Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                    03-08-2024, 10:41 AM

                  ad_right_rmr

                  Collapse

                  News

                  Collapse

                  Topics Statistics Last Post
                  Started by seqadmin, Yesterday, 06:37 PM
                  0 responses
                  8 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, Yesterday, 06:07 PM
                  0 responses
                  8 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 03-22-2024, 10:03 AM
                  0 responses
                  49 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 03-21-2024, 07:32 AM
                  0 responses
                  66 views
                  0 likes
                  Last Post seqadmin  
                  Working...
                  X