<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>

CC BY-SA 3.0
2017-12-14T10:39:06.600Z
So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only.
stack_learner
CC BY-SA 3.0
2017-12-14T10:42:52.913Z
Please show your output to these people.
benn
CC BY-SA 3.0
2017-12-14T10:46:16.243Z
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.
stack_learner
CC BY-SA 3.0
2017-12-14T10:49:19.447Z
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.
benn
CC BY-SA 3.0
2017-12-14T11:00:04.230Z
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
stack_learner
CC BY-SA 3.0
2017-12-14T11:05:08.360Z
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!
benn
CC BY-SA 3.0
2017-12-14T11:23:47.473Z
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`.
Devon Ryan
CC BY-SA 3.0
2017-12-14T12:53:47.253Z
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.
Ian Sudbery
CC BY-SA 3.0
2017-12-14T13:14:30.667Z
Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats.
stack_learner
CC BY-SA 3.0
2017-12-14T13:22:04.320Z
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?
stack_learner
CC BY-SA 3.0
2017-12-14T13:25:09.030Z
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).
benn
CC BY-SA 3.0
2017-12-14T13:36:17.377Z
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.
stack_learner
CC BY-SA 3.0
2017-12-14T13:37:27.340Z
The tmp dir has 6 input.sorted.nnnn.bam files.
stack_learner
CC BY-SA 3.0
2017-12-14T13:45:10.937Z
Please start a new question for that error.
benn
CC BY-SA 3.0
2017-12-14T10:07:31.273Z
2017-12-14T11:21:09.167Z
2017-12-14T11:21:09.167Z
benn
<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>

<rna-seq><sam><samtools>
Read counts from BAM file
# 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.  

---

