Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Read coverage for specific region of the genome

    I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

    To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

    Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

    samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example

  • #2
    Originally posted by Kaas View Post
    I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

    To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

    Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

    samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example
    Hi- You can use coverageBed in bedtools http://bedtools.readthedocs.org/en/latest/:

    Code:
    ## Bed file of intervals where to count reads (i.e. your intron):
    cat b.bed
    NW_003614142.1 35303 35329
    NW_003614142.1 35250-35302
    NW_003614142.1 35329 -35359
    ...
    
    ## Count reads in bed intervals
    coverageBed -abam accepted_hits.bam -b b.bed

    Comment


    • #3
      Hi

      Thank you for the answer. That worked fine... but the results did not match what i found in CLC genomic workbench. Seems that it can only be used for DNA and not RNA data with introns as it probably only use the information in the bam file of start of annealing and length of annealing. Results can be seen here, but the intron has the same coverage as the flanking regions using bedtools.

      Comment


      • #4
        Try samtools view with the -c parameter, which will generate a count summary of your regions.

        Comment


        • #5
          Originally posted by Kaas View Post
          I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

          To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

          Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

          samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example
          You can do this using the featureCounts program included in Subread package (http://subread.sourceforge.net).

          Code:
          # create a tab-delimited annotation file ("annot.txt")
          GeneID	Chr	Start	End	Strand
          1	NW_003614142.1	35303	35329	+
          2	NW_003614142.1	35250	35302	+
          3	NW_003614142.1	35329	35359	+
          
          # count reads overlapping with each region
          featureCounts -a annot.txt -i your_reads.bam -o counts.txt -F 'SAF' -b -f -O

          Comment

          Latest Articles

          Collapse

          • seqadmin
            Best Practices for Single-Cell Sequencing Analysis
            by seqadmin



            While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
            06-06-2024, 07:15 AM
          • 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

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by seqadmin, Today, 07:24 AM
          0 responses
          9 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, Yesterday, 08:58 AM
          0 responses
          11 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 06-12-2024, 02:20 PM
          0 responses
          15 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 06-07-2024, 06:58 AM
          0 responses
          182 views
          0 likes
          Last Post seqadmin  
          Working...
          X