Answers.Body: <p>Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your <code>samtools view -H</code> code. With <code>samtools idxstats file.bam</code> you'll get the reads per chromosome (or per nucleotide sequence in your case).</p>
  
  <p>For use in <code>featureCounts</code> my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).</p>
  
  <p>If your data is properly aligned to hg19 your output should be something like this:</p>
  
  <pre><code>samtools view -H aligned2hg19.bam
  
  @SQ SN:chrM LN:16571
  @SQ SN:chr1 LN:249250621
  @SQ SN:chr2 LN:243199373
  @SQ SN:chr3 LN:198022430
  @SQ SN:chr4 LN:191154276
  @SQ SN:chr5 LN:180915260
  @SQ SN:chr6 LN:171115067
  @SQ SN:chr7 LN:159138663
  @SQ SN:chr8 LN:146364022
  @SQ SN:chr9 LN:141213431
  @SQ SN:chr10    LN:135534747
  @SQ SN:chr11    LN:135006516
  @SQ SN:chr12    LN:133851895
  @SQ SN:chr13    LN:115169878
  @SQ SN:chr14    LN:107349540
  @SQ SN:chr15    LN:102531392
  @SQ SN:chr16    LN:90354753
  @SQ SN:chr17    LN:81195210
  @SQ SN:chr18    LN:78077248
  @SQ SN:chr19    LN:59128983
  @SQ SN:chr20    LN:63025520
  @SQ SN:chr21    LN:48129895
  @SQ SN:chr22    LN:51304566
  @SQ SN:chrX LN:155270560
  @SQ SN:chrY LN:59373566
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T10:39:06.600Z
Answers.Text: So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only.
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T10:42:52.913Z
Answers.Text: Please show your output to these people.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T10:46:16.243Z
Answers.Text: But you mean the bam file I have doesn't have any chr, read count and other details? I was thinking the above posted one are from header and there is another command to get the readcounts and other details.
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T10:49:19.447Z
Answers.Text: It probably has read counts, but not per chromosome but per @SQ sequence (which are RNA sequences). With `samtools idxstats` you'll see the read counts per @SQ.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T11:00:04.230Z
Answers.Text: Ok. I want to use featurecounts in R to get the read counts. I have the bam files and hg19.fasta and hg19.fasta.fai files. May I know how I should give the command to get the counts for all the samples. Thanq
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T11:05:08.360Z
Answers.Text: Maybe you misunderstood, but you need to correctly align to hg19 first in order to do this in featureCounts. From your output it is clear that it is not aligned to hg19, but to a bunch of RNA sequences. Please discuss with your people that your data is not properly aligned to hg19 (show your output, this is proof), or try to align it yourself from fastq. When you have done that, you'll need a proper gtf annotation file from UCSC or ENSEMBL which is necessary for featureCounts. Good luck!
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T11:23:47.473Z
Answers.Text: Your options are to either align to the genome (as @b.nota said) or filter using samtools however you like (e.g., for MAPQ) and then get counts with `samtools idxstats`.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T12:53:47.253Z
Answers.Text: When you are told they are aligned to hg19, I guess they mean tat are aligned to the hg19 version of the transcriptome, not hg19 itself. In this case you done need featureCounts to get read counts (and indeed, it would be hard/impossible to do so). `samtools idxstats` should give you want you are looking for.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:14:30.667Z
Answers.Text: Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats.
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:22:04.320Z
Answers.Text: One question. Do you I need to sort the bam files? samtools idxstats what are the options I should give in the command to get read counts?
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:25:09.030Z
Answers.Text: Yes you'll need to sort bam files first. `samtools idxstats aln.sorted.bam` will give you the results, see [here](http://www.htslib.org/doc/samtools.html).
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:36:17.377Z
Answers.Text: Yes, I saw the manual and tried like this and endup with an error.    
   samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam. This gave me [W::bam_hdr_read] EOF marker is absent. The input is probably truncated; [E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes;
  [bam_sort_core] truncated file. Aborting.
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:37:27.340Z
Answers.Text: The tmp dir has 6 input.sorted.nnnn.bam files.
Answers.UserDisplayName: stack_learner
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T13:45:10.937Z
Answers.Text: Please start a new question for that error.
Answers.UserDisplayName: benn
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T10:07:31.273Z
Answers.LastActivityDate: 2017-12-14T11:21:09.167Z
Answers.LastEditDate: 2017-12-14T11:21:09.167Z
Answers.UserDisplayName: benn
Body: <p>I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:</p>
  
  <pre><code>samtools view -H IonXpress_input.bam
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>@HD VN:1.4  SO:coordinate   GO:none
  @SQ SN:NR_039978    LN:984
  @SQ SN:NM_130786    LN:1766
  @SQ SN:NM_138932    LN:9293
  @SQ SN:NM_033110    LN:2678
  @SQ SN:NM_000014    LN:4678
  @SQ SN:NM_144670    LN:5192
  @SQ SN:NR_040112    LN:1201
  @SQ SN:NM_017436    LN:2106
  @SQ SN:NM_016161    LN:1771
  @SQ SN:NM_015665    LN:1854
  </code></pre>
  
  <p>I guess these are from header. After these I see something like following:</p>
  
  <pre><code>@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4
  </code></pre>
  
  <p>Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?</p>
  
Tags: <rna-seq><sam><samtools>
Title: Read counts from BAM file
ThreadText: # Read counts from BAM file
  
  (Asked by: username_0 on 2017-12-14)
  
  I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:
  
  ```
  samtools view -H IonXpress_input.bam
  
  ```
  
  Output:
  
  ```
  @HD VN:1.4  SO:coordinate   GO:none
  @SQ SN:NR_039978    LN:984
  @SQ SN:NM_130786    LN:1766
  @SQ SN:NM_138932    LN:9293
  @SQ SN:NM_033110    LN:2678
  @SQ SN:NM_000014    LN:4678
  @SQ SN:NM_144670    LN:5192
  @SQ SN:NR_040112    LN:1201
  @SQ SN:NM_017436    LN:2106
  @SQ SN:NM_016161    LN:1771
  @SQ SN:NM_015665    LN:1854
  
  ```
  
  I guess these are from header. After these I see something like following:
  
  ```
  @RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4
  
  ```
  
  Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?
  
  **Tags:** rna-seq, sam, samtools
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2017-12-14. Score: 7 | Accepted answer
  
  Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR\_\* or NM\_\* accession codes with your `samtools view -H` code. With `samtools idxstats file.bam` you'll get the reads per chromosome (or per nucleotide sequence in your case).
  
  For use in `featureCounts` my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).
  
  If your data is properly aligned to hg19 your output should be something like this:
  
  ```
  samtools view -H aligned2hg19.bam
  
  @SQ SN:chrM LN:16571
  @SQ SN:chr1 LN:249250621
  @SQ SN:chr2 LN:243199373
  @SQ SN:chr3 LN:198022430
  @SQ SN:chr4 LN:191154276
  @SQ SN:chr5 LN:180915260
  @SQ SN:chr6 LN:171115067
  @SQ SN:chr7 LN:159138663
  @SQ SN:chr8 LN:146364022
  @SQ SN:chr9 LN:141213431
  @SQ SN:chr10    LN:135534747
  @SQ SN:chr11    LN:135006516
  @SQ SN:chr12    LN:133851895
  @SQ SN:chr13    LN:115169878
  @SQ SN:chr14    LN:107349540
  @SQ SN:chr15    LN:102531392
  @SQ SN:chr16    LN:90354753
  @SQ SN:chr17    LN:81195210
  @SQ SN:chr18    LN:78077248
  @SQ SN:chr19    LN:59128983
  @SQ SN:chr20    LN:63025520
  @SQ SN:chr21    LN:48129895
  @SQ SN:chr22    LN:51304566
  @SQ SN:chrX LN:155270560
  @SQ SN:chrY LN:59373566
  
  ```
  
  > **Comments:**
  > **username_0**: So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only.  
  > **username_1**: Please show your output to these people.  
  > **username_0**: But you mean the bam file I have doesn't have any chr, read count and other details? I was thinking the above posted one are from header and there is another command to get the readcounts and other details.  
  > **username_1**: It probably has read counts, but not per chromosome but per @SQ sequence (which are RNA sequences). With `samtools idxstats` you'll see the read counts per @SQ.  
  > **username_0**: Ok. I want to use featurecounts in R to get the read counts. I have the bam files and hg19.fasta and hg19.fasta.fai files. May I know how I should give the command to get the counts for all the samples. Thanq  
  > **username_1**: Maybe you misunderstood, but you need to correctly align to hg19 first in order to do this in featureCounts. From your output it is clear that it is not aligned to hg19, but to a bunch of RNA sequences. Please discuss with your people that your data is not properly aligned to hg19 (show your output, this is proof), or try to align it yourself from fastq. When you have done that, you'll need a proper gtf annotation file from UCSC or ENSEMBL which is necessary for featureCounts. Good luck!  
  > **username_2**: Your options are to either align to the genome (as @b.nota said) or filter using samtools however you like (e.g., for MAPQ) and then get counts with `samtools idxstats`.  
  > **username_3**: When you are told they are aligned to hg19, I guess they mean tat are aligned to the hg19 version of the transcriptome, not hg19 itself. In this case you done need featureCounts to get read counts (and indeed, it would be hard/impossible to do so). `samtools idxstats` should give you want you are looking for.  
  > **username_0**: Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats.  
  > **username_0**: One question. Do you I need to sort the bam files? samtools idxstats what are the options I should give in the command to get read counts?  
  > **username_1**: Yes you'll need to sort bam files first. `samtools idxstats aln.sorted.bam` will give you the results, see [here](http://www.htslib.org/doc/samtools.html).  
  > **username_0**: Yes, I saw the manual and tried like this and endup with an error.      samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam. This gave me [W::bam_hdr_read] EOF marker is absent. The input is probably truncated; [E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes; [bam_sort_core] truncated file. Aborting.  
  > **username_0**: The tmp dir has 6 input.sorted.nnnn.bam files.  
  > **username_1**: Please start a new question for that error.  
  
  ---
  
  