Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Simon Anders
    replied
    Also, remember that a read that aligns to multiple places will be listed in multiple lines in the SAM file.

    Leave a comment:


  • NitaC
    replied
    Ohhh, good to know. Thank you.

    Leave a comment:


  • maubp
    replied
    Originally posted by NitaC View Post
    Are you saying that after filtering for forward reads, there may still be reads in my results that have not actually mapped to the forward strand but are listed because their mate did?
    Yes in general. This is deliberate so that when coordinate sorted the unmapped read is next to its mapped partner in the BAM file. See "Recommended Practise for the SAM format" in the specification. Also "Bit 0x4 is the only reliable place to tell whether the segment is unmapped".

    Leave a comment:


  • NitaC
    replied
    Hi maubp!

    It is my understanding that 0x4 tells you that the read was unmapped and that 0x10 indicates that the read is on the reverse strand. I don't think either of those flags really address my issue. I've already filtered for unaligned reads.

    And I'm not sure I understand your second response. If I've filtered out unaligned reads and filtered for only the forward strand, ultimately I still need some sort of reference name or gene id to count right? Are you saying that after filtering for forward reads, there may still be reads in my results that have not actually mapped to the forward strand but are listed because their mate did?

    As it turns out, my little pipeline actually matches the results using htseq-count so I'm a little more confident in what I did. And my co-worker's counts were from one sample whereas mine were from all the samples combined. That little tidbit would've saved me from a lot of headaches.

    Leave a comment:


  • maubp
    replied
    You need to read more about the FLAG field in SAM/BAM. Here bits 0x4 and 0x10 will tell you if a read is mapped and to which strand.

    You can't just look at column 3 (RNAME), as that can be a reference name if the read's mate was mapped even if the read itself was not.

    Leave a comment:


  • NitaC
    started a topic procedure for extracting raw read counts

    procedure for extracting raw read counts

    Good morning everyone.

    I am a newb and this is my first time touching any sequencing data. It just so happens that the best way to learn is to dive right in. Incidentally I have hit a bit of a stumbling block and was hoping someone here could shed some light on it.

    I am trying to extract the raw read counts for transcripts mapped only to the forward strand. The mapping was done using Gsnap (I did not do the mapping.) We did paired-end sequencing using Illumina's HiSeq platform. Again, this was my first time working with sam/bam files and it took some time to understand the formats. Here is what I did:

    To get only the forward read counts I took the sam file I was supplied and converted it back to a bam file. Then using the 0x0040 flag I used a command to get just the forward reads and redirected that output to another sam file. Then, looking at column 3 in the sam file, I saw either the name of the transcript or a *. It is my understanding that the "*" is for those reads that did not align. So I wrote a little bash script to filter those out. I was told by a co-worker that using column 3 I could get the raw counts of the mapped transcripts. So, I wrote another little script to do that. I know my script is working correctly, however, my numbers do not match the numbers my co-worker has. In fact, I have counts for some transcripts that my co-worker reported as not detected at all. I know there are tools to automatically get counts, but I want to actually understand these file formats before using some black box of tools.

    This is how one line looks in the sam file:
    [anitabow@kirby2 spikeIn]$ head -n 35 DR38spikeInForward.sam|tail -n 1
    FCC1WG7ACXX:5:1101:3115:1979#GGACCCAT 99 ERCC-00130 851 40 100M = 958 207 CTTGAAGCAGATGATATCCCGATGCTGACGGCTTACAATAAACGTGATCAAAAACTGCCTGATTTTATACCGACCGCCGGAAGGGATCACATTATGGTCA CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIFFHHDGIJIJDEHJIIJJJJIHHIJJHHHHHHHFFFFFDDDDDDDDDDDDDDDDDDDDDDEDDCCD MD:Z:100 NH:i:1 HI:i:1 NM:i:0 SM:i:40 XQ:i:40 X2:i:0

    As you can see, the third column has the name of one transcript. This is what I have been using to extract my counts. So basically, my script will count every time it sees ERCC-00130 for example.

    Thanks in advance for any help.
    Last edited by NitaC; 04-29-2013, 05:57 AM. Reason: added an example of what I was talking about

Latest Articles

Collapse

  • seqadmin
    Exploring the Dynamics of the Tumor Microenvironment
    by seqadmin




    The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
    07-08-2024, 03:19 PM
  • seqadmin
    Exploring Human Diversity Through Large-Scale Omics
    by seqadmin


    In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
    06-25-2024, 06:43 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 06:53 AM
0 responses
12 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-10-2024, 07:30 AM
0 responses
34 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-03-2024, 09:45 AM
0 responses
203 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-03-2024, 08:54 AM
0 responses
213 views
0 likes
Last Post seqadmin  
Working...
X