[Camg]

Yes, I think you already found your problem. You have reads mapping to contigs in your reference file, but these contigs must not have any genes annotated on them and therefore are not present in your .gff file.

You could check this by removing these contigs from your reference file and running your whole pipeline again, or just using an awk script to remove all the SAM lines with alignments to one of these contigs, then run HT-Seq again. That should get rid of a big proportion of the "no_feature" counts.

If you want to get useful information from these reads you're going to have to annotate genes on these contigs, which you could do with something like cufflinks.

Good luck

[Alex.Manuel]

Damn! I just saw that default option in htseq-count is strand-specific...and as I used the RNA TruSeq V2 PE Library prep. I can be as good as certain that it is unstranded... so did HTseq count with the --stranded=no option and get less no_feature counts. Here is the result:

count sum 2464902
no_feature 728598
ambiguous 78167
too_low_aQual 0
not_aligned 0
alignment_not_unique 862188

Before (without --stranded=no)

see above

So the result looks better. And I thought about my question, if it is not possible to get less no feature genes, but I think when the gene is not yet annotated then there is no Gene ID and it wont make sense to include genes without an ID. So I have to live with that...
Can you confirm my thoughts?

Cheers Alex

[Alex.Manuel]

Dear Seqanswer community, dear Simon,

I'm dealing with stickleback RNA-Seq data. Mapping was done with Tophat2. The ref.fasta and the anno.gtf (both from ensembl) were included in the tophat options.
./tophat -p 8 -G stickle_gen_anno.gtf -o 193_ctrl24_CT_Shared_thout stickle_genome 193_ctrl24_CT_1_Shared.fastq 193_ctrl24_CT_2_Shared.fastq

I used samtools for sorting and converting the accepted_hits.bam (from tophat) to sam format and did the HTseq-count analysis with the anno.gtf I already used in tophat.
Coming to my recent problem. Is it normal that so many genes have no feature? So more than genes that are counted?
Here as an example (with one line added to see haw many genes are counted in total:

count sum 1266408
no_feature 2005250
ambiguous 9
too_low_aQual 0
not_aligned 0
alignment_not_unique 862188

So first I checked if the chromosomes displayed in the warning messages of htseq_count
e.g. Warning: Skipping read 'H108:2920J82ACXX:1:1101:3626:48576', because chromosome 'scaffold_893', to which it has been aligned, did not appear in the GFF file.
are not in the anno.gtf file. And indeed they are not. So they have to be just in the ref.fasta file.
But isn't there a possibility to count these genes? If not the loss of information is really big!?

Many thanks for an answer!

Cheers Alex

[rboettcher]

Dear all,

I encounter the same error message ('SAM line does not contain at least 11 tab-delimited fields.') as mentioned before when running htseq-count with this command:

htseq-count merged_align.bam UCSC_annotation_hg19.gtf > htseq_test.txt

However, when I inspect the provided bam file via samtools view it looks like this:

HWI-ST301L:234:C0GECACXX:3:1102:2230:56147 355 chr1 11164 1 101M = 11274 211 GACCGCCCCTTGCTTGCAGCCGGGCACTACAGGACCCGCTTGCTCACGGTGCTGTGCCAGGGCGCCCCCTGCTGGCGACTAGGGCAACTGCAGGGCTCTCT * AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:101 YT:Z:UU NH:i:3 CC:Z:chr15 CP:i:102519906 HI:i:0

The alignments have been done via Tophat2 + bowtie2 (most recent releases).
Does anyone have an idea why this error occurs?

Best regards

[himanshu04]

Thanks a lot for your replies. I had used the -S option funnily though when I used --sam option it did work out fine. Thanks a lot for your feedback I really appreciate it.

[Camg]

Did you use the -S option when running Bowtie? That's what will give you the standard SAM output file. Otherwise I think it gives you the SAM-like file that you currently have. Try running your Bowtie command again with -S and then use that file with HTseq.

[himanshu04]

I created it using Bowtie . I am actually following the same method as described in the seqwiki page. What do you think I should be doing?

[Camg]

I think your "SAM" file has a bit of a strange format. What alignment program did you use to create it?
Take a look at the example of my SAM file a few lines back and you can see some differences right away.

[himanshu04]

htseq-count -m intersection-strict -s no aligned_074262 TAIR10_GTF.gtf > counts_074262
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
478324 GFF lines processed.
Error occured when reading first line of sam file.
Error: ('SAM line does not contain at least 11 tab-delimited fields.', 'line 1 of file aligned_074262')
[Exception type: ValueError, raised in _HTSeq.pyx:1220]

The file looks like
SRR074262.3 HWI-EAS121:2:1:14:1646 + Chr3 328884 AGGGAGGGAGGGAGTGAGATTGNTTCGATCGCCAAT BCCCBCCBACCCBCACBCBBC:!<B@A<?>@:>8AA 0 22:A>N
SRR074262.4 HWI-EAS121:2:1:15:1470 + Chr1 14296071 CTGGGTTTTTGTGTTATTGAGANTCTGAGTTTGAGA BBA>?96>@@>==?A?@A<<=4!.5;7268?@607= 0 22:G>N

[Simon Anders]

I get an error saying unable to read first line of the sam file.

Well, how does the first line of your sam file look like?

[himanshu04]

Dear Simon,
I am using HT-Seq counts for summarizing the data of arabidopsis gene. When i use the command: 'htseq-count -m intersection-strict -s no aligned_074262 TAIR10_GTF.gtf > counts_074262'.
I get an error saying unable to read first line of the sam file. unable to determine 11 delimited tabs in the file. What can be the possible solution for this?
Thanks.

[Simon Anders]

The default value for 'type' is 'exon', not 'CDS'. This is because a GTF files usually have two lines for each exon, one labelled 'exon' and the other 'CDS', the difference being that the former includes untranslated regions and the latter only the coding part of the exon. (For interior exons, both lines contain the same coordinates.) I felt that most users would want htrseq-count to count a read also if it maps to a gene's UTR and this is why the default is 'exon', not 'CDS'.

In your GTF file, all the 'exon' lines seem to be missing, so you hit on the right solution to tell htseq-count to use the 'CDS' lines instead.

[Camg]

So I was messing around with parameters and found that when I changed --type=CDS to match my GTF file the program ran perfectly. I guess I was assuming CDS was the standard feature in column 3 of a GTF file since that's what is shown in the example GTF (link in manual). Also, I'm not sure why not specifying the --type attribute resulted in the error message I was getting.

Anyways, it is working now. Thanks for this great tool.

[Camg]

Error running HTSeq

Hi,

I'm having the same problem as the original post in this thread. When I run HTSeq I'm getting the error:

Warning: Skipping read 'HWI-ST765:7:2307:9917:51869#0', because chromosome 'chr1', to which it has been aligned, did not appear in the GFF file.
Warning: Skipping read 'HWI-ST765:7:2307:9917:76673#0', because chromosome 'chr1', to which it has been aligned, did not appear in the GFF file.

This is the command I used:

python ../../programs/HTSeq-0.5.3p3/HTSeq/scripts/count.py --mode=intersection-nonempty --stranded=no --idattr=transcript_id -o 24.1.htseq 24.1.sorted.sam ../ec_txn.gtf

However, I can see that my SAM file and GTF file are using the same chromosome names so I'm not sure why HTSeq isn't recognizing them.

Here are examples of a few lines from my SAM and GTF files:

HWI-ST765:7:1101:10001:176633#0 147 chr1 97536 255 101M = 97470 -167 ACAAAGATGTCGCACTTGGCATTGATGTACATGCAGACCTCCTGGCAGACCTTGCGCTCGATGTGCTTCTTGAGCTTGTGGAGCTGTCCTGGGCATCGGAG DDDDDDDDDDDDCCDDDDDDDDDDEEFEDDDDDDDCADDFEFFEHHHHHJIJJJIGGIJJIJIHEIJIJJJJJJJJHEJJJJIJJJJIHBHHHFFDB;CC@ NM:i:0 NH:i:1
HWI-ST765:7:1101:10001:176633#0 99 chr1 97470 255 101M = 97536 167 CGGAGGTTTCTGTATCTTCTCCTGTCGTTAACCACGACTCTGTAGAAGTCGCAGTCCCCGACCAGCACAAAGATGTCGCACTTGGCATTGATGTACATGCA BCCFFFBDHHHHHJJJJJJJJJJJJJIFIIJJJJJJIJJJGGIJIJJJHJJIJJIIIJJHHFFDEDDDDDDCDDDCCBDDDDDDDDCDDDDEDEEEEEDDD NM:i:1 NH:i:1
HWI-ST765:7:1101:10002:95570#0 147 chr1 107247 255 101M = 107184 -164 TCCTTCCAAGCTCCAAGCAGCTTCCAATCGAAATGACAGGAATGTTCTTCTTTTTTGCCAGCGCCATGACAAGCTTGGAAATTCTTGGCTCTGCATTCTCG BDDDCDDDDDDDDDDDDDDDDDDDDDEDEEEEFEFFFFFFFHHHHHIIGIIIJIJJJJJJJJJJIJIJIIJJJIJIJIJJJJJJJJJJHHHHHFFFFFCCC NM:i:0 NH:i:1
HWI-ST765:7:1101:10002:95570#0 99 chr1 107184 255 101M = 107247 164 CGGCAACACAGCACCCCTTGCTCCTCACCTTCCCGCTGGAGCTGACATTCTCCACACCGACAATCCTTCCAAGCTCCAAGCAGCTTCCAATCGAAATGACA CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJIIIIJFHIJIJIIIJGGJIIJIGHHBECDDDCDDDDDDDDDCDDDDDDDDDDDDDDDDDDDDDDDDD NM:i:0 NH:i:1

chr1 tophat CDS 4117 4593 . + 0 gene_id "ECU01_0010"; transcript_id "ECU01_0010";
chr1 tophat CDS 4391 5086 . - 0 gene_id "ECU01_0020"; transcript_id "ECU01_0020";
chr1 tophat CDS 5014 7486 . - 0 gene_id "ECU01_0023"; transcript_id "ECU01_0023";
chr1 tophat CDS 7536 8832 . - 0 gene_id "ECU01_0025"; transcript_id "ECU01_0025";
chr1 tophat CDS 11158 11529 . - 0 gene_id "ECU01_0030"; transcript_id "ECU01_0030";

Let me know if you have any ideas on how to solve this issue.
Thanks in advance

[sunkorner]

Thank you very much and also sorry. I am new to RNAseq data analysis and so somehow must have overlooked the details.