Announcement

Collapse

Welcome to the New Seqanswers!

Welcome to the new Seqanswers! We'd love your feedback, please post any you have to this topic: New Seqanswers Feedback.
See more
See less

htseq-count for sam and gff3

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • htseq-count for sam and gff3

    Hi everyone,

    I need some help with running htseq-count.

    I am running htseq-count for a sam file obtained as a result of Illumina reads (not paired) aligned to the genome. I am getting the following Warning for all my reads:

    "Warning: Skipping read 'S5_057841196', because chromosome 'SL2.40ch00', to which it has been aligned, did not appear in the GFF file."

    and of course, all of my reads go in the "alignment_not_unique" pile, so there is no gene count at the end.

    I read a similar problem posted by hibachings2013 in 2010, but unlike his problem, my sam and gff3 files BOTH have the correct chromosome names.

    What am I doing wrong?

    Below are examples of entries in my sam and gff3 files:

    @HD VN:1.0 SO:coordinate
    @SQ SN:SL2.40ch00 LN:21805821
    @SQ SN:SL2.40ch01 LN:90304244
    @SQ SN:SL2.40ch02 LN:49918294
    @SQ SN:SL2.40ch03 LN:64840714
    @SQ SN:SL2.40ch04 LN:64064312
    @SQ SN:SL2.40ch05 LN:65021438
    @SQ SN:SL2.40ch06 LN:46041636
    @SQ SN:SL2.40ch07 LN:65268621
    @SQ SN:SL2.40ch08 LN:63032657
    @SQ SN:SL2.40ch09 LN:67662091
    @SQ SN:SL2.40ch10 LN:64834305
    @SQ SN:SL2.40ch11 LN:53386025
    @SQ SN:SL2.40ch12 LN:65486253
    @PG ID:TopHat VN:1.3.1 CL:/usr/local/bin/tophat -I 5000 --segment-mismatches 1 -o ./Output/tophatSA1_S6/ -G ./Annotation/ITAG2.3_gene_models.gff3 ./Reference/S_lycopersicm_genome ./Input/SA1_S6
    S6_021828409 272 SL2.40ch00 96314 0 49M * 0 0 GACTCTTGAATACAATCTTACAATTTTTCCTCACAAATTGCTACACCCA IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:2 NH:i:7 CC:Z:SL2.40ch02 CP:i:13220492 HI:i:0
    S6_028644958 256 SL2.40ch00 225372 0 49M * 0 0 CCGACACACTAAATAAAAGAACAATATCACATTGCATATCAAACTAATA IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:1 NH:i:11 CC:Z:= CP:i:3434611 HI:i:0
    S6_008579902 272 SL2.40ch00 251575 0 49M * 0 0 CTTATGATTTTAATATGAACACATTTATCACTTTCATCATTCTTCGATC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:2 NH:i:10 CC:Z:SL2.40ch01 CP:i:5753251 HI:i:0
    S6_033284153 272 SL2.40ch00 251575 0 49M * 0 0 CTTATGATTTTAATATGAACACATTTATCACTTTCATCATTCTTCGATC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:2 NH:i:10 CC:Z:SL2.40ch01 CP:i:5753251 HI:i:0
    S6_084652203 272 SL2.40ch00 251582 0 49M * 0 0 TTTTAATATGAACACATTTATCACTTTCATCATTCTTCGATCCATTTGC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:2 NH:i:8 CC:Z:SL2.40ch01 CP:i:5753258 HI:i:0
    S6_015020076 272 SL2.40ch00 251583 0 49M * 0 0 TTTAATATGAACACATTTATCACTTTCATCATTCTTCGATCCATTTGCC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NM:i:2 NH:i:9 CC:Z:SL2.40ch01 CP:i:5753259 HI:i:0

    ##gff-version 3
    ##feature-ontology http://song.cvs.sourceforge.net/*che...?revision=1.93
    ##sequence-region SL2.40ch00 1 21805821
    SL2.40ch00 ITAG_eugene gene 16437 18189 . + . Alias=Solyc00g005000;ID=gene:Solyc00g005000.2;Name=Solyc00g005000.2;from_BOGAS=1;length=1753
    SL2.40ch00 ITAG_eugene mRNA 16437 18189 . + . ID=mRNA:Solyc00g005000.2.1;Name=Solyc00g005000.2.1;Note=Aspartic proteinase nepenthesin I (AHRD V1 **-- A9ZMF9_NEPAL)%3B contains Interpro domain(s) IPR001461 Peptidase A1 ;Ontology_term=GO:0006508;Parent=gene:Solyc00g005000.2;from_BOGAS=1;interpro2go_term=GO:0006508;length=1753;nb_exon=2
    SL2.40ch00 ITAG_eugene exon 16437 17275 . + . ID=exon:Solyc00g005000.2.1.1;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene five_prime_UTR 16437 16479 . + . ID=five_prime_UTR:Solyc00g005000.2.1.0;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene CDS 16480 17275 . + 0 ID=CDS:Solyc00g005000.2.1.1;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene intron 17276 17335 . + . ID=intron:Solyc00g005000.2.1.1;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene exon 17336 18189 . + 0 ID=exon:Solyc00g005000.2.1.2;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene CDS 17336 17940 . + 2 ID=CDS:Solyc00g005000.2.1.2;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene three_prime_UTR 17941 18189 . + . ID=three_prime_UTR:Solyc00g005000.2.1.0;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
    ###
    SL2.40ch00 ITAG_eugene gene 68062 68764 . + . Alias=Solyc00g005020;ID=gene:Solyc00g005020.1;Name=Solyc00g005020.1;from_BOGAS=1;length=703
    SL2.40ch00 ITAG_eugene mRNA 68062 68764 . + . ID=mRNA:Solyc00g005020.1.1;Name=Solyc00g005020.1.1;Note=Unknown Protein (AHRD V1);Parent=gene:Solyc00g005020.1;from_BOGAS=1;length=703;nb_exon=3;eugene_evidence_code=10F0H0E0IEG
    SL2.40ch00 ITAG_eugene exon 68062 68211 . + 0 ID=exon:Solyc00g005020.1.1.1;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene CDS 68062 68211 . + 0 ID=CDS:Solyc00g005020.1.1.1;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene intron 68212 68343 . + . ID=intron:Solyc00g005020.1.1.1;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene exon 68344 68568 . + 0 ID=exon:Solyc00g005020.1.1.2;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene CDS 68344 68568 . + 0 ID=CDS:Solyc00g005020.1.1.2;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene intron 68569 68653 . + . ID=intron:Solyc00g005020.1.1.2;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene exon 68654 68764 . + 0 ID=exon:Solyc00g005020.1.1.3;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    SL2.40ch00 ITAG_eugene CDS 68654 68764 . + 0 ID=CDS:Solyc00g005020.1.1.3;Parent=mRNA:Solyc00g005020.1.1;from_BOGAS=1
    Last edited by sofia17; 09-02-2011, 12:12 PM.

  • #2
    It might be because this is not a GTF file; it does not have 'gene_ID' attribute -- or rather, it has them, but they are named 'Parent'. Hence, you need the option '-i Parent'.

    Comment


    • #3
      Sorry, that was wrong. The GFF file contains a transcript ID in the 'Parent' attribute of 'exon' lines, but we would need the gene ID. You will need to fix this manually.

      Comment


      • #4
        so in the gff3 file I should replace entries such as

        "Parent=mRNA:Solyc00g005000.2.1"

        with

        "ID=gene:Solyc00g005000.2.1"???

        Comment


        • #5
          Nearly. It should be "ID=gene:Solyc00g005000.2". Compare carefully the "mRNA" and the "gene" lines.

          Comment


          • #6
            Within a gene should I do the above replacement for ALL components, i.e. mRNA, exon, intron, CDS, three_prime_UTR, etc...? Or for mRNA lines ONLY?
            Last edited by sofia17; 09-06-2011, 08:29 AM.

            Comment


            • #7
              Only for the "exon" lines, because htseq-count ignores the others.

              Comment


              • #8
                I have done the replacement (please see below that all exons have the correct 'ID=gene:...' field). But it still gives same error:
                Error occured in line 6 of file /Users/vanderlab1/Desktop/try/processed.gff.
                Error: Feature exon:Solyc00g005000.2.1.1 does not contain a 'gene_id' attribute
                Exception type: SystemExit, raised in count.py:55]

                I tried also: htseq-count -i ID=gene <samFile> <gff3File>
                and got the same error

                Error occured in line 6 of file /Users/vanderlab1/Desktop/try/processed.gff.
                Error: Feature exon:Solyc00g005000.2.1.1 does not contain a 'ID=gene' attribute
                Exception type: SystemExit, raised in count.py:55]

                Maybe I use wrongly the -i option?

                ##gff-version 3
                ##feature-ontology http://song.cvs.sourceforge.net/*che...?revision=1.93
                ##sequence-region SL2.40ch00 1 21805821
                SL2.40ch00 ITAG_eugene gene 16437 18189 . + . Alias=Solyc00g005000;ID=gene:Solyc00g005000.2;Name=Solyc00g005000.2;from_BOGAS=1;length=1753
                SL2.40ch00 ITAG_eugene mRNA 16437 18189 . + . ID=mRNA:Solyc00g005000.2.1;Name=Solyc00g005000.2.1;Note=Aspartic proteinase nepenthesin I (AHRD V1 **-- A9ZMF9_NEPAL)%3B contains Interpro domain(s) IPR001461 Peptidase A1 ;Ontology_term=GO:0006508;Parent=gene:Solyc00g005000.2;from_BOGAS=1;interpro2go_term=GO:0006508;length=1753;nb_exon=2
                SL2.40ch00 ITAG_eugene exon 16437 17275 . + . ID=exon:Solyc00g005000.2.1.1;ID=gene:Solyc00g005000.2;from_BOGAS=1
                SL2.40ch00 ITAG_eugene five_prime_UTR 16437 16479 . + . ID=five_prime_UTR:Solyc00g005000.2.1.0;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
                SL2.40ch00 ITAG_eugene CDS 16480 17275 . + 0 ID=CDS:Solyc00g005000.2.1.1;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
                SL2.40ch00 ITAG_eugene intron 17276 17335 . + . ID=intron:Solyc00g005000.2.1.1;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
                SL2.40ch00 ITAG_eugene exon 17336 18189 . + 0 ID=exon:Solyc00g005000.2.1.2;ID=gene:Solyc00g005000.2;from_BOGAS=1
                SL2.40ch00 ITAG_eugene CDS 17336 17940 . + 2 ID=CDS:Solyc00g005000.2.1.2;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
                SL2.40ch00 ITAG_eugene three_prime_UTR 17941 18189 . + . ID=three_prime_UTR:Solyc00g005000.2.1.0;Parent=mRNA:Solyc00g005000.2.1;from_BOGAS=1
                Last edited by sofia17; 09-08-2011, 07:34 AM.

                Comment


                • #9
                  Try '-i ID', not '-i ID=gene'. Only the part before the '=' is the attribute name.

                  Comment


                  • #10
                    Thank you Simon for “holding hands” . It worked with “–i ID”.

                    From your experience, is it normal that out of about 15 million Illumina reads 1.5 million go in the no_feature pile?

                    Comment


                    • #11
                      Originally posted by sofia17 View Post
                      From your experience, is it normal that out of about 15 million Illumina reads 1.5 million go in the no_feature pile?
                      Yes, that does not sound unusual.

                      Comment


                      • #12
                        Dear Simon,

                        I am back in business asking questions.

                        I ran at the R console

                        source("http://www.bioconductor.org/biocLite.R")
                        biocLite("DESeq")

                        to install DESeq but when I try to load the library with

                        library(DESeq)

                        it gives the following err:

                        Error in loadNamespace(i[[1L]], c(lib.loc, .libPaths())) :
                        there is no package called 'annotate'
                        Error: package/namespace load failed for 'DESeq'

                        I am running R 2.13.1 (latest version) on Mac OS 10.6.4. Any suggestions?

                        Comment


                        • #13
                          Oooops. Is working now. Maybe the bioconductor site was down earlier....

                          Comment


                          • #14
                            Problem with htseq-count

                            Dear SImon,
                            I am just starting using HTSeq and I found the same problem as Sofia, but anything of the solutions that you recommended works this time. CAn you help me? This is the command line I use:
                            [email protected]:~/HTSeq-0.5.3p3$ htseq-count -s no s_1_sequence_clipped_CDS_tophat.bam TbruceiTreu927_TriTrypDB-3.32.gff

                            And what I get:
                            Error occured in line 263 of file TbruceiTreu927_TriTrypDB-3.3.2.gff.
                            Error: Feature gene|exon_TB927.5.300b-1 does not contain a 'gene_id' attribute
                            [Exception type: SystemExit, raised in count.py:55]

                            I'd tried to use -i ID but I only obtain this:

                            [email protected]:~/HTSeq-0.5.3p3$ htseq-count -s no -i ID s_1_sequence_clipped_CDS_tophat.bam TbruceiTreu927_TriTrypDB-3.32.gff
                            Error occured in line 29875 of file TbruceiTreu927_TriTrypDB-3.32.gff.
                            Error: Failure parsing GFF attribute line
                            [Exception type: ValueError, raised in __init__.py:171]

                            I am working wih trypanosoma brucei from what I only have the gff data (there is no gtf data available) and when I look for the line this is the information I have in the gff starting the line 263 adn some more lines as example:


                            pidb|BAC26D11 ApiDB exon 53516 54478 . - . ID=apidb|exon_TB927.5.300b-1;Name=exon;description=exon;size=963;Parent=apidb|rna_TB927.5.300b-1
                            apidb|Tb927_04_v4 ApiDB gene 1458621 1459109 . - . ID=apidb|Tb04.24M18.150;Name=Tb04.24M18.150;description=hypothetical+protein%2C+conserved;size=489;web_id=Tb04.24M18.150;locus_tag=Tb04.24M18.150;size=489;Alias=20660476,Tb04.24M18.150
                            apidb|Tb927_04_v4 ApiDB mRNA 1458621 1459109 . - . ID=apidb|rna_Tb04.24M18.150-1;Name=Tb04.24M18.150-1;description=Tb04.24M18.150-1;size=489;Parent=apidb|Tb04.24M18.150;Ontology_term=GO:0044429;Dbxref=ApiDB:Tb04.24M18.150,taxon:185431
                            apidb|Tb927_04_v4 ApiDB CDS 1458621 1459109 . - 0 ID=apidb|cds_Tb04.24M18.150-1;Name=cds;description=.;size=489;Parent=apidb|rna_Tb04.24M18.150-1
                            apidb|Tb927_04_v4 ApiDB exon 1458621 1459109 . - . ID=apidb|exon_Tb04.24M18.150-1;Name=exon;description=exon;size=489;Parent=apidb|rna_Tb04.24M18.150-1
                            apidb|Tb927_04_v4 ApiDB gene 1312951 1313406 . - . ID=apidb|Tb04.3I12.100;Name=Tb04.3I12.100;description=hypothetical+protein;size=456;web_id=Tb04.3I12.100;locus_tag=Tb04.3I12.100;size=456;Alias=Tb04.3I12.100
                            apidb|Tb927_04_v4 ApiDB mRNA 1312951 1313406 . - . ID=apidb|rna_Tb04.3I12.100-1;Name=Tb04.3I12.100-1;description=Tb04.3I12.100-1;size=456;Parent=apidb|Tb04.3I12.100;Dbxref=ApiDB:Tb04.3I12.100,taxon:185431
                            apidb|Tb927_04_v4 ApiDB CDS 1312951 1313406 . - 0 ID=apidb|cds_Tb04.3I12.100-1;Name=cds;description=.;size=456;Parent=apidb|rna_Tb04.3I12.100-1
                            apidb|Tb927_04_v4 ApiDB exon 1312951 1313406 . - . ID=apidb|exon_Tb04.3I12.100-1;Name=exon;description=exon;size=456;Parent=apidb|rna_Tb04.3I12.100-1
                            apidb|Tb927_05_v4 ApiDB gene 345728 346237 . + . ID=apidb|Tb05.28F8.200;Name=Tb05.28F8.200;description=hypothetical+protein;size=510;web_id=Tb05.28F8.200;locus_tag=Tb05.28F8.200;size=510;Alias=Tb05.28F8.200
                            apidb|Tb927_05_v4 ApiDB mRNA 345728 346237 . + . ID=apidb|rna_Tb05.28F8.200-1;Name=Tb05.28F8.200-1;description=Tb05.28F8.200-1;size=510;Parent=apidb|Tb05.28F8.200;Dbxref=ApiDB:Tb05.28F8.200,taxon:185431
                            apidb|Tb927_05_v4 ApiDB CDS 345728 346237 . + 0 ID=apidb|cds_Tb05.28F8.200-1;Name=cds;description=.;size=510;Parent=apidb|rna_Tb05.28F8.200-1
                            apidb|Tb927_05_v4 ApiDB exon 345728 346237 . + . ID=apidb|exon_Tb05.28F8.200-1;Name=exon;description=exon;size=510;Parent=apidb|rna_Tb05.28F8.200-1
                            apidb|Tb927_05_v4 ApiDB gene 235483 235992 . - .


                            What may I change to make htseq-count work?

                            Thank you very much in advance!
                            Sandra

                            Comment


                            • #15
                              Using 'i ID' did solve your first problem, because afterward, the script no longer complained about line 263 of the GFF file but about line 29875. So, have a look there, too.

                              Your next problem will be that htseq-count, at the moment, does not read bam file. You will need to do something like
                              Code:
                              samtools view myalignment.bam | htseq-count [options] - myannotation.gff
                              The lone minus sign tells htseq-count to read the sam file from standard input, i.e., through the pipe from 'samtools view'.

                              Comment

                              Working...
                              X