SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Read counts from SAM file mapped to de novo assembled transcripts using HTSeq-count alan_sm RNA Sequencing 2 06-12-2015 09:54 PM
htseq read counts different in separate runs h_manoj Bioinformatics 0 02-06-2014 11:46 AM
Can CLC genomics read mapping files be used in Bioconductor/R and HTSeq-counts? tdelaney Bioinformatics 1 02-20-2013 10:07 PM
htseq-count with warning for every read to represent all of zero counts in output hibachings2013 RNA Sequencing 10 07-15-2011 11:19 AM
De novo assembly of human genomes with massively parallel short read sequencing dan Literature Watch 0 12-21-2009 05:40 AM

Reply
 
Thread Tools
Old 07-30-2014, 11:29 AM   #1
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default HTSeq, human genomes and low read counts: am I doing anything wrong?

This is my basic workflow:

DATA and ALIGNMENT:

50SE Illumina reads mapped with tophat2 (v 2.0.11; default setting) to a bowtie2 index of hg38 (all chromosomes plus mitochondria, but no alternate chromosomal assemblies)

outfile reports as follows:
Quote:
Reads:
Input : 29397870
Mapped : 28589633 (97.3% of input)
of these: 5179700 (18.1%) have multiple alignments (612 have >20)
97.3% overall read mapping rate.
READ COUNTING

1) Generated a "corrected" gtf file from UCSC (i.e., geneID and transcript IDs are different) using the method described here.
2) converted tophat alignment bam to sam (samtools)
3) counted reads with HTSeq using the general command:
Quote:
htseq-count --stranded=no accepted_hits.sam hg38_UCSC_alt.gtf >HTSeq.count
My HTSeq output reports as follows:

Quote:
no_feature 13107880
ambiguous 5784082
too_low_aQual 0
not_aligned 0
alignment_not_unique 25653910
The result is that I have ~4.5 million unique alignments as counted by HTSeq. That seems low and when I query the bam file,

Quote:
samtools view accepted_hits.bam | grep "NH:i:1" -c
the output reports 25707742 alignments as unique (NH:i:1)

Can anyone offer a simple explanation of what is happening here? The hand-wavy explanation that I have at the moment is that HTSeq count is very conservative, and if a read cannot be assigned unambiguously, and exclusively to a gene, it is not counted. However, the HTSeq output is secretly making me question the integrity of my pipeline and/or my input files.

Thank you for any advice or reassurance.

David

EDIT: and it occurs to me that it might be relevant that the library prep was RiboMinus. Since I'm used to working with poly-A purified samples, my expectations might be too high in this case.

Last edited by ExMachina; 07-30-2014 at 11:42 AM.
ExMachina is offline   Reply With Quote
Old 07-30-2014, 12:07 PM   #2
Xianglinwu
Junior Member
 
Location: Rochester mn

Join Date: Jul 2014
Posts: 1
Default Suresect adapter sequence

what is sureselect adapter sequences? indexing strand of Sureselect is different from indexing strand of Truseq LT, can anyone give some information?
Xianglinwu is offline   Reply With Quote
Old 07-31-2014, 05:26 AM   #3
bioBob
Member
 
Location: Virginia

Join Date: Mar 2011
Posts: 72
Default

Love the location. You should sign your name the same way.

I have not done the 'correction' method you have referenced. I normally download the reference and annotation file from igenomes in tophat ready format. Can you post a snippet from your corrected gtf file?

I believe HTSeq will now output a samfile with a field indicating what the read contributed to , have you looked at that and looked to see where they map?

-o <samout>, --samout=<samout>
write out all SAM alignment records into an output SAM file called <samout>, annotating each line with its assignment to a feature or a special counter (as an optional field with tag ‘XF’)
bioBob is offline   Reply With Quote
Old 07-31-2014, 05:31 AM   #4
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default

One other change I just tried was the HTSeq option "--mode=intersection-nonempty "

I was hoping to get more information, but this option only decreased my "no_feature" count by about 5000. So I'm still left with the question: is it reasonable to have only ~15% of my mapped reads associate with annotated genes?
ExMachina is offline   Reply With Quote
Old 07-31-2014, 05:33 AM   #5
bioBob
Member
 
Location: Virginia

Join Date: Mar 2011
Posts: 72
Default

No, that is not reasonable. Can you post a few lines of your corrected gtf?
bioBob is offline   Reply With Quote
Old 07-31-2014, 05:48 AM   #6
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default

Quote:
Originally Posted by bioBob View Post
I have not done the 'correction' method you have referenced. I normally download the reference and annotation file from igenomes in tophat ready format. Can you post a snippet from your corrected gtf file?
Here's a small chunk from chr1. As you can see some genes have corrected formatting but others (generally non-coding, I think) remain uncorrected:

Quote:
chr1 knownGene exon 69091 70008 . + . gene_id "Q8NH21"; transcript_id "uc001aal.1"; exon_number "1"; exon_id "uc001aal.1.1"; gene_name "Q8NH21";
chr1 knownGene CDS 69091 70005 . + 0 gene_id "Q8NH21"; transcript_id "uc001aal.1"; exon_number "1"; exon_id "uc001aal.1.1"; gene_name "Q8NH21";
chr1 knownGene start_codon 69091 69093 . + 0 gene_id "Q8NH21"; transcript_id "uc001aal.1"; exon_number "1"; exon_id "uc001aal.1.1"; gene_name "Q8NH21";
chr1 knownGene stop_codon 70006 70008 . + 0 gene_id "Q8NH21"; transcript_id "uc001aal.1"; exon_number "1"; exon_id "uc001aal.1.1"; gene_name "Q8NH21";
chr1 knownGene exon 134773 139696 . - . gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "1"; exon_id "uc021oeg.2.1"; gene_name "B4DF06";
chr1 knownGene CDS 138533 139696 . - 0 gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "1"; exon_id "uc021oeg.2.1"; gene_name "B4DF06";
chr1 knownGene exon 139790 139847 . - . gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "2"; exon_id "uc021oeg.2.2"; gene_name "B4DF06";
chr1 knownGene CDS 139790 139792 . - 0 gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "2"; exon_id "uc021oeg.2.2"; gene_name "B4DF06";
chr1 knownGene exon 140075 140566 . - . gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "3"; exon_id "uc021oeg.2.3"; gene_name "B4DF06";
chr1 knownGene start_codon 139790 139792 . - 0 gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "1"; exon_id "uc021oeg.2.1"; gene_name "B4DF06";
chr1 knownGene stop_codon 138530 138532 . - 0 gene_id "B4DF06"; transcript_id "uc021oeg.2"; exon_number "1"; exon_id "uc021oeg.2.1"; gene_name "B4DF06";
chr1 knownGene exon 182393 182746 . + . gene_id "uc031tlc.1"; transcript_id "uc031tlc.1"; exon_number "1"; exon_id "uc031tlc.1.1";
chr1 knownGene exon 183132 183240 . + . gene_id "uc031tlc.1"; transcript_id "uc031tlc.1"; exon_number "2"; exon_id "uc031tlc.1.2";
chr1 knownGene exon 183740 184878 . + . gene_id "uc031tlc.1"; transcript_id "uc031tlc.1"; exon_number "3"; exon_id "uc031tlc.1.3";
I'm also wondering if we're jumping the gun on trying to use the newer hg38 assembly and annotation? As you have alluded, there are more refined input files available for hg19/37.

Currently I am waiting on the imminent ENSEMBL release of 38 to run the same pipeline on an compare results to the UCSC-based results.


Quote:
I believe HTSeq will now output a samfile with a field indicating what the read contributed to , have you looked at that and looked to see where they map?

-o <samout>, --samout=<samout>
write out all SAM alignment records into an output SAM file called <samout>, annotating each line with its assignment to a feature or a special counter (as an optional field with tag ‘XF’)
Wow, that's a great tip. That's exactly what I need to have to see what's going on. Thanks!
ExMachina is offline   Reply With Quote
Old 07-31-2014, 06:56 AM   #7
bioBob
Member
 
Location: Virginia

Join Date: Mar 2011
Posts: 72
Default

If it were me, I would align and count to the previous release just to see what the count assignments look like in terms of the number of reads in no_feature etc.

Another thing, are you sure the chromosome id's are all identical in the fasta and gtf? In the build I have, the chromosomes are labeled as 1, 2 etc instead of chr1, chr2 etc.
bioBob is offline   Reply With Quote
Old 07-31-2014, 07:00 AM   #8
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,080
Default

Give "featureCounts" a try, while you are at it.

As bioBio pointed out, check to make sure your chrom ID's match in your ref/BAM/GTF files.
GenoMax is offline   Reply With Quote
Old 07-31-2014, 07:07 AM   #9
bioBob
Member
 
Location: Virginia

Join Date: Mar 2011
Posts: 72
Default

Oh, and make sure the reads are sorted in the bam file. HTSeq now has a flag for which method the sorting was performed, ie name or position.
bioBob is offline   Reply With Quote
Old 07-31-2014, 07:19 AM   #10
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default

Quote:
If it were me, I would align and count to the previous release just to see what the count assignments look like in terms of the number of reads in no_feature etc.
Good idea and I already tried that
Same general results

Quote:
Another thing, are you sure the chromosome id's are all identical in the fasta and gtf? In the build I have, the chromosomes are labeled as 1, 2 etc instead of chr1, chr2 etc.
The "1, 2..." IDs are ENSEMBL while UCSC uses "chr1, chr2..."

And good point on the sorting bioBob. I checked and all my bam files are sorted ( SO:coordinate)
ExMachina is offline   Reply With Quote
Old 07-31-2014, 07:24 AM   #11
bioBob
Member
 
Location: Virginia

Join Date: Mar 2011
Posts: 72
Default

Right, so you need to set that flag as the default is name. You probably had a lot of messages in the console about mate not found.

-r <order>, --order=<order>¶
For paired-end data, the alignment have to be sorted either by read name or by alignment position. If your data is not sorted, use the samtools sort function of samtools to sort it. Use this option, with name or pos for <order> to indicate how the input data has been sorted. The default is name.
bioBob is offline   Reply With Quote
Old 07-31-2014, 07:39 AM   #12
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default

Quote:
Originally Posted by bioBob View Post
Right, so you need to set that flag as the default is name. You probably had a lot of messages in the console about mate not found.

-r <order>, --order=<order>¶
For paired-end data, the alignment have to be sorted either by read name or by alignment position. If your data is not sorted, use the samtools sort function of samtools to sort it. Use this option, with name or pos for <order> to indicate how the input data has been sorted. The default is name.
Interesting. It looks like "-r" and "--order" are not options in my version (0.5.4). I wonder if that's the problem right there?!

EDIT: actually, I don't think this is a problem for me here, since all I have are SE reads.

Last edited by ExMachina; 07-31-2014 at 12:22 PM.
ExMachina is offline   Reply With Quote
Old 07-31-2014, 12:22 PM   #13
ExMachina
Member
 
Location: AATGCATCTCATGTTATCCTTCTCGAG

Join Date: Mar 2010
Posts: 14
Default

Thanks for the input on this. As an update, here's what I've figured out:

The UCSC gtf file should not be made with the "knownGenes" table but should instead be made from the "refGene" table--the "knownGenes" table has too many overlapping coordinates.

Using the refGene annotations, I now get ~33% of my mapped reads mapping (uniquely) to known genes. I feel a little more confident with this result, especially given the RiboMinus library prep.

More comments are always welcome and I will report back here once the new ENSEMBL annotation is released.

Thanks for the help!

-David
ExMachina is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:36 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO