Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • gsnap produces more reads than fastq

    hello, i am running gsnap with rna seq data against a reference genome:

    gsnap -m 10 -B 5 -t 8 -A sam -d GENOME seqs.fastq > res.sam

    but strangely, the output sam file contains more reads than the original fastq files:

    $ cat seqs.fastq | echo $((`wc -l`/4))
    3776979

    $ cat res.sam | grep -v '^ *@' | wc -l
    6009141



    and I don't understand what that means....

    thank you

  • #2
    I'm not familiar with gsnap but it seems like there could be a few factors, the most likely one being your grep command. When I try to find out how many alignments have been reported I use the command:
    Code:
    cut -f1 out.sam | sort | uniq | wc -l
    That way you don't count nonspecific alignments. That being said, your SAM file may contain multiple alignments for many reads. There should be a flag to control how many alignments per read are reported but if not there are ways around it. I hop this helps!

    Comment


    • #3
      Originally posted by twaddlac View Post
      I'm not familiar with gsnap but it seems like there could be a few factors, the most likely one being your grep command. When I try to find out how many alignments have been reported I use the command:
      Code:
      cut -f1 out.sam | sort | uniq | wc -l
      That way you don't count nonspecific alignments. That being said, your SAM file may contain multiple alignments for many reads. There should be a flag to control how many alignments per read are reported but if not there are ways around it. I hop this helps!
      Besides mapping to multiple locations, the sam file also has a header that can contain many lines and throw off your count.

      Comment


      • #4
        Originally posted by severin View Post
        Besides mapping to multiple locations, the sam file also has a header that can contain many lines and throw off your count.
        I forgot to mention, to omit the header lines you should do
        Code:
        grep -v '^@' out.sam | cut -f1 | sort | uniq | wc -l
        If you convert your SAM file to BAM you will:
        A) reduce the size of the file (BAM < SAM)
        B) view the alignment output without the headers
        Code:
        samtools view out.bam
        To convert the sam file to bam, just do
        Code:
        samtools view -bS out.sam > out.bam

        Comment


        • #5
          If a read is aligned to multiple location, it appears multiple times in the SAm file. Sort the file by read name, cut to the first coulmn (read name), pipe it through 'uniq', then count again.

          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
          7 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, Yesterday, 06:07 PM
          0 responses
          7 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